基于ADMM算法的航空發動機模型預測控制

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

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

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

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

本文通過引入輔助變量將航空發動機MPC中的QP問題轉化為可分結構的優化問題，并基于ADMM算法對其進行求解，與IPM算法進行比較，驗證了ADMM算法的優越性。

1 航空發動機MPC問題

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

在單輸入MPC中，取二次型性能指標，在采樣時刻d之后的預測時域內，航空發動機MPC的目標是求解在預測時域內使性能指標J最小的燃油流量Wfb輸入序列，同時滿足發動機低壓轉子轉速Nl、壓氣機出口靜壓 P3s和低壓渦輪出口總溫T6不超過其最大值。可表示為

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

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

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

利用離散系統狀態方程的響應：

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

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

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

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

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

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

2 交替方向乘子法

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

通過乘子法，引入對偶變量y，構造增廣拉格朗日函數：

式中：ρ＞0為懲罰參數。

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

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

3 航空發動機MPC的ADMM算法實現

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

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

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

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

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

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

式中：σ為懲罰系數。

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

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

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

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

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

根據2個殘差給出收斂判定準則：

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

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

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

1）獲得發動機采樣時刻 d的狀態和狀態方程。

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

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

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

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

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

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

4 仿 真

4.1 控制結構

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

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

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

圖1 航空發動機MPC系統結構Fig.1 Aircraft engine MPC system structure

4.2 仿真結果

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

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

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

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

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

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

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

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

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

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

表1 仿真參數Tab1e 1 Simu1ation parameters

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

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

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

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

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

5 結 論

ADMM算法在航空發動機MPC的滾動優化中具有良好的應用前景。

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

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

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