基于經(jīng)濟(jì)模型預(yù)測(cè)控制的金氰化浸出過程動(dòng)態(tài)實(shí)時(shí)優(yōu)化

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

（1 寧波財(cái)經(jīng)學(xué)院機(jī)械與電氣工程學(xué)院，浙江寧波315175； 2 浙江大學(xué)寧波理工學(xué)院，浙江寧波315100； 3 江南大學(xué)輕工過程先進(jìn)控制教育部重點(diǎn)實(shí)驗(yàn)室，江蘇無錫214122； 4 浙江大學(xué)控制科學(xué)與工程學(xué)院，浙江杭州310027）

引 言

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

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

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

1 金氰化浸出過程

1.1 過程概述

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

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

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

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

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

反應(yīng)速率表示為如下形式

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

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

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

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

綜上，金氰化浸出過程的操作目標(biāo)描述為經(jīng)濟(jì)指標(biāo)最大化問題（成本最?。?，即求解如下最優(yōu)化問題

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

2 經(jīng)濟(jì)模型預(yù)測(cè)控制

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

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

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

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

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

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

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

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

對(duì)控制輸入u，在每個(gè)時(shí)間間隔內(nèi)視為常數(shù)，得到擴(kuò)展向量

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

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

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

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

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

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

（2）計(jì)算卡爾曼增益

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

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

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

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

式中

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

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

3.1 算法描述

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

（1）初始化

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

（2）for k=1 to ∞

①獲得測(cè)量值yk；

3.2 仿真結(jié)果

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

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

圖4 EKF狀態(tài)及參數(shù)估計(jì)效果（θ=[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的運(yùn)行軌跡（標(biāo)稱工況）Fig.5 Trajectories of control inputs QCN(nominal condition)

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

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

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

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

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

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

4 結(jié) 論

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