多回路氣動伺服彈性系統穩定裕度的變結構μ分析方法

劉祥， 孫秦

西北工業大學 航空學院， 西安 710072

多回路氣動伺服彈性系統穩定裕度的變結構μ分析方法

劉祥， 孫秦*

西北工業大學 航空學院， 西安 710072

作為一個多輸入多輸出(MIMO)系統，氣動伺服彈性系統的各控制回路是相互耦合的，但對于各控制回路的穩定裕度目前尚無統一的計算方法。針對MIMO控制系統的穩定裕度計算問題，首先分析了現有的回差陣奇異值方法和μ分析方法，并指出了2種方法各自保守性的來源。在此基礎上，提出了一種變結構μ分析方法，通過迭代調整擾動模型結構來求解穩定裕度，并從理論上證明了算法的單調收斂特性。以某彈性飛機的陣風減緩控制系統為例進行了穩定裕度分析。3種方法結果的對比表明，本文方法能夠有效降低分析結果的保守性。

氣動伺服彈性； 多輸入多輸出系統； 耦合； 穩定裕度； 變結構μ分析

作為涉及柔性飛行器結構、非定常氣動力和飛行控制系統三者相互作用的多學科技術，氣動伺服彈性技術已成為現代飛行器設計的關鍵內容。由于關系到飛行安全，氣動伺服彈性系統的穩定性分析問題一直備受研究者關注[1-2]。在實際工程中，要實現機動控制、載荷減緩、乘坐品質改善及顫振抑制等目的需要多個控制面的參與，這使控制系統具有多個輸入輸出通道，如何有效分析并衡量控制回路的擾動穩定性依然是具有挑戰性的問題。

經典單輸入單輸出(SISO)線性系統的穩定裕度一般通過幅值裕度和相位裕度來度量，且有多種精確計算方法。然而對于各回路之間存在耦合的多輸入多輸出(MIMO)線性系統，若直接采用SISO線性系統的分析方法來求解穩定裕度，得到的結果往往不盡人意[3]。文獻[4-5]通過在控制回路中引入對角擾動矩陣并計算回差陣最小奇異值的方法，得到了MIMO線性系統各控制回路的幅值和相位同時變化時的穩定裕度，并獲得了廣泛的工程應用[6-7]。Doyle[8]通過定義結構奇異值μ(ω)來分析具有結構不確定性的MIMO線性系統的魯棒穩定性。此后，Lind和Brenner[9]提出了基于μ分析的多回路系統穩定性評估方法，可用來確定不同類型的模型不確定性同時存在時氣動伺服彈性系統的穩定飛行范圍。Tsao等[10]在μ方法的基礎上通過構造雙線性變換模型給出了判斷多回路系統是否滿足給定幅值和相位裕度的充分條件。Vartio等[11]采用μ方法檢驗了多回路系統的幅值裕度和相位裕度。針對具有多個實參數擾動的系統，文獻[12-13]根據映射理論提出了一種穩定裕度的精確計算方法，但這種算法的復雜度太高且迭代計算量過大。Al-Shamali等[14]基于頻域下的Nyquist理論定義了臨界擾動半徑與臨界方向等概念，并以此為基礎計算SISO線性系統的Nyquist魯棒穩定裕度，Son等[15]將這一理論推廣到了MIMO線性系統。李信棟和茍興宇[16]利用系統回差陣提出了2種穩定裕度改進方法，并將兩者的結果與回差陣奇異值法的結果綜合得到了保守性相對較低的結果。

針對多回路氣動伺服彈性系統的穩定裕度問題，本文首先分析了現有的回差矩陣奇異值方法和μ分析方法，發現回差陣奇異值方法的保守性來自于方法本身，而μ分析方法的保守性主要是因建模時擾動模型的構造方式所致。為此，本文提出了一種變結構μ分析方法，通過迭代調整擾動模型的結構求解多回路控制系統的穩定裕度，可較好地提高穩定裕度估計精度，并對迭代參數施加約束以保證算法的快速收斂特性。最后，以一彈性飛機為例，采用不同方法計算比較了其陣風減緩控制系統的穩定裕度，驗證了本文方法的優勢。

1 回差陣奇異值法

本節對回差陣奇異值法進行理論推導。系統模型如圖1所示，其中P(s)和K(s)分別為被控對象和控制系統，定義系統回差陣為I+K(s)P(s)。在輸入端引入擾動量測陣

E(s)=diag{klexp(iφl)}l=1,2,…,n

(1)

若閉環系統是漸進穩定的，則回差陣必滿足det(I+KPE)≠0，根據奇異值的性質，則有

(2)

穩定裕度即研究系統在穩定范圍內，量測陣

圖1 系統模型
Fig.1 System model

E(s)中kl和φl同時變化的最大容許值。由E(s)的定義可知kl不為零，故E(s)非奇異。而由原閉環系統穩定可知I+KP非奇異，利用矩陣分離特性可得

I+KPE=

(3)

進一步，若閉環系統穩定則

(4)

(5)

結合式(4)和式(5)可得系統穩定的充分條件為

(6)

應用式(1)，式(6)等價于

(7)

記回差陣最小奇異值為m,根據式(7)，分別令所有φl=0°或所有kl=1，可得到多回路線性系統的幅值裕度GM和相位裕度PM表達式為

(8)

通過上面的分析可看出，保證系統穩定的條件式(6)是一個保守條件。對于某些具體問題，這種近似可能帶來過大的保守性，以致影響到對控制系統的選擇和評估。

2 μ分析方法

本節介紹多回路系統穩定裕度的μ分析方法。μ方法的基本處理途徑是將模型分解為名義模型部分(即確定性部分)和幅值有限的不確定性部分。為此，將式(1)中的擾動量測陣表示為

E(s)=I+Δ=diag{klexp(iφl)}

l=1,2,…,n

(9)

式中：Δ為對角復數不確定性矩陣，即

Δ=diag{δ1,δ2,…,δn}δl∈C;l=1,2,…,n

結合式(9)，圖1可轉換為圖2中的魯棒穩定性分析結構，其中系統傳遞函數矩陣為

圖2 魯棒穩定性分析結構
Fig.2 Robust stability analysis structure

(10)

對于M∈Cn×n，定義其結構奇異值為

μΔ(M)=

(11)

若無Δ1∈Δ使det(I-MΔ1)=0，則μΔ(M)=0。

從式(11)中的定義可以看出，μΔ(M)的倒數是導致反饋系統失穩的攝動量最小值。因此可得到閉環系統魯棒穩定的充要條件為

(12)

結合式(9)，式(12)等價于

l=1,2,…,n

(13)

分別令所有φl=0° 或所有kl=1，可得到幅值裕度GM和相位裕度PM表達式為

(14)

通過以上分析可以看出，式(12)是圖2中模型魯棒穩定的充要條件。但因式(9)中Δ各對角元素的幅值|δl|與E(s)各對角元素的幅值|kl|之間并非單調關系，從M的魯棒穩定性出發估計圖1中系統的穩定裕度勢必會為結果引入保守性。另外，式(14)中幅值裕度與相位裕度間的一一對應關系并不符合客觀情況，這也從另一方面證實了此方法的保守性。

3 變結構μ分析方法

為克服式(9)中由擾動量測陣的構造方法產生的保守性。在其基礎上引入參數γ并重新定義E(s)為

E(s)=γI+Δ=diag{klexp(iφl)}

l=1,2,…,n

(15)

此時圖1的等價模型如圖3所示。

圖3 有攝動的多回路系統
Fig.3 Multiloop system with perturbation

以式(15)為基礎構造圖2中所示的魯棒穩定性分析模型，可得

(16)

對應不同的γ取值，可求得不同的幅值裕度和相位裕度。下面對γ的取值方法進行討論。

當固定γ的取值時，由式(12)和式(15)可知閉環系統穩定的充分必要條件為

l=1,2,…,n

(17)

式中:μΔ,γ(M)為給定γ時的結構奇異值。令式(17)中所有φl=0° 可得

l=1,2,…,n

(18)

(19)

令式(17)中所有kl=1可得

(20)

易推得式(20)有解的充分必要條件為

(21)

這等價于

(22a)

(22b)

即γ在取值時必須滿足式(22a)和式(22b)的約束。

(23)

將式(22a)和式(23)兩側分別相加可得

(24)

參考式(24)，以γ為參變量迭代求解圖3中多回路系統的穩定裕度。令初始值為1，取迭代公式為

(25)

定理1根據式(25)迭代求解穩定裕度時，增益裕度單調遞增且收斂。

證明證明過程分以下2步：首先證明增益上界隨γ的增大而增大，再證明迭代過程中增益裕度的單調收斂性。

1) 對于第j次和第j+1次迭代，由式(17)可知其幅值和相位邊界分別滿足

(26a)

(26b)

這2條邊界必有交點。當γj+1>γj時，易知在交點處式(26)各變量滿足圖4中的三角形關系。

此時，2次迭代過程中的增益上界滿足

|AD|+|CD|-|AB|-|CB|=

|BD|+|CD|-|CB|>0

(27)

即增益上界隨γ的增大而增大。

2) 當γ=γ1時顯然滿足式(22a)和式(22b)。現假設γ=γj時同樣滿足式(22a)和式(22b)，則當γ=γj+1時滿足

圖4 變量關系圖
Fig.4 Variable relation graph

j=1,2,…

(28)

由數學歸納法知γ是單調遞增的，故迭代過程中增益上界也是單調遞增的。參考圖4可得，增益下界滿足

|AD|-|CD|-|AB|+|CB|=

|BD|+|CB|-|CD|>0

(29)

(30)

即增益下界收斂于1。式(30)即可作為迭代過程的收斂判據。

綜上可知，增益裕度單調遞增且收斂。定理1得證。

實際上每次迭代過程得到的穩定區域邊界與系統的真實穩定區域邊界至少有1個交點，而各次迭代得到的穩定區域的并集可以一定程度上逼近實際的穩定區域，因此在迭代過程中減小γ的增量可以得到保守性更低的穩定區域。但為了加速幅值裕度的收斂過程仍依據式(25)更新γ參數。雖然在迭代過程中相位裕度并無單調收斂特性，但可以通過將其取為各次迭代中的最大值以最大程度地降低保守性。

根據以上討論可將迭代過程描述如下：

1) 令迭代序號j=1，γj=1，PM=0°。

2) 根據式(16)求解M，并根據μ分析理論計算結構奇異值μj。

(31)

若PM<φj，則令φj→PM，其中“→”表示賦值運算。

4) 判斷迭代是否收斂。若是，轉至步驟5)；否則根據式(25)更新γ，并令j+1→j，再轉至步驟2)。

5) 按式(32)計算GM，結束循環。

(32)

需要注意的是，本文針對的是多輸入多輸出系統控制回路的穩定裕度計算問題，因此并未具體考慮模型的參數不確定性和未建模不確定性，這似乎限制了本文方法的用途。但實際上并非如此，這可以從以下2個方面解釋：

首先，本文的穩定裕度計算問題是一個相對獨立的問題，對于系統參數變化或建模誤差較小情況下的系統穩定裕度可得到相對其他方法保守性更低的結果。

其次，對于模型的參數不確定性或未建模不確定性不可忽略的情況，可將其此類不確定性加入文中的算法結構，并與式(9)中代表回路增益和幅值變化的復數不確定塊Δ合并為階數更高的不確定塊。然后通過調整算法結構可以進一步分析此類模型的控制回路穩定裕度。但這樣處理的價值和意義有待商榷，因為對于模型不確定性不可忽略的情況，系統的魯棒穩定性分析結果往往是以不確定參數或未建模不確定性的魯棒穩定范圍的形式給出[9]，且已有成熟的分析方法。當然，也可以將本文方法直接應用于名義模型，得到的穩定裕度結果可以作為前述結果的補充，從控制回路的角度反映閉環系統的魯棒穩定性。

4 算 例

4.1 算例模型1

為驗證本文方法的有效性，首先對文獻[16]中的雙輸入雙輸出控制系統進行穩定裕度分析。被控對象和控制器的傳遞函數矩陣分別為

表2給出了分別由本文方法、回差陣奇異值法和μ分析方法計算出的控制回路穩定裕度。易知本文的變結構μ分析方法具有最低的保守性。

表1 增益上界、增益下界和相位裕度的收斂過程Table 1 Convergence process of ， and PM

表2 不同方法對算例1的穩定裕度計算結果

4.2 算例模型2

以文獻[17]中的彈性飛機模型為基礎，對設計出的H∞陣風減緩控制系統進行穩定裕度分析。飛機的氣動力模型如圖5所示，飛行馬赫數和高度分別為0.3和5 000 m。具體模型數據見文獻[18]。

圖5 飛機氣動力模型
Fig.5 Airplane aerodynamic model

在構造全機狀態空間模型時，結構模態取為全機對稱的前7階彈性模態及剛體俯仰和沉浮模態。在構造時域非定常氣動力模型時，采用最小狀態法[17]對由ZAERO軟件計算得到的頻域廣義氣動力進行拉氏域有理擬合，選取4個氣動力滯后根并對其進行非線性優化以減小擬合誤差。滯后根最終取值為[-0.110 -0.297 -0.307 -0.530]。傳感器輸出為重心處俯仰速率及翼尖過載，控制面采用升降舵和副翼，對應的舵機傳遞函數均為

采用離散“1-cos”型陣風模型，垂直陣風速度在飛行軌跡方向的分布如圖6所示。

根據MATLAB的魯棒控制工具箱設計H∞控制器，設計目標為減緩陣風引起的翼根彎矩和俯仰速率響應。因控制器的設計過程與本文主題無關，故在此不予贅述，而僅在附錄A中給出其降階后的狀態空間模型。易知控制系統為雙輸入雙輸出的多回路耦合系統，系統框圖如圖7所示。圖8給出了翼根彎矩和俯仰速率的開環和閉環離散陣風響應。從圖中可看出閉環系統穩定，且翼根彎矩和俯仰速率得到了有效抑制。

采用本文方法分析上述控制回路穩定裕度時，迭代過程在第4次循環達到收斂。表3給出了由本文方法、回差陣奇異值法和μ分析方法計算出的控制回路穩定裕度結果，可以看出本文方法的保守性最低。為保持閉環系統穩定，兩條控制回路的相位不變時可容許的最大增益為8.68 dB；兩條回路的增益不變時可容許的最大相位滯后為63.0°。

圖6 離散陣風
Fig.6 Discrete gust

圖7 彈性飛機陣風減緩控制方案
Fig.7 Gust alleviation control scheme of elastic aircraft

下面進一步通過時域仿真對上述結果進行驗證。根據計算結果可知各回路相位不變時幅值可同時增大2.71倍。圖9(a)分別給出了幅值同時增大2.71倍時翼尖過載及飛機俯仰速率的陣風響應，圖9(b)分別給出了幅值同時增大2.74倍時翼尖過載及飛機俯仰速率的陣風響應。

從圖9可以看出，控制回路增益同時增大2.71 倍時系統是穩定的，而增益繼續增大到2.74倍時系統臨界不穩定。由上述結果可知，利用本文提出的變結構μ分析方法得到的穩定裕度已足夠逼近系統真實穩定裕度。

圖8 翼根彎矩和俯仰速率的離散陣風響應
Fig.8 Wing root bending moment and pitch rate response to discrete gust

表3 不同方法對算例2的穩定裕度計算結果

圖9 k=2.71和k=2.74時翼尖過載及俯仰速率的陣風響應
Fig.9 Wing tip overload and pitch rate response to gust with k=2.71 and k=2.74

5 結 論

1) 建立的變結構μ分析方法能夠方便求解氣動伺服彈性系統的穩定裕度，并具有比已有方法更小的保守性。

2) 建立的變結構μ分析方法在迭代過程具有單調收斂特性，能夠快速求解系統的穩定裕度，且迭代次數并不隨問題規模的增大而增大。

[1] HADDADPOUR H. Aeroservoelastic stability of supersonic slender-body flight vehicles[J]. Journal of Guidance, Control, and Dynamics, 2006, 29(6): 1423-1427.

[2] KARPEL M, MOULIN B, CHEN P C. Extension of the g-method flutter solution to aeroservoelastic stability analysis[J]. Journal of Aircraft， 2005, 42(3): 789-792.

[3] DOYLE J C. Robustness of multiloop linear feedback systems[C]//Proceedings of the IEEE Conference on Decision and Control. Piscataway, NJ: IEEE Press, 1978: 12-18.

[4] LEHTOMAKI N A, SANDELL J N R, ATHANS M. Robustness results in linear-quadratic Gaussian based multivariable control designs[J]. Transactions on Automatic Control, 1981, AC-26(1): 75-93.

[5] MUKHOPADHYAY V, NEWSOM J R. Application of matrix singular value properties for evaluating gain and phase margins of multiloop systems[C]//AIAA Guidance and Control Conference. Reston: AIAA, 1982: 420-428.

[6] BURKEN J J. Flight-determined stability analysis of multiple-input-multiple-output control systems: AIAA-1992-4396[R]. Reston: AIAA, 1992.

[7] 吳斌, 程鵬. 多變量飛控系統的穩定裕度分析[J]. 航空學報, 1998, 19(6): 18-22.

WU B, CHENG P. Stability margin analysis of the multiloop flight control systems[J]. Acta Aeronautica et Astro-nautica Sinica, 1998, 19(6): 18-22 (in Chinese).

[8] DOYLE J. Analysis of feedback systems with structured uncertainties[J]. IEEE Proceedings D: Control Theory and Applications, 1982, 129(6): 242-250.

[9] LIND R, BRENNER M. Analyzing aeroservoelastic stability margins using theμmethod[C]//AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics and Materials Conference. Reston: AIAA, 1998: 1672-1681.

[10] TSAO T T, LEE F C, AUGENSTEIN D. Relationship between robustnessμ-analysis and classical stability margins[C]//Proceedings of the 1998 IEEE Aerospace Conference. Part 1 (of 5). Piscataway, NJ: IEEE Press, 1998: 481-486.

[11] VARTIO E J, SHAW E E, VETTER T. Gust load alleviation flight control system design for a sensor craft vehicle[C]//26th AIAA Applied Aerodynamics Conference. Reston: AIAA, 2008: 1-10.

[12] DE GASTON R R E, SAFONOV M G. Exact calculation of the multiloop stability margin[J]. Transactions on Automatic Control, 1988, 33(2): 156-171.

[13] SIDERIS A, DE GASTON R R E. Multivariable stability margin calculation with uncertain correlated parameters[C]//Proceedings of the 25th IEEE Conference on Decision & Control. Piscataway, NJ: IEEE Press, 1986: 766-771.

[14] AL-SHAMALI S A, JI B, CRISALLE O D, et al. The Nyquist robust sensitivity margin for uncertain closed-loop systems[J]. International Journal of Robust and Nonlinear Control, 2005, 15(14): 619-634.

[15] SON J E, MANI A S, LATCHMAN H A. Robustness analysis for MIMO systems with unstructured uncertainties[C]//2010 8th IEEE International Conference on Control and Automation. Piscataway, NJ: IEEE Press, 2010: 1333-1337.

[16] 李信棟, 茍興宇. 多輸入多輸出線性定常系統穩定裕度的分析與改進[J]. 控制理論與應用， 2014, 31(1): 105-111.

LI X D, GOU X Y. Analysis and improvement of stability margin for multi-input multi-output linear time-invariant systems[J]. Control Theory & Applications， 2014, 31(1): 105-111 (in Chinese).

[17] KARPEL M, MOULIN B, CHEN P C. Dynamic response of aeroservoelastic systems to gust excitation[J]. Journal of Aircraft， 2005, 42(5): 1264-1272.

[18] ZONA Technology, Inc. ZAERO application’s manual Vol.2[M]. Scottsdale: ZONA Technology, Inc. 2008: 121-198.

Variable-structureμsynthesismethodforstabilitymarginofmultiloopaeroservoelasticsystem

LIUXiang,SUNQin*

SchoolofAeronautics,NorthwesternPolytechnicalUniversity,Xi’an710072,China

Aeroservoelasticsystemisamulti-input-multi-output(MIMO)systemwithitscontrolloopscoupledwitheachother,whilecurrentlythereisnounifiedtheoryinregardtothestabilitymarginofthecontrolloops.Thispaperaddressesthisproblembyfirstanalyzingtheexistingreturndifferencematrixmethodandtheμsynthesismethod,andthenproposinganewmethodtoreducetheconservativenessduringanalysis.Thisnewmethod,calledthevariable-structureμsynthesismethod,solvesthestabilitymarginthroughadjustingtheperturbationmodelinseveraliterationsandisprovedtohavefastconvergenceproperty.Finally,differentmethodsareusedtocalculatethestabilitymarginsofthegustalleviationsystemofaflexibleaircraft,andtheresultsshowthattheproposedapproachislessconservativecomparedwiththeexistingmethods.

aeroservoelasticity;multi-input-multi-outputsystem;coupling;stabilitymargin;variable-structureμsynthesis

2016-04-21；Revised2016-05-19；Accepted2016-06-21；Publishedonline2016-06-270836

URL：www.cnki.net/kcms/detail/11.1929.V.20160627.0836.004.html

.E-mailsunqin@nwpu.edu.cn

2016-04-21；退修日期2016-05-19；錄用日期2016-06-21； < class="emphasis_bold">網絡出版時間

時間：2016-06-270836

www.cnki.net/kcms/detail/11.1929.V.20160627.0836.004.html

＊

.E-mailsunqin@nwpu.edu.cn

劉祥， 孫秦． 多回路氣動伺服彈性系統穩定裕度的變結構μ分析方法J． 航空學報,2017,38(4):120350.LIUX,SUNQ.Variable-structureμsynthesismethodforstabilitymarginofmultiloopaeroservoelasticsystemJ.ActaAeronauticaetAstronauticaSinica,2017,38(4):120350.

http://hkxb.buaa.edu.cnhkxb@buaa.edu.cn

10.7527/S1000-6893.2016.0201

V249.1

A

1000-6893(2017)04-120350-09

(責任編輯： 鮑亞平， 張晗)

附錄A:

控制器狀態空間方程為

yc=Ccxc+Dcuc

式中：xc、uc和yc分別為控制器的狀態變量、輸入和輸出；Ac、Bc、Cc和Dc分別為系數矩陣，取值如下：

Ac=

a11=-61.714 4,a12=328.244 5,a21=-328.244 5,a22=-61.714 4,a33=-224.140 3,a34=161.166 2,a43=-161.166 2，a44=-224.140 3,a55=-40.325 9,a56=244.077 8,a65=-244.077 8,a66=-40.325 9,a77=-18.982 0,a78=234.244 1,a87=-234.244 1,a88=-18.982 0,a99=-13.043 2,a9 10=211.119 6,a10 9=-211.119 6,a10 10=-13.043 2,a11 11=-53.624 0,a11 12=65.750 4,a12 11=-65.750 4,a12 12=-53.624 0,a13 13=-62.856 0,a14 14=-17.728 9,a14 15=63.133 9,a15 14=-63.133 9,a15 15=-17.728 9,a16 16=-24.524 5,a16 17=54.873 4,a17 16=-54.873 4,a17 17=-24.524 5,a18 18=-9.759 7,a19 19=-0.005 1,a20 20=-0.936 1