高超聲速滑翔飛行器再入軌跡快速、高精度優化

趙吉松，張建宏，李 爽

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

0 引 言

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

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

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

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

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

1 再入軌跡運動模型

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

(1)

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

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

(2)

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

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

(3)

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

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

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

2 再入軌跡優化問題

2.1 初始條件

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

(4)

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

2.2 路徑約束

再入軌跡優化問題的控制變量為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分別為邊界層外緣(即正激波后)的氣流密度和粘性系數，haw和hw分別為壁面恢復焓和壁面焓，壁面溫度采用輻射平衡壁面溫度。相關參數的計算方法參見文獻[32]。

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

(9)

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

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

2.3 末端約束

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

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

(10)

2.4 目標函數

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

(11)

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

3 快速、高精度軌跡優化技術

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

3.1 稀疏差分算法

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

以狀態方程為例，將其寫成如下向量形式

(12)

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

(13)

(14)

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

(15)

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

3.2 網格細化技術

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

在區間[0, 1]內，由N個均勻分布的初始子區間反復對分得到的廣義二分網格為[21]

-1≤j≤Jmax

(16)

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

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

(17)

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

(18)

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

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

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

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

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

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

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

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

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

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

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

1)賦初值j=-1。

2)求Wj, N和Gold的交集

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

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

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

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

4 仿真結果

本節應用稀疏差分法和多分辨率網格細化技術求解高超聲速滑翔飛行器再入軌跡優化問題，以驗證方法的有效性。軌跡優化的參數取值為：引力常數μ=3.98603195×1014m3/s2；飛行器質量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；普朗特常數Pr=0.715；初始狀態h0=79 248 m，θ0=0°，φ0=0°，v0=7334 m/s，ψ0=0°，γ0=0°；終端狀態約束hf=24 384 m，vf= 762 m/s，γf=-5°。采用SNOPT 7[7]求解由軌跡優化問題離散得到的NLP。采用的計算平臺為MacBook Air(處理器為Intel Core i5-5250U 1.6 GHz，操作系統為Windows 10企業版，編程語言為Matlab R2009a)。

4.1 最大橫向航程軌跡

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

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

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

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

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

圖4 最優狀態變量(74個節點)Fig.4 Optimal state solution(74 points)

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

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

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

4.2 再入落點區

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

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

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

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

5 結 論

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