渦扇發動機的一類大包線系統建模

凌彥聰，周文祥，朱平芳，曾建平

（1.廈門大學航空航天學院，廈門 361005；2.南京航空航天大學能源與動力學院，南京 210016）

非線性和時變是渦扇發動機系統中普遍存在的現象，高精度的非線性時變數學模型是渦扇發動機特性分析和綜合的基礎。早期渦扇發動機建模方法，大都采用線性化模型，未能真實地反映渦扇發動機系統時變和非線性的動態信息。針對這一問題，吳斌等［1］建立線性變參數（Linear parameter?varying，LPV）模型描述渦扇發動機大范圍變化的動態特性，LPV 系統的系統參數可以表示為參數的函數，是基于狀態的增量形式，需要與穩態基線模 型 結 合 使 用。Fu 等［2?3］提 出 一 類 能 夠 描 述 對 象非線性和時變特性的非線性變參數（Nonlinear pa?rameter?varying，NPV）模型。該模型將非線性特性表征為依賴狀態的多項式函數，其狀態是基于對象的物理量，但由于其不需要基線模型則會帶來誤差的積累。因此，在擬合成NPV 模型之前如何高精度地對各個穩態點以及各包線的建模就顯得至關重要。本文采用系統辨識的思想研究渦扇發動機的NPV 系統建模方法。

關于渦扇發動機模型辨識可以分為兩類：一是對試驗數據的辨識修正提高部件法模型的精度。二是選擇輸入輸出數據以及模型類進行辨識建模。美國路易斯研究中心較早地研究系統辨識理論在發動機控制系統設計中的應用，Banazadeh等［4］提出一種掃頻方法，通過發動機頻域響應得到其 頻 率 響 應 函 數。Kaufman 等［5］提 出GE7 燃 氣 輪機建立狀態變量模型的分步辨識方法。Breikin等［6］通過對發動機穩態點的輸入輸出數據，利用最小二乘法辨識狀態空間系統矩陣參數得到發動機多變量的狀態變量模型。Norton 等［7］研究不同的輸入燃油信號對發動機時間序列模型辨識的影響，提出計算發動機時間序列模型參數的方法。As?gari 等［8］與郭政波等［9］提出NARMAX 模型利用正交分解法對發動機起動過程進行系統辨識，建立發動機非線性模型。隨著航空發動機多變量控制規律研究的深入，國內也相應地展開渦扇發動機模型辨識領域的研究，姜銳［10］、吳斌［11］等學者采用正弦燃油信號通過頻率響應的方法得到渦扇發動機狀態依賴的傳遞函數。鄭斐華［12］提出了采用偽隨機信號基于系統辨識的航空發動機建模方法。

基于上述討論，本文借助部件模型、系統辨識等方法，使用混合信號充分刺激部件模型，將某渦扇發動機建模成NPV 系統，該系統能完整地反映渦扇發動機非線性和狀態參數在大范圍內快速變化的時變特性［13］，更加接近實際對象，為渦扇發動機的控制提供一種新途徑。

1 NPV 模型

考慮如下NPV 系統

式中：x∈Rnx、u∈Rnu和y∈Rny分別為系統狀態、系 統 輸 入 和 輸 出，σ(t)∈Rnσ為 時 變 參 數，A(x，σ(t))、B(x，σ(t))、C(x，σ(t))和D(x，σ(t))為關于x和σ(t)的多項式函數。

工程中許多被控對象可采用NPV 系統作為其模型［2?3］。特別地，當系統（1）的系統矩陣只依賴時變參數時，退化為線性變參數系統

當系統矩陣僅為關于狀態量的多項式時，系統（1）退化成多項式非線性系統

因此，NPV 系統可以看作是多項式非線性系統和LPV 系統的推廣。

值得注意的是，文獻［1］將渦扇發動機模型建成如式（2）所示的LPV 系統，其系統矩陣是關于高壓轉子轉速絕對量的多項式，而狀態則是其增量形式，形如

該建模方式只是原系統在工作點的一階近似，系統參數依舊只是調度參數，忽略了渦扇發動機系統的強非線性。本文借助部件級模型和系統辨識等方法，使用混合信號充分激勵部件級模型，把渦扇發動機系統建模成如下多項式模型

繼而建立了渦扇發動機全包線的NPV 模型，系統矩陣是發動機狀態和時變參數的多項式，顯式地體現渦扇發動機的非線性和時變特性，彌補LPV 系統建模的不足。

2 渦扇發動機大包線建模方法

以某型小涵道比雙軸渦扇發動機為研究對象，其非線性數學模型可以描述為［14］

式中：x∈Rnx、u∈Rnu和y∈Rny分別為系統狀態、系統輸入和輸出。

2.1 單個穩態點狀態變量模型

設(xe，ue，ye)為 渦 扇 發 動 機 系 統 的 一 個 穩 態點，該穩態點附近系統（6）可化為

式中：Ae、Be、Ce、De為待定系統矩陣。

本 文 選 擇x=[nLnH]T，u=[Wfb A8]T，y=[nHπ]T，其中，nL為低壓轉子轉速，nH為高壓轉子轉速，Wfb為主燃油量，A8為喉道面積，π為發動機的增壓比。

設采樣周期為T，kT時刻的狀態、輸入和輸出分別記為x(k)、u(k)和y(k)，采用歐拉積分的方法離散化如下

假設發動機非線性模型在穩態點(xe，ye)附近受擾動信號 Δu生成了輸入輸出數據：{(x(k)，y(k))，k=0，1，…，N}，N為采樣個數。

記偏差量數據矩陣為

基于最小二乘法，在部件級模型的狀態量與狀態空間模型的狀態量誤差平方和最小情況下，則系統矩陣參數的解為

由式（9）可確定在穩態點(xe，ue，ye)模型（7）。

2.2 典型工作點的多項式模型

在某一飛行條件下，即取固定高度H和馬赫數Ma，選 取N1個 穩 態 點 (xei，uei，yei)，i=1，2，…，N1。對穩態點(xei，uei，yei)，按2.1 節中方法得到模型參數矩陣為(Aei，Bei，Cei，Dei)。采用多項式擬合方法，建立典型工作點的發動機多項式非線性模型

非線性系統（13）中系統矩陣求解步驟如下：

（1）記

則由實例點構成的集合為

記F是由xei到Aei的函數集合，則

令h(x)=ωx，ω=[ω1…ωN1]T為 多 項式模型的系數向量。

由于渦扇發動機nH更能近似反映發動機的機械載荷和熱載荷變化狀況［16］，以及直接反映動態特性變化，故式（13）中系統矩陣只為nH的函數。

2.3 大包線NPV 模型

在飛行包線內，選取N2個典型工作點(Mai，Hi)，i=1，2，…，N2。按2.2 節中方法建立第i(i=1，2，…，N2)個典型工作點的多項式模型

在 包 線 0 ≤Ma≤1.2，0 ≤H≤8 內，將(A(x)(Mai，Hi)，B(x)(Mai，Hi)，C(x)(Mai，Hi)，D(x)(Mai，Hi))，i=1，2，…，N2，利 用 回 歸 算 法 分 別 擬 合 成 如 下形式

σ(t)=(Ma，H)是關于高度和馬赫數的多項式函數，則得到渦扇發動機大包線模型。

3 某型渦扇發動機建模

本部分基于系統辨識的方法在不同飛行條件下獲取不同高壓轉速下的狀態空間模型，再通過多項式擬合和回歸方法獲取各系統矩陣關于時變參數與狀態變量的多項式描述。為了消除數據特征之間的量綱影響，將調度參數歸一化至[0，1]。

3.1 激勵信號及辨識算法

根據流量連續和功率平衡的共同工作方程建立的部件級模型，對其信號激勵得到輸入輸出數據。采用階躍信號沖擊部件級模型得到的輸入輸出數據無法刻畫渦扇發動機的動態特性，且一步求解最小二乘問題時擬合程度不高。若要達到較好的辨識效果，則需多步系統辨識，計算量較大且動態過程一般。因此，有學者采用偽隨機序列［17］（Pseudo random binary sequence，PRBS）沖擊部件級模型，PRBS 是由周期性數字序列經過濾波處理得到，具有可重復性、自相關性能好的特點，如圖1所示。經過仿真驗證，相比使用階躍信號，采用PRBS 能夠充分激勵部件級模型，得到更好的結果，但其結果模型泛化性能很差。

圖1 偽隨機序列Fig.1 Pseudo random binary sequence

本文結合機器學習中訓練集的思想，提出一種基于階躍信號、正弦信號以及偽隨機信號的混合信號。采用上下界幅度在穩態燃油幅度的±2%左右的混合信號，即一定幅度的M 序列，階躍信號以及一定幅值和頻率的正弦信號，如圖2所示。相比使用單一信號，采用混合信號可以充分激勵部件模型，得到更好的結果，且結果模型泛化性能較好。

圖2 混合信號Fig.2 Mixed signal

辨識采用的最小二乘算法為Levenberg?Mar?quardt（L?M）估計算法［18］。L?M 算法的每一步迭代都可以調整阻尼項以確保誤差下降。當遠離最優解時，算法接近最速下降法。反之，則接近高斯牛頓（Gauss?Newton，G?N）算法，利用二階導數的信息，快速收斂到最優解。

在一定飛行條件和某工作狀態下分別使用階躍信號、偽隨機信號和混合信號激勵部件級模型，利用L?M 算法對部件模型的輸入輸出數據進行動態辨識，驗證激勵信號對于系統辨識的影響。本文以地面條件為例，nH的誤差大小作為評價辨識效果的好壞，結果分別見表1 和圖3～7，其中，圖3、5、7 中兩條曲線分別代表部件級模型的輸出轉速和辨識結果的輸出轉速曲線。圖4 為偽隨機信號激勵下辨識模型與部件級模型階躍響應曲線。圖6為混合信號激勵下辨識模型與部件級模型階躍響應誤差曲線。

圖3 偽隨機信號激勵下的系統辨識結果Fig.3 System identification results under pseudo random signal excitation

圖4 偽隨機信號辨識結果的驗證Fig.4 Validation of identification results of pseudo random signal systems

圖5 混合信號激勵下的系統辨識結果Fig.5 Identification results with mixed signal excitation

圖6 混合信號激勵下的系統辨識誤差Fig.6 System identification error with mixed signal excita?tion

表1 激勵信號的辨識誤差及驗證誤差Table 1 Identification error and validation error of exci?tation signal %

圖7 G-N 算法下的系統辨識結果Fig.7 System identification results under G-N algorithm

由圖3、4 可知，偽隨機信號輸入有很好的動態辨識效果，但在驗證的過程中，該狀態空間模型的泛化性能很差，出現發散現象。

由圖5、6 可知，采用混合信號輸入的辨識結果較好?；旌闲盘柼岣吡俗钚《说臄M合程度，且結果模型具有較好的泛化能力，在驗證過程中誤差不超過0.02%。

由圖5、7 可知，求解渦扇發動機系統辨識最小二乘問題，采用L?M 估計方法的平均誤差為0.01%，采用G?N 算法得到的最優值的平均誤差為0.55%。因此，采用L?M 估計方法求解的精度更高。

3.2 典型工作點的多項式模型驗證

通過混合信號與L?M 算法系統辨識得到地面條件下20 個工作狀態的線性系統，通過多項式擬合建立地面條件下的多項式非線性系統為

非線性系統（15）的矩陣參數擬合結果如圖8～11 所示（篇幅有限，僅列出A(x)擬合結果）。

圖8 系數a11 擬合曲線Fig.8 Fitting curve of a11

為了檢驗多項式非線性系統的精度，將其與非線性動態熱力學模型［19］分別做較大范圍仿真驗證，仿真結果如圖12、13 所示。

由圖12、13 可知，本文建立的多項式非線性系統能夠較精確地反映原系統動態特性，且穩態誤差在0.5%以內。由于篇幅有限，未列出其他飛行條件下的多項式非線性模型。

圖9 系數a12 擬合曲線Fig.9 Fitting curve of a12

圖10 系數a21 擬合曲線Fig.10 Fitting curve of a21

圖12 nH=0.89 至nH=0.95 仿真驗證Fig.12 Simulation verification from nH=0.89 to nH=0.95

圖13 nH=0.95 至nH=0.99 仿真驗證Fig.13 Simulation verification from nH=0.95 to nH=0.99

3.3 大包線NPV 模型驗證

本部分將各飛行條件（表2）的模型通過回歸算法得到大包線NPV 模型如式（16）所示。

表2 大包線H 和MaTable 2 Large envelope H and Ma

式中aij，bij，cij和d2j（i，j=1，2）分別為

為了檢驗NPV 模型的精度，將NPV 系統與非線性動態熱力學模型做仿真驗證，仿真結果如圖14～17 所示。

由圖14～17 可知，NPV 系統與非線性動態熱力學模型誤差較小，不超過1%。

圖14 Ma=0.3,H =3 時nH=0.91～1.04 仿真驗證Fig.14 Simulation verification of nH=0.91—1.04 with Ma=0.3,H =3

圖15 Ma=1.2,H =8 時nH=1.0～1.044 仿真驗證Fig.15 Simulation verification of nH=1.0—1.044 with Ma=1.2,H =8

圖16 Ma=0.7,H =6 時nH=0.973～1.039 仿真驗證Fig.16 Simulation verification of nH=0.973—1.039 with Ma=0.7,H =6

圖17 Ma=0.2,H =2 時nH=0.881～0.946 仿真驗證Fig.17 Simulation verification of nH=0.881—0.946 with Ma=0.2,H =2

4 結 論

本文將一種混合信號作為激勵信號，采用L?M 算法求解最小二乘問題，采用系統辨識方法，建立了渦扇發動機NPV 模型。NPV 系統矩陣參數能夠根據高度、馬赫數以及高壓轉子轉速實時改變，能更好地刻畫渦扇發動機的強非線性、時變特性以及大轉速范圍的動態特性，為渦扇發動機的控制設計提供了一種新途徑，是進一步研究渦扇發動機穩態及過渡態控制的基礎。在該NPV 建模方法的基礎上，下一步將開展增量形式NPV 模型建模研究。