王文虎

中北大學(xué)機(jī)電工程學(xué)院，太原 030051

亞軌道飛行器(Suborbital Launch Vehicle，SLV)是指從地面發(fā)射，在亞軌道空間運(yùn)行作業(yè)，然后重返地面的跨大氣層可重復(fù)使用飛行器[1]。SLV返回段具有飛行高度低、任務(wù)多樣性、飛行環(huán)境多變的特點(diǎn)。

快速、準(zhǔn)確、魯棒的SLV返回軌跡優(yōu)化可以實(shí)現(xiàn)快速任務(wù)規(guī)劃從而極大地降低成本。返回軌跡優(yōu)化本質(zhì)是求解一類高度復(fù)雜的最優(yōu)控制問(wèn)題。常見的數(shù)值方法主要包括間接法和直接法。近年來(lái)，直接法中的偽譜法因其求解精度高、收斂速度快、具有很好的魯棒性而備受關(guān)注[2]。偽譜法主要有：勒讓德偽譜法(Legendre Pseudospectral Method, LPM)、拉道偽譜法(Radau Pseudospectral Method, RPM)和高斯偽譜法(Gauss Pseudospectral Method, GPM)。

國(guó)內(nèi)外關(guān)于GPM及LPM研究較多，而RPM國(guó)內(nèi)未見相關(guān)文獻(xiàn)。Huntington[3]對(duì)上述3種偽譜法在計(jì)算效率和精度方面進(jìn)行了比較，但所用的RPM是向后Radau偽譜法，并且算例較為簡(jiǎn)單，沒有比較各種偽譜法在處理復(fù)雜問(wèn)題時(shí)的能力。目前，國(guó)內(nèi)外還沒有文獻(xiàn)對(duì)向前RPM與GPM進(jìn)行過(guò)比較。

本文針對(duì)亞軌道飛行器，分別采用GPM和向前RPM進(jìn)行了SLV返回軌跡快速優(yōu)化研究，同時(shí)對(duì)2種偽譜法進(jìn)行了比較分析。結(jié)果表明，Gauss偽譜法在處理含控制量約束的問(wèn)題時(shí)比向前Radau偽譜法更具優(yōu)勢(shì)。

1 Gauss偽譜法與向前Radau偽譜法

1.1 Gauss偽譜法與向前Radau偽譜法基本原理

偽譜法是一種正交配點(diǎn)法，其基本思路為：通過(guò)選擇合適的正交配點(diǎn)與節(jié)點(diǎn)，構(gòu)造全局插值多項(xiàng)式來(lái)逼近狀態(tài)和控制變量，配點(diǎn)處狀態(tài)量的導(dǎo)數(shù)可由全局插值多項(xiàng)式求導(dǎo)來(lái)近似，從而將微分方程約束轉(zhuǎn)換為一組代數(shù)約束。性能指標(biāo)中的積分項(xiàng)可由相應(yīng)的Gauss求積公式計(jì)算得到。經(jīng)上述離散化方法，原最優(yōu)控制問(wèn)題轉(zhuǎn)化為非線性規(guī)劃問(wèn)題，而后通過(guò)有效的大規(guī)模稀疏NLP求解器即可求解。

Gauss偽譜法與向前Radau偽譜法基本原理如圖1所示，其主要區(qū)別在于：1)選用的配點(diǎn)不同。GPM選勒讓德-高斯(Legendre-Gauss, LG)點(diǎn)∈(-1,1)，而向前RPM選標(biāo)準(zhǔn)勒讓德-高斯-拉道(Legendre-Gauss-Radau, LGR)點(diǎn)∈[-1,1)；2)近似狀態(tài)量所用的方式不同。GPM用LG點(diǎn)與初始時(shí)刻點(diǎn)為節(jié)點(diǎn)構(gòu)造插值多項(xiàng)式，終端狀態(tài)通過(guò)高斯求積公式得到，而向前RPM用全部離散點(diǎn)為節(jié)點(diǎn)構(gòu)造插值多項(xiàng)式近似所有狀態(tài)量。關(guān)于Gauss偽譜法文獻(xiàn)較多，在此不再贅述。具體GPM離散化方法可參見文獻(xiàn)[4]。

圖1 Gauss偽譜法與向前Radau偽譜法基本原理圖

1.2 向前Radau偽譜法離散化最優(yōu)控制問(wèn)題[5]

1)時(shí)域變換

為將原問(wèn)題由t∈[t0,tf]轉(zhuǎn)化到τ∈[-1,1]，對(duì)時(shí)間變量t作變換：

(1)

時(shí)域變換后的Bolza問(wèn)題如下：

性能指標(biāo)

J=Φ(x(-1),t0,x(1),tf)+

(2)

動(dòng)力學(xué)微分方程約束

(τ),u(τ),τ;t0,tf)

(3)

邊界條件約束

φ(x(-1),t0,x(1),tf)=0

(4)

路徑約束

C(x(τ),u(τ),τ;t0,tf)≤0

(5)

2)配點(diǎn)選取及狀態(tài)與控制變量的近似

選N次與N-1次Legendre正交多項(xiàng)式之和PN(τ)+PN-1(τ)的根即標(biāo)準(zhǔn)LGR點(diǎn)作為配點(diǎn)，配點(diǎn)數(shù)為N。以τN+1=1和N個(gè)標(biāo)準(zhǔn)LGR點(diǎn)為節(jié)點(diǎn)構(gòu)造Lagrange插值多項(xiàng)式來(lái)近似狀態(tài)變量：

(6)

Li(τ)(i=1,…,N+1)為L(zhǎng)agrange基函數(shù)。

以N個(gè)LGR點(diǎn)為節(jié)點(diǎn)構(gòu)造Lagrange插值多項(xiàng)式來(lái)近似控制變量：

(7)

3)動(dòng)力學(xué)微分方程轉(zhuǎn)化

對(duì)(6)式求導(dǎo)，結(jié)合(3)式可以得出：

,Uk,τk;t0,tf)=0

(k=1,…,N)

(8)

其中Xk=X(τk)∈Rn，Uk=U(τk)∈Rm為配點(diǎn)處狀態(tài)和控制量，Dki∈RN×(N+1)為微分矩陣，可通過(guò)式(9)離線確定：

(9)

其中φ(τi)=(1-τi)[PN(τi)+PN-1(τi)]。

4)性能指標(biāo)

性能指標(biāo)中的積分項(xiàng)均通過(guò)高斯拉道求積公式來(lái)近似：

J=Φ(X0,t0,Xf,tf)+

(10)

其中，wk為高斯拉道求積公式權(quán)系數(shù)。

5)邊界條件約束與路徑約束

邊界與路徑約束只需將離散狀態(tài)和控制量代入式(4)，(5)即可。

至此，由式(8)，(10)以及離散化后的式(4)，(5)將原最優(yōu)控制問(wèn)題轉(zhuǎn)化為非線性規(guī)劃問(wèn)題。

2 SLV返回軌跡優(yōu)化

2.1 動(dòng)力學(xué)模型

考慮地球旋轉(zhuǎn)影響，假定飛行器側(cè)滑角為0，SLV無(wú)動(dòng)力返回動(dòng)力學(xué)方程為[6]：

γ

2ω(tanγcosφsinψ-sinφ)-

(11)

從實(shí)際系統(tǒng)考慮，控制舵面偏轉(zhuǎn)限制必然使得攻角α、傾側(cè)角σ存在帶寬和速率限制，在模型中考慮α，σ角速率限制。令

，

(12)

由式(11)與(12)聯(lián)立，可得最優(yōu)控制問(wèn)題狀態(tài)變量X=[rθφVγψασ]T∈R8，而控制變量為U=[uαuσ]T∈R2。uα,uσ稱為“偽控制量”，可以對(duì)其加以限制來(lái)滿足實(shí)際控制能力約束。

為加快優(yōu)化速度、有效地利用自動(dòng)微分技術(shù)，密度ρ(r)、當(dāng)?shù)匾羲賏(r)以及升阻力系數(shù)CL,CD均采用擬合模型。

當(dāng)李莉把交了兩個(gè)月租金的出租屋讓給許峰，自己投奔梅子的時(shí)候，梅子指著她的腦袋痛心疾首：“世界上最傻的女人就是你，許峰那個(gè)白眼狼，一看就是個(gè)攀高枝的。他騙你騙得高段啊，你走得那個(gè)順溜啊?！?/p>

2.2 優(yōu)化性能指標(biāo)與各類約束

1)性能指標(biāo)與終端約束

為描述性能指標(biāo)與終端約束，引入“末端進(jìn)場(chǎng)走廊”(FAC)概念(見圖2)，F(xiàn)AC也可以理解為“進(jìn)場(chǎng)著陸窗口”，是在返回段終端速度、航跡角與航向角、著陸場(chǎng)位置及方向給定情況下，能夠保證飛行器安全返回的終端位置(經(jīng)緯高)的范圍。

圖2 末端進(jìn)場(chǎng)走廊(FAC)[6]

終端約束要求在滿足其它各類約束條件下飛行器能夠準(zhǔn)確的捕獲FAC，即

≤r(tf)≤

V(tf)=Vf

γ(tf)=γf

ψ(tf)=ψf

其中rFAC，θFAC，φFAC分別為FAC中心的地心距、經(jīng)度、緯度，F(xiàn)ACL，F(xiàn)ACW，F(xiàn)ACH分別為三維FAC的長(zhǎng)、寬、高。

由于SLV返回的主要目標(biāo)是能夠安全返回指定的著陸場(chǎng)，因此選擇如下性能指標(biāo)：

w1,2,3為權(quán)重系數(shù)，這里選w1,2,3=1，即使得終端位置與末端進(jìn)場(chǎng)走廊中心位置距離最小。

2)其它約束

3 優(yōu)化算例與結(jié)果分析

3.1 優(yōu)化算例

本文采用X-33總體參數(shù)及氣動(dòng)模型進(jìn)行仿真計(jì)算[6]。SLV質(zhì)量35828kg，參考面積149.3881m2，地球半徑Re取6378.1km。

返回初始狀態(tài)與終端條件約束：

[r0,θ0,φ0,V0,γ0,ψ0]=[6429.1km,-85°,30°,2.6km/s,-1.3°,0°]

[Vf,γf,ψf]=[91.44m/s,-6°,-60°]

[rFAC,θFAC,φFAC]=[6378.7km,-80.7112°,30°]

[FACH/2,FACW/2,FACL/2]=[122m,0.00137°,0.001097°]

返回過(guò)程中狀態(tài)量約束：

≤≤

偽控制量約束：

≤≤

過(guò)載、動(dòng)壓、熱流等路徑約束：

≤≤

3.2 優(yōu)化結(jié)果

分別采用Gauss偽譜法和向前Radau偽譜法對(duì)不采用偽控制量和采用偽控制量2種情況進(jìn)行了優(yōu)化，采用國(guó)外大規(guī)模稀疏NLP軟件包SNOPT求解，優(yōu)化結(jié)果如表1及圖3～7所示。兩種方法均選用30個(gè)配點(diǎn)，優(yōu)化時(shí)間都在10s以內(nèi)，考慮到本文算法在Matlab平臺(tái)下實(shí)現(xiàn)，雅克比矩陣通過(guò)自動(dòng)微分方法計(jì)算，如采用C代碼、解析雅克比矩陣、無(wú)量綱化動(dòng)力學(xué)方程等手段，優(yōu)化速度有望進(jìn)一步提高。優(yōu)化性能指標(biāo)表明終端時(shí)刻飛行器非常接近FAC中心。由圖4可以看出，采用偽控制量后GPM和向前RPM的攻角、傾側(cè)角曲線更為平滑。圖5表明采用GPM，控制量變化率均在要求的±5(°)/s范圍內(nèi)，而采用向前RPM傾側(cè)角變化率在終端時(shí)刻達(dá)到了-5.21738(°)/s，超出了要求的范圍，這是由于標(biāo)準(zhǔn)LGR點(diǎn)不包括終端時(shí)刻，導(dǎo)致向前RPM在處理終端控制約束時(shí)受限。圖6表明返回過(guò)程滿足熱流、法向過(guò)載及動(dòng)壓等路徑約束。兩種方法結(jié)果趨勢(shì)基本吻合，由于配點(diǎn)個(gè)數(shù)較少，也導(dǎo)致了一些差異。

表1 優(yōu)化結(jié)果比較

圖3 狀態(tài)量變化曲線

圖4 控制量變化曲線

圖5 偽控制量變化曲線

圖6 路徑約束變化曲線

圖7 哈密爾敦函數(shù)變化曲線

3.3 可行性與最優(yōu)性驗(yàn)證

對(duì)于采用離散化方法求解連續(xù)最優(yōu)控制問(wèn)題而言，驗(yàn)證解的可行性與最優(yōu)性是非常必要的。

可行性通過(guò)數(shù)值積分結(jié)果與優(yōu)化結(jié)果相比較來(lái)驗(yàn)證，積分時(shí)的控制量由最優(yōu)離散控制量插值得到。從圖3中可以看出GPM及向前RPM優(yōu)化與數(shù)值積分結(jié)果非常吻合，表明所得結(jié)果是可行的。

與一般直接法不同，偽譜法可以根據(jù)協(xié)態(tài)映射原理估算配點(diǎn)處的協(xié)態(tài)信息，從而求得哈密爾敦函數(shù)值以檢驗(yàn)解的最優(yōu)性。對(duì)于存在路徑約束、終端時(shí)刻自由的自治系統(tǒng)最優(yōu)控制問(wèn)題，最優(yōu)軌跡對(duì)應(yīng)的哈密爾敦函數(shù)應(yīng)恒為0。圖7中哈密爾敦函數(shù)接近于0，實(shí)際數(shù)據(jù)在10-11量級(jí)，可以認(rèn)為所得結(jié)果接近最優(yōu)解。

4 結(jié)論

本文利用Gauss偽譜法和向前Radau偽譜法進(jìn)行了亞軌道飛行器返回軌跡優(yōu)化，比較了2種方法在處理復(fù)雜問(wèn)題時(shí)的能力。仿真結(jié)果表明，在滿足控制能力約束、路徑約束等條件下，GPM能夠快速準(zhǔn)確地生成捕獲FAC中心的SLV返回軌跡。這2種方法計(jì)算效率相當(dāng)，但向前RPM在處理存在控制量約束問(wèn)題時(shí)效果不佳，不能滿足終端控制能力約束。

參 考 文 獻(xiàn)

[1] Martin J C and LAW G W.Suborbital Reusable Launch Vehicles and Applicable Market[R].The Aerospace Corporation Report, 2002.

[2] Fahroo F, Ross I M.Advances in Pseudospectral Methods for Optimal Control[C].AIAA Guidance, Navigation and Control Conference and Exhibit.Honolulu, Hawaii, 2008:1-23.AIAA 2008-7309.

[3] Huntington G T, Benson D, Rao A V.A Comparison of Accuracy and Computational Efficiency of Three Pseudospectral Methods[C].AIAA Guidance, Navigation, and Control Conference and Exhibit.Hilton Head, South Carolina, 2007:1-25.AIAA 2007-6405.

[4] Huntington G T.Advancement and Analysis of a Gauss Pseudospectral Transcription for Optimal Control Problems[D].Cambridge, MA: Massachusetts Institute of technology, 2007.

[5] Garg D, Patterson M A, Darby C L, et al.Direct Trajectory Optimization and Costate Estimation of General Optimal Control Problems Using a Radau Pseudospectral Method[C].AIAA Guidance, Navigation, and Control Conference.Chicago, Illinois, 2009:1-29.AIAA 2009-5989.

[6] Singh B, Bhattacharya R.Optimal Guidance of Hypersonic Vehicles Using B-Splines and Galerkin Projection[C].AIAA Guidance, Navigation, and Control Conference and Exhibit.Honolulu, Hawaii, 2008:1-11.AIAA 2008-7263.