李 瑾,王永海,趙雅心,趙亮博

(1. 北京航天長征飛行器研究所,北京 100076;2. 中國運載火箭技術(shù)研究院,北京 100076)

1 引言

升力式飛行器返回大氣層過程中,可充分利用其氣動性能和機動能力,完成大范圍、長距離、高馬赫數(shù)下的飛行,可對敵防御系統(tǒng)進行有效突破。在飛行過程中,其會面臨復(fù)雜的外部氣動環(huán)境變化:非線性的駐點熱流條件、飛行器過載以及動壓限制等,同時也需考慮飛行路徑約束(如禁飛區(qū))。因此針對飛行器再入段的軌跡優(yōu)化近年來成為研究的熱點,并已經(jīng)形成諸多有效的方法,如直接法、間接法、精確解法和現(xiàn)代啟發(fā)式算法等。由于偽譜法[1]和凸優(yōu)化方法作為直接法中兩個比較高效的方法,其得到廣泛使用。目前主要針對提高優(yōu)化效率、縮短規(guī)劃時間對上述方法展開大量研究和改進,以實現(xiàn)可應(yīng)用于未來彈上實時軌跡優(yōu)化的解決方案[2]。

針對偽譜法,其主要包括Gauss、Legendre和Redau三類偽譜法,并且目前大部分解決方案是對上述三種方法的適應(yīng)性改進[3]。但偽譜法的本質(zhì)是一種離散和轉(zhuǎn)化方法,其最終計算性能指標(biāo)需與求解方案相結(jié)合,如大多數(shù)情況下偽譜法與序列二次規(guī)劃(SQP)一同使用以完成優(yōu)化計算,考慮到SQP問題中存在的對初值靈敏度、收斂速度等方面的不足,其仍需要進一步對算法進行改進[3]。而現(xiàn)代啟發(fā)式算法如遺傳算法、粒子群算法等雖然在迭代速度方面有著獨特的優(yōu)勢,但對再入過程的復(fù)雜約束條件會出現(xiàn)精度不高,計算后期容易陷入局部極小值等問題。

而凸優(yōu)化作為近年來逐漸興起的優(yōu)化算法,其對于初值不敏感,并且當(dāng)問題滿足凸優(yōu)化的模型時,具有局部最優(yōu)解為全局最優(yōu)解的優(yōu)勢,可有效避免陷入局部極小值的問題。并且隨著計算方法的提升,大規(guī)模凸優(yōu)化問題能夠在短時間內(nèi)求得最優(yōu)解[4],該方法也逐漸展開應(yīng)用。其率先在宇航領(lǐng)域展開應(yīng)用,美國JPL實驗室的Acikmese[5]等人率先將凸優(yōu)化方法應(yīng)用于火星著陸器動力下降段制導(dǎo)任務(wù)中,并取得很好的效果;劉新福和Lu等采用凸優(yōu)化方法對飛行器不同場景下軌跡優(yōu)化方法進行了詳細(xì)的分析。在可重復(fù)使用火箭回收[6]、空間軌道轉(zhuǎn)移[7]、地外天體著陸以及可重復(fù)使用飛行器軌跡優(yōu)化場景中多次應(yīng)用該方法[8-10],通過對非凸問題進行轉(zhuǎn)化,得到凸優(yōu)化計算的模型,高效快速完成求解計算。周祥[11]等提出指令反解方法嵌套在凸優(yōu)化模型中,使得計算精度得到提高,完成再入三維軌跡剖面的快速規(guī)劃。王振博[12]采用序列凸優(yōu)化方法,通過增加傾側(cè)角變化率作為新的控制輸入,對高非線性條件下的行星再入最優(yōu)控制問題進行了快速求解。楊奔[13]針對凸優(yōu)化方法解決再入飛行器軌跡優(yōu)化過程中出現(xiàn)的“偽不可行”問題采用b樣條曲線離散控制量并基于回溯直線搜索方法,提高數(shù)值仿真的穩(wěn)定性和光滑性。王勁博[14]、馬林[15]等結(jié)合序列凸優(yōu)化和無損凸化方法展開了系列研究并給出相應(yīng)理論證明,應(yīng)用于可重復(fù)使用火箭軌跡優(yōu)化和動力軟著陸問題中。

本文以升力式飛行器三自由度歸一化再入運動模型為基礎(chǔ),以能量e為自變量,減少運動方程數(shù)目,同時避免終端時間自由問題。采用基于可變信賴域的序列凸化和松弛方法對約束條件進行凸化處理。通過自動調(diào)節(jié)信賴域范圍和權(quán)值函數(shù)大小,以減少迭代中的振蕩問題。最后將凸問題離散化處理后,進行數(shù)值仿真計算。通過分析優(yōu)化結(jié)果,驗證算法的可行性和計算效率。

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

2.1 再入運動模型

升力式飛行器再入段沒有發(fā)動機推力作用,僅在重力和空氣動力作用力下完成飛行。為簡化優(yōu)化數(shù)目,選取能量為自變量,考慮地球自轉(zhuǎn)角速度的影響,建立再入段的三自由度無量綱運動模型如下

(1)

(2)

上式中需要注意的,速度V無量綱,密度ρ為有量綱的氣動密度,是關(guān)于r的函數(shù)ρ=ρ0exp(-(r-1)R0/hs),hs表示高度系數(shù),參考面積Sref和質(zhì)量m均為有量綱的物理量。氣動升力系數(shù)CL和阻力系數(shù)CD本文假設(shè)只與飛行馬赫數(shù)以及攻角有關(guān)。

由地球自轉(zhuǎn)引起的哥式加速度項C1、C2以及牽連加速度項C3、C4如下所示。

(3)

2.2 約束條件

終端等式約束為高度和經(jīng)緯度,飛行航跡角給定在一個范圍內(nèi)。同理航向角由于滑翔段后還將有其它飛行階段(如能量管理階段或者是末制導(dǎo)階段),其也需限制在一定范圍內(nèi)。不失一般性,這里假設(shè)終端約束條件如下所示

(4)

其飛行過程中,需滿足熱流、動壓和過載要求,其無量綱形式如下所示

(5)

除上述約束外,控制量傾側(cè)角約束如下

(6)

2.3 優(yōu)化問題描述

綜合考慮上述升力式飛行器運動模型和相關(guān)約束條件,考慮時間最短優(yōu)化指標(biāo),其軌跡優(yōu)化問題(P1)可以表示為

subto(1) (4) (5) (6)

(7)

3 凸化與離散

由于升力式飛行器軌跡優(yōu)化問題是一高度非線性的最優(yōu)控制問題,其中的性能指標(biāo)、狀態(tài)方程以及約束條件大部分難以直接利用凸優(yōu)化方法進行求解。因此需對其進行凸化處理,本文主要采用序列凸優(yōu)化方法(SCP)和松弛方法進行凸化,通過求解一系列凸優(yōu)化的子問題,從而達到逐步逼近原有問題的解。

3.1 運動方程和控制量凸化

根據(jù)上文選取確定的狀態(tài)量和控制量,將運動方程組進行線性化處理,改寫為以下形式:

(8)

其中

b(xk)=f0(xk)-A(xk)xk+fΩ(xk)

A中各項系數(shù)可通過推導(dǎo)計算求得。式中的xk為參考軌跡的參考點的狀態(tài)量。并注意矩陣A中對矢徑r求導(dǎo)的相關(guān)項中,阻力加速度D對r的導(dǎo)數(shù),Dr=-DR0/hs。由于線性化展開發(fā)生在xk處,需引入信賴域約束以確保其線性化的有效性,即

|x-xk|≤δ

(9)

針對控制量u的約束條件(6)中,第一個不等式約束滿足凸化要求,但第二條約束條件是非凸的二次等式約束,擬對其采用松弛方法將其轉(zhuǎn)化為不等式約束

(10)

3.2 約束和指標(biāo)凸化

針對過程約束,將其轉(zhuǎn)化為關(guān)于矢徑r的對數(shù)形式的約束,并進行整合,如下所示

(11)

將上述約束進行整理改寫為

(12)

針對飛行過程中,還將有飛行路徑約束(如禁飛區(qū))和其它約束,這部分將在后文中進行添加以驗證算法能力。

針對性能指標(biāo)凸化問題,通過輔助變量將待積函數(shù)轉(zhuǎn)化為線性形式,考慮優(yōu)化初期經(jīng)緯度等式約束較為嚴(yán)格,難以得到全局最優(yōu)解,以罰函數(shù)形式將其加入性能指標(biāo)當(dāng)中。同時為確?？刂屏考s束松弛帶來的最優(yōu)解與原問題保持一致,增加正則項,因此性能指標(biāo)如下所示

(13)

其中,Δθ,Δφ為待優(yōu)化變量,cθ,cφ為懲罰項因子,迭代過程進行適應(yīng)性調(diào)整。

3.3 最優(yōu)問題離散化

為了進行后續(xù)數(shù)值求解,需要將上述最優(yōu)問題進行離散化處理,將能量區(qū)間[e0,ef]進行等距劃分,得到離散的自變量E={e0,e1,e2,…,eN},其中N為離散參數(shù)。對應(yīng)的離散狀態(tài)量和離散控制量分別如下所示

X={x0,x1,x2,…,xN},U={u0,u1,u2,…,uN}

(14)

其中,xi=[ri,θi,φi,γi,ψi]T,ui=[u1,u2]T。根據(jù)上述線性化的方程以及梯形積分公式有

i=1,2,…,N

(15)

合并之后得

(16)

Hi-1,Ki,Gi相關(guān)參數(shù)由式(15)易知。將上式整理成矩陣形式Mz=Y,其選取優(yōu)化變量z為z=[x0,x1,x2,…,xN,u0,u1,u2,…,uN,Δθ,Δφ],M矩陣和Y向量如下所示

同樣,其性能指標(biāo)也需要離散化處理,主要為積分項進行離散化處理。并將其轉(zhuǎn)化為J=Lz的形式。如下所示

(17)

通過上述離散化處理,凸問題P1的離散形式如P2所示,P2問題已經(jīng)轉(zhuǎn)化為標(biāo)準(zhǔn)的具有線性性能指標(biāo)函數(shù),可以進行凸優(yōu)化求解。

(18)

3.4 改進信賴域算法

針對序列凸優(yōu)化算法,其旨在通過序列線性化方式將非凸問題轉(zhuǎn)化,但該過程中需要滿足線性化要求的信賴域指標(biāo)。固定的信賴域會在前期優(yōu)化中會出現(xiàn)較大的振蕩現(xiàn)象。為優(yōu)化該過程,設(shè)計了基于可變信賴域和權(quán)值函數(shù)的方法如下所示,使迭代過程中更加接近優(yōu)化值。

δ=-Kδ0

K=-ae-bk

(19)

其中,k為迭代次數(shù),a,b分別為設(shè)置參數(shù),可根據(jù)需求選取合適數(shù)值,本文中設(shè)置a=1.19,b=-0.36。同理性能指標(biāo)中的懲罰函數(shù)因子cθ和cφ也按照上述方法進行設(shè)置。

程序流程圖如下所示。

圖1 優(yōu)化流程圖

4 數(shù)值仿真

4.1 問題描述

本文中,擬采用美國某一可重復(fù)使用飛行器(RLV)的飛行模型[16]。其質(zhì)量為m=104305kg,參考面積為S=391.22m2,其氣動模型擬合函數(shù)如下所示

(20)

圖2 氣動力系數(shù)與速度之間的關(guān)系

攻角方案采用給定的分段函數(shù)式剖面,其具體形式如下所示

(21)

相關(guān)設(shè)定參數(shù)如下表1所示。

表1 初始相關(guān)參數(shù)

仿真中為確保序列線性化的有效性和最終結(jié)果的準(zhǔn)確性,下面給出信賴域和收斂范圍。

(22)

根據(jù)上述的迭代流程,由于該方法求解凸優(yōu)化問題對于初始值不敏感。結(jié)合相關(guān)研究,確定第一步參與計算的軌跡狀態(tài)量為初始狀態(tài)量與終端狀態(tài)的線性插值得到。

4.2 計算結(jié)果分析

本文的計算仿真采用MATLAB平臺下的Mosek求解器進行求解凸優(yōu)化問題,其可大大提高計算的效率。運行計算機搭載IntelCorei5-6500 3.20GHz處理器。每次優(yōu)化CPU的平均時間為0.174s。根據(jù)飛行結(jié)果,其飛行時間為1374.05秒,飛行距離為5824.06千米,下面對部分飛行參數(shù)結(jié)果進行展示分析。

圖3中為飛行高度的變化情況,結(jié)合下圖8的飛行軌跡圖可知其優(yōu)化曲線滿足約束要求,并且在逐步迭代過程中曲線逐漸趨于最優(yōu)解。分析下述飛行終點與設(shè)置值大小如下表2所示,可知終端高度和航向角滿足等式約束,而經(jīng)緯度為滿足優(yōu)化初期結(jié)果正確性,于性能指標(biāo)中增加了懲罰因子,使得終端緯度與目標(biāo)值有一定誤差。速度傾角為給定的區(qū)間范圍,其數(shù)值滿足終端條件。

表2 規(guī)劃終端結(jié)果表

圖3 迭代飛行軌跡變化情況

由飛行控制量傾側(cè)角的變化情況可知,其數(shù)值變化主要集中于10度上下的區(qū)間范圍,由于飛行過載的限制,整個飛行軌跡為一個大半徑的圓弧,其形式與傾側(cè)角變化過程正好對應(yīng)。除在接近終端約束過程中需要調(diào)節(jié)速度傾角和終端位置坐標(biāo),因此傾側(cè)角會出現(xiàn)一個幅度較大的變化,其滿足飛行要求。如圖4。

圖4 傾側(cè)角變化情況

由圖5,6所示,速度傾角在滑翔飛行過程中始終保持在0度附近,滿足滑翔飛行要求;而航向角由于初始值和終端值的設(shè)定,其需要完成一定范圍的角度變化。并且同飛行高度變化曲線和控制量變化結(jié)合來看,其速度傾角在滑翔飛行段由于高度變化存在小幅度的跳動現(xiàn)象,并且在飛行終端,為滿足終端硬性約束條件,傾側(cè)角和速度傾角同步發(fā)生一定的振蕩。飛行高度與速度的變化曲線,如圖8所示,其過程滿足速度變化要求。

圖5 速度傾角變化情況

圖6 飛行航向角變化情況

圖7 飛行高度隨速度的變化情況

圖8 地面軌跡變化曲線

考慮禁飛區(qū)的路徑約束如下所示

(θ-θz)2+(φ-φz)2≥Rz

(23)

其中,θz、φz分別為禁飛區(qū)所在圓的經(jīng)緯度,Rz為半徑,易知該約束不滿足凸優(yōu)化要求,再此對其進行松弛凸化處理。

2(φk-φz)φk

(24)

本文設(shè)立兩處禁飛區(qū)限制,其具體位置如下表3所示。其禁飛區(qū)下的飛行軌跡如圖9 所示?？梢娊?jīng)過迭代計算得到的飛行軌跡可以有效地躲避禁飛區(qū),并且仍能滿足相關(guān)約束限制。

表3 禁飛區(qū)位置

圖9 飛行軌跡(含禁飛區(qū))經(jīng)度緯度及高度變化規(guī)律

通過仿真計算,如下圖10可得在增加可變信賴域之后,曲線的振蕩幅度明顯變小,可以有效避免在迭代過程中出現(xiàn)大幅度抖振而帶來的仿真優(yōu)化失敗。而且在權(quán)值函數(shù)的作用下,優(yōu)化終點的參數(shù)更加趨向于終端約束。

圖10 飛行軌跡局部圖

5 結(jié)論

本文以升力式飛行器為研究對象,針對其再入段進行軌跡快速優(yōu)化計算,并通過可變信賴域的改進序列凸優(yōu)化算法,有效減少優(yōu)化過程中的抖振現(xiàn)象,使得優(yōu)化結(jié)果更加準(zhǔn)確。并且該算法相較其它方法,其計算耗時短,對初值的要求不高,可在未知飛行軌跡情況下進行大致估計而不影響最終結(jié)果,因此其未來將具備很好的應(yīng)用價值。