單睿斌，李秋紅，何鳳林，馮海龍，管庭筠

（南京航空航天大學 能源與動力學院 江蘇省航空動力系統(tǒng)重點實驗室，南京210016）

模型預測控制（Model Predictive Control，MPC）具有顯式處理航空發(fā)動機工作中約束的能力，可以簡化航空發(fā)動機控制系統(tǒng)的結(jié)構(gòu)，在發(fā)動機控制領域獲得了廣泛的關注[1-4]。MPC通過在線優(yōu)化的手段獲得最優(yōu)控制量，因此高效的在線優(yōu)化算法可以提升MPC的實時性，有助于其在航空發(fā)動機上的應用。常用于MPC的優(yōu)化算法有二次規(guī)劃（Quadratic Programming，QP）[5]和序列二次 規(guī) 劃 （Sequential Quadratic Programming，SQP）[6]。SQP方法在每次迭代中需要計算一個QP子問題，計算量大，實時性較差；求解QP問題的方法包括內(nèi)點法（Interior Point Method，IPM）和有效集方法（Active Set Method）。有效集方法每次迭代需要判定有效集操作，且無法利用MPC問題稀疏性求解，適合控制變量和約束數(shù)量較少的情況[7]；IPM在每次迭代中需要求解 Karush-Kuhn-Tucker（KKT）系統(tǒng)，可以利用 MPC的稀疏性加速求解，但是在迭代過程中KKT系統(tǒng)會改變，每次迭代需要進行矩陣分解或者求逆操作[8]，在約束多或預測時域較長的情況下計算耗時較多。

交替方向乘子法（Alternating Direction Method of Multipliers，ADMM）于 1975年提出，在20世紀80年代獲得廣泛討論[9]，其將一個大的、具有可分結(jié)構(gòu)的凸優(yōu)化問題分解為幾個小的優(yōu)化問題，降低了求解的難度和算法的復雜性。近10年來，原本用來求解變分不等式的 ADMM在優(yōu)化計算中被廣泛采用[10-12]。

本文通過引入輔助變量將航空發(fā)動機MPC中的QP問題轉(zhuǎn)化為可分結(jié)構(gòu)的優(yōu)化問題，并基于ADMM算法對其進行求解，與IPM算法進行比較，驗證了ADMM算法的優(yōu)越性。

1 航空發(fā)動機MPC問題

航空發(fā)動機控制系統(tǒng)具有穩(wěn)態(tài)控制、加減速控制[13]、超限保護控制[14-18]、抗執(zhí)行機構(gòu)積分飽和等功能。傳統(tǒng)的發(fā)動機控制系統(tǒng)各功能模塊通過min／max選擇、切換等邏輯綜合為一個整體，控制結(jié)構(gòu)復雜。在MPC中，控制量通過求解一個約束最優(yōu)化問題來獲得，只需要一個優(yōu)化控制器，就可以實現(xiàn)上述所有控制功能，且避免了設計過程中的保守性，因而得到了廣泛的關注。

在單輸入MPC中，取二次型性能指標，在采樣時刻d之后的預測時域內(nèi)，航空發(fā)動機MPC的目標是求解在預測時域內(nèi)使性能指標J最小的燃油流量Wfb輸入序列，同時滿足發(fā)動機低壓轉(zhuǎn)子轉(zhuǎn)速Nl、壓氣機出口靜壓 P3s和低壓渦輪出口總溫T6不超過其最大值?？杀硎緸?/p>

式中：ref為高壓轉(zhuǎn)子轉(zhuǎn)速指令；Nh為高壓轉(zhuǎn)子轉(zhuǎn)速；λ為權重系數(shù)；np為預測時域的長度。

求解式（1）需要獲得各個輸出量和 Wfb之間的解析表達關系。在采樣時刻d將發(fā)動機表示為一個線性系統(tǒng)為

式中：N=[Nl，Nh]T為狀態(tài)量；ye=Nh為被控制量；yc=[Nl，P3s，T6]T為需要限制的量；A、B、C和D為系統(tǒng)的系數(shù)矩陣；Cc、Dc為限制量對應的輸出矩陣。

利用離散系統(tǒng)狀態(tài)方程的響應：

可以獲得輸出量在預測時域np內(nèi)的預測值，其預測方程為

式（1）可以改寫為一個帶有線性不等式約束的二次規(guī)劃問題：

式中：Ymax為由 Nlmax、P3smax和 T6max組成的適當維數(shù)的與對應的限制序列；Umax為由Wfbmax構(gòu)成的輸入限制序列。

通過求解QP問題式（7），可以獲得滿足約束條件下，使性能指標最小的最優(yōu)控制輸入序列

為了實現(xiàn)閉環(huán)非線性控制，將序列的第一個值Wfb（d+1）作為發(fā)動機所需燃油輸入量，然后在d+1采樣時刻，根據(jù)發(fā)動機的狀態(tài)變化重復進行線性化，構(gòu)造QP問題，求解輸入序列。這樣在每個采樣周期求解的是一個線性MPC問題，但在整個時域內(nèi)是一個非線性MPC問題。

MPC通過迭代進行優(yōu)化搜索來獲得最優(yōu)的控制序列，如果迭代過程耗時較多將限制算法的應用。ADMM算法具有對偶上升法的可分解性以及乘子法的全局收斂性，近年來得到了大量關注，由于其所需計算量小、算法結(jié)構(gòu)簡單，本文將其用于求解航空發(fā)動機MPC中的優(yōu)化問題。

2 交替方向乘子法

ADMM算法可用于求解形如的約束規(guī)劃問題。式中：x∈Rn，z∈Rm為待優(yōu)化變量；f（x）+g（z）為待優(yōu)化的目標函數(shù)；L∈Rp×n，M∈Rp×m，c∈Rp，Lx+Mz=c為問題需滿足的線性等式約束。

通過乘子法，引入對偶變量y，構(gòu)造增廣拉格朗日函數(shù)：

式中：ρ＞0為懲罰參數(shù)。

ADMM迭代過程與對偶上升法和乘子法相似，包含3部分：原始變量x的迭代更新、原始變量z的迭代更新和對偶變量y的更新過程，更新策略為[19]

可以看出，在 ADMM算法中，原始變量 x和原始變量z的迭代更新是一個交替進行的過程，每個求解過程只需求解部分變量，降低了求解規(guī)模。

3 航空發(fā)動機MPC的ADMM算法實現(xiàn)

將式（7）所描述航空發(fā)動機MPC中的QP問題簡寫為

式中：P=HTH+λI為正定對稱矩陣；

為應用ADMM算法，首先添加輔助變量，將式（11）改寫為等式約束的形式：

由于式（12）所描述的QP問題不具有可分形式，不能直接采用ADMM算法進行求解。為此將x和z視作為一個整體變量（x，z），引入輔助變量約束）=（x，z），式（12）可以寫為

式中：f2（x，z）為集合的指示函數(shù)集合的指示函數(shù)。

引入對偶變量 （w，y），根據(jù)式（10），可以得到求解式（13）的ADMM迭代過程為[20]

式中：σ為懲罰系數(shù)。

式（14）表示輔助變量的更新過程，是一個等式約束的QP問題，其一階最優(yōu)性條件（KKT條件）為

式中：v（k+1）為拉格朗日乘子。消去，可以獲得

式（18）的系數(shù)矩陣稱作KKT矩陣，該矩陣是一個列滿秩的對稱矩陣，因此方程具有唯一解。觀察式（18），KKT矩陣的結(jié)構(gòu)只和式（11）與參數(shù)ρ相關，在迭代過程中是不變的。~x（k+1）與 v（k+1）可以通過式（20）獲得：

為了加速算法收斂，根據(jù)文獻[12]，加入松弛因子 α∈[1，2]，迭代過程式（22）和式（16）變成

根據(jù)QP問題的一階最優(yōu)性條件，定義QP問題式（12）的原始殘差rprim和對偶殘差rdual：

根據(jù)2個殘差給出收斂判定準則：

一般取εrel=10-3，εabs根據(jù)需要精度選擇。

使用ADMM算法還需要一個初始值，由于MPC問題的特殊性，可以將d-1時刻的最優(yōu)解（x*（d-1），y*（d-1））作為 d時刻的算法初始值，即

因此可以獲得航空發(fā)動機MPC求解控制輸入Wfb的流程：

1）獲得發(fā)動機采樣時刻 d的狀態(tài)和狀態(tài)方程。

2）以二次型性能指標和約束，根據(jù)狀態(tài)方程建立問題式（7）。

3）將問題式（7）改寫為適合用 ADMM算法的形式（13）。

4）根據(jù)式（29），將 d-1時刻得到最優(yōu)解作為d時刻問題求解的初始值。

5）根據(jù)迭代過程式（19）、式（20）、式（23）和式（24）計算并更新各個變量。

6）根據(jù)原始殘差rprim和對偶殘差rdual判斷迭代過程是否滿足精度要求。若滿足，則終止迭代，將獲得的最優(yōu)解的第一項作為所需Wfb送入發(fā)動機，進入步驟7）；若沒有，進入步驟4），直至達到最大迭代次數(shù)。

7）進入采樣時刻d+1，重復步驟1）。

4 仿 真

4.1 控制結(jié)構(gòu)

圖1為航空發(fā)動機MPC系統(tǒng)結(jié)構(gòu)，主要包括模型預測控制器和發(fā)動機模型。MPC根據(jù)輸入的參考指令ref，使用第3節(jié)所述流程求解控制量Wfb；實時非線性模型用于獲得預測模型，并且通過模型的輸出與真實測量值之差作為反饋校正。實時非線性模型為平衡流形展開模型[21]，離線得到，其輸入為 Wfb，狀態(tài)量為 Nl和 Nh，輸出為 Nl、Nh、P3s和T6，調(diào)度參數(shù)為發(fā)動機進口溫度 T2和Nh；基于發(fā)動機進口條件和當前工作轉(zhuǎn)速，以此平衡流形展開模型獲得預測模型。使用部件級模型作為真實發(fā)動機，以ADMM算法對控制量進行迭代尋優(yōu)。

平衡流形展開模型通過實驗數(shù)據(jù)離線計算所得，其輸出與發(fā)動機傳感器輸出存在一定的誤差，使得預測模型也存在誤差，為了滿足控制精度的要求，需要對預測方程的輸出進行修正。

在 d時刻，預測模型的輸出 Nl、Nh、P3s和 T6為 ym（d），真實發(fā)動機傳感器輸出為 ys（d），二者的誤差為 e（d）=ys（d）-ym（d），則未來時刻模型的輸出可以校正為

圖1 航空發(fā)動機MPC系統(tǒng)結(jié)構(gòu)Fig.1 Aircraft engine MPC system structure

4.2 仿真結(jié)果

使用MATLAB進行仿真，MPC算法由 MATLAB實現(xiàn)，模擬真實發(fā)動機的部件級模型使用VC++開發(fā)，并且通過MEX方法在MATLAB中調(diào)用。

指令跟蹤的目的是讓發(fā)動機Nh跟隨參考指令ref的變化，MPC的性能指標J的第一部分即是指令跟蹤的誤差（見式（1））。通過在高度 H=0 km，Ma=0，在發(fā)動機慢車以上給定發(fā)動機Nh指令，使用圖1所示的控制結(jié)構(gòu)，進行Nh閉環(huán)控制。

ADMM算法最大迭代次數(shù)為500次，收斂準則為εabs=1×10-5；控制器預測時域長度 np=2，權重系數(shù)λ=700。Nh響應過程、相對跟蹤誤差、Wfb、T6和P3s的響應過程如圖2所示，數(shù)據(jù)均經(jīng)過歸一化處理。

圖2（a）為Nh在給定指令下的響應過程，從圖中可以看出，在仿真時間里，使用ADMM算法的航空發(fā)動機MPC均能使 Nh跟隨指令變化。圖2（b）為 Nh與 ref的相對誤差，從圖中可以看出，相對誤差在Nh穩(wěn)定時接近0，滿足控制所需精度，表明ADMM算法作為在線優(yōu)化算法可以滿足控制器所需精度。

圖2（c）為控制器計算所得的 Wfb輸入曲線，圖 2（d）和圖 2（e）為 T6和 P3s的響應曲線。在仿真中，T6的響應速度最快，跟隨 Wfb的輸入變化；P3s的響應速度較 Nh快，比 T6略慢一些；Nh的響應速度最慢，這是由于轉(zhuǎn)子動力學特性是所建立的部件級模型的主要動態(tài)特性，其慣性較大；由于建模中忽略了熱慣性，T6緊跟Wfb變化；P3s受轉(zhuǎn)子轉(zhuǎn)速和燃氣流量的共同作用，因此響應速度介于二者之間。觀察到在仿真過程中，其在溫度和壓力較低區(qū)域有明顯的超調(diào)。因為發(fā)動機是一個非線性系統(tǒng)，其動態(tài)特性隨轉(zhuǎn)速變化而變化，因此在仿真權重系數(shù)不變的情況下，低轉(zhuǎn)速區(qū)有超調(diào)，而高轉(zhuǎn)速沒有。

為了驗證基于ADMM算法的MPC對于發(fā)動機約束的處理能力，將T6限制適當調(diào)小，使發(fā)動機在仿真中能觸及 T6限制。圖3（a）為 T6將限制調(diào)小到0.87時 Nh階躍的響應，圖3（b）為 T6響應?？梢钥闯?，由于限制的存在，Nh不能跟蹤給定的指令，而T6達到限制值，并且沒有明顯的超調(diào)。仿真表明了基于ADMM算法的MPC良好的約束處理能力，采用MPC可以取消超限保護控制回路。

為了驗證ADMM算法在實時性方面的優(yōu)勢，使用文獻[22]中所述IPM算法作為對比算法，在H=0 km，Ma=0，慢車以上轉(zhuǎn)速，給定不同的階躍幅值ΔNh和預測時域np，進行6 s的階躍響應仿真，記錄使用2種算法平均單步仿真耗時并進行對比。

圖2 基于ADMM算法的航空發(fā)動機MPC響應Fig.2 Response of aircraft engine MPC based on ADMM

表1中給出了10個仿真過程中轉(zhuǎn)速階躍幅值ΔNh及預測時域設置。表2是在這些設置情況下使用IPM完成一次控制序列優(yōu)化計算所需平均時間和使用ADMM所需時間。

圖3 T6觸及限制下MPC仿真Fig.3 MPC simulation with T6 hitting limit

表1 仿真參數(shù)Tab1e 1 Simu1ation parameters

從表2中可以看出，無論ΔNh是較大還是較小的動態(tài)響應，使用ADMM的仿真耗時均低于使用IPM的仿真。這是由于ADMM算法在迭代過程中，在給定罰參數(shù)不變的情況下，KKT系數(shù)矩陣不變，因此在每次迭代過程中不必重復計算KKT矩陣的逆矩陣或分解，減少了計算量；而IPM算法在每次迭代過程中需要調(diào)整KKT矩陣，在每次迭代求解中都需要進行KKT矩陣求逆或者矩陣分解操作。假設求解需要M次迭代，矩陣分解耗時為tf，KKT系統(tǒng)求解回代耗時 tb，ADMM算法求解需要時間為 tf+Mtb，而 IPM算法耗時為M（tf+tb），M為求解迭代次數(shù)。

預測時域np可以表示MPC滾動優(yōu)化問題的規(guī)模大小。np越大，則待求解的燃油序列維數(shù)越大，需求解的QP問題中的不等式約束矩陣和系數(shù)矩陣的維數(shù)也越大，增加了求解的計算量。

預測時域的增大可以改善航空發(fā)動機MPC的性能。圖4為在H=0km，Ma=0，慢車以上轉(zhuǎn)速，給定階躍幅值 ΔNh=0.18，改變 MPC的預測時域np，得到的完成一次控制序列優(yōu)化所需平均時間t與預測時域 np的關系。隨著 np的增加，QP的規(guī)模增加，構(gòu)成的KKT系數(shù)矩陣的維數(shù)也增加，求解 KKT方程的耗時增加。仿真結(jié)果表明，使用IPM算法的仿真耗時增加隨np增加更為劇烈，而使用ADMM算法的仿真耗時增加幅度較小，證實了本文分析的合理性，表明了ADMM算法在航空發(fā)動機MPC中有比IPM更好的實時性。

表2 兩種算法完成一次控制序列優(yōu)化所需時間對比Tab1e 2 Time consumption comparison of two methods in finishing one-time contro1 sequence optimization

圖4 控制序列優(yōu)化所需平均時間與預測時域的關系Fig.4 Average time consumption for control sequence optimization vs predictive horizon

5 結(jié) 論

ADMM算法在航空發(fā)動機MPC的滾動優(yōu)化中具有良好的應用前景。

1）在指令跟蹤控制中，Nh與指令的相對跟蹤誤差均為0；表明ADMM可以滿足發(fā)動機指令跟蹤和約束管理的精度要求。

2）ADMM算法每次迭代過程只需計算一次矩陣分解，與IPM相比，計算量大大降低，實時性得到了很大的提高。

3）ADMM算法相對IPM有更好的計算耗時穩(wěn)定性。在10組仿真測試中，預測時域從2增加到14，使用IPM算法的單次仿真耗時從10.67 ms增加到224.1 ms，擴大了20多倍；而ADMM算法從3.67 ms增加到15.53 ms，擴大不足5倍。