姚寅偉 李華濱

北京航天自動(dòng)控制研究所，北京 100854

高超聲速再入飛行器利用氣動(dòng)力在大氣層邊緣長距離滑翔，全程機(jī)動(dòng)飛行。再入飛行過程中受到如氣動(dòng)、熱流等嚴(yán)格的過程約束；同時(shí)，再入初始狀態(tài)根據(jù)任務(wù)的不同會(huì)有較大變化，而且飛行過程中會(huì)根據(jù)信息鏈更改滑翔結(jié)束點(diǎn)的使用狀態(tài)。因此，針對(duì)此類飛行器需要在飛行過程中進(jìn)行在線軌跡規(guī)劃。這就要求優(yōu)化算法既要滿足一定的計(jì)算精度，又要滿足快速性的要求。

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

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

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

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

1.1 無動(dòng)力再入運(yùn)動(dòng)模型

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

e=1/r-V2/2

(1)

它隨時(shí)間單調(diào)變化，符合作為自變量的要求，而且只要已知再入初始和終端的狀態(tài)，運(yùn)動(dòng)方程就可在固定自變量區(qū)間[e0,ef]內(nèi)進(jìn)行積分計(jì)算，滿足軌跡優(yōu)化算法的要求。e對(duì)無量綱時(shí)間τ求導(dǎo)得：

(2)

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

θ·

(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分別為飛行器參考面積和質(zhì)量；CL，CD是與α相關(guān)的升力系數(shù)和阻力系數(shù)。

1.2 再入軌跡約束條件

(1)飛行過程約束

根據(jù)飛行器結(jié)構(gòu)和外形的特點(diǎn)，飛行過程中需要滿足：熱流約束、動(dòng)壓約束、過載約束、給定攻角下的平衡滑翔條件，即：

≤

(10)

(11)

n=Lcosα+Dsinα≤nmax

(12)

(13)

(2)終端約束

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

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)

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

(3)控制量約束

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

|γ(e)|≤γmax

|α(e)|≤αmax

(16)

1.3 軌跡優(yōu)化性能指標(biāo)

選取飛行過程總的氣動(dòng)加熱為優(yōu)化指標(biāo)，目標(biāo)函數(shù)表達(dá)式為對(duì)駐點(diǎn)熱流密度的積分，即：

(17)

2 高斯偽譜軌跡優(yōu)化方法

2.1 基本原理

高斯偽譜法將狀態(tài)變量和控制變量在一系列高斯點(diǎn)上離散，并以這些離散點(diǎn)為節(jié)點(diǎn)構(gòu)造拉格朗日插值多項(xiàng)式來逼近狀態(tài)和控制變量。通過對(duì)全局插值多項(xiàng)式求導(dǎo)來近似狀態(tài)變量的導(dǎo)數(shù)，從而將微分方程約束轉(zhuǎn)換為一組代數(shù)約束。性能指標(biāo)中的積分項(xiàng)由高斯積分計(jì)算。終端狀態(tài)也由初始狀態(tài)和右函數(shù)的積分獲得。經(jīng)上述變換，可將最優(yōu)控制問題轉(zhuǎn)化為具有一系列代數(shù)約束的參數(shù)優(yōu)化問題，即非線性規(guī)劃問題(NLP)。

2.2 連續(xù)最優(yōu)控制問題的離散

(1)區(qū)間變換

高斯偽譜法是在[-1,1]區(qū)間內(nèi)通過Lagrange多項(xiàng)式擬合狀態(tài)量和控制量，因此需要對(duì)原積分區(qū)間[e0,ef]進(jìn)行區(qū)間變換，用自變量τ代替e：

(18)

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

1)根據(jù)方程(17)得到目標(biāo)函數(shù)：

τ

(19)

2)根據(jù)方程(3)～(8)，得到狀態(tài)方程：

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

(20)

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

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

(21)

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

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

(22)

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

高斯偽譜法對(duì)軌跡進(jìn)行離散處理，選取N個(gè)Legendre正交多項(xiàng)式PN(τ)的零點(diǎn)(Legendre-Gauss點(diǎn))τk(k=1,2,…,N)和起始點(diǎn)τ0=-1，終端點(diǎn)τf=1作為節(jié)點(diǎn)，構(gòu)成Lagrange插值多項(xiàng)式，并以此為基函數(shù)構(gòu)造狀態(tài)變量和控制變量的近似表達(dá)式，即：

(23)

(24)

其中Lagrange插值基函數(shù)：

(25)

(26)

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

末端節(jié)點(diǎn)因?yàn)橐獫M足微分方程約束，因此不能用上述插值多項(xiàng)式近似表達(dá)，按下式計(jì)算：

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

(27)

其中，ωk為高斯權(quán)重，計(jì)算公式為[6]：

(28)

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

對(duì)式(23)進(jìn)行求導(dǎo)：

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

(29)

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

τk)

(30)

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

(31)

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

(4)性能指標(biāo)函數(shù)的近似

將Lagrange型性能指標(biāo)函數(shù)中的積分項(xiàng)用Gauss積分近似，得：

τk)

(32)

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

g(xf,τf)=0

(33)

和飛行過程約束

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

(34)

從而將原最優(yōu)控制問題轉(zhuǎn)化為一個(gè)一般的非線性規(guī)劃問題(NLP)，即：

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

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

(35)

其中，x為設(shè)計(jì)變量，包括6個(gè)狀態(tài)變量和傾側(cè)角、攻角2個(gè)控制變量，p為邊界條件個(gè)數(shù)，l為過程約束個(gè)數(shù)。

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

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

3 初值軌跡計(jì)算方法

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

3.1 初始下降段

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

δ

(36)

(37)

3.2 擬平衡滑翔段

(1)縱向軌跡計(jì)算方法

由于初始下降段結(jié)束后，飛行器達(dá)到平衡滑翔狀態(tài)，因此采用基于擬平衡滑翔條件的方法計(jì)算滑翔段縱向再入軌跡。

忽略地球旋轉(zhuǎn)的條件下，不考慮橫向運(yùn)動(dòng)，剩余航程對(duì)能量微分的表達(dá)式為：

θ·

(38)

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

(39)

假設(shè)平衡滑翔起點(diǎn)和終點(diǎn)對(duì)應(yīng)的剩余航程分別為stogo0和stogo1，將傾側(cè)角的大小設(shè)計(jì)為速度的分段線性函數(shù)：

(40)

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

(2)橫向軌跡計(jì)算方法

縱向軌跡計(jì)算得到了傾側(cè)角的大小，橫向控制則通過改變傾側(cè)角的符號(hào)來完成。對(duì)于航程遠(yuǎn)，橫向機(jī)動(dòng)大的再入飛行器，單次或者兩次反號(hào)可能無法實(shí)現(xiàn)軌跡的精確控制，因此，本文采用航向角誤差走廊控制傾側(cè)角符號(hào)的變化，即航向角誤差超出設(shè)定值就改變傾側(cè)角符號(hào)。

4 數(shù)值仿真實(shí)例

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

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

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

表1 終端狀態(tài)允許偏差

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

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

4.2 Gauss偽譜優(yōu)化再入軌跡數(shù)值仿真

1)初始狀態(tài)：

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

2)終端狀態(tài)：

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

以再入全程總吸熱量最小作為性能指標(biāo)，如式(17)所示。Gauss節(jié)點(diǎn)數(shù)N=30。初值軌跡采用表2中的4號(hào)軌跡。終端狀態(tài)收斂條件同表1。

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

圖2 經(jīng)度-緯度曲線

圖3 能量-傾側(cè)角曲線

圖4 能量-攻角曲線

表3 過程約束

5 結(jié)論

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

參 考 文 獻(xiàn)

[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] 雍恩米.高超聲速滑翔式再入飛行器軌跡優(yōu)化與制導(dǎo)方法研究[D].長沙：國防科學(xué)技術(shù)大學(xué)，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] 胡正東.天基對(duì)地打擊武器軌道規(guī)劃與制導(dǎo)技術(shù)研究[D].長沙：國防科學(xué)技術(shù)大學(xué)博士論文，2009，120-123.

[6] 宗群，田栢苓，竇立謙.基于Gauss偽譜法的臨近空間飛行器上升段軌跡優(yōu)化[J].宇航學(xué)報(bào)，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.