趙吉松，張建宏，李 爽

(1. 南京航空航天大學(xué)航天學(xué)院，南京 210016；2. 北京空天技術(shù)研究所，北京 100074)

0 引 言

高超聲速滑翔飛行器在空天返回、遠(yuǎn)程快速打擊等領(lǐng)域具有重要的應(yīng)用背景。由于自身的氣動布局以及惡劣飛行環(huán)境的影響，高超聲速飛行器動力學(xué)模型具有高度非線性、強(qiáng)耦合及參數(shù)不確定等特點(diǎn)，因而其制導(dǎo)控制系統(tǒng)設(shè)計(jì)具有挑戰(zhàn)性[1-5]。軌跡優(yōu)化對于制導(dǎo)控制系統(tǒng)設(shè)計(jì)乃至總體設(shè)計(jì)具有重要意義。高超聲速滑翔軌跡優(yōu)化屬于最優(yōu)控制問題，其求解方法主要分為間接法和直接法[6]。間接法借助變分法或最大值原理，把泛函優(yōu)化轉(zhuǎn)化為兩點(diǎn)邊值問題求解；直接法通過對控制變量和/或狀態(tài)變量進(jìn)行離散把泛函優(yōu)化轉(zhuǎn)化為非線性規(guī)劃(Nonlinear programming, NLP)求解。直接法中的配點(diǎn)法[6]由于不需要推導(dǎo)最優(yōu)性必要條件，并且對初值的敏感性較低，容易收斂，近年來得到廣泛研究。

高超聲速滑翔飛行器再入軌跡優(yōu)化問題[6]是復(fù)雜的多約束、非線性最優(yōu)控制問題，優(yōu)化計(jì)算量通常比較大。目前，以序列二次規(guī)劃(Sequential quadratic programming，SQP)算法為代表的NLP求解器需要提供NLP的目標(biāo)函數(shù)和約束對自變量的一階偏導(dǎo)數(shù)甚至二階偏導(dǎo)數(shù)。其中，基于一階偏導(dǎo)數(shù)的NLP求解器只需要提供一階偏導(dǎo)數(shù)，然后采用擬牛頓法構(gòu)造近似的二階偏導(dǎo)數(shù)，比如SNOPT[7]；基于二階偏導(dǎo)數(shù)的NLP求解器除了需要提供一階偏導(dǎo)數(shù)，還需要提供二階偏導(dǎo)數(shù)，比如IPOPT[8]。文獻(xiàn)[9-12]研究發(fā)現(xiàn)，與直接計(jì)算NLP的偏導(dǎo)數(shù)相比，利用NLP的偏導(dǎo)數(shù)的稀疏性并且將其轉(zhuǎn)化為軌跡優(yōu)化問題對自變量的偏導(dǎo)數(shù)能夠顯著提高軌跡優(yōu)化的效率。考慮到實(shí)際飛行器的氣動數(shù)據(jù)通常為離散形式，本文研究采用文獻(xiàn)[11]給出的稀疏差分方法，通過計(jì)算軌跡優(yōu)化問題的一階偏導(dǎo)數(shù)高效計(jì)算NLP的一階偏導(dǎo)數(shù)，從而提高軌跡優(yōu)化的效率。

由于受到熱流、過載和動壓等多種約束，高超聲速飛行器再入軌跡通常是不光滑的。對于非光滑軌跡，如果需要高精度求解軌跡優(yōu)化問題，均勻加密不是一種較好的方案。因?yàn)檫@樣處理會導(dǎo)致計(jì)算量顯著增加，需要更多的計(jì)算資源，并且隨著節(jié)點(diǎn)數(shù)目的增加，算法的收斂性通常變差[13]。這種情況下有必要引入網(wǎng)格細(xì)化技術(shù)對節(jié)點(diǎn)進(jìn)行局部加密或者調(diào)整，在軌跡變化平緩區(qū)域采用稀疏網(wǎng)格，在軌跡變化劇烈區(qū)域采用稠密網(wǎng)格。目前，國內(nèi)外研究者在網(wǎng)格細(xì)化方面已做了不少工作[13-28]。對于局部配點(diǎn)法，由于離散格式本身對于離散節(jié)點(diǎn)分布沒有限制，因而局部配點(diǎn)法的離散節(jié)點(diǎn)可以根據(jù)需要任意布置。適用于局部配點(diǎn)法的網(wǎng)格細(xì)化方法主要有小波分析法[13-14]、離散誤差分析法[15]、多分辨率技術(shù)[16-21]、密度函數(shù)法[22]等。全局配點(diǎn)法(又稱偽譜法或者正交配點(diǎn)法)的離散節(jié)點(diǎn)是正交多項(xiàng)式的根，其節(jié)點(diǎn)分布特點(diǎn)是中間稀兩端密。全局配點(diǎn)法的優(yōu)勢是對于光滑問題具有非常高的代數(shù)精度，但是其離散節(jié)點(diǎn)不可以隨意布置。因此，全局配點(diǎn)法在網(wǎng)格細(xì)化方面沒有局部配點(diǎn)法靈活。全局配點(diǎn)法的網(wǎng)格細(xì)化方法通常是通過添加結(jié)點(diǎn)(段與段之間的連接點(diǎn)稱為結(jié)點(diǎn))將軌跡分段優(yōu)化，利用結(jié)點(diǎn)附近的離散節(jié)點(diǎn)較密的特點(diǎn)達(dá)到加密網(wǎng)格的效果，比如結(jié)點(diǎn)偽譜法[23]以及各種自適應(yīng)偽譜法[25-28]。分段優(yōu)化的主要問題是通常很難確定需要分多少段以及在哪分段。若通過優(yōu)化迭代確定這些信息會導(dǎo)致分段過多、采用的離散節(jié)點(diǎn)較多從而增加計(jì)算量等弊端[23-28]。為了提高網(wǎng)格選取的靈活性，文獻(xiàn)[29]提出了一種適應(yīng)于任意網(wǎng)格的偽譜法，但是隨著非高斯網(wǎng)格區(qū)域的增加，方法的精度顯著降低。

在眾多網(wǎng)格細(xì)化算法中，多分辨率技術(shù)[16-21]具有吸引力。多分辨率技術(shù)的思想是分析離散節(jié)點(diǎn)的插值誤差，如果某個(gè)節(jié)點(diǎn)的插值誤差較大，則在該節(jié)點(diǎn)附近細(xì)化網(wǎng)格，否則不進(jìn)行細(xì)化。多分辨率技術(shù)原理簡單、魯棒性好，并且采用的離散節(jié)點(diǎn)較少。文獻(xiàn)[16-17]給出的多分辨率方法基于二進(jìn)網(wǎng)格，其初始網(wǎng)格節(jié)點(diǎn)數(shù)必須是2j+1(j為整數(shù))，應(yīng)用受到限制。文獻(xiàn)[19-20]定義了一種節(jié)點(diǎn)數(shù)不受限制的廣義二分網(wǎng)格，但是給出的網(wǎng)格細(xì)化算法要求初始網(wǎng)格節(jié)點(diǎn)數(shù)目必須是最粗糙網(wǎng)格的二倍，使用仍然不夠方便。最近，文獻(xiàn)[21]提出了一種新的廣義二分網(wǎng)格，克服了這一缺陷。將這種新型網(wǎng)格與多分辨率技術(shù)相結(jié)合可構(gòu)造出相應(yīng)的網(wǎng)格細(xì)化算法。

本文研究采用稀疏差分法和基于新型廣義二分網(wǎng)格的網(wǎng)格細(xì)化算法快速、高精度求解高超聲速滑翔飛行器再入軌跡優(yōu)化問題，并且通過優(yōu)化落點(diǎn)區(qū)邊界軌跡研究了再入飛行器的落點(diǎn)區(qū)范圍。

1 再入軌跡運(yùn)動模型

將飛行器簡化為質(zhì)點(diǎn)，記狀態(tài)向量x=[r,θ,φ,v,ψ,γ]T分別表示飛行器的矢徑、經(jīng)度、緯度、速度、航向角和航跡角，那么描述飛行器質(zhì)心運(yùn)動的微分方程組為(忽略地球自轉(zhuǎn)的影響)[6]

(1)

式中：g為飛行器受到的地球引力加速度，g=μ/r2;μ為地球引力常數(shù)。

空氣動力產(chǎn)生的加速度沿飛行軌跡的切向、法向和側(cè)向分量as，an，aw分別為

(2)

式中：σ為速度傾側(cè)角；m為飛行器的質(zhì)量。

升力和阻力的表達(dá)式分別為

(3)

式中：ρ為大氣密度，采用美國標(biāo)準(zhǔn)大氣模型(1976年)[30]插值計(jì)算；A為參考面積；CL和CD分別為升、阻力系數(shù)，Ma和α分別為馬赫數(shù)和攻角。

圖1給出升力系數(shù)和阻力系數(shù)隨馬赫數(shù)和攻角的變化特性[31]。與解析形式的氣動數(shù)據(jù)相比，離散形式的氣動數(shù)據(jù)更為常見，更具有代表性。

圖1 氣動力系數(shù)隨馬赫數(shù)和攻角變化規(guī)律Fig.1 Variation of aerodynamic coefficients with Mach number and angle of attack

2 再入軌跡優(yōu)化問題

2.1 初始條件

高超聲速滑翔飛行器的再入初始條件為

(4)

式中：h為飛行高度，h=r-Re，Re為地球半徑。

2.2 路徑約束

再入軌跡優(yōu)化問題的控制變量為u=[α(t),σ(t)]T，其變化范圍受到限制，即

(5)

式中：αmin和αmax分別為攻角的下邊界和上邊界；σmin和σmax分別為傾側(cè)角的下邊界和上邊界。

為了使得高超聲速滑翔飛行器安全返回，其再入飛行過程的氣動加熱受到限制，即

(qs)3D≤qU

(6)

式中：(qs)3D為三維駐點(diǎn)熱流；qU為熱流上限。

對于主曲率半徑分別為R1和R2(R1≤R2) 的任意形狀三維駐點(diǎn)，其熱流表達(dá)式為[32]

(7)

式中：(qs)AXI為軸對稱駐點(diǎn)熱流，k=R1/R2。

軸對稱駐點(diǎn)的熱流表達(dá)式如下

(8)

式中：ρe和μe分別為邊界層外緣(即正激波后)的氣流密度和粘性系數(shù)，haw和hw分別為壁面恢復(fù)焓和壁面焓，壁面溫度采用輻射平衡壁面溫度。相關(guān)參數(shù)的計(jì)算方法參見文獻(xiàn)[32]。

參數(shù)due/dx為駐點(diǎn)速度梯度，表達(dá)式如下

(9)

式中：u∞為來流速度大小;ρ∞為來流速度大小;ρs為正激波后氣流等熵滯止密度(駐點(diǎn)密度)。

與其他軌跡優(yōu)化研究經(jīng)常采用的氣動加熱計(jì)算方法相比，式(8)和(9)給出的氣動加熱計(jì)算方法考慮了更多的影響因素，因而具有更高的精度，計(jì)算的駐點(diǎn)熱流與風(fēng)洞實(shí)驗(yàn)和高可信度CFD方法的結(jié)果比較一致[32]。此外，該方法能夠考慮馬赫數(shù)沿飛行軌跡的變化(以及由此導(dǎo)致的相關(guān)變量的變化)，因而沿整個(gè)飛行軌跡都具有較高的預(yù)測精度。

2.3 末端約束

再入軌跡的末端約束參照航天飛機(jī)的末端能量管理窗口選取，具體取值如下

h(tf)=hf,v(tf)=vf,γ(tf)=γf

(10)

2.4 目標(biāo)函數(shù)

橫向航程是衡量高超聲速飛行器再入機(jī)動能力的重要指標(biāo)。本文取最大橫向航程為優(yōu)化指標(biāo)。由于θ0=0°，φ0=0°，因而目標(biāo)函數(shù)可寫為

(11)

綜上所述，高超聲速滑翔飛行器再入軌跡優(yōu)化問題描述為：確定控制變量u=[α(t),σ(t)]T，使得目標(biāo)函數(shù)(11)最小化，并且滿足狀態(tài)方程(1)，初始條件(4)，路徑約束(5)和(6)以及終端約束(10)。

3 快速、高精度軌跡優(yōu)化技術(shù)

本文采用局部配點(diǎn)法[6]求解高超聲速滑翔再入軌跡優(yōu)化問題。局部配點(diǎn)法通過在一組離散節(jié)點(diǎn)上對軌跡優(yōu)化問題離散化，將其轉(zhuǎn)化為NLP問題。目前的NLP求解器(SNOPT、IPOPT等)在求解NLP時(shí)需要反復(fù)計(jì)算其偏導(dǎo)數(shù)。因此，提高NLP偏導(dǎo)數(shù)的計(jì)算效率對于提高軌跡優(yōu)化效率具有重要意義。稀疏差分算法[9-11]通過分析NLP偏導(dǎo)數(shù)的稀疏特性，將NLP偏導(dǎo)數(shù)分解為原始軌跡優(yōu)化問題的偏導(dǎo)數(shù)，從而能夠顯著減小計(jì)算量，并且不影響優(yōu)化精度。此外，離散節(jié)點(diǎn)的數(shù)量直接決定了NLP的規(guī)模，與軌跡優(yōu)化的計(jì)算量密切相關(guān)。多分辨率網(wǎng)格細(xì)化技術(shù)[16-17]能夠采用較少的節(jié)點(diǎn)取得較高的優(yōu)化精度。本文將稀疏差分法與多分辨率網(wǎng)格細(xì)化技術(shù)相結(jié)合，快速、高精度求解軌跡優(yōu)化問題。

3.1 稀疏差分算法

采用基于梯度的算法(比如SQP算法)求解NLP時(shí)，需要提供其目標(biāo)函數(shù)和約束對優(yōu)化變量的一階甚至二階偏導(dǎo)數(shù)才能提高求解效率。為了高效計(jì)算偏導(dǎo)數(shù)，本文采用一種稀疏差分算法[11]。

以狀態(tài)方程為例，將其寫成如下向量形式

(12)

式中：x為n維狀態(tài)向量;f為n維狀態(tài)方程函數(shù)向量;u為m維控制向量;t為時(shí)間變量。

(13)

(14)

根據(jù)矩陣鏈?zhǔn)角髮?dǎo)法則可得到狀態(tài)方程離散格式約束對全部自變量的偏導(dǎo)數(shù)如下

(15)

可見，經(jīng)過上述推導(dǎo)，狀態(tài)方程離散殘差約束(13)對自變量的一階偏導(dǎo)數(shù)可以分解為原始軌跡優(yōu)化問題的狀態(tài)方程對其自變量的偏導(dǎo)數(shù)。將方程(15)與狀態(tài)方程(1)相結(jié)合可以得到離散殘差約束的偏導(dǎo)數(shù)矩陣的稀疏型(判斷矩陣的每一項(xiàng)是否為零)和非零項(xiàng)的高效計(jì)算方法。采用類似的方法可以得到NLP的目標(biāo)函數(shù)和其他約束對自變量偏導(dǎo)數(shù)的高效計(jì)算方法，具體參見文獻(xiàn)[11]。由于原始軌跡優(yōu)化問題的約束和變量的數(shù)量與離散后的NLP相比大幅度減少，因而這樣處理可以顯著提高偏導(dǎo)數(shù)的計(jì)算效率，從而提高軌跡優(yōu)化問題的求解效率。

3.2 網(wǎng)格細(xì)化技術(shù)

傳統(tǒng)的多分辨率技術(shù)基于二進(jìn)制網(wǎng)格，其初始節(jié)點(diǎn)數(shù)目必須為2j+1，應(yīng)用不方便。文獻(xiàn)[19-20]定義了一種節(jié)點(diǎn)數(shù)不受限制的廣義二分網(wǎng)格，但是要求軌跡優(yōu)化的初始網(wǎng)格節(jié)點(diǎn)數(shù)必須是最粗糙網(wǎng)格的2倍，使用仍然不夠方便。文獻(xiàn)[21]定義了一種新型廣義二分網(wǎng)格能夠克服這一缺陷。本文將這種新型廣義二分網(wǎng)格[21]與多分辨率技術(shù)[16]相結(jié)合構(gòu)造網(wǎng)格細(xì)化技術(shù)，細(xì)化再入滑翔軌跡的離散節(jié)點(diǎn)。

在區(qū)間[0, 1]內(nèi)，由N個(gè)均勻分布的初始子區(qū)間反復(fù)對分得到的廣義二分網(wǎng)格為[21]

-1≤j≤Jmax

(16)

式中：整數(shù)j表示網(wǎng)格分辨率，整數(shù)k表示節(jié)點(diǎn)位置，整數(shù)Jmax表示網(wǎng)格最高分辨率。由于j的最小值為-1，因而整數(shù)N必須為偶數(shù)。采用Wj, N表示屬于Vj+1, N但不屬于Vj, N的節(jié)點(diǎn)，即

0≤k≤2jN-1}， 0≤j≤Jmax-1

(17)

因此，節(jié)點(diǎn)τj+1, k∈Vj+1, N滿足如下關(guān)系

(18)

廣義二分網(wǎng)格的子空間Vj, N滿足嵌套關(guān)系，即V0, N?V1, N?…?VJmax, N，并且當(dāng)Jmax→∞時(shí)，VJmax, N= [0, 1]。子空間Wj, N滿足Wj, N∩Wl, N=?(j≠l)。由定義和可知，二分網(wǎng)格是將區(qū)間[0, 1]不斷對分獲得，并且當(dāng)k≥j時(shí)滿足Wk, N∩Vj, N=?。圖2給出一組廣義二分網(wǎng)格的示例(N= 6,Jmax= 5)。其中V0, N=V-1, N∪W-1, N通常作為初始網(wǎng)格。

圖2 廣義二分網(wǎng)格的示例(N=6, Jmax= 5)Fig.2 Generalized dyadic mesh with N=6 and Jmax= 5

對于一組離散最優(yōu)解，多分辨率技術(shù)[16]的基本思想是在每個(gè)離散節(jié)點(diǎn)處通過其周圍節(jié)點(diǎn)插值計(jì)算該節(jié)點(diǎn)的函數(shù)值。如果插值結(jié)果與離散最優(yōu)解差別較大，則需要在該節(jié)點(diǎn)附近細(xì)化網(wǎng)格；否則不進(jìn)行細(xì)化。為了應(yīng)用多分辨率技術(shù)，需要結(jié)合軌跡優(yōu)化的精度要求，選取初始網(wǎng)格參數(shù)N，對應(yīng)于控制變量的細(xì)化誤差限ε，最高分辨率Jmax和最大細(xì)化迭代次數(shù)Imax(Imax≥Jmax+1，Imax的默認(rèn)值為Jmax+1)。賦初值I=1，Gold=V0, N，估計(jì)NLP的優(yōu)化初值，將其記為Xold。那么，網(wǎng)格細(xì)化方法的流程如下：

1)以Gold為初始網(wǎng)格、Xold為初值求解NLP。若I≥Imax，結(jié)束算法；否則，轉(zhuǎn)入下一步。

2)執(zhí)行步驟(1)～(4)，對Gold進(jìn)行細(xì)化。

(1)令Φold={uj,k:τj,k∈Gold}，并將其改寫為Φold={φl(τj,k):l=1,…,m,τj,k∈Gold}。

(2)初始化中間網(wǎng)格Gint=V0, N以及相應(yīng)的函數(shù)值Φint={φl(τ0,k)∈Φold, 0 ≤k≤N,l=1, …,m}。

(3)多分辨率細(xì)化算法。

(4)本次細(xì)化得到的非均勻網(wǎng)格為Gnew=Gint,對應(yīng)的函數(shù)值為Φnew=Φint。

3)將變量I加1。如果細(xì)化后的網(wǎng)格Gnew與Gold相同，則停止細(xì)化；否則，將步驟1)求解出的NLP最優(yōu)解插值到新網(wǎng)格Gnew，并以此作為新的優(yōu)化初值Xold，令Gold=Gnew，返回步驟1)。

在上文網(wǎng)格細(xì)化方法的流程中，步驟2)對應(yīng)的多分辨率細(xì)化算法是關(guān)鍵，具體算法如下：

1)賦初值j=-1。

2)求Wj, N和Gold的交集

3)賦初值i=1，執(zhí)行如下步驟(1)～(6)。

(5)將新增加節(jié)點(diǎn)對應(yīng)的函數(shù)值添加到Φint。若新增節(jié)點(diǎn)對應(yīng)的函數(shù)值未知，則利用網(wǎng)格Gold和函數(shù)Φold通過p階ENO插值計(jì)算。

4)將變量j增加1。若j≤Jmax，則返回步驟2)；否則，轉(zhuǎn)入下一步。

5)結(jié)束本次迭代的網(wǎng)格細(xì)化算法。

4 仿真結(jié)果

本節(jié)應(yīng)用稀疏差分法和多分辨率網(wǎng)格細(xì)化技術(shù)求解高超聲速滑翔飛行器再入軌跡優(yōu)化問題，以驗(yàn)證方法的有效性。軌跡優(yōu)化的參數(shù)取值為：引力常數(shù)μ=3.98603195×1014m3/s2；飛行器質(zhì)量m=17000 kg，氣動力參考面積A=71.4 m2；控制量邊界αmin=-10°，αmax= 35°，σmin=-90°，σmax= 90°；熱流上限qU= 4.54×105W/m2；駐點(diǎn)主曲率半徑R1=R2= 0.98 m；普朗特常數(shù)Pr=0.715；初始狀態(tài)h0=79 248 m，θ0=0°，φ0=0°，v0=7334 m/s，ψ0=0°，γ0=0°；終端狀態(tài)約束hf=24 384 m，vf= 762 m/s，γf=-5°。采用SNOPT 7[7]求解由軌跡優(yōu)化問題離散得到的NLP。采用的計(jì)算平臺為MacBook Air(處理器為Intel Core i5-5250U 1.6 GHz，操作系統(tǒng)為Windows 10企業(yè)版，編程語言為Matlab R2009a)。

4.1 最大橫向航程軌跡

最大橫向航程可用于評估高超聲速滑翔飛行器的再入機(jī)動能力。采用本文第3節(jié)的方法研究該問題。網(wǎng)格細(xì)化參數(shù)取值為：N=20，ε=0.1°，Jmax= 3，p= 2，N1= 1，N2= 1。仿真結(jié)果表明，本文方法耗時(shí)約5 s即可優(yōu)化出最大橫向航程再入軌跡，單側(cè)最大橫向航程為25.85°(對應(yīng)2874.5 km)。

圖3給出優(yōu)化的攻角和速度傾側(cè)角隨時(shí)間的變化曲線。圖中圓圈為離散最優(yōu)解，實(shí)線由ENO插值得出。網(wǎng)格細(xì)化算法共迭代4次，從321個(gè)均勻離散節(jié)點(diǎn)中選取74個(gè)節(jié)點(diǎn)描述最優(yōu)攻角和速度傾側(cè)角。由圖3可知，最優(yōu)攻角和速度傾側(cè)角并非光滑變化，而本文方法能夠在攻角和傾側(cè)角變化較劇烈的區(qū)域加密網(wǎng)格，在其變化平緩的區(qū)域采用稀疏網(wǎng)格，整體效果是采用較少的離散節(jié)點(diǎn)高精度求解了最優(yōu)攻角和速度傾側(cè)角。如果采用相同數(shù)目的均勻離散節(jié)點(diǎn)，顯然無法取得這種效果。

圖4給出最優(yōu)狀態(tài)變量的變化曲線。圖4中圓圈為離散最優(yōu)解，細(xì)實(shí)線為采用數(shù)值積分方法(四階Runge-Kutta方法)沿最優(yōu)控制變量積分動力學(xué)微分方程得到的結(jié)果?？梢?，數(shù)值積分結(jié)果與離散最優(yōu)解非常一致，證實(shí)了方法的精度。

圖3 最優(yōu)控制變量(74個(gè)節(jié)點(diǎn))Fig.3 Optimal control solution(74 points)

圖5給出駐點(diǎn)熱流沿最優(yōu)軌跡的變化曲線。從圖5中的局部放大圖可以看出，優(yōu)化的駐點(diǎn)熱流嚴(yán)格滿足路徑約束，再次展示了方法的精度。

圖4 最優(yōu)狀態(tài)變量(74個(gè)節(jié)點(diǎn))Fig.4 Optimal state solution(74 points)

圖5 駐點(diǎn)熱流變化曲線Fig.5 Time history of stagnation heating rate

為與前文了對比，圖6給出采用開源最優(yōu)控制程序GPOPS 4.1(其原理是hp自適應(yīng)偽譜法[24])求解該問題得到的最優(yōu)攻角和速度傾側(cè)角。其中，GPOPS 4.1的網(wǎng)格細(xì)化迭代次數(shù)與本文方法相同，其余參數(shù)設(shè)置取GPOPS 4.1的默認(rèn)值。仿真結(jié)果表明，GPOPS 4.1采用了228個(gè)離散節(jié)點(diǎn)描述最優(yōu)解，將近是本文方法的3.1倍。分析圖6可知，GPOPS 4.1除了在控制變量變化劇烈區(qū)域進(jìn)行了細(xì)化，還在控制變量變化平緩區(qū)域進(jìn)行了細(xì)化，而后者不是必須的，因此導(dǎo)致離散節(jié)點(diǎn)過多。可見，本文方法在精度相當(dāng)?shù)那闆r下，需要的離散節(jié)點(diǎn)數(shù)量較少。

圖6 hp自適應(yīng)偽譜法求解的最優(yōu)控制變量(228個(gè)節(jié)點(diǎn))Fig.6 Optimal control solution from hp adaptive pseudospectral method (228 points)

4.2 再入落點(diǎn)區(qū)

高超聲速滑翔飛行器的再入落點(diǎn)區(qū)對于飛行器返回機(jī)場的設(shè)計(jì)、應(yīng)急機(jī)場的選取，以及高超聲速再入滑翔武器的攻擊區(qū)分析等具有重要意義。本文首先采用軌跡優(yōu)化方法得到最大縱向航程軌跡和最小縱向航程軌跡，然后在不同的縱向航程約束情況下分別優(yōu)化最大橫向航程，得到了再入落點(diǎn)區(qū)的邊界(軌跡優(yōu)化的初始條件、路徑約束和終端約束與最大橫向航程問題相同)。圖7給出右側(cè)再入落點(diǎn)區(qū)的若干條邊界軌跡(由于不考慮地球自轉(zhuǎn)，左側(cè)再入落點(diǎn)區(qū)與右側(cè)再入落點(diǎn)區(qū)呈對稱分布)。為了便于展示，圖中還給出了再入軌跡在地球表面的投影以及落點(diǎn)區(qū)的邊界?？梢?，高超聲速滑翔再入飛行器能夠在利用氣動力在大范圍內(nèi)機(jī)動飛行。

圖7 再入飛行器右半側(cè)落點(diǎn)區(qū)邊界軌跡示例Fig.7 Examples of trajectories at the boundary of reentry vehicle’s right-half landing footprint

將對應(yīng)于不同縱向航程的落點(diǎn)區(qū)邊界軌跡的末端位置相連即可得到落點(diǎn)區(qū)邊界，如圖8所示。圖中的符號“+”表示飛行器的初始再入點(diǎn)位置，即(θ0,φ0)=(0°, 0°)?？梢?，再入飛行器的落點(diǎn)區(qū)接近于橢圓形，縱向范圍θ=12.62°～102.81°(對應(yīng)1401.1～11 432.3 km)，橫向范圍φ=-25.85°～25.85°(對應(yīng)-2874.5～2874.5 km)。高超聲速滑翔飛行器能夠在滿足軌跡約束的情況下到達(dá)落點(diǎn)區(qū)內(nèi)的任意位置。由于本文研究中沒有考慮地球自轉(zhuǎn)的影響，因而落點(diǎn)區(qū)呈左右對稱分布。若考慮地球自轉(zhuǎn)，那么落點(diǎn)區(qū)的形狀會發(fā)生畸變，不再對稱分布；但是這種變化相對于落點(diǎn)區(qū)范圍而言通常比較小。

圖8 再入飛行器落點(diǎn)區(qū)邊界Fig.8 Boundary of entry vehicle’s landing footprint

5 結(jié) 論

本文研究了采用稀疏差分法和網(wǎng)格細(xì)化技術(shù)快速、高精度求解高超聲速滑翔飛行器三維再入軌跡優(yōu)化問題。建立了與實(shí)際工程問題比較一致的三維再入軌跡優(yōu)化模型(含有熱流約束、控制量邊界約束和末端能量窗口約束等，采用離散大氣模型、離散氣動數(shù)據(jù)和較為詳細(xì)的駐點(diǎn)熱流模型)。從兩個(gè)方面研究了提高軌跡優(yōu)化的效率:1)采用了一種高效的稀疏差分法分析由軌跡優(yōu)化問題離散得到的NLP一階偏導(dǎo)數(shù)的稀疏特性并給出NLP偏導(dǎo)數(shù)非零項(xiàng)的計(jì)算方法，從而提高了NLP的求解效率;2)采用了一種新型基于廣義二分網(wǎng)格的網(wǎng)格細(xì)化算法對離散節(jié)點(diǎn)進(jìn)行細(xì)化，以較少的節(jié)點(diǎn)高精度描述最優(yōu)控制變量，從而在保證精度的前提下進(jìn)一步減小軌跡優(yōu)化的計(jì)算量。仿真結(jié)果表明，所述方法在桌面計(jì)算機(jī)上耗時(shí)約5 s即可優(yōu)化出一條嚴(yán)格滿足各種約束的三維最優(yōu)再入軌跡，并且數(shù)值積分結(jié)果與優(yōu)化的離散解非常一致。此外，所述方法能夠快速分析高超聲速滑翔飛行器的再入落點(diǎn)區(qū)，進(jìn)一步證實(shí)了方法的有效性和應(yīng)用潛力。