趙吉松，張建宏，李 爽

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

0 引 言

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

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

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

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

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

1 再入軌跡運動模型

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

(1)

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

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

(2)

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

升力和阻力的表達式分別為

(3)

式中：ρ為大氣密度，采用美國標準大氣模型(1976年)[30]插值計算；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分別為傾側角的下邊界和上邊界。

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

(qs)3D≤qU

(6)

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

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

(7)

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

軸對稱駐點的熱流表達式如下

(8)

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

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

(9)

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

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

2.3 末端約束

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

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

(10)

2.4 目標函數(shù)

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

(11)

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

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

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

3.1 稀疏差分算法

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

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

(12)

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

(13)

(14)

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

(15)

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

3.2 網(wǎng)格細化技術

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

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

-1≤j≤Jmax

(16)

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

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

(17)

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

(18)

廣義二分網(wǎng)格的子空間Vj, N滿足嵌套關系，即V0, N?V1, N?…?VJmax, N，并且當Jmax→∞時，VJmax, N= [0, 1]。子空間Wj, N滿足Wj, N∩Wl, N=?(j≠l)。由定義和可知，二分網(wǎng)格是將區(qū)間[0, 1]不斷對分獲得，并且當k≥j時滿足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)解，多分辨率技術[16]的基本思想是在每個離散節(jié)點處通過其周圍節(jié)點插值計算該節(jié)點的函數(shù)值。如果插值結果與離散最優(yōu)解差別較大，則需要在該節(jié)點附近細化網(wǎng)格；否則不進行細化。為了應用多分辨率技術，需要結合軌跡優(yōu)化的精度要求，選取初始網(wǎng)格參數(shù)N，對應于控制變量的細化誤差限ε，最高分辨率Jmax和最大細化迭代次數(shù)Imax(Imax≥Jmax+1，Imax的默認值為Jmax+1)。賦初值I=1，Gold=V0, N，估計NLP的優(yōu)化初值，將其記為Xold。那么，網(wǎng)格細化方法的流程如下：

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

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

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

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

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

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

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

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

1)賦初值j=-1。

2)求Wj, N和Gold的交集

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

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

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

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

4 仿真結果

本節(jié)應用稀疏差分法和多分辨率網(wǎng)格細化技術求解高超聲速滑翔飛行器再入軌跡優(yōu)化問題，以驗證方法的有效性。軌跡優(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；駐點主曲率半徑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。采用的計算平臺為MacBook Air(處理器為Intel Core i5-5250U 1.6 GHz，操作系統(tǒng)為Windows 10企業(yè)版，編程語言為Matlab R2009a)。

4.1 最大橫向航程軌跡

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

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

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

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

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

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

圖5 駐點熱流變化曲線Fig.5 Time history of stagnation heating rate

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

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

4.2 再入落點區(qū)

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

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

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

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

5 結 論

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