侯忻宜，李惠峰

(北京航空航天大學(xué)宇航學(xué)院，北京 100191)

0 引言

吸氣式高速飛行器一般由渦輪或沖壓發(fā)動(dòng)機(jī)提供動(dòng)力，可充分利用大氣層內(nèi)的氧氣進(jìn)行高效率巡航飛行。吸氣式高速飛行器要經(jīng)歷長時(shí)間遠(yuǎn)距離的高速飛行，因此需要合理地設(shè)計(jì)優(yōu)化軌跡，以滿足不同任務(wù)形式的需要。在線自主軌跡優(yōu)化有利于提高飛行器的自適應(yīng)和自主決策能力。

目前針對(duì)吸氣式高速飛行器的軌跡規(guī)劃研究分別集中在爬升段和巡航段。楊志紅等針對(duì)吸氣式高速飛行器爬升軌跡規(guī)劃問題，采用高斯偽譜法離散，并用序列二次規(guī)劃方法求解問題。Zhou等、龔春林等、趙長見等針對(duì)吸氣式推進(jìn)和火箭推進(jìn)的組合動(dòng)力運(yùn)載器的上升段軌跡規(guī)劃問題進(jìn)行了研究。Liu等采用粒子群和序列二次規(guī)劃算法求解最優(yōu)穩(wěn)態(tài)巡航條件。Kang等采用遺傳算法和直接打靶法解決了一種考慮燃油消耗和熱流率最小的多目標(biāo)周期巡航軌跡優(yōu)化問題。Gao 等通過擬合特定周期函數(shù)的參數(shù)，將偽譜法用于周期巡航軌跡規(guī)劃問題中。相關(guān)研究主要集中在縱向穩(wěn)態(tài)巡航和周期巡航這種固定飛行模式的優(yōu)化問題中，但是當(dāng)軌跡規(guī)劃任務(wù)涉及目標(biāo)到達(dá)和禁飛區(qū)規(guī)避時(shí)，這種固定的飛行形式限制了飛行器性能的發(fā)揮，并且難以在線應(yīng)用。對(duì)于無動(dòng)力再入的情況，Liu等、Wang等、Wang等采用基于凸優(yōu)化的方法對(duì)此進(jìn)行了研究，但是并未考慮吸氣式飛行器的推力與狀態(tài)量之間的相互影響。

本文針對(duì)飛行器在三維空間中的運(yùn)動(dòng)和目標(biāo)到達(dá)任務(wù)，將原最優(yōu)控制問題凸化和離散化，轉(zhuǎn)化為一系列凸問題進(jìn)行求解，設(shè)計(jì)算法的外環(huán)迭代策略。

數(shù)值仿真實(shí)驗(yàn)表明，該方法具有可行性，并且相比于非線性規(guī)劃算法，在性能指標(biāo)近似相等的情況下可以大大縮短優(yōu)化時(shí)間，證明了其針對(duì)吸氣式高速飛行器的軌跡優(yōu)化問題的最優(yōu)性和在線應(yīng)用的能力。

1 優(yōu)化問題模型

1.1 任務(wù)描述

對(duì)于采用沖壓發(fā)動(dòng)機(jī)作為動(dòng)力的吸氣式高速飛行器，一般需要經(jīng)歷機(jī)載發(fā)射和火箭助推段，助推段結(jié)束后關(guān)機(jī)分離，飛行器進(jìn)入超燃沖壓發(fā)動(dòng)機(jī)工作窗口。傳統(tǒng)的軌跡設(shè)計(jì)思路是引導(dǎo)飛行器爬升到固定的高度，然后進(jìn)行穩(wěn)態(tài)巡航飛行。這種軌跡形式需要進(jìn)行一定的約束，并且巡航段的飛行模式較為固定，性能難以達(dá)到最優(yōu)。針對(duì)這種情況，設(shè)計(jì)了如圖1所示的吸氣式高速飛行器的爬升-巡航飛行任務(wù)剖面，將飛行器的高度和速度約束放寬，在發(fā)動(dòng)機(jī)工作窗口內(nèi)尋找最優(yōu)解，有利于最大限度發(fā)揮飛行器的能力，得到最優(yōu)飛行軌跡。優(yōu)化階段從發(fā)動(dòng)機(jī)點(diǎn)火開始，到達(dá)目標(biāo)上空時(shí)結(jié)束。

圖1 吸氣式高速飛行器任務(wù)剖面圖Fig.1 Mission profile of the air breathing high velocity vehicle

1.2 數(shù)學(xué)模型建立

本文研究飛行器的質(zhì)心運(yùn)動(dòng)模型，將地球視為均質(zhì)圓球，考慮地球曲率的影響，建立無量綱化后的三自由度運(yùn)動(dòng)方程如下

sinsinsin)

2(tansincos-sin)-

(1)

升力和阻力用海平面重力加速度無量綱化，具體表達(dá)式如下

(2)

式中，升力系數(shù)和阻力系數(shù)是關(guān)于攻角和馬赫數(shù)的函數(shù)，為飛行器的參考面積，取值為1 m。為大氣密度，本文采用“準(zhǔn)等溫”大氣模型，為海平面的大氣密度，為基準(zhǔn)高度，取值為7.11 km，大氣密度可表示為高度的函數(shù)

(3)

沖壓發(fā)動(dòng)機(jī)的推力和燃料秒耗量一般與飛行器的高度、馬赫數(shù)、攻角和燃油當(dāng)量比有關(guān)，可以表示為=(,,,)，=(,,,)，根據(jù)飛行器和發(fā)動(dòng)機(jī)模型進(jìn)行插值或者擬合得到。

1.3 最優(yōu)控制問題建模

1.2節(jié)給出了飛行器的運(yùn)動(dòng)方程以及相應(yīng)的狀態(tài)量和控制量。本節(jié)將建立相應(yīng)的最優(yōu)控制問題模型。飛行器的狀態(tài)變量選取為=[,,,,,,]，控制變量選取為=[,,]。

1)性能指標(biāo)：在飛行器總航程固定的情況下，優(yōu)化指標(biāo)設(shè)置為飛行器到目標(biāo)時(shí)的終端質(zhì)量最大

Minimize:=-

(4)

2)運(yùn)動(dòng)方程：1.2節(jié)式(1)所示的運(yùn)動(dòng)方程，用狀態(tài)量和控制量寫為緊湊形式

(5)

約束條件：優(yōu)化的初始狀態(tài)為助推段結(jié)束時(shí)的狀態(tài)，此時(shí)飛行器進(jìn)入沖壓發(fā)動(dòng)機(jī)開機(jī)窗口。初始狀態(tài)=[,,,,,,]中的參數(shù)也相應(yīng)進(jìn)行無量綱化處理，約束如下

()=

(6)

式中,=+表示時(shí)刻飛行器的地心距，和是飛行器初始時(shí)刻所處的經(jīng)度和緯度，，和分別為時(shí)刻的速度大小、飛行路徑角和航向角，為飛行器在時(shí)刻的質(zhì)量。

飛行器在巡航段末端需要到達(dá)目標(biāo)點(diǎn)，用經(jīng)緯度表示為(,)，終端時(shí)刻處的經(jīng)度和緯度為()和()，因此終端時(shí)刻的約束條件為

()=,()=

(7)

受沖壓發(fā)動(dòng)機(jī)開機(jī)條件限制，飛行器的攻角需要保持在一定范圍內(nèi)，和分別為攻角能夠取得的最小值和最大值。發(fā)動(dòng)機(jī)的燃油當(dāng)量比∈[,]，本文中的范圍取為[02，1]。和分別為地心距的下限和上限。和為維持發(fā)動(dòng)機(jī)工作所需要的最小速度和最大速度。為飛行器初始質(zhì)量，表示飛行器的凈質(zhì)量。狀態(tài)量和控制量的不等式約束形式如(8)所示

≤≤，≤≤，
≤≤，≤≤，
≤≤

(8)

一般來講，飛行器在巡航飛行段可能會(huì)遇到禁飛區(qū)，用經(jīng)度和緯度表示禁飛區(qū)的圓心位置，表示半徑，禁飛區(qū)不等式約束如下

(-)+(-)≥

(9)

式中,和分別為禁飛區(qū)圓心的經(jīng)度和緯度，為禁飛區(qū)半徑，以弧度為單位。

綜上所述，飛行器巡航段軌跡優(yōu)化問題可建模為如下形式的最優(yōu)控制問題Problem1。

Problem1：

Minimize:(4)

Subject to:(5)(6)(7)(8)(9)

2 軌跡優(yōu)化算法

2.1 最優(yōu)控制問題線性化

(10)

式中

(11)

在性能指標(biāo)中加入控制量鄰近項(xiàng)

(12)

=+

(13)

(14)

可以看出，式(12)經(jīng)過松弛后，轉(zhuǎn)換為如式(13)和(14)所示的線性形式的性能指標(biāo)和不等式約束，相應(yīng)地得到Problem2。

Problem2:

Minimize:(13)

Subject to:(10)(6)(7)(8)(11)(14)

2.2 離散化

在最優(yōu)控制問題Problem2中，狀態(tài)變量和控制變量均為關(guān)于時(shí)間的連續(xù)函數(shù)，采用一階保持器離散方法將狀態(tài)量和控制量進(jìn)行離散化處理。

取離散點(diǎn)數(shù)目為，初始時(shí)間為，終端時(shí)間為，間隔為Δ=(-)(-1)，相應(yīng)的離散時(shí)間點(diǎn)為=+(-1)Δ,=1,2,…,。將控制輸入?yún)?shù)化為樣本點(diǎn)之間的線性函數(shù)

(15)

經(jīng)過狀態(tài)轉(zhuǎn)移離散后的系統(tǒng)為

+1=+1+2+1+,
=1,2,…,-1

(16)

其中系數(shù)矩陣為

(17)

式(10)的線性微分方程轉(zhuǎn)化為式(16)所示的線性等式約束，其中待優(yōu)化的變量只有和。

對(duì)于存在禁飛區(qū)的情況，將式(11)表示的約束離散化，得到線性不等式約束

(18)

然后，對(duì)式(14)也進(jìn)行離散化處理

(19)

離散后的=[1,2,…,]，式(19)可進(jìn)一步寫為

(20)

將式(20)轉(zhuǎn)化為標(biāo)準(zhǔn)的二階錐約束形式

(21)

Problem3是一個(gè)標(biāo)準(zhǔn)的二階錐規(guī)劃問題，它屬于迭代算法中的一個(gè)子問題，優(yōu)化變量包括狀態(tài)量={,,,,,,}和控制量={,,},=1,2,…,以及參數(shù)，=1,2,3。

Problem3:

Minimize:-+++

Subject to:

=,=,=

≤≤,≤≤,
≤≤，≤≤,
≤≤,=1,2,…,

(22)

2.3 序列凸優(yōu)化迭代算法

序列凸優(yōu)化算法偽代碼見表1，選取離散形式的兩次迭代地心距之差最大值作為誤差，當(dāng)兩次迭代結(jié)果之間的誤差小于規(guī)定值時(shí)，可認(rèn)為軌跡收斂。

表1 算法偽代碼Tab.1 Algorithm pseudocode

序列凸優(yōu)化算法的流程如圖2所示。

圖2 算法流程圖Fig.2 Algorithm flow chart

3 數(shù)值仿真算例

針對(duì)某型吸氣式高速飛行器進(jìn)行巡航段的軌跡優(yōu)化數(shù)值仿真實(shí)驗(yàn)，優(yōu)化階段包括沖壓發(fā)動(dòng)機(jī)從開機(jī)到關(guān)機(jī)的整個(gè)飛行過程，飛行時(shí)間為920 s。

飛行器的初始高度為20 km，速度為1 200 m/s，質(zhì)量為5 000 kg，初始經(jīng)緯度和飛行路徑角分別為0°，航向角為30°。終端經(jīng)度和緯度分別為9.5°和5.5°。飛行過程中高度的取值范圍為20～30 km，速度范圍為1 200～1 800 m/s，攻角范圍為0°～5°，燃油當(dāng)量比范圍為0.2～1，飛行器凈質(zhì)量為4 000 kg。

仿真中采用的參數(shù)如表2所示。經(jīng)緯度初始參考選為初始位置和終端位置兩點(diǎn)連線。質(zhì)量的初始參考設(shè)計(jì)為從飛行器總質(zhì)量到凈質(zhì)量的線性變化形式，燃油當(dāng)量比的初始參考也采用從到的線性變化形式。其他狀態(tài)量和控制量的初始參考無特殊要求，一般為初始時(shí)刻的常值即可。

表2 仿真參數(shù)值Tab.2 Simulation parameters

采用MATLAB軟件進(jìn)行數(shù)值仿真實(shí)驗(yàn)，實(shí)驗(yàn)用筆記本電腦處理器型號(hào)為i5-10210U，主頻為1.60 GHZ，每次迭代的二階錐規(guī)劃問題采用ECOS求解器進(jìn)行求解，當(dāng)最大位置誤差小于時(shí)迭代結(jié)束。軌跡優(yōu)化算法經(jīng)過19次迭代后，相鄰兩次迭代的位置誤差小于規(guī)定值，算法終止，仿真耗時(shí)0.065 4 s，滿足在線軌跡優(yōu)化的要求。使用GPOPS工具箱的優(yōu)化結(jié)果作為最優(yōu)參考，優(yōu)化時(shí)間為12.713 4 s，共進(jìn)行了3次網(wǎng)格自適應(yīng)迭代，網(wǎng)格劃分精度為10。

圖3和圖4分別為狀態(tài)量和控制量曲線的收斂結(jié)果與最優(yōu)參考的對(duì)比，可以看出飛行器完成了目標(biāo)到達(dá)任務(wù)，軌跡形式與預(yù)期一致。在20 km高處，發(fā)動(dòng)機(jī)點(diǎn)火，飛行器開始爬升，在26 km附近進(jìn)行波動(dòng)飛行，受攻角條件和發(fā)動(dòng)機(jī)推力約束限制，飛行高度在24～30 km范圍內(nèi)變化，接近目標(biāo)上空時(shí)高度下降至23 km。優(yōu)化結(jié)果曲線和GPOPS結(jié)果曲線變化趨勢一致，收斂至近似局部最優(yōu)解，仿真結(jié)果證明了算法的最優(yōu)性。

(a) 高度-時(shí)間曲線

(b) 速度-時(shí)間曲線

(c) 質(zhì)量-時(shí)間曲線圖3 狀態(tài)量對(duì)比圖Fig.3 Comparison of state

(a)攻角-時(shí)間曲線

(b) 燃油當(dāng)量比-時(shí)間曲線

(c)推力-時(shí)間曲線圖4 控制量對(duì)比圖Fig.4 Comparison of control

從圖3中可以看出，飛行器的高度和速度均滿足約束，處于沖壓發(fā)動(dòng)機(jī)開機(jī)窗口范圍內(nèi)。質(zhì)量變化曲線圖顯示，終端燃料剩余142 kg，GPOPS優(yōu)化結(jié)果為144 kg，相差只有2 kg，二者優(yōu)化得到的性能指標(biāo)差異約為1%。

從圖4可以看出，飛行器的攻角約束滿足，在2°和5°之間變化。燃油當(dāng)量比隨時(shí)間增加逐漸從1變?yōu)?.2，與初始參考的變化趨勢一致，相應(yīng)的推力變化范圍為0～30 kN，隨節(jié)流系數(shù)的增加，推力總體呈下降趨勢。

圖5為終端質(zhì)量性能指標(biāo)隨迭代次數(shù)的變化曲線圖，在經(jīng)過19次迭代后，終端質(zhì)量達(dá)到4 142 kg，此時(shí)迭代結(jié)束。圖6顯示了相鄰兩次迭代的位置和速度的誤差絕對(duì)值的最大值隨迭代次數(shù)變化的曲線。位置和速度誤差在迭代初期變化較大，后期逐漸減小，最后一次迭代和上一次迭代之間的差值分別為7.6 m和0.11 m/s，證明軌跡收斂。圖7顯示了相鄰兩次迭代的攻角和燃油當(dāng)量比的誤差絕對(duì)值的最大值隨迭代次數(shù)變化的曲線。攻角誤差曲線在前期有小范圍的波動(dòng)，11次迭代后小于10。燃油當(dāng)量比誤差曲線在2～12次迭代之間略有波動(dòng)，與其他變量相比收斂難度更大，但12次迭代后趨于收斂。仿真結(jié)果表明軌跡收斂效果良好。

圖5 性能指標(biāo)隨迭代次數(shù)變化圖Fig.5 Value of the performance index at each step

圖6 位置和速度誤差隨迭代次數(shù)變化圖Fig.6 Position and velocity error at each step

圖7 攻角和燃油當(dāng)量比誤差隨迭代次數(shù)變化圖Fig.7 AoA and throttling coefficient error at each step

作為對(duì)比，采用GPOPS進(jìn)行多階段軌跡優(yōu)化，將不同工況下采用爬升+穩(wěn)態(tài)巡航方式的優(yōu)化結(jié)果列在表3中。表格中的數(shù)據(jù)分別表示不同巡航高度和速度下，飛行器到達(dá)目標(biāo)時(shí)剩余燃料的質(zhì)量。表3中巡航高度為26 km時(shí)，最大的剩余燃料質(zhì)量為65 kg，相比于圖2的終端剩余燃料142 kg，爬升+穩(wěn)態(tài)巡航的飛行方式多消耗燃料約7.7%。仿真結(jié)果證明，本文提出的爬升巡航優(yōu)化軌跡比穩(wěn)態(tài)巡航軌跡更省燃料，飛行器具有更大的能力應(yīng)對(duì)突發(fā)狀況或進(jìn)行后續(xù)機(jī)動(dòng)。

表3 穩(wěn)態(tài)巡航優(yōu)化結(jié)果Tab.3 Steady-state cruise optimization results

對(duì)于存在禁飛區(qū)約束的情況，算法步驟同上。本文設(shè)置3個(gè)禁飛區(qū)，在經(jīng)緯度圖中表示禁飛區(qū)中心和半徑，如表4所示。

表4 禁飛區(qū)中心和半徑Tab.4 Center and radius of the no-fly zone

圖8為優(yōu)化得到的經(jīng)緯度圖，圖中黑色虛線表示初始參考軌跡，紅色點(diǎn)劃線表示算法迭代兩次后得到的軌跡，藍(lán)色實(shí)線表示最終收斂結(jié)果。仿真結(jié)果表明，飛行器能夠完成禁飛區(qū)規(guī)避任務(wù)，到達(dá)目標(biāo)點(diǎn)，算法求解時(shí)間0.088 5 s，滿足巡航段在線軌跡優(yōu)化要求。

圖8 經(jīng)度-緯度圖Fig.8 Longitude-latitude

4 結(jié)論

1)對(duì)吸氣式高速飛行器進(jìn)行了爬升巡航段在線軌跡優(yōu)化，將爬升段和巡航段作為整體進(jìn)行優(yōu)化，擴(kuò)大了飛行器的高度和速度約束范圍。爬升巡航整體優(yōu)化得到的軌跡燃料消耗更小，飛行器具有更大的能力應(yīng)對(duì)突發(fā)狀況或進(jìn)行后續(xù)機(jī)動(dòng)，避免了傳統(tǒng)的定高定速巡航的軌跡形式具有的局限性。

2)設(shè)計(jì)基于凸優(yōu)化的序列優(yōu)化算法在線求解軌跡優(yōu)化問題，將原最優(yōu)控制問題在參考軌跡處進(jìn)行凸化和離散化后，轉(zhuǎn)化為標(biāo)準(zhǔn)的二階錐規(guī)劃問題，基于內(nèi)點(diǎn)法進(jìn)行求解。數(shù)值仿真實(shí)驗(yàn)表明，該算法的求解時(shí)間滿足在線軌跡優(yōu)化的要求。算法求解得到的軌跡與最優(yōu)參考接近，性能指標(biāo)相差僅為1%，驗(yàn)證了算法的最優(yōu)性。