基于平衡理論的大柔性飛行器模型降階

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

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

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

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

1 大柔性飛行器模型

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

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

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

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

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

圖1 “禿鷹”飛行器結構及坐標系圖Fig. 1 Vehicle structure and coordinate systems for vulture

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

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

2 模型降階

本文主要采用基于平衡實現的平衡殘差降階方法，來簡化系統模型。

2.1 平衡實現和平衡殘差法

平衡實現理論主要由Moore[6]提出。 考慮如下線性時不變系統：

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

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

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

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

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

則模型降階為

2.2 大柔性飛行器模型降階

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

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

圖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 狀態空間Newmark 法求解方程

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

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

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

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

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

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

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

則時刻的速度微分方程為

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

由于本文所研究模型已經是一階線性模型， 故狀態向量即為，則

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

4 仿真分析

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

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

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

圖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

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

5 結論

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

[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].現代應用物理， 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.