郭敏文 王大軼

1. 北京控制工程研究所,北京 100190 2. 空間智能控制技術(shù)重點(diǎn)實(shí)驗(yàn)室,北京 100190

探月返回器將以接近第二宇宙速度高速再入地球大氣層。返回器再入的過(guò)程考慮為配平攻角飛行，即忽略了攻角和側(cè)滑角的控制，僅僅通過(guò)調(diào)整傾側(cè)角來(lái)控制飛行軌跡。要實(shí)現(xiàn)大航程、強(qiáng)約束的再入飛行，對(duì)小升阻比探月返回器的制導(dǎo)控制提出了全新的要求。

再入制導(dǎo)算法主要分為2類(lèi)[1]：標(biāo)準(zhǔn)軌道法和預(yù)測(cè)校正法[2]，本文主要研究標(biāo)準(zhǔn)軌道法的參考軌跡優(yōu)化設(shè)計(jì)。從六七十年代至今，不少學(xué)者針對(duì)跳躍式再入軌跡設(shè)計(jì)問(wèn)題提出了解析方法，這些方法主要可分為2類(lèi)，一類(lèi)主要用匹配漸近展開(kāi)方法描述跳躍式再入軌跡。其中，Kuo和Vinh[3]利用改進(jìn)的匹配漸近展開(kāi)方法得到跳躍式再入軌跡的較高精度的二階近似表達(dá)式，Kluever[4]則基于該方法解析得到跳躍式再入?yún)⒖架壽E，并分段給出了跟蹤參考軌跡的傾側(cè)角指令。另一類(lèi)為參考軌跡阻力剖面的設(shè)計(jì)方法。Garcia-Llama[5]將跳躍式再入軌跡第一次再入段(初始再入點(diǎn)到跳出點(diǎn))的阻力剖面設(shè)計(jì)為四階多項(xiàng)式，并應(yīng)用反饋線性化方法跟蹤其設(shè)計(jì)的阻力剖面。其阻力剖面的設(shè)計(jì)求解思想有別于文獻(xiàn)[6]將再入軌跡阻力剖面設(shè)計(jì)為二次函數(shù)的做法。

近年來(lái)，許多學(xué)者針對(duì)軌跡規(guī)劃問(wèn)題提出了多種數(shù)值解法，且取得了重要進(jìn)展。這些方法主要可分為：直接法[7]和間接法[8-9]。間接法將最優(yōu)控制問(wèn)題最終轉(zhuǎn)化為兩點(diǎn)邊值問(wèn)題，該問(wèn)題求解過(guò)程對(duì)協(xié)態(tài)變量的初值高度敏感，收斂性差，同時(shí)軌跡規(guī)劃過(guò)程中存在的各種約束使其求解過(guò)程變得較為復(fù)雜。相對(duì)于間接法，直接法在收斂的魯棒性和解決實(shí)際問(wèn)題的適用性上更具有優(yōu)勢(shì)。最近發(fā)展起來(lái)的一種求解最優(yōu)控制問(wèn)題的改進(jìn)Gauss偽譜方法[10]是一種基于全局和局部插值多項(xiàng)式混合的直接配點(diǎn)法，其優(yōu)點(diǎn)在于可用盡量少的節(jié)點(diǎn)獲得較高的擬合精度。

應(yīng)該指出，軌跡優(yōu)化的理論已經(jīng)比較成熟，但因跳躍式再入軌跡的優(yōu)化設(shè)計(jì)有其新的特點(diǎn)，如強(qiáng)約束、多約束等，有必要對(duì)這方面進(jìn)行深入地研究。目前該研究方向的文獻(xiàn)相對(duì)較少。南英[11]對(duì)單次再入飛行、二次再入軌跡、多次再入軌跡在過(guò)載和熱流方面進(jìn)行了比較研究，并得出二次再入飛行為最優(yōu)返回軌跡的結(jié)論。針對(duì)跳躍式軌跡的優(yōu)化設(shè)計(jì)，Istratie[12]也提出了使多種不同性能指標(biāo)最小的優(yōu)化方案，但均未給出完整的最優(yōu)跳躍式再入軌跡，只是針對(duì)第一次再入段的飛行軌跡進(jìn)行了優(yōu)化設(shè)計(jì)。

為滿(mǎn)足過(guò)載和熱流等過(guò)程量的約束條件，本文利用改進(jìn)的Gauss偽譜方法[10]對(duì)完整的跳躍式再入軌跡進(jìn)行優(yōu)化設(shè)計(jì)。利用該方法計(jì)算比較了返回器在各相同初始再入條件下不發(fā)生跳躍和發(fā)生跳躍時(shí)的最遠(yuǎn)飛行距離，并對(duì)返回器在不發(fā)生跳躍時(shí)的飛行能力進(jìn)行了分析。最后通過(guò)仿真表明該方法能很好地實(shí)現(xiàn)小升力體月球返回的跳躍式再入軌跡優(yōu)化，有效滿(mǎn)足過(guò)載及駐點(diǎn)熱流的約束，并分析得出對(duì)于小升阻比返回器，當(dāng)再入航程要求大于3500km時(shí)，需采用跳躍式軌跡再入的結(jié)論。

1 探月返回器跳躍式再入軌跡優(yōu)化問(wèn)題描述

1.1 再入動(dòng)力學(xué)方程

整個(gè)飛行過(guò)程如圖1所示，圖中r0和rf分別為初始和末端時(shí)刻的地心距。假設(shè)探月返回器處于配平飛行狀態(tài)，且地球是一個(gè)均勻球體，不考慮地球扁率、地球公轉(zhuǎn)及地球自轉(zhuǎn)，大氣模型采用美國(guó)1976年標(biāo)準(zhǔn)大氣模型，返回器再入的無(wú)量綱質(zhì)心動(dòng)力學(xué)方程[1]為

(1)

圖1 探月返回器跳躍式再入軌跡示意圖

1.2 約束條件

因探月返回器高速再入的特點(diǎn)，過(guò)載和熱流約束變得十分苛刻。為了保證再入過(guò)程的安全，需要嚴(yán)格滿(mǎn)足以下約束條件：

(1)氣動(dòng)加熱約束

為減小氣動(dòng)加熱，要求駐點(diǎn)熱流不超過(guò)給定的最大值，即

(2)

(2)過(guò)載約束

為了減小再入時(shí)的過(guò)載，要求瞬時(shí)過(guò)載小于最大允許過(guò)載，即

(3)

其中nmax為最大允許過(guò)載值。

(3)終端狀態(tài)約束

考慮經(jīng)度、緯度、高度滿(mǎn)足終端約束條件：

(4)

(4)控制量約束

考慮到實(shí)際傾側(cè)角機(jī)動(dòng)不可能瞬時(shí)完成，需要對(duì)傾側(cè)角進(jìn)行限幅，即：

|σ|≤σmax

(5)

1.3 目標(biāo)函數(shù)

在高速再入過(guò)程中，為了減小防熱系統(tǒng)的重量，性能指標(biāo)一般取為再入過(guò)程的熱流最?。?/p>

(6)

2 跳躍式軌跡航程分析

(1)跳躍式軌跡擴(kuò)大航程能力分析

當(dāng)?shù)厍蜃赞D(zhuǎn)和正對(duì)著陸點(diǎn)的航向偏移量可以忽略時(shí)，相對(duì)著陸點(diǎn)的航程變量s可以由下式得到：

(7)

將式(7)增加為系統(tǒng)描述的第7個(gè)方程，即將縱程s考慮為一個(gè)新的變量x7，其中rf為再入開(kāi)傘點(diǎn)的高度。

由狀態(tài)方程(1)和(7)可看出，其中狀態(tài)變量r，v，γ和變量s的變化規(guī)律與另3項(xiàng)狀態(tài)變量解耦，這里通過(guò)單獨(dú)分析這4個(gè)量的運(yùn)動(dòng)規(guī)律，優(yōu)化計(jì)算各不同初始再入角情況下發(fā)生跳躍和不發(fā)生跳躍時(shí)的最大航程，以對(duì)跳躍式再入軌跡在實(shí)現(xiàn)長(zhǎng)縱程飛行方面的優(yōu)勢(shì)進(jìn)行比較分析。因此優(yōu)化性能指標(biāo)取為新變量x7的末狀態(tài)量，即

J=-x7f

(8)

這里，發(fā)生跳躍和不發(fā)生跳躍時(shí)的狀態(tài)約束如下，其中不發(fā)生跳躍的飛行路徑角保持小于0，即無(wú)上升飛行階段，軌跡不發(fā)生跳躍。

發(fā)生跳躍：

狀態(tài)量r和γ在整個(gè)飛行過(guò)程中被約束為：

-90°≤γ≤10°

rf≤r≤r0

不發(fā)生跳躍：

狀態(tài)量r和γ在整個(gè)飛行過(guò)程中被約束為：

-90°≤γ≤0°

rf≤r≤r0

(2)不發(fā)生跳躍時(shí)飛行能力分析

這里通過(guò)考慮升阻比偏差及大氣密度偏差對(duì)再入過(guò)程的影響，得到不同情況下返回器不發(fā)生跳躍時(shí)的最大航程和過(guò)載峰值，并對(duì)這些表征其飛行能力的數(shù)據(jù)進(jìn)行分析，以得到需采用跳躍式軌跡再入的航程要求。

氣動(dòng)系數(shù)偏差采用如下表達(dá)式

(9)

由于高層大氣變化極為復(fù)雜，從人造地球衛(wèi)星上天以來(lái)，各國(guó)學(xué)者一直致力于開(kāi)發(fā)高精度的標(biāo)準(zhǔn)大氣模型及其擾動(dòng)模型，采用了多個(gè)指標(biāo)(考慮經(jīng)緯度、季風(fēng)、大氣成分以及太陽(yáng)活動(dòng)等多種因素的影響)來(lái)描述大氣的擾動(dòng)。期間，形成了諸多不同的大氣攝動(dòng)模型，簡(jiǎn)單的有固定值偏差模型，復(fù)雜的如美國(guó)開(kāi)發(fā)的GRAM(Global Reference Atmosphere Model)模型。

本文取大氣攝動(dòng)模型為固定值偏差模型來(lái)分析不發(fā)生跳躍時(shí)再入軌跡的最大航程，其變化范圍取為[-30%, +30%]。

3 Gauss 偽譜法基本原理

高斯偽譜法將再入動(dòng)力學(xué)微分方程轉(zhuǎn)化為代數(shù)約束方程，將制導(dǎo)問(wèn)題轉(zhuǎn)化為不需要積分彈道的最優(yōu)規(guī)劃問(wèn)題。偽譜方法的一個(gè)顯著特征為譜收斂，即收斂速度大于N-m，其中N是配點(diǎn)個(gè)數(shù)，m是任意有限數(shù)值。

偽譜法的求解步驟：

1)不失一般性，將最優(yōu)控制問(wèn)題轉(zhuǎn)化為Bolza形式，即最小性能指標(biāo)函數(shù)：

(10)

滿(mǎn)足動(dòng)力學(xué)約束、等式邊值條件約束、不等式路徑約束如下：

(11)

這里t和τ之間的映射關(guān)系為：

(12)

2)基于Gauss偽譜方法的Bolza問(wèn)題離散化

(13)

(14)

(15)

對(duì)狀態(tài)變量表達(dá)式進(jìn)行微分得到，

(16)

(17)

其中k=1,…,N;i=0,…,N。從而動(dòng)力學(xué)微分方程約束轉(zhuǎn)化為如下的代數(shù)約束：

(18)

其中Xk≡X(τk)∈Rn，Uk≡U(τk)∈Rm。離散化后的動(dòng)力學(xué)方程約束只在LG點(diǎn)被滿(mǎn)足，而離散狀態(tài)初值為X0=X(-1)，末端時(shí)刻狀態(tài)值可通過(guò)Gauss積分近似得到

同樣，性能指標(biāo)也可通過(guò)高斯積分近似得到

(20)

其中wk為高斯權(quán)值。同時(shí)將式(11)中的邊值約束以及路徑約束分別離散為

以上性能指標(biāo)函數(shù)(20)及代數(shù)約束方程(18),(19)和(21)定義了一個(gè)非線性規(guī)劃問(wèn)題，該問(wèn)題的求解過(guò)程詳見(jiàn)參考文獻(xiàn)[13]，其解就是上述連續(xù)Bolza問(wèn)題的近似解。

4 改進(jìn)算法

以上所描述的Gauss偽譜法是全局配點(diǎn)法，若要提高擬合精度，則需要增加擬合多項(xiàng)式的階次。然而隨著多項(xiàng)式階次的增加，配點(diǎn)數(shù)也相應(yīng)增加，這使得收斂速度減慢，非線性規(guī)劃問(wèn)題的計(jì)算變得難以處理。

這里采用改進(jìn)的Gauss偽譜法，將全局配點(diǎn)法和局部配點(diǎn)法混合使用來(lái)提高擬合精度，即并不只是一味地提高擬合多項(xiàng)式的階次，而是考慮將軌跡分成S個(gè)小段，然后分別對(duì)各小段選擇合適階次的多項(xiàng)式進(jìn)行擬合。每段編號(hào)設(shè)為s∈[1,…,S]，假定第s小段用Ns階的多項(xiàng)式進(jìn)行擬合，其狀態(tài)擬合多項(xiàng)式Xs(τ)如式(13)所示。為達(dá)到提高擬合精度的目的，有2個(gè)方法可以考慮：將每一特定小段再繼續(xù)分段，或者增加該段擬合多項(xiàng)式的階次。而算法的關(guān)鍵在于如何對(duì)2種策略折衷選擇。

算法的基本思路描述如下：

(22)

(23)

則認(rèn)為小段s擬合偏差一致，此時(shí)通過(guò)增加該小段多項(xiàng)式階次Ns來(lái)提高擬合精度；否則，如果存在任何一個(gè)局部偏差最大點(diǎn)不滿(mǎn)足不等式(23)的約束，則將s小段在該點(diǎn)處進(jìn)行分段。這里，λ值的選取體現(xiàn)了局部和全局配點(diǎn)策略的一個(gè)權(quán)衡。當(dāng)λ值較小時(shí)，局部偏差最大點(diǎn)與偏差均值的比值則較易超過(guò)λ值而在某點(diǎn)處再進(jìn)行分段；當(dāng)λ值較大時(shí)，算法類(lèi)似采用了全局配點(diǎn)方法。本文仿真時(shí)采用文獻(xiàn)[9]中所提到的GPOPS優(yōu)化軟件，其具體分段數(shù)取決于擬合精度及λ值的大小，λ值一般取為3.5。

當(dāng)確定某一特定段內(nèi)偏差一致時(shí)，需要采用增加擬合多項(xiàng)式階次的策略。假設(shè)前2次配點(diǎn)數(shù)分別為Nk-1，Nk，且擬合誤差為ο(10-mk-1)，ο(10-mk)，而允許誤差為ο(ε)。假定mk>mk-1，則估計(jì)應(yīng)增加的配點(diǎn)數(shù)為Nk+1

(24)

該算法能自適應(yīng)地更新軌跡的分段數(shù)和各段擬合的配點(diǎn)數(shù)，從而在設(shè)計(jì)盡量少的配點(diǎn)數(shù)的情況下，達(dá)到需要的擬合精度。

5 仿真分析

本文以探月返回器為仿真對(duì)象，質(zhì)量9500kg，最大橫截面積23.8m2。

5.1 跳躍式再入軌跡優(yōu)化

具體仿真參數(shù)設(shè)置見(jiàn)表1，控制量σ約束在[-70°, 70°]之間變化，升阻比0.35，升力系數(shù)CL為0.44，阻力系數(shù)CD為1.25。軌跡狀態(tài)量、過(guò)載、熱流和控制量變化曲線如圖2～5所示。由圖3和4可見(jiàn)，設(shè)計(jì)的再入軌跡滿(mǎn)足過(guò)載和熱流約束。

表1 初末狀態(tài)設(shè)置和過(guò)程約束

圖2 狀態(tài)量變化軌跡

圖3 過(guò)載變化軌跡

圖4 駐點(diǎn)熱流變化軌跡

圖5 控制變量變化軌跡

5.2 跳躍式再入軌跡航程仿真分析

根據(jù)文獻(xiàn)[14]的再入走廊分析結(jié)果，這里在[-6.4°, -5.4°]之間均勻選擇6個(gè)初始再入角值進(jìn)行仿真，分別得到發(fā)生跳躍和不發(fā)生跳躍情況下的最大航程，見(jiàn)表2。為了清晰比對(duì)，這里給出了初始再入角為-6.0°時(shí)發(fā)生跳躍和不發(fā)生跳躍2種情況下的再入軌跡比較，如圖6所示。同時(shí)給出了初始再入角分別為-5.6°，-6.0°，-6.4°時(shí)，跳躍高度h為120km的飛行情況比較，如圖7所示，可知再入角越大(再入角絕對(duì)值)，飛行的最大距離越小，且過(guò)載峰值越大。

由表2中數(shù)據(jù)可見(jiàn)，再入軌跡在不發(fā)生跳躍時(shí)的最大航程僅有3000多公里，由圖6中的曲線對(duì)比發(fā)現(xiàn)當(dāng)跳躍高度為120km時(shí)再入軌跡最大航程遠(yuǎn)遠(yuǎn)大于不發(fā)生跳躍時(shí)的最大再入航程。為得出長(zhǎng)距離飛行是否必須采用跳躍式再入的結(jié)論，下面進(jìn)一步仿真分析在不同升阻比和存在不同程度大氣密度偏差時(shí)不發(fā)生跳躍的最大再入航程。

表2 兩種情況下最大航程比較

圖6 初始再入角為-6.0°時(shí)飛行情況比較

圖7 跳躍高度為120km各不同再入角時(shí)飛行情況比較

利用Gauss偽譜法計(jì)算得到不同升阻比、不同初始再入角情況下的最大航程和過(guò)載峰值，并將數(shù)據(jù)繪制曲線如圖8和9所示。由于航程隨著再入角的減小而增大，這里選擇再入角為-5.6°，仿真得到不同升阻比和不同大氣密度下的最大航程和過(guò)載峰值，并將數(shù)據(jù)繪制曲線如圖10和11所示。

由圖8～11可看出：再入飛行的最大航程，隨著初始再入角的增大而減小，隨著升阻比的增大而增加，隨著大氣密度的變大而減??；過(guò)載峰值則隨著初始再入角的增大而增大，隨著密度的增加，減小較快。而過(guò)載峰值與升阻比之間的關(guān)系與初始再入角的大小有關(guān)：當(dāng)初始再入角較大時(shí)，隨著升阻比的增大，過(guò)載峰值減小較快；當(dāng)初始再入角較小時(shí)，過(guò)載峰值反而隨著升阻比的增大而緩慢增大。

綜上可以得到最后結(jié)論，不發(fā)生跳躍的再入軌跡的最大航程不會(huì)超過(guò)3500km，則當(dāng)再入航程要求大于3500km時(shí)需采用跳躍式軌跡再入。

圖8 最大航程隨升阻比和再入角變化情況

圖9 過(guò)載峰值隨升阻比和再入角變化情況

圖10 最大航程隨大氣密度和升阻比變化情況

圖11 過(guò)載峰值隨大氣密度和升阻比變化情況

6 結(jié)論

本文首先基于改進(jìn)的Gauss偽譜法對(duì)探月返回器高速再入任務(wù)進(jìn)行了軌跡優(yōu)化設(shè)計(jì)，得到了滿(mǎn)足過(guò)載和熱流等強(qiáng)約束且總吸熱量最小的跳躍式再入軌跡。然后利用該方法優(yōu)化計(jì)算了在不同升阻比和存在不同程度大氣密度偏差的情況下，再入過(guò)程發(fā)生跳躍和不發(fā)生跳躍時(shí)的最大航程及過(guò)載峰值。通過(guò)比較可知，跳躍式再入軌跡能在較大程度上擴(kuò)大小升力體再入航程，并得出結(jié)論：當(dāng)再入航程要求大于約3500km時(shí)，需要采用跳躍式軌跡再入。

參 考 文 獻(xiàn)

[1] Wingrove R C. Survey of Atmosphere Re-Entry Guidance and Control Methods[J]. AIAA Journal, 1963, 1(9): 2019-2029.

[2] 李惠峰, 張蕊.探月飛船預(yù)測(cè)-校正再入制導(dǎo)律設(shè)計(jì)[J].空間控制技術(shù)與應(yīng)用, 2009, 35(1): 19-24. (LI Huifeng, ZHANG Rui. Design of Predictor-Corrector Reentry Guidance Law for Lunar Mission Spacecraft. Aerospace Control and Application, 2009, 35(1): 19-24.)

[3] Kuo Z S, Vinh N X. Improved Mathced Asymptotic Solutions for Three-Dimensional Atmospheric Skip Trajectories[J]. Journal of Spacecraft and Rockets, 1997, 34(4): 496-502.

[4] Kluever C A. Entry Guidance Using Analytical Atmospheric Skip Trajectories[J]. Journal of Guidance, Control, and Dynamics, 2008, 31(5): 1531-1534.

[5] García-Llama E. Analytic Development of a Reference Trajectory for Skip Entry[J].Journal of Guidance, Control, and Dynamics, 2011, 34(1): 311-317.

[6] Harpold J C, Graves C A. Shuttle Entry Guidance[J]. Journal of Astronautical Sciences, 1979, 37(3): 239-268.

[7] 姚寅偉, 李華濱.基于Guass偽譜法的高超聲速飛行器多約束三維再入軌跡優(yōu)化[J]. 航天控制, 2012, 30(2): 33-38. (YAO Yinwei, LI Huabin. The Generation of Three-dimensional Contrained Entry Trajectories for Hypersonic Vehicle Based on the Gauss Pseudospectral Method[J]. Aerospace Control, 2012, 30(2): 33-38.)

[8] Fahroo F, Ross I M. Trajetory Optimization by Indirect Spetral-olloation Methods[C].AIAA Astrodynamics Specialist Conference, Denver, CO, 2000.

[9] 高艾, 崔平遠(yuǎn), 崔祜濤.深空著陸器隨機(jī)優(yōu)化制導(dǎo)律設(shè)計(jì)[J].宇航學(xué)報(bào), 2011, 32(9): 1884-1889. (GAO Ai, CUI Pingyuan, CUI Hutao. Stochastically Optimized Guidance Law Design of the Deep Space Lander[J]. Journal of Astronautics, 2011, 32(9): 1884-1889.)

[10] Darby C L, Rao A V. A State Approximation- Based Mesh Refinement Algorithm for Solving Optimal Control Problems Using Pseudospectral Methods[C].AIAA Guidance, Navigation, and Control Conference, Chicago, Illinois, 2009.

[11] 南英, 陸宇平, 龔平.登月返回地球再入軌跡的優(yōu)化設(shè)計(jì)[J].宇航學(xué)報(bào), 2009, 30(5): 1842-1847. (NAN Ying, LU Yuping, GONG Ping. Optimal Reentry Trajectory Design for Mooncraft Returning to the Earth[J]. Journal of Astronautics, 2009, 30(5): 1842-1847.)

[12] Istratie V. Optimal Skip Entry into Atmophere with Minimum Heat and Constraints[C].AIAA Atmospheric Mechanics Conference and Exhibit, Providence, Rhode Island, 2004.

[13] Benson D A, et al. Direct Trajectory Optimization and Costate Estimation via an Orthogonal Collocation Method[J].Journal of Guidance, Control, and Dynamics, 2006, 29(6): 1435-1440.

[14] Souza S D, Sarigul-Klijn N. An Analytical Approach to Skip Earth Entry Guidance of a Low L/D Vehicle[C]. 46th AIAA Aerospace Sciences Meeting and Exhibit, Reno, Nevada, 2008.