劉 旭，李 爽*，曹 濤，宋 婷

(1.南京航空航天大學(xué)航天學(xué)院，南京 210016; 2. 上海航天控制技術(shù)研究所，上海 201109)

1 引言

中國載人航天工程已經(jīng)掌握近地軌道載人飛行器的返回與回收技術(shù)[1]，隨著嫦娥三號(hào)和嫦娥四號(hào)分別在月球正面和背面成功著陸[2-3]，載人探月工程已經(jīng)成為中國航天未來發(fā)展方向之一。

相比近地軌道再入返回，載人探月飛行器的返回軌跡往往是偏心率較大的橢圓軌道或者雙曲線軌道，因此再入速度較大，接近第二宇宙速度。若沿用近地再入任務(wù)的直接彈道再入方式，熱流和過載峰值都將超過飛行器和航天員所能承受的極限。此外，受測控能力限制，且著陸場的選取約束要求陸地定點(diǎn)著陸，因而采用跳躍式再入是目前較為理想的月地高速再入方案[4]。跳躍再入軌跡一般劃分為初次再入段、躍升段、開普勒段、二次再入段和著陸段。再入過程中飛行器采用配平迎角飛行，通過升阻比控制或傾側(cè)角控制實(shí)現(xiàn)再入軌跡調(diào)整，從而對(duì)過載、熱流、落點(diǎn)散布等指標(biāo)進(jìn)行控制[5]，因此飛行器的升阻比和傾側(cè)角控制能力直接決定了落點(diǎn)的可調(diào)整范圍和精度。

典型的跳躍再入飛行器包括美國的Apollo飛船和獵戶座飛船[6-8]以及中國探月工程三期再入返回試驗(yàn)器(CE-5T1)[9]。其中，Apollo飛船第一次進(jìn)入大氣層的速度在10.7 km/s，升力控制躍起高度為80 km左右，具有2200～4600 km的航程機(jī)動(dòng)的能力，側(cè)向航程修正能力優(yōu)于370 km，再入最大過載處于4～7g0，落點(diǎn)精度約為50 km；獵戶座飛船外形與Apollo飛船類似，升阻比約為0.25～0.27，全程采用反作用控制系統(tǒng)(Reaction Control System, RCS)模式控制跳躍式再入返回，但其搭載的新一代人-機(jī)協(xié)同控制策略可以讓返回艙軌道平面可調(diào)，不再需要將著陸場設(shè)置在航天器軌道平面內(nèi)，而且采用了精度更高的制導(dǎo)控制算法，再入段的峰值過載、熱流情況較阿波羅飛船均有較大改善，落點(diǎn)控制精度可達(dá)到5 km內(nèi)；CE-5T1的預(yù)期再入航程介于5600～7100 km 之間，實(shí)際航程超過6000 km[9]。

跳躍再入軌跡規(guī)劃本質(zhì)上是復(fù)雜的多約束、非線性最優(yōu)控制問題，常用的求解方法包括間接法和直接法。間接法基于拉格朗日乘子法和龐特里亞金極值原理，將多約束最優(yōu)控制問題轉(zhuǎn)化為多點(diǎn)邊值問題，然后采用優(yōu)化算法求解邊值問題獲得最優(yōu)解。早期Shi[10]基于Chapman方程和匹配漸近展開法，分別求解了二維和三維場景下的無約束大氣進(jìn)入軌跡，但僅能得到精度有限的一階近似解；隨后，Kuo等[11]引入修正的查普曼(Chapman)變量，采用改進(jìn)的匹配漸近展開法獲得了地球大氣進(jìn)入軌跡的一階和二階近似解；近來，吳旭忠等[12]基于間接法提出一種結(jié)合自動(dòng)判斷約束結(jié)構(gòu)和初始化對(duì)應(yīng)打靶方程技術(shù)的同倫方法，增大了跳躍式再入軌跡優(yōu)化問題的收斂半徑。然而數(shù)值優(yōu)化算法對(duì)協(xié)態(tài)變量初值非常敏感，同時(shí)協(xié)態(tài)變量沒有明確的物理意義，其初值難以猜測，導(dǎo)致間接法算法收斂性較差。同時(shí)，間接法還存在必須推導(dǎo)一階最優(yōu)條件、收斂半徑小、路徑約束問題中控制切換結(jié)構(gòu)難以預(yù)知的缺點(diǎn)[13]，因此應(yīng)用間接法進(jìn)行跳躍式再入軌跡優(yōu)化具有明顯的局限性。

相比之下，直接法不需要推導(dǎo)一階最優(yōu)條件，適用于多段、多約束軌跡優(yōu)化問題，且不需要猜測協(xié)態(tài)變量初值，對(duì)狀態(tài)變量初值也不敏感，具有良好的通用性和收斂性[13]。直接法中的局部配點(diǎn)法是一種通用性較強(qiáng)的軌跡優(yōu)化方法，通過在一組離散節(jié)點(diǎn)上對(duì)最優(yōu)控制問題進(jìn)行離散，將軌跡優(yōu)化問題轉(zhuǎn)化為非線性規(guī)劃(Nonlinear Programming, NLP)問題，并采用序列二次規(guī)劃(Sequential Quadratic Programming, SQP)算法進(jìn)行求解[14-15]。局部配點(diǎn)法在國外已經(jīng)到得到廣泛的應(yīng)用，軌跡優(yōu)化軟件OTIS[16]就是基于局部配點(diǎn)法開發(fā)的。

本文采用基于稀疏差分[17-18]和多分辨率技術(shù)的改進(jìn)局部配點(diǎn)法[19-21]開展載人高速再入總體方案研究，通過求解定點(diǎn)著陸問題、落點(diǎn)區(qū)域問題、再入?yún)?shù)敏感度問題和軌跡重規(guī)劃問題的最優(yōu)再入軌跡，論證載人飛行器在不同工況下的可行再入方案，確定再入軌跡對(duì)各項(xiàng)再入?yún)?shù)的敏感程度，為后續(xù)載人高速再入任務(wù)提供參考。

2 再入問題描述

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

對(duì)于再入軌跡優(yōu)化問題，一般遵循非旋轉(zhuǎn)均質(zhì)地球假設(shè)，且可忽略地球自轉(zhuǎn)影響，此時(shí)描述大氣再入飛行器質(zhì)心運(yùn)動(dòng)的微分方程為式(1)～式(6)：

(1)

(2)

(3)

(4)

(5)

(6)

其中，r為地心到飛行器質(zhì)心的徑向距離，V為飛行器的速度值，θ為經(jīng)度，φ為緯度，γ為航跡傾角，ψ為航跡偏角(指向正東方向?yàn)?°，逆時(shí)針為正)，σ為傾側(cè)角，地球引力g=μ/r2，m為飛行器的質(zhì)量，L和D分別為升力和阻力，表達(dá)式如式(7)所示：

L=ρV2SCL/2,D=ρV2SCD/2

(7)

其中，CL和CD分別為升力系數(shù)和阻力系數(shù)，S為航天器的參考面積。ρ為大氣密度，表達(dá)式如式(8)所示：

ρ=ρ0e-h/H

(8)

其中，ρ0為天體表面大氣密度，h為飛行器高度，H為大氣模型的參考高度。

2.2 約束條件

飛行器再入過程中應(yīng)當(dāng)滿足式(9)～式(11)給定的熱流密度、動(dòng)壓和法向過載約束：

(9)

(10)

(11)

飛行器二次再入段的初、末狀態(tài)應(yīng)當(dāng)滿足給定的邊界條件，且飛行過程中各狀態(tài)量均不能超過一定范圍，定義狀態(tài)的初、末約束E和過程約束C如式(12)所示：

E:x(t0)=x0,x(tf)=xf;C:xmin≤x≤xmax

(12)

其中，狀態(tài)量x=[r,θ,φ,V,γ,ψ]T，x0和xf表示給定的初始和末端狀態(tài)，xmin和xmax表示各狀態(tài)量的的下限與上限。飛行任務(wù)一般要求滿足全部初始約束，末端狀態(tài)滿足部分即可，例如定點(diǎn)再入時(shí)滿足末端狀態(tài)高度、速度和經(jīng)緯度約束即可，末端航跡傾角和航跡偏角不做具體要求。

再入過程中控制量為傾側(cè)角，σ的大小和變化速率均有一定變化范圍，定義控制約束如式(13)所示：

(13)

若給定相應(yīng)的目標(biāo)函數(shù)J=M+L(M為終端項(xiàng)，L為積分項(xiàng))，則載人多約束跳躍再入軌跡優(yōu)化問題可描述為：確定控制量傾側(cè)角σ，使得目標(biāo)函數(shù)最小化，并滿足動(dòng)力學(xué)方程式(1)～(6)和約束條件式(9)～(13)。

3 改進(jìn)的局部配點(diǎn)法

為了提高局部配點(diǎn)法的魯棒性和優(yōu)化效率，本文采用稀疏差分和基于廣義二分網(wǎng)格的多分辨率技術(shù)，從NLP的稀疏性和離散節(jié)點(diǎn)數(shù)量兩方面提高計(jì)算效率和精度。稀疏差分法采用原始軌跡優(yōu)化問題的一階偏導(dǎo)數(shù)表示NLP的一階偏導(dǎo)數(shù)，構(gòu)建稀疏偏導(dǎo)數(shù)矩陣，從而可以顯著減小計(jì)算量[18]；多分辨率方法基于新型廣義二分網(wǎng)格技術(shù)對(duì)離散網(wǎng)格節(jié)點(diǎn)進(jìn)行調(diào)整，能以少量節(jié)點(diǎn)高精度地描述最優(yōu)控制變量，從而減小軌跡優(yōu)化的規(guī)模和計(jì)算時(shí)間[19]。本文改進(jìn)局部配點(diǎn)法的流程如圖1所示，在傳統(tǒng)配點(diǎn)法的基礎(chǔ)上，采用了稀疏差分技術(shù)進(jìn)行偏導(dǎo)數(shù)計(jì)算，同時(shí)利用多分辨率技術(shù)進(jìn)行網(wǎng)格細(xì)化，從而在保證算法精度的同時(shí)減小計(jì)算量。

圖1 改進(jìn)局部配點(diǎn)法的算法流程圖Fig.1 Flowchart of modified local collocation method

3.1 稀疏差分

軌跡優(yōu)化問題中，狀態(tài)量x、過程約束C和目標(biāo)函數(shù)J中積分項(xiàng)L都定義在整個(gè)時(shí)域區(qū)間，因而可以將這3項(xiàng)對(duì)自變量的偏導(dǎo)數(shù)定義在一起，如式(14)所示：

(14)

其中，G1的每一項(xiàng)仍為矩陣，以?f/?xT為例，如式(15)所示：

(15)

通常情況下，G1是稀疏矩陣。為了描述G1的稀疏型，定義struct函數(shù)如式(16)、(17)所示：

(16)

記：

S1=struct(G1)

(17)

式中，struct(G1)表示對(duì)G1的每個(gè)元素進(jìn)行struct運(yùn)算。S1表示G1的稀疏型。為了得到S1，不需要計(jì)算G1的每個(gè)元素的具體值，只需要判斷是否為0。

類似地，可以定義邊界約束E和目標(biāo)函數(shù)J中Mayer項(xiàng)M對(duì)自變量的依賴關(guān)系矩陣和稀疏型如式(18)、(19)所示：

(18)

S1=struct(G1)

(19)

上述離散格式將同一個(gè)節(jié)點(diǎn)處的變量或約束記為一個(gè)列向量，這種記法與數(shù)值積分格式的形式一致，但是不利于推導(dǎo)偏導(dǎo)數(shù)矩陣的稀疏特性。為此，本文定義一種新的變量記法，將變量或約束的同一個(gè)分量在不同節(jié)點(diǎn)的值記為一個(gè)新的向量。以狀態(tài)方程離散殘差ξ為例，見式(20)：

ξ:,j=(ξ0,j,ξ1,j,…,ξN-1,j)T,(j=1,2,…,6)

(20)

利用這種新的變量記法，可將離散化的NLP問題重新描述為求解優(yōu)化變量z，使得目標(biāo)函數(shù)J(z)最小化，并且滿足約束條件如式(21)所示：

Fmin≤F(z)≤Fmax

(21)

式中，優(yōu)化變量z和約束函數(shù)F(z)的定義如下式(22)所示:

(22)

其中，t0，tf和Δt=tf-t0分別為軌跡優(yōu)化初始時(shí)刻、末端時(shí)刻和優(yōu)化時(shí)間。

3.2 廣義二分網(wǎng)格

傳統(tǒng)二進(jìn)網(wǎng)格的初始節(jié)點(diǎn)數(shù)目必須是2j+1(j為非負(fù)整數(shù))，使用不方便。為此，本文引入一種廣義二分網(wǎng)格。

在區(qū)間[0, 1]內(nèi)，由N個(gè)均勻分布的初始子區(qū)間對(duì)分得到的二分網(wǎng)格為式(23)：

(23)

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

(24)

因此，τj+1, k∈Vj+1, N滿足如式(25)所示關(guān)系:

(25)

二分網(wǎng)格的子空間Vj, N滿足嵌套關(guān)系，即V0,N?V1,N?…?VJmax,N，并且當(dāng)Jmax→∞時(shí)，VJmax,N=[0, 1]。子空間Wj, N滿足Wj, N∩Wl, N=? (j≠l)。由定義式(23)和(25)可知，二分網(wǎng)格是將區(qū)間[0, 1]不斷細(xì)分獲得，并且當(dāng)k≥j時(shí)滿足Wk,N∩Vj,N=?。

3.3 網(wǎng)格細(xì)化

網(wǎng)格細(xì)化采用基本無振蕩(Essentially Non-Oscillatory, ENO)插值算法離散節(jié)點(diǎn)插值，得到NLP問題的插值解，并根據(jù)插值解的變化率增刪節(jié)點(diǎn)。本文采用牛頓插值法構(gòu)造ENO插值，假設(shè)當(dāng)前插值節(jié)點(diǎn)為τi,τi+1, …,τi+n，則新增加節(jié)點(diǎn)τnew的選取方式如式(26)所示：

(26)

其中，f為待插值函數(shù)，對(duì)應(yīng)本文的動(dòng)力學(xué)公式(1)～(6)，f(τi-1,…,τi+n-1)和f(τi+1,…,τi+n+1)表示f的差商，定義如式(27)所示：

(27)

在進(jìn)行多分辨率優(yōu)化之前，結(jié)合軌跡優(yōu)化問題的精度要求，選取初始節(jié)點(diǎn)參數(shù)N，網(wǎng)格細(xì)化誤差ε，最高分辨率Jmax，最大迭代次數(shù)Imax(Imax≥Jmax/2)。采用四階Runge-Kutta法將軌跡優(yōu)化問題轉(zhuǎn)化NLP。設(shè)置I=1，Gold=V0, N，估計(jì)NLP優(yōu)化變量初值，將其記為Xold。那么，改進(jìn)的多分辨率方法的流程如下：

1)以Gold為初始網(wǎng)格、Xold為初值求解NLP。

2)網(wǎng)格細(xì)化方法如下(步驟①～④)：

①令Φold={uj,k:τj,k∈Gold}，并將其改寫為

Φold={φl(τj,k):l=1,…,m,τj,k∈Gold}。

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

③細(xì)化算法。

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

3)置I=I+1。如果細(xì)化后網(wǎng)格保持不變，那么結(jié)束細(xì)化；否則，將步驟1)解出的NLP最優(yōu)解插值到新網(wǎng)格Gnew，并以此作為新的優(yōu)化初值Xold, 令Gold=Gnew, 返回步驟1)，繼續(xù)細(xì)化。

以上細(xì)化算法是關(guān)鍵，具體如下：

1)求Wj,N和Gold的交集。

2)設(shè)置i=1，執(zhí)行如下步驟(①～⑥)。

⑤將新增加節(jié)點(diǎn)對(duì)應(yīng)的函數(shù)值添加到Φint。若新增節(jié)點(diǎn)對(duì)應(yīng)的函數(shù)值未知，則利用Gold和Φold通過p階ENO插值計(jì)算。

3)將j增加1。若j≤Jmax，返回步驟1)，否則轉(zhuǎn)入下一步。

4)結(jié)束細(xì)化算法。

4 再入軌跡優(yōu)化

載人探月飛行器再入時(shí)刻的各項(xiàng)參數(shù)取值如下：飛行器質(zhì)量m為6800 kg，參考面積S為12.88 m2，升力系數(shù)CL和阻力系數(shù)CD分別為0.3892和1.3479，kQ為1.9027×10-4，地球海平面大氣密度ρ0為1.225 kg/m3，參考高度H為7200 m，地球引力系數(shù)為3.986×1014m3/s2，地球半徑取6378 km，地球表面重力加速度g0取9.81 m/s2，其余飛行器系統(tǒng)參數(shù)和再入條件如表1所示，其中初次再入段和二次再入段的法向過載峰值分別為6倍g0和4倍g0。

表 1 數(shù)值仿真參數(shù)

此外，改進(jìn)局部配點(diǎn)法采用基于SQP算法的SNOPT7求解器進(jìn)行求解，運(yùn)行環(huán)境為搭載Intel Core i5-4200 M，4 GB內(nèi)存的個(gè)人計(jì)算機(jī)，編譯環(huán)境為MATLAB R2015b。

4.1 定點(diǎn)再入

定點(diǎn)再入要求載人探月飛行器在指定區(qū)域開傘，同時(shí)確保氣動(dòng)加熱和過載處于飛行器和航天員所能承受的安全范圍，對(duì)軌跡優(yōu)化算法的精度要求較高。假設(shè)目標(biāo)開傘點(diǎn)的經(jīng)緯度坐標(biāo)為θf=100°，φf=-4°，則定義目標(biāo)函數(shù)如式(28)所示：

J=(θ-θf)2+(φ-φf)2

(28)

其中，初末條件、過程約束和控制約束見表1。采用改進(jìn)局部配點(diǎn)法可在5 s內(nèi)獲得滿足全部約束條件的定點(diǎn)再入軌跡。圖2～圖7分別給出定點(diǎn)再入軌跡的三維再入軌跡曲線、速度曲線、傾側(cè)角曲線以及過程約束曲線。

圖2 定點(diǎn)再入三維軌跡Fig.2 Three-dimensional path of pin-point reentry

圖3 定點(diǎn)再入時(shí)間-速度曲線Fig.3 Time-speed profile of pin-point reentry

圖4 定點(diǎn)再入時(shí)間-傾側(cè)角曲線Fig.4 Time-bank profile of pin-point Reentry

圖5 定點(diǎn)再入時(shí)間-熱流曲線Fig.5 Time-heat flux profile of pin-point reentry

圖6 定點(diǎn)再入時(shí)間-動(dòng)壓曲線Fig.6 Time- pressure profile of pin-point reentry

圖7 定點(diǎn)再入時(shí)間-過載曲線Fig.7 Time-load profile of pin-point reentry

采用改進(jìn)的局部配點(diǎn)法進(jìn)行軌跡優(yōu)化耗時(shí)約2.126 s，而且從優(yōu)化結(jié)果可知再入軌跡具有典型的跳躍再入特征，最優(yōu)軌跡的末端狀態(tài)、過程約束等參數(shù)如表2所示。其中，末端高度10.6 km，末段經(jīng)緯度為100°、 -4°，末端速度308 m/s，末端航向角為-52.6°，末端航跡角為-2.4°，滿足定點(diǎn)再入約束；且最大熱流密度為239 W/cm2，動(dòng)壓為16.8 kPa，過載為4.5g0/ 3.3g0。同時(shí)，最優(yōu)軌跡的傾側(cè)角曲線光滑，便于制導(dǎo)系統(tǒng)跟蹤，表明算法具有良好的求解精度。此外由定點(diǎn)再入軌跡曲線可知，算法在軌跡變化平緩區(qū)域采用稀疏網(wǎng)格，在軌跡變化劇烈區(qū)域采用稠密網(wǎng)格，證明多分辨率技術(shù)可以有效控制離散節(jié)點(diǎn)的數(shù)量和分布，減小計(jì)算量。

表 2 最優(yōu)軌跡末端狀態(tài)和過程約束

4.2 落點(diǎn)區(qū)域

飛行器的再入落點(diǎn)區(qū)域?qū)τ谠u(píng)估飛行器機(jī)動(dòng)能力和著陸場地的規(guī)劃等具有重要意義。本文采用改進(jìn)局部配點(diǎn)法求解最大落點(diǎn)區(qū)域問題，即在不同的縱向航程約束下以最大橫向航程為目標(biāo)函數(shù)，過程約束同表1保持一致，得到再入落點(diǎn)區(qū)域的邊界，其中縱向航程和橫向航程均為球面上的一段大圓弧。由于不考慮地球自轉(zhuǎn)，正緯度方向上的再入落點(diǎn)區(qū)與負(fù)緯度方向上的再入落點(diǎn)區(qū)呈對(duì)稱分布，圖8給出負(fù)緯度方向上的再入落點(diǎn)區(qū)內(nèi)的若干條邊界軌跡。

圖8 落點(diǎn)區(qū)域右半側(cè)邊界軌跡Fig.8 Boundary paths of parachute deployment footprint

為了便于展示，圖9給出了再入軌跡在經(jīng)緯度平面的投影，將對(duì)應(yīng)于不同縱向航程的落點(diǎn)區(qū)域邊界軌跡的末端位置相連即可得到落點(diǎn)區(qū)域邊界。圖中紅色圓圈表示返回艙的初始再入點(diǎn)對(duì)應(yīng)的經(jīng)緯度位置，即(θ0,φ0)=(0°, 0°)?？梢娫偃腼w行器的落點(diǎn)區(qū)域接近于橢圓形，縱向范圍約為θ=18° ～ 185°(對(duì)應(yīng)2001.5 ～ 20 571.1 km)，橫向范圍φ=-5.488° ～ 5.488°(對(duì)應(yīng)-610.2 ～ 610.2 km)，飛行器能在滿足約束的情況下到達(dá)落點(diǎn)區(qū)內(nèi)的任意位置。若考慮地球自轉(zhuǎn)，那么落點(diǎn)區(qū)域的形狀會(huì)發(fā)生改變，但總體趨勢一致，呈現(xiàn)出橫向航程(或緯度絕對(duì)值)隨縱向航程(或緯度)增大而先增大、后減小的特征。

圖9 落點(diǎn)區(qū)域邊界Fig.9 Boundaries of parachute deployment footprint

4.3 再入?yún)?shù)敏感度

實(shí)際任務(wù)中由于受到各項(xiàng)干擾和誤差，載人探月飛行器的再入界面(或再入狀態(tài))將不可避免地偏離標(biāo)準(zhǔn)值，而飛行器的控制能力有限，因此確定當(dāng)前控制能力所能承受的最大初始偏差對(duì)任務(wù)設(shè)計(jì)具有重要意義。本文參考表1飛行器的再入控制能力，仿真條件參考定點(diǎn)再入問題，開傘點(diǎn)的空間位置保持為θf=100°、φf=-4°，采用改進(jìn)局部配點(diǎn)法研究再入軌跡對(duì)初始參數(shù)的敏感性，探討各初始參數(shù)的允許變化范圍。經(jīng)過單項(xiàng)拉偏仿真，表3給出了再入界面的全部參數(shù)的單項(xiàng)最大拉偏值以及對(duì)應(yīng)的最大過載和開傘點(diǎn)位置誤差。其中初始參數(shù)的基準(zhǔn)值為表1中各個(gè)參數(shù)的初始值，某項(xiàng)參數(shù)的正偏差或者負(fù)偏差是指當(dāng)其他參數(shù)取基準(zhǔn)值時(shí)，該參數(shù)允許的正向或負(fù)向最大拉偏值。拉偏值的選擇是通過多次軌跡優(yōu)化得到的，如以0.1 km為步長不斷增大/減小高度，直到初次再入或二次再入的最大過載超過6倍和4倍g0，從而確定高度的拉偏上下極限。同理可得其余各狀態(tài)量的拉偏范圍。

表3 再入界面變化范圍

表3表明即使存在較大高度偏差或經(jīng)緯度偏差，改進(jìn)局部配點(diǎn)法仍能優(yōu)化出滿足定點(diǎn)再入約束的再入軌跡，這表明算法對(duì)于初始參數(shù)偏差具有一定魯棒性，但過載峰值可能超出安全范圍，而算法質(zhì)量偏差對(duì)再入任務(wù)影響不大。此外，當(dāng)飛行器控制能力給定時(shí)，定點(diǎn)再入軌跡對(duì)于初始速度和航跡傾角最為敏感，速度變化范圍為0.21～0.53 km/s，而航跡傾角的最大變化范圍為-0.2°～+0.4°，過大的速度和航跡傾角偏差將無法保證再入的過程約束和落點(diǎn)精度，從而威脅航天員的生命，因此為保障任務(wù)安全，必須對(duì)初始再入速度的大小和方向進(jìn)行有效控制。

4.4 在線軌跡重規(guī)劃

再入過程中可能發(fā)生機(jī)構(gòu)故障或存在較大制導(dǎo)誤差，本文充分利用跳躍再入的特點(diǎn)，基于改進(jìn)局部配點(diǎn)法設(shè)計(jì)了在線軌跡重規(guī)劃方案：當(dāng)初次再入的制導(dǎo)誤差大于給定閾值，進(jìn)行二次再入界面預(yù)測，并在線重規(guī)劃二次再入軌跡。在線軌跡重規(guī)劃方案圖如圖10所示，制導(dǎo)控制流程如圖11所示。

圖10 在線軌跡重規(guī)劃示意圖Fig.10 Schematic diagram of trajectory re-planning

圖11 在線軌跡重規(guī)劃流程圖Fig.11 Flowchart of trajectory re-planning

從任務(wù)背景方面來看，跳躍再入軌跡的特殊性在于飛行器由于氣動(dòng)力作用會(huì)跳出大氣層，進(jìn)行一段慣性飛行。假設(shè)大氣層高度為120 km，以圖12所示的定點(diǎn)再入軌跡為例，飛行器在大約t=450 s時(shí)跳出大氣層，在大約t=1120 s時(shí)二次再入大氣層，時(shí)間相差長達(dá)670 s。假設(shè)在跳出大氣層時(shí)(A點(diǎn))開始對(duì)二次再入軌跡進(jìn)行優(yōu)化，只要在飛行器進(jìn)入大氣層之前(B點(diǎn))得到新的最優(yōu)軌跡，即可用于制導(dǎo)算法進(jìn)行跟蹤，并不要求實(shí)時(shí)計(jì)算出最優(yōu)軌跡。而飛行器從A點(diǎn)到B點(diǎn)的軌跡是大氣層外慣性軌跡，因此可以準(zhǔn)確預(yù)測。

圖12 定點(diǎn)再入時(shí)間-高度曲線Fig.12 Time-altitude profile of pin-point reentry

從在線優(yōu)化效率方面來看，當(dāng)前離線軌跡優(yōu)化耗時(shí)約5 s，在線優(yōu)化時(shí)還可以利用離線優(yōu)化的結(jié)果作為初值，加快收斂。目前軌跡優(yōu)化算法是基于MATLAB語言，將其改為能夠在嵌入式系統(tǒng)運(yùn)行的C語言從原理上是完全可行的，并且SNOPT求解器也有C語言版本可供使用。同時(shí)，參考獵戶座搭載的處理器性能(型號(hào)為PowerPC 750FX，運(yùn)行速度為900 MHz，2003年發(fā)布的iBook G3個(gè)人筆記本搭載的處理器)可知，載人飛行器的計(jì)算能力已有很大提高。綜上判斷改進(jìn)局部配點(diǎn)法在計(jì)算效率方面能夠滿足在線需求。

軌跡重規(guī)劃仿真時(shí)假設(shè)二次再入時(shí)的狀態(tài)變量為：r=6497.96 km，θ=69.8128°，φ=-3.1671°，V=7608.2 m/s，ψ=-1.8423°，γ=-11.3147°。以該初始狀態(tài)，采用軌跡優(yōu)化技術(shù)重新規(guī)劃一條定點(diǎn)著陸軌跡(θf=100°，φf=-4°)，圖13～圖16給出二次再入軌跡的三維再入軌跡以及過程約束曲線，重規(guī)劃軌跡的末端狀態(tài)和過程約束如表4所示。

圖13 軌跡重規(guī)劃的三維軌跡Fig.13 Three-dimensional path of trajectory re-planning

圖14 軌跡重規(guī)劃的時(shí)間-熱流曲線Fig.14 Time-heat flux profile of trajectory re-planning

圖15 軌跡重規(guī)劃的時(shí)間-動(dòng)壓曲線Fig.15 Time -pressure profile of trajectory re-planning

圖16 軌跡重規(guī)劃的時(shí)間-過載曲線Fig.16 Time-load profile of trajectory re-planning

表 4 重規(guī)劃軌跡的末端狀態(tài)和過程約束

4.5 算法對(duì)比

考慮定點(diǎn)再入問題，采用成熟的自適應(yīng)偽譜法工具包GPOPS-2[22]進(jìn)行軌跡優(yōu)化，該方法同樣采用基于SQP算法的SNOPT7求解器進(jìn)行NLP問題求解，其中網(wǎng)格細(xì)化采用基于離散誤差分析方法，計(jì)算NLP一階偏導(dǎo)數(shù)采用有限差分算法。仿真中各個(gè)狀態(tài)變量的初始值和取值范圍參考表1，目標(biāo)再入點(diǎn)的經(jīng)緯度和4.1節(jié)保持一致，分別為θf=100°，φf=-4°，末端高度為10 km，末端速度為0.1671 km/s，采用自適應(yīng)偽譜法的定點(diǎn)再入仿真結(jié)果如圖17～圖21所示。

圖17 定點(diǎn)再入三維軌跡(GPOPS-2)Fig.17 Three-dimensional path of pin-point reentry (GPOPS-2)

圖18 定點(diǎn)再入時(shí)間-傾側(cè)角曲線(GPOPS-2)Fig.18 Time-bank profile of pin-point reentry (GPOPS-2)

圖19 定點(diǎn)再入時(shí)間-熱流曲線(GPOPS-2)Fig.19 Time-heat flux profile of pin-point reentry (GPOPS-2)

圖20 定點(diǎn)再入時(shí)間-動(dòng)壓曲線(GPOPS-2)Fig.20 Time-pressure profile of pin-point reentry (GPOPS-2)

圖21 定點(diǎn)再入時(shí)間-過載曲線(GPOPS-2)Fig.21 Time-pressure profile of pin-point reentry (GPOPS-2)

采用GPOPS-2進(jìn)行軌跡規(guī)劃共耗時(shí)11.78 s，但峰值熱流達(dá)到300 W/cm2，超過了給定的最大熱流值260 W/cm2，且傾側(cè)角剖面變化較為劇烈。對(duì)比4.1節(jié)中改進(jìn)局部配點(diǎn)法的定點(diǎn)再入性能可知，在優(yōu)化精度相當(dāng)?shù)那闆r下，改進(jìn)局部配點(diǎn)法的計(jì)算效率比采用全局配點(diǎn)的自適應(yīng)偽譜法計(jì)算效率提高5倍，同時(shí)還能輸出更為平滑的控制量，具有更強(qiáng)的過程約束和控制約束處理能力。

5 結(jié)論

1)跳躍式再入最大縱向航程可達(dá)到20 571.1 km，將顯著減小熱流和過載峰值，可滿足中國測控和著陸場對(duì)航程的要求；

2)定點(diǎn)再入任務(wù)對(duì)初始速度大小和方向的控制精度較高，尤其是初始航跡傾角的誤差應(yīng)小于0.2°；

3)根據(jù)跳躍再入軌跡特征提出的在線軌跡重規(guī)劃方案提高了載人任務(wù)的魯棒性和安全性，將為載人登月再入任務(wù)提供重要技術(shù)支撐。