基于Gauss偽譜法的高超聲速飛行器多約束三維再入軌跡優化

姚寅偉 李華濱

北京航天自動控制研究所，北京 100854

高超聲速再入飛行器利用氣動力在大氣層邊緣長距離滑翔，全程機動飛行。再入飛行過程中受到如氣動、熱流等嚴格的過程約束；同時，再入初始狀態根據任務的不同會有較大變化，而且飛行過程中會根據信息鏈更改滑翔結束點的使用狀態。因此，針對此類飛行器需要在飛行過程中進行在線軌跡規劃。這就要求優化算法既要滿足一定的計算精度，又要滿足快速性的要求。

近幾年發展起來的偽譜法(Pseudospectral Method)，以其求解精度高和收斂速度快的特點，在復雜的最優控制問題中得以廣泛的應用。根據離散點的選取方式不同，常見的偽譜方法有Legendre偽譜法、Gauss偽譜法。其中，通過Gauss偽譜法轉化的NLP(Nonlinear Programming Problem非線性規劃問題)求解的KKT(Karush-Kuhn-Tucker最優化條件)條件和連續最優控制問題的一階必要條件相同，這一點優于Legendre偽譜法[1]。影響偽譜法計算精度和收斂速度的主要因素有： 1)參數初值； 2)離散節點的數量。因為偽譜法能夠以較少的節點數量獲得較精確的解[2]，所以影響偽譜方法的主要因素就是參數初值。采用隨機算法求解NLP，如遺傳算法，粒子群算法等，這類方法不依賴參數初值的選取，但是計算時間長，無法滿足快速性的要求。文獻[3]提出一種從可行解到最優解的串行優化策略，以較少的Gauss節點進行優化得到的可行解作為初值[3]，然后采用具有超一次收斂性的SQP方法求解NLP，但初始節點的選取仍然要人為選取，不具有在線自主規劃的能力。

在現有的快速軌跡規劃方法中，文獻[4]中所述方法能夠快速計算一條滿足各種過程約束和狀態約束的三維再入軌跡，但其橫向軌跡設計部分僅通過一次改變傾側角的符號進行控制，不太適用于大范圍橫向機動情況。

因此，本文在現有研究成果的基礎上，采用擬平衡滑翔條件結合橫向誤差走廊的方法計算滿足各種約束的三維再入軌跡，作為高斯偽譜法的初值，彌補了算法對初值的敏感性。然后采用“Gauss偽譜+SQP”的優化方法，進一步求解滿足各種約束、同時使性能指標最小的再入軌跡，最后對通過數值積分的結果和優化的插值計算結果進行比較，分析了此方法的有效性。

1 再入軌跡優化問題

1.1 無動力再入運動模型

最優控制問題的數學描述一般要求給定自變量的初始值和終端值，對于軌跡優化來說，初始時刻是已知的，但終端時刻無法確定，因此引入反值能量e代替時間作為自變量，其表達式為：

e=1/r-V2/2

(1)

它隨時間單調變化，符合作為自變量的要求，而且只要已知再入初始和終端的狀態，運動方程就可在固定自變量區間[e0,ef]內進行積分計算，滿足軌跡優化算法的要求。e對無量綱時間τ求導得：

(2)

再入飛行過程中，假設地球模型為均質圓球，考慮地球旋轉，側滑角為0，采用式(1)所述的能量e代替時間作自變量，可得再入飛行器無量綱三自由度運動方程：

θ·

(3)

(4)

(5)

(6)

(7)

(8)

Kθ3= 2ωeVcosBsinA

KA3= 2ωeV(tanθcosAcosB-sinB)

L=RmρV2SrefCL/2m

D=RmρV2SrefCD/2m

(9)

ρ為大氣密度；Sref，m分別為飛行器參考面積和質量；CL，CD是與α相關的升力系數和阻力系數。

1.2 再入軌跡約束條件

(1)飛行過程約束

根據飛行器結構和外形的特點，飛行過程中需要滿足：熱流約束、動壓約束、過載約束、給定攻角下的平衡滑翔條件，即：

≤

(10)

(11)

n=Lcosα+Dsinα≤nmax

(12)

(13)

(2)終端約束

再入滑翔段的終端狀態同時也是末段制導的初始狀態，根據末段制導的要求，滑翔段終端狀態需要滿足一定約束，即位置約束(包括地心距rf、經度λf、緯度Bf)和速度約束(包括速度大小Vf，彈道傾角θf，航向角Af)，其表達式為：

r(e)f=rf,λ(e)f=λf,B(e)f=Bf

(14)

|V(e)f-Vf|≤ΔV,|θ(e)f-θf|≤Δθ

|A(e)f-Af|≤ΔA

(15)

“Δ”表示各終端狀態允許的誤差范圍，根據飛行任務的不同，部分狀態約束條件可忽略(如：彈道傾角、航向角等)。

(3)控制量約束

飛行過程中，為防止控制量過大，需要對其限幅，即：

|γ(e)|≤γmax

|α(e)|≤αmax

(16)

1.3 軌跡優化性能指標

選取飛行過程總的氣動加熱為優化指標，目標函數表達式為對駐點熱流密度的積分，即：

(17)

2 高斯偽譜軌跡優化方法

2.1 基本原理

高斯偽譜法將狀態變量和控制變量在一系列高斯點上離散，并以這些離散點為節點構造拉格朗日插值多項式來逼近狀態和控制變量。通過對全局插值多項式求導來近似狀態變量的導數，從而將微分方程約束轉換為一組代數約束。性能指標中的積分項由高斯積分計算。終端狀態也由初始狀態和右函數的積分獲得。經上述變換，可將最優控制問題轉化為具有一系列代數約束的參數優化問題，即非線性規劃問題(NLP)。

2.2 連續最優控制問題的離散

(1)區間變換

高斯偽譜法是在[-1,1]區間內通過Lagrange多項式擬合狀態量和控制量，因此需要對原積分區間[e0,ef]進行區間變換，用自變量τ代替e：

(18)

變換后，得到新的區間τ∈[τ0,τf]=[-1,1]。軌跡優化問題可描述為Lagrange型的最優控制問題：

1)根據方程(17)得到目標函數：

τ

(19)

2)根據方程(3)～(8)，得到狀態方程：

τ(τ),U(τ),τ)

(20)

3)根據方程(14)和(15)，得到邊界條件

g(X(τf),τf)=0

(21)

4)根據方程(10)～(13)飛行過程約束

h[X(τ),U(τ),τ]≤0

(22)

(2)狀態變量和控制變量的近似

高斯偽譜法對軌跡進行離散處理，選取N個Legendre正交多項式PN(τ)的零點(Legendre-Gauss點)τk(k=1,2,…,N)和起始點τ0=-1，終端點τf=1作為節點，構成Lagrange插值多項式，并以此為基函數構造狀態變量和控制變量的近似表達式，即：

(23)

(24)

其中Lagrange插值基函數：

(25)

(26)

x(τ)是節點τ=τi(i=0,1,…,N)處的狀態量,u(τ)是節點τ=τk(k=1,2,…,N)處的控制量。由式(23)～(26)可以得到，在節點處，實際狀態量(控制量)與近似狀態量(控制量)相等，即X(τi)=x(τi)，U(τk)=u(τk)。

末端節點因為要滿足微分方程約束，因此不能用上述插值多項式近似表達，按下式計算：

ωkF(x(τk),u(τk),τk)

(27)

其中，ωk為高斯權重，計算公式為[6]：

(28)

(3)微分方程約束的轉化

對式(23)進行求導：

τk)≈(τ)x(τi)

(29)

其中微分矩陣D∈RN×(N+1)可以離線確定，即：

τk)

(30)

利用式(29)代替離散點處微分方程的左邊式，從而將微分方程約束轉換成代數約束：

(31)

其中k=1,2,…,N。

(4)性能指標函數的近似

將Lagrange型性能指標函數中的積分項用Gauss積分近似，得：

τk)

(32)

經過上述變換，基于高斯偽譜法的最優控制問題可以描述為：已知初始點的狀態x0，求離散點上的狀態變量xi和控制變量ui(i=1,2,…,N)，使性能指標(32)最小，并且滿足微分方程約束(31)和終端狀態約束(27)，以及原最優控制問題的邊界條件

g(xf,τf)=0

(33)

和飛行過程約束

h[xk,uk,τk]≤0(k=1,2,…,N)

(34)

從而將原最優控制問題轉化為一個一般的非線性規劃問題(NLP)，即：

s.t.gj(x)≥0,j=1,2,…,p

hj(x)≥0,j=1,2,…,l

(35)

其中，x為設計變量，包括6個狀態變量和傾側角、攻角2個控制變量，p為邊界條件個數，l為過程約束個數。

2.3 高斯偽譜法求解NLP問題的方法

目前，求解NLP的方法中，SQP算法是公認的最有效的算法之一。它具有整體收斂性和局部超一次收斂性[7]。因此，本文采用SQP算法對上述NLP進行求解。

3 初值軌跡計算方法

在實際應用過程中，Gauss偽譜法對待優化參數初值較敏感，初值的好壞直接影響算法的收斂速度和計算精度，即初值越接近最優解，算法收斂越快。而接近最優解的必要條件是滿足各種約束，因此，本文不考慮性能指標，只考慮滿足飛行過程約束和終端狀態，快速計算一條三維再入軌跡，在其上選擇一系列離散點作為Gauss偽譜法待優化參數的初值。

3.1 初始下降段

初始下降段在高度-速度(r-V)剖面內不滿足平衡滑翔約束。因此要軌跡快速平滑地過渡到平衡滑翔狀態，即：

δ

(36)

(37)

3.2 擬平衡滑翔段

(1)縱向軌跡計算方法

由于初始下降段結束后，飛行器達到平衡滑翔狀態，因此采用基于擬平衡滑翔條件的方法計算滑翔段縱向再入軌跡。

忽略地球旋轉的條件下，不考慮橫向運動，剩余航程對能量微分的表達式為：

θ·

(38)

將式(38)與式(6)相除，滿足平衡滑翔條件(13)時，θ≈0，可近似得到：

(39)

假設平衡滑翔起點和終點對應的剩余航程分別為stogo0和stogo1，將傾側角的大小設計為速度的分段線性函數：

(40)

其中，V0和V1分別表示滑翔段初始點和終點的速度；γ0和γ1分別表示滑翔段初始點和終點的傾側角；Vmid=0.75V1表示中間速度，γmid表示Vmid所對應的傾側角的大小，為待優化變量。將式(39)在區間[stogo0,stogo1]內積分得到的末速度隨傾側角γmid單調變化[3]，因此通過割線法搜索滿足剩余航程要求的γmid。因為γ1是根據平衡滑翔條件和滑翔段終點的地心距和速度求出，所以當滑翔終點速度滿足終端約束后，滑翔終點地心距亦滿足約束。

(2)橫向軌跡計算方法

縱向軌跡計算得到了傾側角的大小，橫向控制則通過改變傾側角的符號來完成。對于航程遠，橫向機動大的再入飛行器，單次或者兩次反號可能無法實現軌跡的精確控制，因此，本文采用航向角誤差走廊控制傾側角符號的變化，即航向角誤差超出設定值就改變傾側角符號。

4 數值仿真實例

數值仿真程序包括初值軌跡生成和高斯偽譜優化，均在C++環境下編寫，計算機CPU為Inter(R)Core(TM)2 Duo 3.0GHz。

4.1 軌跡初值對Gauss偽譜法收斂速度的影響

表1為高斯偽譜法的滑翔飛行段終端狀態允許偏差，也是算法收斂的退出條件。表2給出了對應不同軌跡初值時，Gauss偽譜法在滿足表1的收斂條件下的迭代次數。

表1 終端狀態允許偏差

表2 Gauss偽譜法在不同軌跡初值下收斂的迭代次數

從表2中可以看出在橫向狀態(經度、緯度、航向角)偏差基本相等的情況下，初值軌跡的縱向狀態(高度、速度)偏差越小，算法迭代次數越少。其中初值軌跡4是采用本文第3部分所述的方法計算得到，相對于其他3個終端狀態偏差較大的軌跡初值而言，從收斂速度上有明顯優勢。

4.2 Gauss偽譜優化再入軌跡數值仿真

1)初始狀態：

[h0,λ0,B0,V0,θ0,A0]=[55km,30°,5°,5000m/s,-1°,90°]

2)終端狀態：

[hf,λf,Bf,Vf,θf,Af]=[32km,55.24°,4.53°,2000m/s,0°,92.14°]

以再入全程總吸熱量最小作為性能指標，如式(17)所示。Gauss節點數N=30。初值軌跡采用表2中的4號軌跡。終端狀態收斂條件同表1。

圖1～圖4給出了數值仿真的結果，圖例中Initial表示初值軌跡；GP表示節點，是Gauss偽譜法的設計變量；GPM表示通過Gauss偽譜法優化出的節點插值得到結果；ODE表示將插值得到的控制量代入原運動方程積分得到的結果。圖1表示高度曲線，初值軌跡，即表1中的4號軌跡，接近最終優化后的軌跡，因此，算法迭代次數較少，整個優化程序在第4部分所述計算機環境下運行時間小于10s，說明該初值算法與Gauss偽譜法相結合的方法在優化快速性上具有一定優勢，可行有效。圖2表示經緯度曲線，即地軌跡。由圖1和圖2可以看出優化后的結果和積分得到的結果基本一致，說明高斯偽譜法對實際運動模型的近似具有較高的精度。圖3和圖4分別表示控制量傾側角和攻角。表3列出了優化算法對于飛行過程的滿足情況，可以看到：所有約束均小于允許值。而且只要算法收斂，即說明算法滿足表1所示的終端狀態約束，因此，Gauss偽譜法優化得到的結果滿足所有約束。

圖2 經度-緯度曲線

圖3 能量-傾側角曲線

圖4 能量-攻角曲線

表3 過程約束

5 結論

研究了臨近空間高超聲速飛行器的三維再入軌跡的快速優化方法。采用以能量作為自變量的運動方程，結合“初值軌跡生成+高斯偽譜法+SQP求解NLP”的優化方法，相對于一般不含初值生成的直接法，能以較少的迭代次數收斂，提高了算法的計算速度。通過對優化得到的結果和數值積分得到的結果進行比較分析，驗證了算法的實際可飛性。本文提出的方法可作為快速軌跡優化方法的理論儲備，同時可以作為下一步制導問題的基礎。

參 考 文 獻

[1] David Benson.A Gauss Pseudospectral Transcraption for Optimal Control [D].Doctor Dissertation of Massachusetts Institute of Technology, 2005: 109-111.

[2] Huntington G T.Advancement and Analysis of a Gauss Pseudospectral Transcription [D].Cambridge, MA: Massachusetts Institute of Technology, 2007.

[3] 雍恩米.高超聲速滑翔式再入飛行器軌跡優化與制導方法研究[D].長沙：國防科學技術大學，2008：49-50, 58-60.

[4] Zuojun Shen, Ping Lu.On-board Generation of Three-dimensional Constrained Entry Trajectories [J].AIAA Guidance, Navigation, and Control Conference and Exhibit.Monterey, California, 2002-4455.

[5] 胡正東.天基對地打擊武器軌道規劃與制導技術研究[D].長沙：國防科學技術大學博士論文，2009，120-123.

[6] 宗群，田栢苓，竇立謙.基于Gauss偽譜法的臨近空間飛行器上升段軌跡優化[J].宇航學報，2010,31(7):1776-1780.(Zong Qun, Tian Bai-ling, Dou Li-qian.Ascent Phase Trajectory Optimization for Near Space Vehicle Based on Gauss Pseudospectral Methods[J].Journal of Astronautics, 2010,31(7):1776-1780 (in Chinese).)

[7] Philip E G.An SQP Algorithm for Large Scale Constrained Optimization [J].SIAM Review, 2005, 47(1):99-131.