火星大氣進(jìn)入段縱向可達(dá)區(qū)生成的解析同倫法

趙澤端，崔平遠(yuǎn)，朱圣英

(1. 北京理工大學(xué)宇航學(xué)院，北京 100081；2. 深空自主導(dǎo)航與控制工信部重點(diǎn)實(shí)驗(yàn)室，北京 100081；3.飛行器動力學(xué)與控制教育部重點(diǎn)實(shí)驗(yàn)室，北京 100081)

0 引 言

目前的火星著陸任務(wù)從進(jìn)入到最后著陸過程均采用與1976年海盜號(Viking)任務(wù)類似的著陸框架，著陸階段分為大氣進(jìn)入段、傘降段、動力下降段和著陸段[1]?；鹦谴髿膺M(jìn)入段指的是探測器從進(jìn)入火星大氣(距離火星表面約125 km)到降落傘展開的飛行階段。在進(jìn)入段，探測器需利用火星大氣消耗掉自身的絕大部分能量，將速度從進(jìn)入點(diǎn)約6 km/s減到開傘點(diǎn)約500 m/s，高度從進(jìn)入點(diǎn)約125 km 降至開傘點(diǎn)約10 km左右[2-3]。

可達(dá)區(qū)指探測器在同一初始狀態(tài)及任務(wù)參數(shù)約束下的所有可達(dá)終端狀態(tài)，與可控集一樣反應(yīng)了探測器的飛行能力[4-6]?？蛇_(dá)區(qū)分析對于任務(wù)設(shè)計(jì)、著陸點(diǎn)選取及風(fēng)險評估等均具有重要意義[7-8]。針對中高升阻比飛行器，Lu等[9]基于準(zhǔn)平衡滑翔假設(shè)及簡化的動力學(xué)方程，提出一種解析可達(dá)區(qū)求解方法，該方法可最終將可達(dá)區(qū)求解轉(zhuǎn)換成區(qū)間范圍已知的單參數(shù)尋優(yōu)問題，保證了求解的效率和可靠性。但準(zhǔn)平衡滑翔假設(shè)不適用于目前火星任務(wù)采用的小升阻比飛行器[10]。針對小升阻比飛行器，Saraf等[11]提出根據(jù)滿足路徑約束的最大最小阻力加速度剖面來得到最小最大航程，并通過對最大最小阻力加速度插值的方式得到介于最小最大航程之間的軌跡。該方法可以得到一個保守的可達(dá)區(qū)，是探測器實(shí)際可達(dá)區(qū)內(nèi)的子集。與此類似，Li等[12]通過將過程約束描述在r-V平面內(nèi)，可以快速獲得飛行器的可達(dá)區(qū)內(nèi)邊界(靠近初始位置)。Benito等[4]通過將可能的可達(dá)區(qū)域劃分成網(wǎng)格，對每一個網(wǎng)格點(diǎn)利用直接優(yōu)化方法檢測其可行性，最終得到可達(dá)區(qū)及每一個可達(dá)點(diǎn)對應(yīng)的最大開傘高度。該方法需要對每一個網(wǎng)格點(diǎn)單獨(dú)求取最優(yōu)化問題，計(jì)算量大，求解過程的收斂性及可靠性難以保證。Helsen等[13]通過假設(shè)一個比實(shí)際可達(dá)區(qū)更大的邊界，并利用直接優(yōu)化方法單獨(dú)求取離每一個邊界點(diǎn)最近的點(diǎn)，這些點(diǎn)組成飛行器的可達(dá)區(qū)。其思想與Lu等[9]采用的最優(yōu)接近思想相同，求取方法與Benito等[4]的方法相似，屬于直接優(yōu)化方法。

直接法通過將狀態(tài)和控制量離散，將原最優(yōu)問題轉(zhuǎn)換為非線性規(guī)劃(NLP)問題，借助商業(yè)軟件或工具箱求解，具有通用性高、問題建模相對容易的特點(diǎn)。但直接法的求解依賴于初值猜測，當(dāng)初值猜測不理想時，求解過程可能會陷入局部最優(yōu)解，不滿足最優(yōu)必要條件，也可能不收斂。間接法的理論基礎(chǔ)是龐特里亞金極小值原理，通過將原優(yōu)化問題轉(zhuǎn)換為一個等價的兩點(diǎn)邊值問題求解，對最優(yōu)必要條件的滿足較好。但由于協(xié)狀態(tài)缺乏物理意義，且問題對協(xié)狀態(tài)初值比較敏感，故間接法存在協(xié)態(tài)初值猜測困難的問題。

對于可以橫縱向解耦的小升阻比火星進(jìn)入探測器，基于其縱向動力學(xué)方程的最優(yōu)控制問題一般均屬于bang-bang控制，初值收斂域小且對初值敏感。Bertrand等[14]針對bang-bang最優(yōu)控制問題，應(yīng)用同倫法分析不同的初始輔助優(yōu)化問題對于解的收斂性和求解速度的影響。同倫法在求解小推力星際轉(zhuǎn)移軌道最優(yōu)問題方面已經(jīng)有了深入的研究[15]。在火星大氣進(jìn)入制導(dǎo)方面，Zheng等[16]利用同倫法求解了最大開傘高度問題。

本文利用解析同倫法思想，在縱向動力學(xué)背景下，考慮數(shù)值延拓的效率和可靠性，首先為無路徑約束的最大/最小開傘高度、最大/最小航程問題分別創(chuàng)建合適的輔助優(yōu)化問題。并從最大/最小開傘點(diǎn)高度問題出發(fā)，延拓出縱向可達(dá)區(qū)。由于利用解析同倫法思想，從已知最優(yōu)解的最優(yōu)輔助問題出發(fā)延拓到所關(guān)心的最優(yōu)問題，避免了初值猜測困難的問題；通過四個基本的最優(yōu)問題延拓出縱向可達(dá)區(qū)，避免了對可達(dá)區(qū)內(nèi)部點(diǎn)可行性的驗(yàn)證。

1 火星大氣進(jìn)入段動力學(xué)模型

由于火星大氣進(jìn)入段持續(xù)時間較短，可以忽略火星自轉(zhuǎn)的影響，在火星固連系下建立三自由度動力學(xué)模型如下

(1)

(2)

式中：ρ為大氣密度，CD為阻力系數(shù)，Sref為防熱底的參考面積，m為探測器質(zhì)量。此處需強(qiáng)調(diào)的是，除系統(tǒng)狀態(tài)變量外，其他變量若未特別說明外均未做無量綱化處理。

由于目前的探測器在火星大氣進(jìn)入過程中均采用平衡攻角，為簡化分析過程，可近似將升阻比當(dāng)成常數(shù)，故可得到無量綱的升力加速度為

L=D(L/D)

(3)

式中：L/D為升阻比。在本文的研究中，火星大氣密度模型采用如下的指數(shù)模型

ρ=ρ0e-h/hs

(4)

式中：ρ0為參考大氣密度，h為探測器當(dāng)前高度，hs為參考高度。

由于目前所有火星著陸任務(wù)的探測器進(jìn)入級均采用半錐角為70°的球錐形防熱大底，該氣動外型具有較小升阻比(升阻比小于0.3)，進(jìn)入段攻角保持在平衡攻角附近，通過調(diào)整傾側(cè)角的大小和方向來控制探測器的運(yùn)動，故可在制導(dǎo)相關(guān)設(shè)計(jì)中將縱向和側(cè)向解耦分析。建立縱向動力學(xué)模型如下

(5)

式中：s為運(yùn)動軌跡長度在水平方向上的投影，即本文中的航程，用來近似替代水平位置。

考慮到火星大氣進(jìn)入段的擾動和不確定性，軌跡分析和設(shè)計(jì)時需留有足夠的控制余量σ∈[σmin,σmax]，其中π>σmax>σmin>0。本文控制量取為傾側(cè)角余弦u=cosσ，有如下的上下界約束

cosσmax≤u≤cosσmin

(6)

本文僅研究無路徑約束下的最大/最小航程、最大/最小開傘高度及縱向可達(dá)區(qū)問題，有路徑約束的相關(guān)問題可在本文的基礎(chǔ)上，利用精確罰函數(shù)法和同倫法求解[14]。

2 路徑自由的縱向可達(dá)區(qū)邊界求解

本文所定義的縱向可達(dá)區(qū)主要是指同一初始條件下，探測器可以飛達(dá)的航程-高度范圍。進(jìn)入段的終端對應(yīng)降落傘的展開條件，本文采取速度作為開傘觸發(fā)條件，可達(dá)區(qū)分析時需滿足的初始和終端條件如下

(7)

V(τf)=Vf

(8)

對于可達(dá)高度-航程剖面而言，其邊界指的是最大/最小航程和最大/最小開傘高度。

2.1 最大/最小航程問題

2.1.1最大/最小航程原問題

最大航程對應(yīng)探測器可飛達(dá)的最遠(yuǎn)距離，對應(yīng)優(yōu)化問題的性能指標(biāo)為

Js,max=-s(τf)

(9)

最小航程對應(yīng)探測器可飛達(dá)的最近距離，對應(yīng)優(yōu)化問題的性能指標(biāo)為

Js,min=s(τf)

(10)

則路徑自由的最大/最小航程問題可描述如下：

(11)

同時滿足式(5)～式(8)。式(11)中，Js取-s(τf)或s(τf)時，分別對應(yīng)最大航程問題Ps,max或最小航程問題Ps,min。

Hs=λTf

(12)

由歐拉-拉格朗日方程可得：

(13)

由橫截條件及式(7)和式(8)可得：

(14)

由于終端時間自由且哈密頓函數(shù)中不顯含時間，故

Hs(τ)≡0

(15)

由于控制變量u在哈密頓函數(shù)中是一階形式，根據(jù)極小值原理，最優(yōu)控制將是bang-bang的形式，切換條件為λγL/V，具體的切換過程為

(16)

由分析可知，λγL/V只可能在某些單點(diǎn)為零，故問題Ps并非嚴(yán)格意義上的奇異。從簡化問題的角度考慮，可將λγL/V=0時的控制并入λγL/V>0的情況中。

綜上，問題Ps可轉(zhuǎn)換為如下的兩點(diǎn)邊值問題：

問題Ts：尋找由協(xié)狀態(tài)初值和終端飛行時間組成的變量

(17)

滿足式(7)、(8)和(14)描述的邊值條件，式(5)和(13)描述的微分方程條件，式(15)中的哈密頓函數(shù)條件及式(16)中的控制切換條件。

Ts,max和Ts,min分別對應(yīng)最大/最小高度問題兩點(diǎn)邊值問題。由分析可知，問題Ts的控制為bang-bang形式，問題的求解將會對初值猜測特別敏感。由于協(xié)狀態(tài)初值沒有實(shí)際的物理意義，初值猜測比較困難。本文將構(gòu)造最優(yōu)解已知的輔助優(yōu)化問題，通過解析同倫法求解該問題。

2.1.2最大/最小航程輔助優(yōu)化問題

考慮到最大/最小航程問題的bang-bang最優(yōu)控制特點(diǎn)，若輔助優(yōu)化問題與原問題的最優(yōu)解相差太遠(yuǎn)，則存在延拓失敗的危險，即不能從輔助問題最優(yōu)解最終得到原問題最優(yōu)解。即使可以得到，也存在延拓過程過長，迭代計(jì)算量過大的缺點(diǎn)。故通過對原優(yōu)化問題的初步分析得出更貼近原問題的輔助優(yōu)化問題，可加快同倫迭代過程，相應(yīng)減少計(jì)算量，并在一定程度上保證求解過程的穩(wěn)定性。

對于火星大氣進(jìn)入段探測器，當(dāng)控制量u取最大值，即傾側(cè)角σ取最小值時，向上的升力將使得探測器在更高的高度飛行，飛行過程中較小的阻力將使得自身能量消耗速度更慢，探測器將能飛得更遠(yuǎn)，故取最大航程的輔助優(yōu)化問題性能指標(biāo)為

(18)

當(dāng)探測器保持在最大傾側(cè)角附近時，由于升力方向向下，探測器的高度降低較快，飛行過程中阻力較大，則其能量消耗得較快，飛得越近。故構(gòu)造最小航程輔助優(yōu)化問題性能指標(biāo)為

(19)

最大/最小航程的輔助優(yōu)化問題可統(tǒng)一表達(dá)為：

(20)

(21)

(22)

(23)

λV(τf)=0

(24)

根據(jù)式(22)和式(24)的協(xié)狀態(tài)終端值、式(13)的協(xié)狀態(tài)微分方程以及正向積分得到的最優(yōu)狀態(tài)和飛行時間逆向積分即可得到輔助問題最優(yōu)軌跡對應(yīng)的協(xié)狀態(tài)值及最優(yōu)控制量。

2.1.3最大/最小航程問題同倫求解過程

構(gòu)建如下的優(yōu)化問題：

(25)

該問題的哈密頓函數(shù)為

(26)

(27)

控制量u仍以一階形式出現(xiàn)，故該問題的最優(yōu)解仍是bang-bang形式。經(jīng)分析可得，最優(yōu)解具有如下的形式

(28)

其中，對λγL/V+ε-1=0情形的分析與問題Ps中的分析過程相似。下面不再贅述。

協(xié)狀態(tài)的微分方程與問題Ps相同。橫截條件如下：

(29)

2.2 最大/最小開傘高度問題

2.2.1最大/最小開傘高度原問題

由于火星大氣相對稀薄，在進(jìn)入段探測器飛行高度較高，大氣不能為其提供充足的減速能力。為了能為傘降段及后續(xù)的著陸段提供足夠的時間和空間，同時也為了能探測火星南極海拔更高的古大陸，提高探測器進(jìn)入段末端的開傘高度一直是進(jìn)入段制導(dǎo)控制研究的重點(diǎn)。對于最大開傘高度問題，取性能指標(biāo)如下

Jh,max=-r(τf)

(30)

實(shí)際工程設(shè)計(jì)中，可能的最低開傘高度影響著后續(xù)故障預(yù)案的設(shè)計(jì)，以及標(biāo)稱航程在可達(dá)區(qū)內(nèi)的選擇，具有十分重要的意義。對于最小開傘高度問題，取性能指標(biāo)如下

Jh,min=r(τf)

(31)

則路徑自由的最大/最小化開傘高度優(yōu)化問題可描述如下

(32)

同時滿足：式(5)～式(8)。式(32)中，Jh取-r(τf)或r(τf)時分別對應(yīng)最大高度問題Ph,max或最小開傘高度問題Ph,min。

(33)

問題Ph的求解可轉(zhuǎn)換為求解兩點(diǎn)邊值問題Th：尋找由協(xié)狀態(tài)初值和終端飛行時間組成的變量z，滿足式(7)、(8)和(33)描述的邊值條件，式(5)和(13)描述的微分方程條件，式(15)中的哈密頓函數(shù)條件及式(16)中的控制切換條件。

2.2.2最大/最小開傘高度輔助優(yōu)化問題

(34)

同時滿足式(5)～式(8)。

2.2.3最大/最小開傘高度問題同倫求解過程

(35)

(36)

(37)

由于積分項(xiàng)中不顯含狀態(tài)變量，故該問題的協(xié)狀態(tài)微分方程與問題Ps相同，橫截條件為

(38)

由于式(36)中，控制量u以二階的形式出現(xiàn)，由二次函數(shù)在u∈[umin,umax]的最小值分析可知，該問題的最優(yōu)控制結(jié)構(gòu)如下

(39)

3 縱向可達(dá)區(qū)求解

為了從第2節(jié)的子優(yōu)化問題得到縱向可達(dá)區(qū)，取性能指標(biāo)為Jsk，構(gòu)建如下優(yōu)化模型。

(40)

同時滿足式(5)～式(8)以及式(41)

s(τf)=sk

(41)

當(dāng)Jsk=-r(τf)時，對應(yīng)求解可達(dá)航程及其最大高度邊界；當(dāng)Jsk=r(τf)對應(yīng)求解可達(dá)航程范圍及最低高度邊界。

相比于問題Ph，問題RA多了一個終端航程約束。問題RA的哈密頓函數(shù)與問題Ph相同，且對于時間保持恒為零值，故二者的協(xié)狀態(tài)微分方程相同。問題RA的橫截條件為

(42)

構(gòu)造兩個單調(diào)航程序列：{s0,…,si,si+1,…,sI}與{s0,…,sj,sj+1,…,sJ}，滿足s0為問題Ph對應(yīng)的航程，sI為問題Ps,min對應(yīng)的航程，sJ為問題Ps,max對應(yīng)的航程，分別從問題Ph的最優(yōu)解出發(fā)，按照序列依次求解問題RA，其中第k次的最優(yōu)解zk作為第k+1次問題求解初值，最終可得到縱向可達(dá)區(qū)。

4 仿真校驗(yàn)

本節(jié)探測器的物理參數(shù)、初始條件及開傘條件部分根據(jù)文獻(xiàn)[15]選取。其中質(zhì)量m為2800 kg，參考面積Sref為15.9 m2，升阻比L/D為0.24，阻力系數(shù)CD為1.45，初始條件及開傘條件具體如表1所示。

表1 探測器初始條件及開傘條件Table 1 Initial and terminal conditions of the vehicle

為了給制導(dǎo)控制系統(tǒng)留有控制余量，取最大傾側(cè)角σmax為120°，最小傾側(cè)角σmin為30°。在同倫延拓過程中，以下4個優(yōu)化子問題的ε均從0開始，其中ε=0代表輔助優(yōu)化問題的解。具體仿真結(jié)果如下。

4.1 最大航程問題

相比于σ≡σmin的進(jìn)入軌跡，航程在ε由0變化到1時增加了約150 m，出現(xiàn)這種現(xiàn)象的主要原因在于飛行時間的增加，相比于ε由0.1變化到1的絕大部分過程，最大航程問題對應(yīng)的末端航程僅增加了3 m，如圖2所示。故在任務(wù)的快速分析設(shè)計(jì)中，可用σ≡σmin的進(jìn)入軌跡近似代替最大航程軌跡來分析最大航程。

圖1 問題中初始協(xié)狀態(tài)值隨參數(shù)ε的變化Fig.1 Variation of initialcostates with ε in

圖2 問題中航程隨參數(shù)ε的變化關(guān)系Fig.2 Variation of downrange withε in

4.2 最小航程問題

最小航程軌跡比傾側(cè)角保持為σ≡σmax的進(jìn)入軌跡末端航程高135 m，如圖4所示，這種現(xiàn)象有可能是輔助問題求解時采取0.5 s步長，積分精度較低引起的。在任務(wù)的快速分析設(shè)計(jì)中，可用σ≡σmax的進(jìn)入軌跡近似代替最小航程軌跡。

圖3 問題中初始協(xié)狀態(tài)值隨參數(shù)ε的變化Fig.3 Variation of initialcostates with ε in

圖4 問題中航程隨參數(shù)ε的變化關(guān)系Fig.4 Variation of downrange with ε in

4.3 最大開傘高度問題

圖5 問題中初始協(xié)狀態(tài)值隨參數(shù)ε的變化Fig.5 Variation of initialcostates with ε in

圖6 問題中開傘高度隨參數(shù)ε的變化關(guān)系Fig.6 Variation of terminal altitude with ε in

4.4 最小開傘高度問題

圖7 問題中初始協(xié)狀態(tài)值隨參數(shù)ε的變化Fig.7 Variation of initialcostates with ε in

圖8 問題中開傘高度隨參數(shù)ε的變化關(guān)系Fig.8 Variation of terminal altitude with ε in

4.5 縱向可達(dá)區(qū)

另外選取兩個不同的初始航跡角：γ0=-14°和γ0=-13°。在分別求取4個子優(yōu)化問題的基礎(chǔ)上分別求解問題RA。其中γ0=-14°和γ0=-13°的子問題既可以分別從輔助優(yōu)化問題出發(fā)求解，也可以從γ0=-12.5°的4個子問題最優(yōu)解出發(fā)，構(gòu)造{γk1,γk2,…,γkN}的序列，采用同倫法求解，其中γk1=-12.5°,γkN為-14°或-13°，也可以是其他任務(wù)設(shè)計(jì)時感興趣的初始路徑角。3個不同初始路徑角下，無路徑約束的縱向可達(dá)區(qū)如圖9所示。從圖9可以看出，進(jìn)入角更陡時，航程范圍更窄。初始進(jìn)入角對航程的影響更大，對于可達(dá)的最大/小開傘高度影響較小。

圖9 不同初始航跡角下的縱向可達(dá)區(qū)Fig.9 Reachable area with different γ(τ0)

5 結(jié) 論

針對火星大氣進(jìn)入段飛行器的飛行能力分析問題，提出一種基于解析同倫法的縱向可達(dá)區(qū)(開傘點(diǎn)高度-航程剖面)生成方法。在縱向動力學(xué)背景下，考慮數(shù)值延拓的效率和可靠性，首先為無路徑約束的最大/最小開傘高度、最大/最小航程問題分別構(gòu)造最優(yōu)解已知的初始輔助優(yōu)化問題，并以此為出發(fā)點(diǎn)最終求取原最優(yōu)問題的解。在求解4個子優(yōu)化問題的基礎(chǔ)上，通過構(gòu)造合適的同倫參數(shù)，延拓出縱向可達(dá)區(qū)，避免了對可達(dá)區(qū)內(nèi)部點(diǎn)可行性的驗(yàn)證。仿真結(jié)果表明，該方法有效避免了最優(yōu)問題協(xié)態(tài)初值猜測困難的問題，同倫求解過程穩(wěn)定?？梢钥焖偾笕∪我獬跏悸窂浇窍碌哪┒俗畲?最小航程問題、最大/最小開傘高度問題以及縱向可達(dá)區(qū)問題。未來可在本文研究的基礎(chǔ)上，結(jié)合精確罰函數(shù)方法，進(jìn)一步求得存在路徑約束時的縱向可達(dá)區(qū)。