模型預測控制器權重參數整定非線性規劃法

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

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

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

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

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

1 模型預測控制算法

1.1 被控對象模型

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

(1)

式中：xk∈Rn、uk∈Rm、yk∈Rq和dk∈Rl分別為第k個時間采樣點的狀態、輸入、輸出和可測擾動變量維度；A0∈Rn×n、B0∈Rn×m、C0∈Rq×n和Bd0∈Rn×l分別為狀態、輸入、輸出和可測擾動變量的轉化矩陣，其中n、m、q和l均為向量或矩陣的行數或列數；w和v分別為狀態和可測輸出的不確定性噪音。此處假設方程(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集成到狀態維度Δxk中，得到新的狀態空間模型：

(3)

1.3 狀態觀測器

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

(4)

(5)

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

1.4 模型預測控制器

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

(6)

2 模型預測控制器權重參數整定

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

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

(7)

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

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

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

將方程(6)所示的MPC子優化問題轉化為標準二次規劃形式：

(8)

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

(9)

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

(10)

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

subject to:

Qmin≤Q≤Qmax，

Rmin≤R≤Rmax，

(11)

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

3 仿真案例

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

3.1 單輸入單輸出模型

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

(12)

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

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

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

表1 SISO系統整定參數

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

圖1 SISO系統設定點跟蹤MPC動態響應Fig. 1 Dynamic responses of MPC controller for SISO system under setpoint changes

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

(13)

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

表2 SISO系統輸出動態響應IAE值

3.2 連續攪拌反應器

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

圖2 連續攪拌反應器示意圖Fig .2 Schematic diagram of CSTR

(14)

(15)

式中：V為反應器體積；T為反應器內液體溫度；k0為反應速率常數；E為反應活化能；Cp為反應器內液體比熱容；ρL為反應器內液體密度；ΔH為反應熱；Rg為氣體常數。

反應物A的溫度T0設為可測擾動量，因此，反應器的操作目標為通過調節反應物摩爾濃度CA0和熱流速率QR使反應器內反應物濃度CA和溫度T維持在設定點。故反應器操作過程中，其控制變量、操作變量及其可測擾動可表示為方程(16)～(18)，其中下標s表示每個變量的穩態值(見表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模型參數

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

表4 變量范圍

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

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

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

表5 CSTR系統整定參數

圖3 CSTR系統T0升高5 K時MPC動態響應Fig. 3 Dynamic responses of MPC controller for CSTR system under 5 K increase of T0

表6 原料溫度升高5 K時CSTR系統輸出響應IAE值

圖4 CSTR系統CA設定點跟蹤MPC動態響應Fig. 4 Dynamic responses of MPC controller for CSTR system under the changes of CA setpoint

表7 CA設定點跟蹤CSTR系統輸出動態響應IAE值

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

圖5 MPC權重參數整定多目標優化Pareto解Fig. 5 Pareto solution of multi-objective optimization for MPC weight parameter tuning

4 結 論

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

1)基于單層結構的非線性規劃整定方法能夠快速實現MPC權重參數的整定，極大地降低權重參數整定優化問題的求解時間。對于SISO系統和CSTR系統，非線性規劃權重參數整定優化問題的求解時間僅為5～10 s和10～90 s，而多目標優化的整定方法則分別為1.0 h和1.5 h。

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