馮澤民，李 喬，譚陸西,2c，董立春,2c

(1. 重慶科技學院 安全工程學院，重慶 401331；2. 重慶大學 a.化學化工學院；b.化工過程強化與反應(yīng)國家地方聯(lián)合工程實驗室；c.低品位能源利用技術(shù)及系統(tǒng)教育部重點實驗室，重慶 400044)

模型預(yù)測控制(model predictive control，MPC)由于能夠很好地處理多變量約束控制問題，自20世紀70年代以來已被廣泛應(yīng)用于石油化工生產(chǎn)過程中[1-2]。MPC是一種基于模型的控制算法，模型是否能夠準確預(yù)測過程的動態(tài)特征是MPC實現(xiàn)良好控制性能的內(nèi)在關(guān)鍵因素，MPC控制器相關(guān)參數(shù)的整定是影響控制性能的外在主要因素。

MPC中多變量之間的耦合，使采樣時間、預(yù)測時域、控制時域、輸入和輸出變量的權(quán)重參數(shù)整定極具挑戰(zhàn)性。與過程動態(tài)特征相關(guān)的采樣時間、預(yù)測時域和控制時域一般可通過分析過程動態(tài)響應(yīng)特征得到，而輸入和輸出變量的權(quán)重參數(shù)之間由于存在相互耦合，很難通過分析的方法獲得最優(yōu)值，因此，如何整定權(quán)重參數(shù)對MPC取得良好控制性能具有重要意義。目前主流的MPC權(quán)重參數(shù)整定方法可分為兩類[3]：1)將過程模型等價為一階、二階模型，應(yīng)用動態(tài)響應(yīng)性能分析的方法獲得權(quán)重參數(shù)[4-7]；2)將權(quán)重參數(shù)整定轉(zhuǎn)化為優(yōu)化問題，以優(yōu)化控制性能為目標而獲得最優(yōu)權(quán)重參數(shù)[8-12]。前者較適于單輸入單輸出MPC的權(quán)重參數(shù)整定，而后者由于更易于實施，且適用于多變量系統(tǒng)，不需要操作者具備良好的控制理論知識而備受關(guān)注。Yamashita等[12-13]以最小化每個輸出變量的參考軌跡跟蹤積分方差為優(yōu)化目標，應(yīng)用多目標優(yōu)化的方法獲得輸出變量和輸入變量變化速率的最優(yōu)權(quán)重值。Giraldo等[14]提出了基于雙層優(yōu)化的廣義預(yù)測控制器權(quán)重參數(shù)整定方法。Gutiérrez-Urquídez等[15]將多目標優(yōu)化應(yīng)用于拉蓋爾多項式參數(shù)化的MPC權(quán)重參數(shù)的整定。Santamaría等[16]及Vallerio等[17]則將多目標優(yōu)化應(yīng)用于非線性MPC權(quán)重參數(shù)的整定。

多目標優(yōu)化權(quán)重參數(shù)整定方法通常為雙層結(jié)構(gòu)，其中內(nèi)層為整定時域內(nèi)每個采樣時間點的MPC子優(yōu)化問題的求解，而外層為整定時域內(nèi)以權(quán)重參數(shù)為決策變量的控制性能目標函數(shù)的優(yōu)化。這種雙層結(jié)構(gòu)的多目標優(yōu)化問題的求解不僅存在內(nèi)外層優(yōu)化問題之間的迭代，而且需要大量的隨機樣本點來評估目標函數(shù)以獲得最優(yōu)Pareto解，計算結(jié)構(gòu)復雜，時間成本較高。Tran等[18]將無約束MPC優(yōu)化問題的解析解引入MPC權(quán)重參數(shù)整定優(yōu)化問題中，提出了基于單層非線性規(guī)劃的MPC權(quán)重參數(shù)整定方法，可快速獲得最優(yōu)權(quán)重參數(shù)，減小計算時間成本，然而該方法卻不適用于有約束MPC權(quán)重參數(shù)的整定。為此，筆者提出了一種針對有約束MPC權(quán)重參數(shù)整定的非線性規(guī)劃方法，該方法將MPC子優(yōu)化問題轉(zhuǎn)化為MPC權(quán)重參數(shù)整定優(yōu)化問題的一階最優(yōu)KKT(Karush-Kuhn-Tucker)條件，可快速獲得MPC的最優(yōu)權(quán)重參數(shù)，與其他多目標優(yōu)化方法相比可顯著降低計算耗時成本。

1 模型預(yù)測控制算法

1.1 被控對象模型

線性模型由于其魯棒性能較強、計算負荷較小，被廣泛應(yīng)用于過程工業(yè)中MPC的實施，本研究的對象是方程(1)所示的離散線性時不變狀態(tài)空間模型。

(1)

式中：xk∈Rn、uk∈Rm、yk∈Rq和dk∈Rl分別為第k個時間采樣點的狀態(tài)、輸入、輸出和可測擾動變量維度；A0∈Rn×n、B0∈Rn×m、C0∈Rq×n和Bd0∈Rn×l分別為狀態(tài)、輸入、輸出和可測擾動變量的轉(zhuǎn)化矩陣，其中n、m、q和l均為向量或矩陣的行數(shù)或列數(shù)；w和v分別為狀態(tài)和可測輸出的不確定性噪音。此處假設(shè)方程(1)的矩陣對(A0，B0)和(A0，C0)分別可控和可觀測，且w和v服從高斯分布。

1.2 增廣模型

在方程(1)中引入差分算子后可得：

(2)

式中：Δxk+1=xk+1-xk，Δyk+1=yk+1-yk，Δuk=uk-uk-1，Δdk=dk-dk-1。

將輸出維度yk集成到狀態(tài)維度Δxk中，得到新的狀態(tài)空間模型：

(3)

1.3 狀態(tài)觀測器

方程(1)所示的被控對象模型通常是通過系統(tǒng)辨識或非線性機理模型線性化得到，其狀態(tài)變量通常不可測，故需在每個采樣時間點以當前測得的輸出值和工藝擾動量估計當前的實際狀態(tài)。以方程(4)所示的狀態(tài)觀測器估計當前的狀態(tài)：

(4)

(5)

式中：Qw和Rv分別為方程(3)的狀態(tài)與輸出維度的噪音方差矩陣，本研究中分別取Rv=Iq×q和Qw=CTRvC。

1.4 模型預(yù)測控制器

基于以上方程，MPC控制器的優(yōu)化問題可描述為：

(6)

2 模型預(yù)測控制器權(quán)重參數(shù)整定

2.1 雙層多目標優(yōu)化整定方法(MOO)

在MPC控制器的參數(shù)整定過程中，本研究中僅考慮了MPC控制器輸出變量權(quán)重參數(shù)Q和輸入變量變化速率權(quán)重參數(shù)R的最優(yōu)選擇，預(yù)測時域、控制時域以及采樣時間均通過分析過程的動力學特征得到。Q越大則MPC控制器能夠更快速地跟蹤參考軌跡設(shè)定點；R越大則MPC動態(tài)響應(yīng)越趨于平滑，但不能快速地跟蹤參考軌跡設(shè)定點。為此，將整定時域內(nèi)的MPC控制器跟蹤參考軌跡的平均積分方差(Ψ1)和輸入變量變化速率的二次方的平均值(Ψ2)作為優(yōu)化目標，構(gòu)成如下所示的權(quán)重參數(shù)整定多目標優(yōu)化問題[19-20]：

(7)

方程(7)所示的多目標優(yōu)化問題應(yīng)用Shama等[21]開發(fā)的基于MS Excel的多目標優(yōu)化求解器(EMOO)求解，該求解器采用NSGA-II算法[22]，能夠很好地處理含有約束和整型變量的多目標優(yōu)化問題，已被廣泛應(yīng)用于換熱器設(shè)計[23]及化工過程優(yōu)化[24-26]中。多目標優(yōu)化得到的Pareto解為一系列最優(yōu)解的集合，因此需進一步選擇一組最佳的解作為MPC最優(yōu)權(quán)重參數(shù)。Wang等[27]比較了典型的10種用于從Pareto解中選取最優(yōu)解的方法，表明灰色關(guān)聯(lián)法(gray relational analysis，GRA)具有較好的選擇性能，且該方法可直接選出最優(yōu)解，而不需要進一步設(shè)定每個目標的權(quán)重值。因此，用GRA方法從多目標優(yōu)化所得的Pareto解中選取最優(yōu)解。

2.2 單層非線性規(guī)劃整定方法(NLP)

方程(7)所示的多目標優(yōu)化問題在求解過程中不僅存在內(nèi)層MPC子優(yōu)化問題和外層權(quán)重參數(shù)整定優(yōu)化問題之間的迭代，而且需要大量的子代數(shù)和種群數(shù)量作為評估目標方程的隨機樣本點以求得最優(yōu)Pareto解，因此MOO權(quán)重參數(shù)整定方法通常計算復雜，耗時較長，不能快速求解權(quán)重參數(shù)整定優(yōu)化問題。然而在方程(7)中，有約束的MPC子優(yōu)化問題可轉(zhuǎn)化為外層MPC權(quán)重參數(shù)整定優(yōu)化問題的一階KKT最優(yōu)條件，從而將方程(7)所示的雙層結(jié)構(gòu)MPC權(quán)重參數(shù)整定優(yōu)化問題轉(zhuǎn)化為單層非線性規(guī)劃問題，運用基于梯度的非線性規(guī)劃方法進行求解，可極大地縮短計算時間。

將方程(6)所示的MPC子優(yōu)化問題轉(zhuǎn)化為標準二次規(guī)劃形式：

(8)

式中：Δuk=[Δuk,Δuk+1,…,Δuk+Nc-1]T為tk時刻MPC子優(yōu)化問題的控制動作變化速率維度；H和g分別為方程(6)所示的MPC子優(yōu)化問題的海森矩陣和雅克比矩陣；Cons為常數(shù)項；Ab和b分別為線性不等式約束的變換矩陣及其值。故方程(6)所示的MPC子優(yōu)化問題的拉格朗日方程可表示為：

(9)

式中：λi和si分別為第i個線性約束的拉格朗日乘子和松弛因子；nc為MPC子優(yōu)化問題線性約束不等式的個數(shù)；Abi和bi分別為第i個線性約束表達式的系數(shù)和約束值。故方程(8)的最優(yōu)KKT條件可表示為：

(10)

方程(7)所示的優(yōu)化目標Ψ1和Ψ2可通過權(quán)重化集成為單目標，并將方程(10)作為MPC權(quán)重參數(shù)整定優(yōu)化問題的約束條件，可得方程(11)所示的MPC權(quán)重參數(shù)整定單層非線性規(guī)劃法。

subject to:

Qmin≤Q≤Qmax，

Rmin≤R≤Rmax，

(11)

式中：γ∈[0,1]為目標函數(shù)的權(quán)重因子，其值越大，優(yōu)化所得的權(quán)重參數(shù)使得MPC能夠更快地跟蹤參考軌跡設(shè)定值；Δu0,j為控制動作序列Δuj中的第一個控制動作，由轉(zhuǎn)換矩陣M得到。

3 仿真案例

在Matlab 2019a中編譯MPC控制算法，并通過Matlab中的Spreadlink工具箱實現(xiàn)Matlab與在MS Excel中建立的EMOO多目標優(yōu)化求解器之間的數(shù)據(jù)傳輸。在Matlab 2019a中應(yīng)用Casadi算法網(wǎng)執(zhí)行非線性規(guī)劃權(quán)重參數(shù)整定算法，并用ipopt求解器求解方程(11)所示的非線性規(guī)劃問題。仿真案例的測試在Simulink 2019a中實施。本節(jié)分別以一個單輸入單輸出(SISO)二階時間延時傳遞函數(shù)和連續(xù)攪拌反應(yīng)器(CSTR)模型闡述所提出的MPC權(quán)重參數(shù)整定算法的有效性。

3.1 單輸入單輸出模型

方程(12)給出了一個SISO二階時間延時傳遞函數(shù)模型[12]。

(12)

式中：G為被控對象模型；s為復變量。

該傳遞函數(shù)為無約束系統(tǒng)，首先以采樣時間ts=0.5 min離散化為方程(1)所示的狀態(tài)空間模型，并取Np和Nc分別為60和3個采樣點，故對于MPC控制器還需確定權(quán)重參數(shù)Q和R。為了簡化計算，對于SISO系統(tǒng)，取輸出權(quán)重Q=1，則優(yōu)化過程中只需確定相應(yīng)的R即可。故可通過優(yōu)化整定時域內(nèi)MPC控制器對參考軌跡的跟蹤性能而獲得最優(yōu)權(quán)重參數(shù)值，選取整定時域為100采樣點，起始點u0=0，y0=0，ysp=0；當t=2.5 min時，設(shè)定ysp=2；當t=25.0 min時，設(shè)定ysp=-2。

在MOO整定方法中，將遺傳算法的代數(shù)和種群數(shù)分別設(shè)定為100和50，而在NLP整定方法中將權(quán)重因子γ分別設(shè)定為1.00、0.75和0.50以考察γ對MPC整定結(jié)果的影響。表1為優(yōu)化所得的R值。由于應(yīng)用GRA方法從Pareto解中選取的最優(yōu)解更傾向于使積分跟蹤均方差Ψ1最小化，因此MOO整定方法所得的最優(yōu)R值與γ=1.00時NLP整定方法所得的最優(yōu)解相近。

表1 SISO系統(tǒng)整定參數(shù)

由圖1描述的以不同整定方法所得的最優(yōu)權(quán)重R作為MPC權(quán)重參數(shù)時SISO系統(tǒng)的動態(tài)響應(yīng)可以看出，隨著權(quán)重因子γ逐漸增大，MPC趨于更快速地跟蹤參考軌跡設(shè)定點(圖1(a))，但控制動作的變化幅度也同時增大；反之，MPC輸出響應(yīng)更趨于平緩。因此，在NLP整定方法中，可通過調(diào)節(jié)γ的大小優(yōu)化MPC動態(tài)響應(yīng)特性。當γ=1.00時，NLP和MOO兩種整定方法所得的MPC動態(tài)響應(yīng)性能相近。

圖1 SISO系統(tǒng)設(shè)定點跟蹤MPC動態(tài)響應(yīng)Fig. 1 Dynamic responses of MPC controller for SISO system under setpoint changes

不同整定方法所得的MPC的動態(tài)響應(yīng)性能由方程(13)所示的積分絕對偏差(integral absolute error，IAE)予以定量評估。

(13)

表2給出了應(yīng)用不同整定方法所得的最優(yōu)權(quán)重R時MPC輸出動態(tài)響應(yīng)的IAE值?？梢钥闯?，應(yīng)用NLP整定方法時隨著γ值減小，MPC輸出響應(yīng)的IAE逐漸變大；而當γ=1.00時，應(yīng)用NLP和MOO整定方法所得的MPC的輸出動態(tài)響應(yīng)的IAE值非常接近。NLP整定方法的計算用時僅需5～10 s，而MOO整定方法則需約1 h，由此可見，NLP整定方法能夠快速地整定MPC權(quán)重參數(shù)，且整定性能與MOO整定方法相近。

表2 SISO系統(tǒng)輸出動態(tài)響應(yīng)IAE值

3.2 連續(xù)攪拌反應(yīng)器

一連續(xù)攪拌反應(yīng)器如圖2所示，反應(yīng)物A以體積流量F、摩爾濃度CA0、溫度T0進入反應(yīng)器內(nèi)，通過液相二階可逆反應(yīng)轉(zhuǎn)化為產(chǎn)物B，反應(yīng)產(chǎn)生的熱量經(jīng)夾套冷劑以熱流速率QR移走。反應(yīng)器可描述為方程(14)和(15)所示的數(shù)學模型，其參數(shù)值見表3[28]。

圖2 連續(xù)攪拌反應(yīng)器示意圖Fig .2 Schematic diagram of CSTR

(14)

(15)

式中：V為反應(yīng)器體積；T為反應(yīng)器內(nèi)液體溫度；k0為反應(yīng)速率常數(shù)；E為反應(yīng)活化能；Cp為反應(yīng)器內(nèi)液體比熱容；ρL為反應(yīng)器內(nèi)液體密度；ΔH為反應(yīng)熱；Rg為氣體常數(shù)。

反應(yīng)物A的溫度T0設(shè)為可測擾動量，因此，反應(yīng)器的操作目標為通過調(diào)節(jié)反應(yīng)物摩爾濃度CA0和熱流速率QR使反應(yīng)器內(nèi)反應(yīng)物濃度CA和溫度T維持在設(shè)定點。故反應(yīng)器操作過程中，其控制變量、操作變量及其可測擾動可表示為方程(16)～(18)，其中下標s表示每個變量的穩(wěn)態(tài)值(見表3)。

y=[y1,y2]T=[CA-CAs,T-Ts]T；

(16)

u=[u1,u2]T=[CA0-CA0s,QR-QRs]T；

(17)

d=[T0-T0s]。

(18)

表3 CSTR模型參數(shù)

操作變量、控制變量和可測擾動的范圍見表4。

表4 變量范圍

方程(14)和(15)在穩(wěn)態(tài)操作點經(jīng)線性化并以采樣時間為0.02 h離散化后可得方程(1)所示的線性動態(tài)模型，其參數(shù)分別為：

為了方便權(quán)重參數(shù)整定優(yōu)化問題的求解，模型參數(shù)需根據(jù)變量各自的范圍區(qū)間轉(zhuǎn)化為無維度值，并取Np和Nc分別為10和3個采樣點，整定時域選取100個采樣點。在MOO整定方法中，遺傳算法的代數(shù)和種群數(shù)分別設(shè)定為100和50，NLP整定方法中γ分別設(shè)定為1.00、0.75和0.50，以考察γ對于MPC整定結(jié)果的影響。該系統(tǒng)為抗擾動干擾的MPC權(quán)重參數(shù)整定，因此，在整定過程中設(shè)定，當t=0.02 h時，可測擾動原料溫度升高5 K，即方程(1)中d=5(轉(zhuǎn)化后的無維度值)。

表5給出了優(yōu)化所得的權(quán)重參數(shù)Q和R的值。圖3描述了原料溫度T0升高5 K時，MPC控制器輸入與輸出變量動態(tài)響應(yīng)?？梢钥闯?，NLP和MOO兩種整定方法所得的MPC控制器都具有很好的抗可測溫度擾動的能力，能夠快速平滑地回到設(shè)定點。表6給出了原料溫度T0升高5 K后CA和T的動態(tài)響應(yīng)IAE值。可以看出，雖然MOO整定方法所得的MPC控制器CA的IAE值較NLP整定方法所得的IAE值小，但MOO整定方法所得的MPC控制器T的IAE值較γ=1.00時NLP整定方法所得的IAE值大。此外，表6中γ=1.00時NLP整定方法所得的MPC輸出響應(yīng)IAE值總和遠較MOO整定方法所得的IAE值總和小。因此，NLP整定方法所得的MPC控制性能更好，或與MOO整定方法所得的MPC控制器性能相近。

表5 CSTR系統(tǒng)整定參數(shù)

圖3 CSTR系統(tǒng)T0升高5 K時MPC動態(tài)響應(yīng)Fig. 3 Dynamic responses of MPC controller for CSTR system under 5 K increase of T0

表6 原料溫度升高5 K時CSTR系統(tǒng)輸出響應(yīng)IAE值

圖4 CSTR系統(tǒng)CA設(shè)定點跟蹤MPC動態(tài)響應(yīng)Fig. 4 Dynamic responses of MPC controller for CSTR system under the changes of CA setpoint

表7 CA設(shè)定點跟蹤CSTR系統(tǒng)輸出動態(tài)響應(yīng)IAE值

綜上所述， NLP整定方法所得的MPC控制性能優(yōu)于或接近于MOO整定方法所得的MPC，且在整定過程中，NLP整定方法求解優(yōu)化問題所需的計算時間僅為10～90 s，而MOO整定方法所需的計算時間則約為1.5 h。由圖5所示的MOO整定方法所得的Pareto解可以看到GRA方法選取的最優(yōu)點(藍色五角星)趨向于最小化參考軌跡跟蹤目標Ψ1，這與方程(11)中γ=1.00時的優(yōu)化意義一致，從另一層面表明了兩種整定方法的一致性。

圖5 MPC權(quán)重參數(shù)整定多目標優(yōu)化Pareto解Fig. 5 Pareto solution of multi-objective optimization for MPC weight parameter tuning

4 結(jié) 論

權(quán)重參數(shù)的整定是MPC控制器取得良好控制性能的重要因素。針對基于雙層結(jié)構(gòu)多目標優(yōu)化的MPC權(quán)重參數(shù)整定方法存在求解過程較慢、耗時較長的問題，建立了一種非線性規(guī)劃的權(quán)重參數(shù)整定方法。該方法將MPC權(quán)重參數(shù)整定中每個時間采樣點的MPC子優(yōu)化問題等價為外層MPC權(quán)重參數(shù)整定優(yōu)化問題的最優(yōu)KKT(Karush-Kuhn-Tucker)條件，將MPC權(quán)重參數(shù)整定的雙層多目標優(yōu)化問題轉(zhuǎn)化為單層非線性規(guī)劃問題，并應(yīng)用一個SISO時間延時傳遞函數(shù)和CSTR反應(yīng)器仿真案例評估該算法的有效性。結(jié)果表明：

1)基于單層結(jié)構(gòu)的非線性規(guī)劃整定方法能夠快速實現(xiàn)MPC權(quán)重參數(shù)的整定，極大地降低權(quán)重參數(shù)整定優(yōu)化問題的求解時間。對于SISO系統(tǒng)和CSTR系統(tǒng)，非線性規(guī)劃權(quán)重參數(shù)整定優(yōu)化問題的求解時間僅為5～10 s和10～90 s，而多目標優(yōu)化的整定方法則分別為1.0 h和1.5 h。

2)基于單層非線性規(guī)劃整定方法的MPC控制器控制性能好于或接近基于多目標優(yōu)化整定方法的MPC控制器控制性能。