基于ADMM算法的航空發(fā)動機(jī)模型預(yù)測控制

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

（南京航空航天大學(xué) 能源與動力學(xué)院 江蘇省航空動力系統(tǒng)重點(diǎn)實(shí)驗(yàn)室，南京210016）

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2 交替方向乘子法

ADMM算法可用于求解形如的約束規(guī)劃問題。式中：x∈Rn，z∈Rm為待優(yōu)化變量；f（x）+g（z）為待優(yōu)化的目標(biāo)函數(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的迭代更新是一個(gè)交替進(jìn)行的過程，每個(gè)求解過程只需求解部分變量，降低了求解規(guī)模。

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

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

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

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

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

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

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

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

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

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

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

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

根據(jù)QP問題的一階最優(yōu)性條件，定義QP問題式（12）的原始?xì)埐顁prim和對偶?xì)埐顁dual：

根據(jù)2個(gè)殘差給出收斂判定準(zhǔn)則：

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

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

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

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

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

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

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

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

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

7）進(jìn)入采樣時(shí)刻d+1，重復(fù)步驟1）。

4 仿 真

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

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

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

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

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

4.2 仿真結(jié)果

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

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

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

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

圖2（c）為控制器計(jì)算所得的 Wfb輸入曲線，圖 2（d）和圖 2（e）為 T6和 P3s的響應(yīng)曲線。在仿真中，T6的響應(yīng)速度最快，跟隨 Wfb的輸入變化；P3s的響應(yīng)速度較 Nh快，比 T6略慢一些；Nh的響應(yīng)速度最慢，這是由于轉(zhuǎn)子動力學(xué)特性是所建立的部件級模型的主要?jiǎng)討B(tài)特性，其慣性較大；由于建模中忽略了熱慣性，T6緊跟Wfb變化；P3s受轉(zhuǎn)子轉(zhuǎn)速和燃?xì)饬髁康墓餐饔茫虼隧憫?yīng)速度介于二者之間。觀察到在仿真過程中，其在溫度和壓力較低區(qū)域有明顯的超調(diào)。因?yàn)榘l(fā)動機(jī)是一個(gè)非線性系統(tǒng)，其動態(tài)特性隨轉(zhuǎn)速變化而變化，因此在仿真權(quán)重系數(shù)不變的情況下，低轉(zhuǎn)速區(qū)有超調(diào)，而高轉(zhuǎn)速沒有。

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

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

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

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

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

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

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

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

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

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

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

5 結(jié) 論

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

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

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

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