基于經濟模型預測控制的金氰化浸出過程動態(tài)實時優(yōu)化

關宏偉，葉凌箭，沈非凡，顧德，宋執(zhí)環(huán)

（1 寧波財經學院機械與電氣工程學院，浙江寧波315175； 2 浙江大學寧波理工學院，浙江寧波315100； 3 江南大學輕工過程先進控制教育部重點實驗室，江蘇無錫214122； 4 浙江大學控制科學與工程學院，浙江杭州310027）

引 言

濕法冶金工業(yè)中，氰化浸出是從金礦石中提取金元素的一種常用工藝，具有結構簡單、效率高等諸多優(yōu)點[1-2]。傳統(tǒng)的金氰化浸出過程以人工操作為主，隨著自動化技術的發(fā)展，近年來對金氰化浸出過程的建模及自動控制方法陸續(xù)開展了研究。文獻[2-4]基于采集的工業(yè)過程數(shù)據研究了金氰化浸出過程的動力學模型。在相關工作的基礎上，文獻[5]使用浸出過程的數(shù)學模型開展了仿真研究。文獻[6]進一步討論了浸出槽體積的工藝優(yōu)化設計問題。這些方法有助于減少金氰化浸出過程的人工干預、提高運行效率，進而實現(xiàn)過程的安全、穩(wěn)定操作。

市場經濟背景下，通過操作優(yōu)化提高金氰化浸出過程的經濟效益顯得尤為重要。離線優(yōu)化方法基于系統(tǒng)標稱模型，可預先計算出最優(yōu)操作輸入，在一定程度上提高收率或節(jié)能降耗[7]。但真實過程存在模型不準確、不確定擾動等因素。在線實時優(yōu)化[8-11]是改進離線優(yōu)化缺點的手段之一。以預測控制為代表的先進控制方法，通過改善關鍵質量指標的控制效果，實現(xiàn)卡邊優(yōu)化[12]。這類技術主要針對的是約束變量，一般來說，在線實時優(yōu)化需要基于過程的非線性模型，進行擾動估計、反復優(yōu)化計算等計算量較大的步驟[13-15]。對此，張俊等[16-18]和葉凌箭等[19-20]分別基于修正項自適應（modifier adaptation）[21]和 自 優(yōu) 化 控 制（self-optimizing control）[22-26]研究了金氰化浸出過程的實時優(yōu)化問題，取得了有益效果。但這些研究對金氰化浸出過程的優(yōu)化均為穩(wěn)態(tài)優(yōu)化，即只追求系統(tǒng)處于穩(wěn)態(tài)時的操作最優(yōu)性。金氰化過程的時間常數(shù)大，穩(wěn)態(tài)等待時間可長達數(shù)十小時[27]，若只考慮穩(wěn)態(tài)優(yōu)化仍可能引起較大的經濟損失。如果能在暫態(tài)過程中實行動態(tài)優(yōu)化策略，則可以帶來更大的經濟收益。近年來得到了廣泛研究的經濟模型預測控制（EMPC）[28-29]，將經濟指標作為滾動優(yōu)化的目標函數(shù)，在每個采樣點直接優(yōu)化未來時間窗的操作輸入序列。和傳統(tǒng)的分層遞階結構中的模型預測控制相比，EMPC 同時集成了優(yōu)化和控制目標，具有更大的經濟潛力。在此背景下，本文基于EMPC 研究金氰化浸出過程的動態(tài)優(yōu)化問題，以改進現(xiàn)有穩(wěn)態(tài)優(yōu)化方法的不足。

模型準確性對EMPC 的性能發(fā)揮著關鍵作用，而金氰化浸出過程在運行時受到測量噪聲、不確定參數(shù)等大量不確定因素的影響。現(xiàn)有的實時優(yōu)化方法中，如傳統(tǒng)的參數(shù)自適應二步法[30]，首先利用測量值對不確定參數(shù)進行估計，再執(zhí)行優(yōu)化運算。但如果只利用穩(wěn)態(tài)測量值，同樣只具有穩(wěn)態(tài)優(yōu)化效果。本文在動態(tài)過程中對系統(tǒng)狀態(tài)和模型參數(shù)進行估計，并在此基礎上進行動態(tài)優(yōu)化，對參數(shù)估計的實時性提出了更高要求。文獻[31-32]使用Tikhonov 正則化方法估計動力學反應速率，并對相關參數(shù)進行辨識。在測量噪聲和不確定參數(shù)的共同影響下，本文提出使用擴展卡爾曼濾波（EKF）對系統(tǒng)狀態(tài)和模型參數(shù)進行動態(tài)估計。在卡爾曼濾波框架下，可以通過構造增廣系統(tǒng)，將不確定參數(shù)視為擴展狀態(tài)進行同時估計，提高了方法的適用性。和其他估計方法相比，基于卡爾曼濾波的方法在統(tǒng)計意義上具有最優(yōu)性，保證了狀態(tài)和參數(shù)估計效果。另一方面，此類方法在滾動時域上進行迭代估計，在線計算量不大，具有工程可行性。

1 金氰化浸出過程

1.1 過程概述

金濕法冶金包括氰化浸出、鋅粉置換和濃縮洗滌等工序。本文考慮其中第一道工序，即金氰化浸出工藝，該環(huán)節(jié)對整個煉金工藝中具有關鍵影響。金氰化浸出的典型結構如圖1 所示：緩沖箱中盛放經過調漿后的礦漿，用輸送泵打入n 級串聯(lián)的氣力浸出槽（后文取n=4），在溢流作用下礦漿從前一級浸出槽不斷流入后續(xù)浸出槽。在每一級浸出槽中均通入浸出劑（氰化鈉溶液），同時通入壓縮空氣提供反應所需的溶解氧，并產生氣力攪拌作用加速反應。產生的化學反應為

圖1 金氰化浸出過程Fig.1 Gold cyanidation leaching process

1.2 過程模型及優(yōu)化目標

基于物料守恒方程對金氰化浸出過程建模，包括如下固液相中的金元素、氰根離子守恒關系

式中，Cs、Cl和CCN分別是固相金品位、液相金濃度和液相氰離子濃度；Qs、Ql和QCN是礦漿固、液相流量和浸出劑氰化鈉的流量；Ms和Ml是停留在浸出槽中的固、液相質量；下角標i對應第i個浸出槽。

反應速率表示為如下形式

式中，rAu和rCN是金和氰離子的反應速率；k1～k6是反應速率參數(shù)；C∞是理想情況下的固相殘留金品位，是礦石平均粒徑d 的函數(shù)；Co是溶解氧濃度。此外，有如下關系

式中，Cw為礦漿中的固相濃度；V 是浸出槽體積；ρs和ρl分別是礦漿固、液相密度。

如果供氧成本可以忽略，則風機功率可調節(jié)為充分大，使溶液中的溶解氧濃度最大化以加快浸出反應。此外，假定浸出過程的礦石處理量恒定（Qs、Ql為常數(shù)），則金氰化浸出過程的操作變量為每個浸出槽的浸出劑流量。金氰化浸出過程的模型參數(shù)標稱值列于表1。

表1 金氰化浸出過程模型參數(shù)Table 1 Model parameters for gold cyanidation leaching process

綜上，金氰化浸出過程的操作目標描述為經濟指標最大化問題（成本最小），即求解如下最優(yōu)化問題

式中，J 為成本函數(shù)，由浸出劑成本、殘液處理費用和未浸出的金元素損失三部分構成，PCN、PCNd和PAu代表各自的單價；a 和a*=0.9 分別代表浸出率和最小允許的浸出率；QCN,max=10 kg/h 是能提供的最大浸出劑流量（輸送能力）。

2 經濟模型預測控制

2.1 離散系統(tǒng)優(yōu)化問題

現(xiàn)有文獻報道的金氰化浸出過程的操作優(yōu)化方法，大都只考慮穩(wěn)態(tài)優(yōu)化，即將狀態(tài)微分方程式（2）～式（4）中的微分方程設為0，再將式（9）轉換為小規(guī)模NLP 問題求解。為了進一步在動態(tài)過程中優(yōu)化經濟性能，經濟模型預測控制（EMPC）[28]延續(xù)了傳統(tǒng)MPC“預測模型—反饋校正—滾動優(yōu)化”的基本思路，而將經濟指標納入到MPC 的目標函數(shù)中，直接優(yōu)化操作輸入。

在k 時刻，考慮時間窗口[k,…,k+N]內的非線性動態(tài)優(yōu)化問題

其中，JN為時間窗內的總經濟性能指標，由終端指標φ(xk+N)和中間指標ψ(xi,ui)組成，x 和u 為狀態(tài)和控制變量，f和G分別為系統(tǒng)的動態(tài)非線性模型和約束條件，X、U和Xf為可行域。

式（10）表述為離散系統(tǒng)的動態(tài)優(yōu)化問題，金氰化浸出過程的模型方程式（2）～式（8）定義為微分代數(shù)方程（DAE）。因此，首先應將連續(xù)系統(tǒng)離散化，轉化為有限維問題。

圖2 直接配置法（n=3）Fig.2 Direct collocation method(n=3)

本文使用直接配置法（direct collocation）[31]進行離散化。如圖2 所示，在每一個時間間隔內選擇n個配置點，將配置點上的狀態(tài)變量x和代數(shù)變量z近似為多項式方程。時間窗[k,…,k+N]內，所有配置點上的狀態(tài)和代數(shù)變量擴充為

為保證狀態(tài)變量在時間間隔點的連續(xù)性，引入約束條件

對控制輸入u，在每個時間間隔內視為常數(shù)，得到擴展向量

經上述離散化后，可將浸出過程動態(tài)優(yōu)化問題轉化為有限維大規(guī)模NLP 問題，可使用典型的NLP求解器求解，如IPOPT等[32]。

2.2 基于EKF的狀態(tài)及不確定參數(shù)估計

上文所述的EMPC，假設了以下條件：（1）k 時刻，系統(tǒng)狀態(tài)xk能準確測量；（2）所有模型參數(shù)已知，滾動求解式（11）時沿時間軸變化的只有xk。而過程測量值包含大量隨機噪聲，并且系統(tǒng)受到各種擾動影響，參數(shù)具有不確定性和時變性，難以滿足前述條件。

擴展卡爾曼濾波（EKF）[33]是一種有效的非線性系統(tǒng)狀態(tài)估計方法，考慮如下非線性系統(tǒng)

式中，f 和h 分別是狀態(tài)和測量方程，過程擾動w ～N(0,Q) 和測量噪聲v ～N(0,R)服從高斯分布。在第k次迭代中，EKF執(zhí)行以下步驟[33]

（1）估計狀態(tài)和誤差協(xié)方差

（2）計算卡爾曼增益

（3）更新狀態(tài)和誤差協(xié)方差

（4）更新雅可比矩陣（線性化）

參數(shù)估計：對EKF 進行適當拓展，可用于實現(xiàn)參數(shù)估計[33]。引入增廣狀態(tài)x'=[xTθT]T，θ是待估計的不確定參數(shù)，并假設

式中，wθ,k-1～N(0,Qθ)是虛擬噪聲，可設為很小的數(shù)。增廣的非線性系統(tǒng)為

式中

增廣系統(tǒng)式（22）在形式上和式（14）保持一致，因此可以應用EKF 估計增廣狀態(tài)x'=[xTθT]T，實現(xiàn)原始狀態(tài)x和未知參數(shù)θ的同步估計。

3 金氰化浸出過程的動態(tài)實時優(yōu)化

3.1 算法描述

基于經濟模型預測控制EMPC和EKF參數(shù)估計（EMPC+EKF）的金氰化浸出過程動態(tài)實時優(yōu)化方法，其算法流程如下：

（1）初始化

①選擇采樣時間和窗口長度N，基于浸出過程的微分代數(shù)方程式（2）～式（8），應用直接配置法將其離散化，將EMPC 描述為如式（10）所示的NLP問題；

（2）for k=1 to ∞

①獲得測量值yk；

3.2 仿真結果

狀態(tài)估計：為測試EKF 的狀態(tài)估計性能，設置固相金品位Cs、液相金濃度Cl和液相氰離子濃度CCN，測量噪聲分別為N(0,0.12)、N(0,0.52)和N(0,52)，采樣時間1 h。此時，假設控制輸入QCN=[0.37 0.45 0.48 0]T保持不變，系統(tǒng)初態(tài)隨機分布，仿真時長200 h。圖3顯示了標稱工況下的狀態(tài)估計效果，最上層4 個子圖依次為第1～4 個浸出槽的固相金品位Cs的估計結果，中間4 個子圖為第1～4 個浸出槽的液相金濃度Cl，最下層4 個子圖為液相氰離子濃度CCN。如圖3 所示，EKF 能在大量噪聲下對真實變量進行有效估計，經歷30 步左右，各變量的估計值大致收斂至真實值。

圖3 EKF狀態(tài)變量估計效果（標稱工況）Fig.3 Estimation performances for state variables using EKF(nominal condition)

圖4 EKF狀態(tài)及參數(shù)估計效果（θ=[0.0012 1.917 3.866×10-9]T）Fig.4 Estimation performances for state variables and parameters using EKF (θ=[0.0012 1.917 3.866×10-9]T)

圖5 操作輸入QCN的運行軌跡（標稱工況）Fig.5 Trajectories of control inputs QCN(nominal condition)

參數(shù)估計：假設反應速率參數(shù)θ=[k1k2k5]T為未知參數(shù)，其真實值為θ=[0.0012 1.917 3.866×10-9]T。設k=0 時刻參數(shù)的初始估計為標稱值θ0=[0.0011 2.13 3.68×10-9]T，圖4 顯示了EKF 對狀態(tài)變量及未知參數(shù)的估計效果。如圖所示，狀態(tài)變量x（第1～3 層子圖）和未知參數(shù)θ（最下層3個子圖）都得到了準確估計，其中參數(shù)θ 的收斂速度較狀態(tài)變量x 更慢，一個可能的原因是θ自身非測量變量。

EMPC 動態(tài)優(yōu)化：首先考慮θ 已知且為標稱值時，對比直接應用EMPC（帶噪聲的測量值）和EMPC+EKF（濾波后的測量值）的優(yōu)化方法。如圖5所示，由于測量噪聲的影響，直接應用EMPC 得到的輸入優(yōu)化軌跡波動劇烈，難以收斂。EMPC+EKF 方法得到了平穩(wěn)的結果，并且4 個輸入變量在經歷若干步后均收斂到穩(wěn)態(tài)最優(yōu)值。圖6顯示了經濟成本J 的變化軌跡，雖然兩種方案均呈現(xiàn)出下降的趨勢，EMPC+EKF 給出了更理想的結果。此外，黑色虛線是直接將輸入變量設置并保持為穩(wěn)態(tài)最優(yōu)值（假設已知）的經濟指標。可以看到，EMPC+EKF 方案在初期（0～5 h）的經濟性能劣于直接將輸入變量設置為穩(wěn)態(tài)最優(yōu)值，這是由于EKF 的狀態(tài)估計仍有較大偏差，隨著狀態(tài)估計效果的提高，EMPC+EKF 的經濟性能逐漸高于穩(wěn)態(tài)優(yōu)化（15～40 h），表明了動態(tài)優(yōu)化的優(yōu)越性。隨著系統(tǒng)趨于穩(wěn)定，EMPC+EKF 最終收斂到穩(wěn)態(tài)最優(yōu)點。

圖6 成本函數(shù)J的軌跡（標稱工況）Fig.6 Trajectories of cost function J(nominal condition)

圖7 操作輸入QCN的運行軌跡（參數(shù)θ未知）Fig.7 Trajectories of control inputs QCN(unknown parameters θ)

考慮模型參數(shù)真實值為θ= [0.0013 1.917 3.866×10-9]T（假設未知），同步執(zhí)行參數(shù)估計和EMPC。如圖7所示的運行結果，即使系統(tǒng)受未知參數(shù)影響，操作輸入QCN也能在EMPC和EKF協(xié)同作用下逐漸逼近真實最優(yōu)值。此時需要約150 h 操作輸入才全部收斂，慢于單獨使用EMPC 或EKF，這是由于兩者交叉作用所致。經濟指標J 的運行軌跡如圖8 所示，可以看到J 最終趨向于真實的穩(wěn)態(tài)最優(yōu)值。注意，該最優(yōu)值在參數(shù)不確定時是未知的，而本文方法仍然能實現(xiàn)最優(yōu)操作。

圖8 成本函數(shù)J的軌跡（參數(shù)θ未知）Fig.8 Trajectories of cost function J(unknown parameters θ)

4 結 論

本文提出了一種基于經濟模型預測控制（EMPC）的金氰化浸出過程動態(tài)實時優(yōu)化方法。以經濟指標作為EMPC 的性能指標，在動態(tài)過程中求解相應的非線性規(guī)劃問題來滾動優(yōu)化操作輸入，實現(xiàn)優(yōu)化與控制一體化。為解決金氰化浸出過程中普遍存在的測量噪聲和參數(shù)可變等不確定性問題，使用了擴展卡爾曼濾波對系統(tǒng)狀態(tài)和不確定參數(shù)進行同步估計，使EMPC 模型和真實過程相吻合。仿真結果表明，本文方法一方面能通過參數(shù)估計確認出浸出過程的真實最優(yōu)點，另一方面，基于EMPC的動態(tài)優(yōu)化策略能充分提高系統(tǒng)的經濟性能，具有良好的應用前景。