孫明敏，陸宇平，徐 亮，沈華勛

（南京航空航天大學 自動化學院，江蘇 南京210016）

近年來，高空長航時（HALE）飛行器因其超長續(xù)航能力日益受到重視。這類飛行器不但可以執(zhí)行軍事上的機載情報、監(jiān)視與偵察等任務， 而且可用于民用上的網(wǎng)絡通信以及一般大氣研究[1-2]，具有廣闊的應用前景和研究價值。 然而由于任務需要，理想HALE 飛行器具有機翼展弦比大、機身輕而薄的特點，大柔性飛行器（VFA）就是此類飛行器最典型的代表。VFA 的大柔性結(jié)構(gòu)可能導致飛行器在常規(guī)飛行條件下發(fā)生大的變形、結(jié)構(gòu)低頻顫振與剛體運動耦合，以及氣動失速等問題， 具有與剛體飛行器顯著不同的氣動彈性和飛行動力學特性[3-4]，這意味著飛行器的模型十分復雜，維數(shù)較高，且存在較多非線性，給控制設計帶來困難。因此，在控制設計之前，有必要研究VFA 模型的模型降階，便于對此類復雜系統(tǒng)進行理論分析，減少數(shù)據(jù)運算量。

本文采用基于平衡實現(xiàn)理論的平衡殘差降階方法，針對某高階大柔性飛行器線性化模型，進行了模型降階，并利用基于狀態(tài)空間Newmark 法給出了詳細的降階模型仿真分析。

1 大柔性飛行器模型

本文研究對象是基于美國國防部高級研究計劃局（DARPA）“禿鷹”計劃所設計的大柔性飛行器，簡稱“禿鷹”模型[5]。 “禿鷹”模型的結(jié)構(gòu)及坐標系如圖1 所示，該飛行器擁有15 個向前的推力發(fā)動機和6 個升降舵，其中兩個升降舵位于機翼兩端，剩余4 個升降舵位于中間尾翼的尾部。 需注意，翼根、吊艙和發(fā)動機坐標系與慣性坐標系一致。

“禿鷹”模型是一個高保真線性動力學模型，由真實測量數(shù)據(jù)生成，在本文中基于平衡點將數(shù)據(jù)進行線性化處理，其中一個平衡點為：

模型有n=707 個狀態(tài)量，m=381 個輸入量和p=404 個輸出量，在一個線性化單點飛行條件下，其開環(huán)線性時不變狀態(tài)空間形式如下：

該模型描述了慣性空間6 自由度(x,y,z,Φ,θ,Ψ)及其速度(u,v,w,p,q,r)的動力學方程，特別之處在于模型中引入了空氣來流wi、機體彈性位置xflex和彈性速度vflex：

381 個輸入量中有33 個控制輸入 （15 個發(fā)動機推力、6個升降舵偏轉(zhuǎn)角δ、6 個偏轉(zhuǎn)角速度δ˙和6 個偏轉(zhuǎn)角加速度δ¨），剩余348 個輸入表示陣風擾動：

圖1 “禿鷹”飛行器結(jié)構(gòu)及坐標系圖Fig. 1 Vehicle structure and coordinate systems for vulture

404 個輸出量中，翼根傳感器測量飛行器重心處的18 個變量，每個吊艙傳感器測量兩個變量：局部迎角和側(cè)滑角。 每 個節(jié)點傳感器測量18 個局部狀態(tài)：

傳統(tǒng)控制技術(shù)對高階“禿鷹”模型不再適用，因此必須對其進行狀態(tài)空間模型降階，將較大系統(tǒng)轉(zhuǎn)化為一個近似的較小系統(tǒng)，使得降階模型仍能保持原模型的一些性能，便于控制器的設計。

2 模型降階

本文主要采用基于平衡實現(xiàn)的平衡殘差降階方法，來簡化系統(tǒng)模型。

2.1 平衡實現(xiàn)和平衡殘差法

平衡實現(xiàn)理論主要由Moore[6]提出。 考慮如下線性時不變系統(tǒng)：

定義系統(tǒng)可控性和可觀測性Gram 矩陣Wc，Wo分別為:

如果系統(tǒng)是穩(wěn)定可控的，則Wc滿秩；如果系統(tǒng)是穩(wěn)定可觀的，則Wo滿秩。 進而可知，線性可控性和可觀測性Gram 矩陣Wc，Wo是下面Lyapunov 方程的唯一正定解：

對于系統(tǒng){A,B,C,D}，如果存在非奇異變換矩陣T，將原系統(tǒng)變換為等價系統(tǒng)，其中相應可控性和可觀測性Gram 矩陣均為對角陣且相等，即

式中σ1≥σ2≥…≥σn≥0，則稱}系統(tǒng)為原系統(tǒng)的平衡實現(xiàn)，為系統(tǒng)的Hankel 陣奇異值。

Hankel 陣奇異值提供了狀態(tài)重要性的測度。 即大奇異值對應的狀態(tài)受控制輸入影響最大， 而輸出也受該狀態(tài)變化的影響最大。 因此，對應最大奇異值的狀態(tài)影響系統(tǒng)輸入輸出行為最大。如果存在n＞r＞0，滿足σr+1，σr+2，…，σn遠小于σ1，σ2，σr則可將σr以后對應的狀態(tài)變量進行剩余殘化， 僅保留σr之前的狀態(tài)。 具體方法是將平衡系統(tǒng)進行矩陣分塊：

則模型降階為

2.2 大柔性飛行器模型降階

對于本文研究對象“禿鷹”模型，定義對象動力學模型(2)為，不考慮陣風擾動，則，取，對應于前33 個輸入量的數(shù)據(jù)。設定截斷值，則忽略奇異值小于的部分，保留大于的部分，對應的降階后模型階數(shù)為。需要注意的是，對原模型進行模型降階前，應該先分離不穩(wěn)定的模態(tài)，降階后再把該部分加上。

圖2 分別給出了原模型和降階后模型的傳遞函數(shù)矩陣最大和最小奇異值。 由圖可知，80 階系統(tǒng)在頻率低于100 Hz 時與原模型曲線基本一致，充分保持了原系統(tǒng)特性。

圖2 原模型和降階模型傳函矩陣最大、最小奇異值對比Fig. 2 Comparison of the maximum and minimum singular values of the transfer function matrix between the original model and the reduced－order model

3 狀態(tài)空間Newmark 法求解方程

本文研究系統(tǒng)為一階線性方程組， 由于該系統(tǒng)方程維數(shù)過大，而Matlab 中常用的四階龍格－卡塔等算法精度雖高，但求解過程太過復雜，時間太長，在求解該線性系統(tǒng)模型時1 小時只能完成1 s 時間仿真運算。 且在運算中最好采取變步長算法，若設置定步長大于0.000 1 s 時，在運算一段時間后會出現(xiàn)計算錯誤無法求解情況。 因此需要選擇合適的一階微分求解方法進行數(shù)值運算， 這里采用基于平均速度的狀態(tài)空間Newmark 求解法[7]，具體方法如下。

通過引入狀態(tài)空間向量，一個一般動力學系統(tǒng)的二階微分方程可以寫成以下一階狀態(tài)空間方程形式：

其中，設定{q}=[{x˙}，{x}]T為狀態(tài)向量，{x}為位移向量，[A]為系統(tǒng)運動矩陣，{F}為力向量，[A]和{F}的具體形式如下：

[M]，[C]，[K]分別為系統(tǒng)的慣性、阻尼和剛度矩陣，{f(t)}為外力向量。 如圖3 所示，在時間tn和tn+1（或者步長和）的間隔△t 間定義新的變量τ(0≤τ≤△t，定義內(nèi)平均速度{q˙(t)}為：

令τ=0 時，{q(τ)}={q}n則

根據(jù)以上方程，令τ=△t 或t=tn+1，有

圖3 狀態(tài)空間Newmark 法的平均速度示意圖Fig. 3 Constant average velocity for the state－space newmark method

則時刻的速度微分方程為

將原始一階方程代入上式，可以得到tn+1時刻速度的解

由于本文所研究模型已經(jīng)是一階線性模型， 故狀態(tài)向量即為，則

由以上表達式可知第時刻的速度可以由第時刻的速度求得，這種基于平均速度的Newmark 方法比傳統(tǒng)方法更簡單直觀， 具有很高的可靠性和準確性， 且更容易在軟件程序中實現(xiàn)。 實驗證明此方法求解方程精度較高，與使用四階龍格－庫塔法求解結(jié)果幾乎沒有差別， 但是運算時間要比龍格－庫塔法快至少100 倍，非常適用于高階微分方程的求解。

4 仿真分析

下面給定相同輸入條件，對80 階降階系統(tǒng)和原系統(tǒng)輸出曲線進行對比分析。

假設起始時刻飛行器處于靜穩(wěn)定狀態(tài), 令均勻分布于主翼的15 個發(fā)動機在0～10 s 同時產(chǎn)生10lbs 的推力，仿真時間30 s， 則飛行器受到均勻推力作用時的輸出曲線如圖4～6 所示，其中翼根處測量的是飛行器質(zhì)心的狀態(tài)量，吊艙3 和測量點11 的位置見圖1。

施加向前推力主要影響飛行器縱向運動特性， 水平速度增大，則升力增大，垂直速度隨之變大，高度上升（翼根坐標系中向下為正）；飛行器首先產(chǎn)生抬頭力矩，隨著升力增大產(chǎn)生負的俯仰力矩使得迎角減小。由于推力均勻施加在主翼上，橫側(cè)向運動主要是由機翼各模塊相互耦合作用的， 受推力直接影響不大。

圖4 飛行器翼根處速度及角度變化曲線Fig. 4 Velocities and angles measured at the wing root

圖5 吊艙3 處迎角變化曲線Fig. 5 Angle of attack measured at the pod 3

圖6 測量點11 處速度變化曲線Fig. 6 Velocities measured at the point 11

由仿真曲線可知，在推力作用下，降階模型的縱向和垂直方向速度、角度、迎角等特性變化與原模型基本一致，且符合理論分析結(jié)果。仿真曲線中側(cè)向運動幅度較小，且受耦合作用影響，振蕩頻率較大；對比原模型和降階模型側(cè)向運動可知，雖然兩者差別較大， 但是相對整個飛行器模型影響很小。 可見，降階模型能夠保持原飛行器系統(tǒng)的運動特性，效果較好，滿足精度要求。

5 結(jié)論

本文在忽略陣風擾動的情況下，利用基于平衡實現(xiàn)的平衡殘差方法，對某高保真大柔性飛行器模型進行了模型降階，同時針對此高階模型， 運用基于平均速度的狀態(tài)空間Newmark 法求解方程[8]，給出原系統(tǒng)和降階系統(tǒng)的輸出響應分析。仿真結(jié)果表明，采用平衡殘差法能得到近似度較高的降階模型，對比原系統(tǒng)模型，降階模型能夠很好地反映飛行器系統(tǒng)的縱向和橫側(cè)向響應特性，且誤差很小，驗證了上述降階方法的可行性。 此外，仿真結(jié)果與理論分析的定性比較，從一定程度上證明了模型的正確性。 本文的結(jié)果可為下一步基于大柔性飛行器控制系統(tǒng)的設計提供參考。

[1] Shearer C M,Cesnik C E S.Nonlinear flight dynamics of very flexible aircraft[J].Journal of Aircraft,2007,5(44):1528-1545.

[2] Whitson S. The Proteus, Giving Shape to Form Unknown [J].Private Pilot, 1988-12, 33(12): 44-50.

[3] 謝長川, 吳志剛, 楊超. 大展弦比柔性機翼的氣動彈性分析[J]. 北京航空航天大學學報, 2005, 29(12): 1087-1090.XIE Chang-chuan,WU Zhi -gang,YANG Chao.Aeroelastic analysis of large aspect ratio wing [J].Journal of Beijing University of Aeronautics and Astronautics, 2005, 29(12): 1087-1090.

[4] Shearer C M. Coupled Nonlinear Aeroelasticity and Flight Dynamics of Fully Flexible Aircraft[D]. Michigan: University of Michigan,2008.

[5] Gadient R,Wise K,Lavretsky E.Very flexible aircraft control challenge problem[C]//Guidance,Navigation, and Control and Co -located Conferences, Minneapolis, Minnesota, AIAA,2012,4793:1-6.

[6] Moor B C. Principal component analysis in linear systems:Controllability, observability, and model reduction [J]. IEEE Transactions on Computer -Aided Design of Integrated Circuits and Systems, 1998,17(8): 645-654.

[7] Byung Ok Kim, An SungLee. A transient response analysis in the state-space applying the average velocity concept[J].Journal of Sound and Vibration, 2005,281: 1023-1035.

[8] 田傳艷， 胡軍照， 楊黔龍，等.滑翔航彈半實物實時仿真軟件的設計[J].現(xiàn)代應用物理， 2014（3）：239-244.TIAN Chuan-yan，HU Jun-zhao，YANG Qian-long，et al.A hardware-in-the-loop real time simulation software for the gliding missile[J].Modern Applied Physics ， 2014（3）：239-244.