Gauss偽譜法的再入可達域計算方法*

王　濤,張洪波,李永遠,湯國建

(1.國防科技大學(xué) 航天科學(xué)與工程學(xué)院， 湖南 長沙　410073； 2.中國運載火箭技術(shù)研究院 研發(fā)中心， 北京　100076)

?

Gauss偽譜法的再入可達域計算方法*

王濤1,張洪波1,李永遠2,湯國建1

(1.國防科技大學(xué) 航天科學(xué)與工程學(xué)院， 湖南 長沙410073； 2.中國運載火箭技術(shù)研究院 研發(fā)中心， 北京100076)

摘要：為了加快優(yōu)化速度和提高優(yōu)化質(zhì)量，提出一種基于Gauss偽譜法的再入可達域計算方法。鑒于再入時一般采用固定的攻角剖面，將攻角作為狀態(tài)變量，僅對傾側(cè)角進行單變量尋優(yōu)。優(yōu)化過程中，再入縱程被視為終端約束，以獲取不同縱程下的最大橫程，將速度傾角視為過程約束，以消除彈道的跳躍現(xiàn)象。通過仿真，求解出了通用航空飛行器的再入可達域，結(jié)果與間接法的理論證明一致。

關(guān)鍵詞：再入；可達域；Gauss偽譜法；攻角剖面；橫程

飛行器再入可達域是飛行器在滿足約束條件下再入落點的集合[1]。研究再入可達域?qū)υ偃肴蝿?wù)規(guī)劃及應(yīng)急迫降著陸點選擇具有重要意義。計算可達域的關(guān)鍵是確定落點區(qū)域的邊界，可將其描述為一個最優(yōu)控制問題。根據(jù)優(yōu)化方法的不同，前人的研究可以分為兩類，即直接優(yōu)化方法和間接優(yōu)化方法。

Ngo[2]以最大橫程為優(yōu)化指標，利用間接法求解傾側(cè)角的表達式，由于計算復(fù)雜，仿真中采用近似的傾側(cè)角表達式。Lu[1]假設(shè)可達域的周圍存在一系列虛擬目標點，通過使飛行器落點盡可能接近虛擬目標點，從而得到可達域的輪廓。算法基于擬平衡滑翔條件(Quasi-EquilibriumGlideCondition，QEGC)，利用最優(yōu)控制原理巧妙地將可達域計算問題轉(zhuǎn)化為單個參數(shù)求解問題。Li[3]，Zhang[4]等對此方法展開進一步研究，提出了更實用的虛擬目標選取方法。然而，采用間接法計算可達域時一般基于QEGC，實際飛行中不一定滿足此條件。

Saraf[5]基于航天飛機的制導(dǎo)模式提出了一種解析的可達域計算方法，通過規(guī)劃阻力剖面，得到近似的可達域。Chang[6]對該方法進行了改進，對末端傾側(cè)角進行了修正，從而使可達域的計算更加精確。雖然該方法計算效率較高，但無法證明可達域的最優(yōu)性。

文獻[7-9]利用直接法求解可達域。其中Fahroo[7-8]利用拉格朗日偽譜法求解再入可達域，由于約束條件處理不當，得到的軌跡并不光滑，可跟蹤性較差。文獻[10]和文獻[11]分別針對超低再入情況和沿赤道再入情況進行研究，采用偽譜法優(yōu)化可達域，然而，這兩種情況下地球自轉(zhuǎn)的影響較小，不具有一般性。

近年來，偽譜法在軌跡優(yōu)化領(lǐng)域得到了廣泛的應(yīng)用[12]。Gauss偽譜法作為偽譜法的一種，在處理約束方面有明顯優(yōu)勢。

1問題描述

1.1一般再入問題

1.1.1動力學(xué)方程

假設(shè)地球為旋轉(zhuǎn)圓球，描述飛行器運動狀態(tài)的變量包括地心距r、經(jīng)度λ、地心緯度φ、速度大小V、當?shù)厮俣葍A角θ和航跡偏航角ψ。以時間t為自變量的高超聲速滑翔飛行器運動方程為：

(1)

式中：σ為傾側(cè)角；g為地球引力加速度；而L，D分別為升力和阻力加速度大小。則有：

(3)

1.1.2約束條件

(4)

其中：k為常數(shù)，與飛行器相關(guān)；g0為海平面引力系數(shù)。

對于一般再入過程，控制量通常取攻角和傾側(cè)角，其約束為：

(5)

為了滿足能量管理段(TerminalAreaEnergyManagement，TAEM)的需要，應(yīng)對再入終端進行一定的約束?？紤]到本文著重研究再入可達域，對終端位置不作要求，終端只需滿足：

(6)

式中，μ為地球引力場數(shù)。

1.2P坐標系

如圖1所示，OE和I分別為地心和飛行器質(zhì)心位置。根據(jù)再入起點及方位角確定一個再入大圓弧平面，作為新的“赤道”平面，并定義一個地心坐標系OE-XpYpZp：Xp軸沿飛行起始點地心矢徑方向；Yp軸在再入大圓弧平面內(nèi)垂直于Xp軸指向目標方向；Zp軸與Xp，Yp軸構(gòu)成右手系，指向新的極點P，其經(jīng)緯度為(λp,φp)。在P坐標系中，起始點的經(jīng)度為零，經(jīng)度和緯度描述了再入縱程和再入橫程。利用這些特性，可極大地簡化彈道規(guī)劃算法。

圖1　P坐標系示意圖Fig.1　P coordinate system

圖1中，N為北極點，飛行器滑翔起始點I的坐標為(λ0,φ0)，方位角為ψ0。因此一般坐標系與P坐標系的轉(zhuǎn)換關(guān)系為[13]：

(7)

此外，在P坐標系中，哥氏加速度和牽連加速度也有相應(yīng)的變化，具體見文獻[13]。

1.3可達域問題

1.3.1優(yōu)化指標

再入可達域問題可描述為求解不同縱程條件下的最大橫程。在P坐標系內(nèi)，落點的橫程可用緯度表征：

J=min{cos(φpf)}

(8)

式中：φpf為P坐標系中落點的緯度，如圖2所示。

圖2　固定縱程條件下的最大橫程Fig.2　Maximum cross-range trajectories with fixed down-range

圖2中，Lset為設(shè)定的縱程，可以把Lset視為終端約束，通過改變Lset的值來得到不同縱程條件下的最大橫程彈道。λpf為落點在P坐標系中的經(jīng)度。為了提高計算效率，計算零緯度以上的最大橫程時，可以把σ的尋優(yōu)范圍縮小至[-90, 0]；針對零緯度以下的最大橫程，可以把σ的尋優(yōu)范圍縮小至[0, 90]。

1.3.2攻角剖面的處理

已知標稱攻角剖面為能量的線性分段函數(shù)。

(9)

其中，最大飛行攻角αmax=20°，最大升阻比對應(yīng)攻角αmaxL/D=10°，E1和E2為攻角曲線分段節(jié)點處的能量，可根據(jù)飛行器防熱及航程需求來確定。攻角剖面的一階微分為：

(10)

優(yōu)化過程中，把攻角當作狀態(tài)變量來處理，這樣優(yōu)化變量就只剩下傾側(cè)角，從而大大提高了尋優(yōu)的效率。式(10)與動力學(xué)方程組(1)構(gòu)成一組新的狀態(tài)方程，狀態(tài)變量α的初始值為αmax。

1.3.3速度傾角的處理

飛行器再入過程中，不希望出現(xiàn)跳躍現(xiàn)象。然而在使用Gauss偽譜法尋優(yōu)時，再入軌跡往往出現(xiàn)多次躍起現(xiàn)象，為了消除這種現(xiàn)象，對速度傾角進行以下限制：

θ≤0

(11)

式(11)與方程組(4)構(gòu)成了新的約束方程組。

2Gauss偽譜法

Gauss偽譜法以Legendre多項式的根為離散點，將狀態(tài)變量和控制變量同時離散化，從而將彈道優(yōu)化的最優(yōu)控制問題轉(zhuǎn)換為非線性規(guī)劃問題進行求解。采用Gauss偽譜法求得的非線性規(guī)劃問題的解，是間接法一階最優(yōu)性的必要條件。

設(shè)最優(yōu)控制問題的時間區(qū)間為[t0, tf]，采用Gauss偽譜法則需將時間區(qū)間轉(zhuǎn)換到[-1, 1]，因此對時間變量t作變換：

(12)

Gauss偽譜法的離散點為 N 階Legendre多項式的根，Legendre多項式為：

(13)

Legendre-Gauss點分布在區(qū)間(-1,1)，增加τ0=-1，得到區(qū)間[-1,1)，共N+1個插值點。以Lagrange插值多項式為基函數(shù)描述狀態(tài)變量和控制變量：

(14)

其中，

狀態(tài)變量的一階微分可通過對式(14)求導(dǎo)來近似，將動力學(xué)微分方程約束轉(zhuǎn)換為代數(shù)約束。

(15)

所以動力學(xué)方程滿足：

(16)

式中，n=1,…,N。

終端狀態(tài)Xf可通過拉格朗日積分得到：

U(τn),τn;t0,tf]

(17)

Gauss偽譜法中的性能指標函數(shù)為：

(18)

將連續(xù)系統(tǒng)離散化后，最優(yōu)控制問題轉(zhuǎn)化為非線性規(guī)劃問題，利用序列二次規(guī)劃算法進行尋優(yōu)，可得到需要的再入軌跡。

3仿真分析

表1　再入初始條件

3.1求解最大橫程

在固定攻角剖面的條件下，以最大橫程為優(yōu)化指標，不對再入縱程進行約束，采用Gauss偽譜法進行尋優(yōu)，得到最大橫程的彈道，如圖3、圖4所示。優(yōu)化過程中，由于對傾側(cè)角進行了限制，σ∈[-90, 0]，所以只得到單側(cè)向最大橫程的彈道。從圖4中可以看出，攻角作為狀態(tài)變量，嚴格符合式(9)。由于傾側(cè)角是唯一的優(yōu)化變量，僅需要十幾秒就找到了良好的優(yōu)化結(jié)果。

圖3　最大橫程彈道的地面軌跡Fig.3　Ground trajectory of the maximum cross-range

圖4　最大橫程彈道的控制量Fig.4　Control variables of the maximum cross-range

3.2求解可達域

圖5　可達域彈道簇Fig.5　Trajectory cluster

在固定攻角剖面的情況下，計算飛行器再入可達域，如圖5～8所示。圖5給出了不同縱程條件下的最大橫程彈道簇，彈道簇的落點所圍成的區(qū)域便是再入可達域，由于優(yōu)化過程中沒有考慮QEGC，各彈道均存在跳躍現(xiàn)象。圖6給出了彈道簇的地面軌跡，虛線為可達域的邊界，由于地球自轉(zhuǎn)的影響，可達域呈現(xiàn)非對稱的形狀。圖7給出了彈道簇的最優(yōu)傾側(cè)角變化規(guī)律：隨著時間的變化，傾側(cè)角最終趨向于零，這與文獻[1]的結(jié)果一致，文獻[1]中已給出了證明。圖8給出了彈道簇的熱流變化，所有彈道的熱流都被限制在允許的范圍內(nèi)， 圖7中指出了熱流約束對傾側(cè)角的影響。

圖6　可達域彈道簇的地面軌跡Fig.6　Ground trajectory cluster

圖7　可達域彈道簇的傾側(cè)角Fig.7　Bank angle of trajectory cluster

圖8　可達域彈道簇的熱流Fig.8　Heat flux of trajectory cluster

3.3防止彈道跳躍的可達域

圖9　可達域彈道簇的高度變化(without θ<0)Fig.9　Altitude variation of trajectory cluster (without θ<0)

圖10　可達域彈道簇的傾側(cè)角(with θ<0)Fig.10　Bank angle curve cluster (with θ<0)

圖11　可達域彈道簇的地面軌跡(with θ<0)Fig.11　Ground trajectory cluster (with θ<0)

圖12　可達域彈道簇的熱流(with θ<0)Fig.12　Heat flux of trajectory cluster (with θ<0)

限制再入速度傾角小于零，重新計算可達域，得到圖9～12。圖9給出了彈道簇(withoutθ<0)和彈道簇(withθ<0)的高度變化曲線，通過比較可知，對速度傾角約束之后，成功消除了跳躍現(xiàn)象。圖10給出了彈道簇(withθ<0)的傾側(cè)角，當飛行器進入平衡滑翔段后，為了防止跳躍，傾側(cè)角需適當增大，從而導(dǎo)致圖中的空白區(qū)域。圖11給出了彈道簇的地軌跡，與圖6相比，彈道的最大縱程明顯縮短。圖12給出了熱流的變化，熱流被限制在允許的范圍內(nèi)。

圖13給出了兩種情況下的可達域。通過比較可知，在對速度傾角進行限制后，可達域明顯縮小。

圖13　可達域的比較Fig.13　Comparison of entry footprints

4結(jié)論

采用Gauss偽譜法求解固定攻角剖面的再入可達域。求解過程中，攻角被視為狀態(tài)變量，再入縱程被視為終端約束，速度傾角被視為過程約束。經(jīng)過這些處理之后，以再入橫程為優(yōu)化指標，對傾側(cè)角進行單變量尋優(yōu)，得到了固定攻角剖面下的可達域。仿真結(jié)果與文獻[1]中間接方法的結(jié)論一致，即彈道的傾側(cè)角趨向于零。與文獻[1]中

所提方法相比，不需要大量的推導(dǎo)計算，除了Gauss偽譜法的離散誤差外，沒有對動力學(xué)模型進行任何近似。

參考文獻(References)[1]Lu P, Xue S B. Rapid generation of accurate entry landing footprints[J]. Journal of Guidance, Control, and Dynamics, 2010, 33(3): 756-765.

[2]NgoAD,DomanBD.Footprintdeterminationforreusablelaunchvehiclesexperiencingcontroleffectorfailures[C]//ProceedingsofAIAAGuidance,Navigation,andControlConferenceandExhibit,California, 2002.

[3]LiHF,ZhangR,LiZY,etal.Footprintproblemwithangleofattackoptimizationforhighliftingreentryvehicle[J].ChineseJournalofAeronautics, 2012, 25(2): 243-251.

[4]ZhangR,LiZY,ShiLN.Ageneralfootprintgenerationapproachforliftingreentryvehicle[C]//Proceedingsof28thInternationalCongressoftheAeronauticalSciences, 2012.

[5]SarafA,LeavittJA,MeaseKD.Landingfootprintcomputationforentryvehicles[C]//ProceedingsofAIAAGuidance,Navigation,andControlConferenceandExhibit,RhodeIsland, 2004.

[6]常松濤, 楊業(yè), 王永驥， 等. 計算再入飛行器可達區(qū)域的快速算法[J]. 華中科技大學(xué)學(xué)報：自然科學(xué)版, 2012, 40(7): 1-5.CHANGSongtao,YANGYe,WANGYongji,etal.Arapidalgorithmforgeneratinglandingfootprintforentryvehicles[J].JournalofHuazhongUniversityofScienceandTechnology:NatureScience, 2012, 40(7): 1-5.(inChinese)[7]FahrooF,DomanBD,NgoAD.Footprintgenerationforreusablelaunchvehiclesusingadirectpseudospectralmethod[C]//ProceedingsoftheAmericanControlConference, 2003: 2163-2168.

[8]FahrooF,DomanBD.Adirectmethodforapproachandlandingtrajectoryreshapingwithfailureeffectestimation[C]//ProceedingsofAIAAGuidance,Navigation,andControlConferenceandExhibit, 2004.

[9]BollinoK,RossIM,DomanDB.Sixdegreesoffreedomtrajectoryoptimizationforreusablelaunchvehiclefootprintdetermination[C]//ProceedingsofAAS/AIAAAstrodynamicsSpecialistsConference,Springfield, 2007.

[10]解永鋒, 唐碩. 亞軌道飛行器再入可達域快速計算方法[J]. 飛行力學(xué), 2011, 29(4): 72-75.

XIEYongfeng,TANGShuo.Rapidcalculationofentryfootprintofsuborbitallaunchvehicle[J].FlightDynamics, 2011, 29(4): 72-75.(inChinese)

[11]李柯, 聶萬勝, 馮必鳴. 助推-滑翔飛行器可達區(qū)域影響因素研究[J]. 現(xiàn)代防御技術(shù). 2013, 41(3): 42-47.

LIKe,NIEWansheng,FENGBiming.Affectingfactoroffootprintforboost-glidevehicle[J].ModernDefenceTechnology, 2013, 41(3): 42-47.(inChinese)

[12]HuntingtonGT,RaoAV.OptimalconfigurationofspacecraftformationviaaGausspseudospectralmethod[C]//ProceedingsofAASSpaceflightMechanicsMeeting, 2005.[13]XIEYu,LiuLH,LiuJ,etal.Rapidgenerationofentrytrajectorieswithwaypointandno-flyzoneconstraints[J].ActaAstronautica, 2012, 77(4): 167-181.

Landing footprint generation of entry vehicle based onGauss pseudospectral method

WANG Tao1, ZHANG Hongbo1, LI Yongyuan2, TANG Guojian1

(1.CollegeofAerospaceScienceandEngineering,NationalUniversityofDefenseTechnology,Changsha410073,China;2.Research&DevelopmentCenter,ChinaAcademyofLaunchVehicleTechnology,Beijing100076,China)

Abstract：Toreducecalculationtimeandimproveoptimizationresult,aGausspseudospectralmethodwasproposed.Fortheattackanglewasregardedasastatevariableandonlybankanglewasoptimized,whichwasreasonablesinceattackangleprofilewasoftenpredesignedinengineering.Down-rangeofvehiclewasregardedaseventconstrainttoobtainthemaximumcross-range,whiletheflightpathanglewasregardedaspathconstrainttoeliminateskipsonthetrajectory.Throughsimulation,thefootprintofCAV(commonaerovehicle)isgenerated,whichisinaccordancewiththetheoreticalanalysisoftheindirectoptimizationmethod.

Keywords：entry;footprint;Gausspseudospectralmethod;attackangleprofile;cross-range

doi:10.11887/j.cn.201603013

收稿日期：2015-03-13

基金項目：航天科技創(chuàng)新基金資助項目(CASC201306)

作者簡介：王濤(1989—)，男，河南開封人，博士研究生，E-mail：wangtao_smile@126.com； 湯國建(通信作者)，男，教授，博士，博士生導(dǎo)師，E-mail：tangguojian@nudt.edu.cn

中圖分類號：V448.2

文獻標志碼：A

文章編號：1001-2486(2016)03-075-06

http://journal.nudt.edu.cn