劉 旭，葉 松，林子瑞，黃翔宇，李 爽

(1.南京航空航天大學(xué)航天學(xué)院，南京 211106；2.北京航天自動(dòng)控制研究所，北京 100854；3.北京控制工程研究所，北京 100094)

0 引 言

火星大氣進(jìn)入是指探測(cè)器進(jìn)入火星大氣層邊界到超聲速減速傘展開(kāi)的這一階段，時(shí)間持續(xù)約4 min，期間探測(cè)器的速度將從5～7 km/s減速到馬赫數(shù)2.2以下。進(jìn)入過(guò)程中，探測(cè)器應(yīng)滿足熱流密度、動(dòng)壓、過(guò)載等約束以保證順利開(kāi)傘，同時(shí)也應(yīng)滿足開(kāi)傘點(diǎn)的高度和經(jīng)緯度，從而確保后續(xù)飛行任務(wù)的時(shí)間和空間余量。以“火星科學(xué)實(shí)驗(yàn)室”為例，任務(wù)要求減速傘展開(kāi)時(shí)速度為1.4～2.6，動(dòng)壓為250～850 Pa，開(kāi)傘點(diǎn)高度不低于5 km，過(guò)載峰值低于13，熱流峰值低于100 W·cm。因此，多約束條件下末端高度最大化問(wèn)題是火星大氣進(jìn)入軌跡優(yōu)化的重點(diǎn)之一。

火星大氣進(jìn)入軌跡優(yōu)化主要采用間接法和直接法。間接法方面，鄭藝裕將無(wú)約束火星進(jìn)入末端高度最大化問(wèn)題的解作為初始猜值，并采用精確罰函數(shù)方法將路徑約束處理為等式約束后增廣到目標(biāo)函數(shù)，進(jìn)而將路徑約束問(wèn)題轉(zhuǎn)化為無(wú)約束問(wèn)題求解。Long等使用間接法確定火星進(jìn)入末端高度最大化問(wèn)題的最優(yōu)傾側(cè)角剖面為Bang-Bang結(jié)構(gòu)，并且至少切換兩次。直接法則直接離散狀態(tài)量和控制量，將最優(yōu)控制問(wèn)題離散為非線性規(guī)劃問(wèn)題進(jìn)行求解，求解方法包括直接配點(diǎn)法、正交配點(diǎn)法(即偽譜法)、凸優(yōu)化法等。偽譜法包括Legendre偽譜法、Chebyshev偽譜法、Gauss偽譜法以及Radau偽譜法等。直接配點(diǎn)法方面，Zhao等采用廣義二分網(wǎng)格布置配點(diǎn)，結(jié)合網(wǎng)格細(xì)化和稀疏差分技術(shù)，將火星進(jìn)入軌跡優(yōu)化問(wèn)題轉(zhuǎn)化為稀疏非線性規(guī)劃問(wèn)題后快速求解。

目前，凸優(yōu)化方法因其具有多項(xiàng)式計(jì)算復(fù)雜度和理論全局最優(yōu)性，因而在航空航天領(lǐng)域得到廣泛關(guān)注，應(yīng)用場(chǎng)景包括無(wú)人機(jī)編隊(duì)飛行、行星著陸、火箭發(fā)射和回收、運(yùn)載器制導(dǎo)等。在大氣進(jìn)入問(wèn)題方面，Liu等將傾側(cè)角的正、余弦值作為控制量，并引入二階錐約束，從而將非凸問(wèn)題轉(zhuǎn)化為凸問(wèn)題，并采用序列凸優(yōu)化方法求解。Wang等則引入傾側(cè)角速度作為控制量，實(shí)現(xiàn)控制量和狀態(tài)量解耦，從而消除傾側(cè)角高頻振蕩。Zhou等則引入多分辨率技術(shù)，對(duì)均勻節(jié)點(diǎn)進(jìn)行網(wǎng)格細(xì)化，提高了凸優(yōu)化收斂解的數(shù)值精度。此外，Sagliano等結(jié)合Flipped Radau偽譜法和凸優(yōu)化方法，采用偽譜離散格式提高問(wèn)題離散精度，有效地提高了算法的收斂速度。Wang等則基于Radau偽譜凸優(yōu)化方法，進(jìn)一步給出根據(jù)線性化誤差動(dòng)態(tài)調(diào)整信賴域大小的算法，提高了序列凸優(yōu)化方法的收斂性。

以上軌跡優(yōu)化方法均為確定性優(yōu)化方法，而面對(duì)真實(shí)環(huán)境中的大量不確定性，確定性優(yōu)化方法的局限性日益顯現(xiàn)。為解決不確定條件下的軌跡優(yōu)化問(wèn)題，NASA已經(jīng)在航空航天領(lǐng)域開(kāi)展了不確定性量化研究，該方法主要環(huán)節(jié)包括不確定性建模、不確定性分析和不確定性優(yōu)化。不確定性建模時(shí)，根據(jù)概率信息的完備情況，可以考慮隨機(jī)不確定性和認(rèn)知不確定性兩類，其中隨機(jī)不確定性的建?？梢圆捎酶怕收摲椒?，而認(rèn)知不確定性由于無(wú)法準(zhǔn)確給出概率密度模型或高階統(tǒng)計(jì)矩的概率信息，因而只能用上下界表示。不確定性分析的主要目的是估算系統(tǒng)輸出不確定性分布的期望或方差，常見(jiàn)的方法主要有蒙特卡洛法、敏感度矩陣法、協(xié)方差分析法和廣義多項(xiàng)式混沌法、隨機(jī)配點(diǎn)法、響應(yīng)面法、拉丁超立方采樣法等。而不確定性優(yōu)化的研究主要集中于魯棒優(yōu)化和可靠性優(yōu)化方法，其中魯棒優(yōu)化要求在最差工況下仍能滿足任務(wù)約束并最小化目標(biāo)函數(shù)，而可靠性優(yōu)化則給出最小化目標(biāo)函數(shù)時(shí)違反任務(wù)約束的概率。在大氣進(jìn)入軌跡不確定性量化方面，Halder等，Prabhakar等較早開(kāi)展了相關(guān)研究，國(guó)內(nèi)崔平遠(yuǎn)、李海陽(yáng)、李爽等人的課題組也開(kāi)展了相關(guān)工作。

大氣進(jìn)入問(wèn)題建模時(shí)通常以時(shí)間或能量為自變量，但以時(shí)間為自變量時(shí)，需要給定或者優(yōu)化末端時(shí)間，問(wèn)題求解難度較大。而以能量為自變量時(shí)，不僅忽略中心天體的自轉(zhuǎn)角速度，導(dǎo)致模型精度損失，同時(shí)要求已知末端高度和速度來(lái)確定末端能量，因此不適用于末端高度最大化問(wèn)題。盡管李俊等提出以弧長(zhǎng)為自變量建立大氣進(jìn)入模型，但仍需要對(duì)末端弧長(zhǎng)進(jìn)行尋優(yōu)。本文提出以縱向航程角為自變量建立火星進(jìn)入模型，將末端時(shí)間自由問(wèn)題轉(zhuǎn)化為末端縱向航程角固定問(wèn)題，并結(jié)合Legendre偽譜離散格式和序列凸優(yōu)化方法，將火星大氣進(jìn)入末端高度最大化問(wèn)題轉(zhuǎn)化為凸優(yōu)化問(wèn)題后求解。

1 問(wèn)題描述

1.1 火星大氣進(jìn)入模型

假設(shè)火星大氣相對(duì)火星表面靜止，且火星為均勻球體，則彈道升力式火星探測(cè)器在火星大氣中的無(wú)動(dòng)力飛行過(guò)程可以用一組微分方程表示：

(1)

(2)

(3)

sincoscos)

(4)

(5)

(6)

(7)

式中：、、和分別表示火星探測(cè)器的升力系數(shù)、阻力系數(shù)、參考面積和質(zhì)量，表示火星大氣密度，表達(dá)式為：

=e-

(8)

式中：為火星表面大氣密度;為飛行高度;為標(biāo)準(zhǔn)大氣模型的參考高度;e為自然對(duì)數(shù)的底數(shù)。

火星大氣進(jìn)入問(wèn)題的典型約束條件包括過(guò)程約束和邊界約束，過(guò)程約束為熱流密度、動(dòng)壓和過(guò)載約束：

(9)

()=,()=,≤≤

(10)

式中：狀態(tài)量=[,,,,,],和表示給定的初始和末端狀態(tài);和表示各狀態(tài)量的變化范圍上下限。飛行任務(wù)一般要求滿足全部初始約束，末端約束滿足部分即可。

此外,傾側(cè)角的大小和角速度也有一定變化范圍,定義控制約束如下:

(11)

=-

(12)

1.2 新的自變量

式(1)～式(6)為常用的大氣進(jìn)入動(dòng)力學(xué)模型，自變量為無(wú)量綱時(shí)間，但在求解火星進(jìn)入末端高度最大化問(wèn)題時(shí)需要同時(shí)對(duì)末端時(shí)間進(jìn)行優(yōu)化，求解難度較大。此外，常用的模型還包括以負(fù)的比機(jī)械能(簡(jiǎn)稱能量)=1-2為自變量的模型，但根據(jù)其定義可知需要已知末端高度和速度才能確定末端能量，而在火星進(jìn)入問(wèn)題中，末端高度和速度在一定區(qū)間內(nèi)即可，而非固定值，因此該模型不適合末端高度最大化問(wèn)題。同時(shí)基于能量的模型需要忽略火星自轉(zhuǎn)才能獲得簡(jiǎn)潔的動(dòng)力學(xué)模型，即假設(shè)dd≈，因此動(dòng)力學(xué)模型精度有所下降。

圖1 縱向航程和縱向航程角示意圖

縱向航程角的初始值為0，末端值則可以根據(jù)初始和末端時(shí)刻的經(jīng)緯度確定：

cos=sinsin+coscoscos(-)

(13)

同時(shí)，縱向航程角始終單調(diào)遞增，相對(duì)時(shí)間的微分表達(dá)式也很簡(jiǎn)潔，因此本文選擇作為新的自變量建立火星進(jìn)入動(dòng)力學(xué)模型。相比基于時(shí)間的模型，本文建模方法可以將末端時(shí)間自由問(wèn)題轉(zhuǎn)化為末端縱向航程角固定問(wèn)題，避免對(duì)末端時(shí)間進(jìn)行優(yōu)化。同時(shí)，相比基于能量的模型，不需要事先已知末端高度和速度，也不需要忽略火星自轉(zhuǎn)項(xiàng)，模型精度得以保證。以縱向航程角為自變量的狀態(tài)微分方程為：

(14)

1.3 狀態(tài)量擴(kuò)充

注意到式(14)中不包含飛行時(shí)間，因此本文將無(wú)量綱時(shí)間擴(kuò)充為狀態(tài)變量，所以本文方法也可以通過(guò)定義時(shí)間最優(yōu)目標(biāo)函數(shù)=來(lái)求解時(shí)間最優(yōu)問(wèn)題。引入飛行時(shí)間后，狀態(tài)量擴(kuò)展為：

=[,,,,,,]

(15)

此外，選擇傾側(cè)角的導(dǎo)數(shù)作為控制變量可以抑制傾側(cè)角的高頻振蕩，則狀態(tài)微分方程可改寫(xiě)為控制仿射形式：

′=()+()+()

(16)

式中：=[,,,,,,,],=dd，而()=[,,…,],()=[,,…,]和()=[ 1, 2,…, 8]的表達(dá)式分別為：

(17)

(18)

(19)

因此，火星進(jìn)入末端高度最大化問(wèn)題P可以表述為：

(20)

s.t.(9),(10),(11),(16)

(21)

2 偽譜凸優(yōu)化方法

Legendre偽譜凸優(yōu)化(Legendre pseudospectral convex programming, LPCP)方法結(jié)合偽譜法和凸優(yōu)化方法，將狀態(tài)量和控制量在(Legendre-Gauss-Lobatto,LGL)點(diǎn)處離散，并通過(guò)Lagrange插值多項(xiàng)式逼近狀態(tài)量和控制量，從而將狀態(tài)微分方程和目標(biāo)函數(shù)中的積分運(yùn)算轉(zhuǎn)化為代數(shù)運(yùn)算，再結(jié)合線性化方法，將最優(yōu)控制問(wèn)題轉(zhuǎn)化為凸優(yōu)化問(wèn)題，最后通過(guò)求解凸優(yōu)化問(wèn)題得到原始最優(yōu)控制問(wèn)題的近似解。本文直接給出LPCP方法的一般步驟，有關(guān)Legendre偽譜法和凸優(yōu)化方法的內(nèi)容詳見(jiàn)文獻(xiàn)[9]和文獻(xiàn)[35]。

2.1 一階線性化

凸優(yōu)化問(wèn)題是指目標(biāo)函數(shù)和約束條件都為凸函數(shù)的最優(yōu)化問(wèn)題，而火星進(jìn)入末端高度最大化問(wèn)題的非凸性來(lái)源于狀態(tài)微分方程和過(guò)程約束。因此，本文采用一階泰勒展開(kāi)方法將這兩項(xiàng)約束線性化。

首先將狀態(tài)微分方程(16)在參考軌跡處線性化：

(22)

式中：()中元素(,=1,2,…,8)的表達(dá)式為：

(23)

(24)

(25)

(26)

(27)

其余元素均為0，此外：

=-,=

(28)

=-,=

(29)

同理，式(9)中的過(guò)程約束可線性化為不等式約束：

(30)

(31)

(32)

(33)

至此，問(wèn)題P可以近似為連續(xù)凸優(yōu)化問(wèn)題P：

(34)

s.t.(10),(11),(22),(30)

(35)

2.2 偽譜離散化

一般的序列凸優(yōu)化(Sequential convex programming, SCP)方法采用均勻節(jié)點(diǎn)離散問(wèn)題，并通過(guò)梯形積分處理狀態(tài)微分方程和目標(biāo)函數(shù)中的積分項(xiàng)。為保證離散精度，通常需要布置較多離散點(diǎn)，因此問(wèn)題規(guī)模較大。而偽譜法將問(wèn)題在一系列全局正交配點(diǎn)上離散，并通過(guò)Gauss積分處理狀態(tài)微分方程和目標(biāo)函數(shù)的積分項(xiàng)，可以用較少的離散節(jié)點(diǎn)保證離散精度。本文采用Legendre偽譜離散格式將問(wèn)題在LGL節(jié)點(diǎn)處離散。

(36)

式中：()為L(zhǎng)agrange插值基函數(shù),且有：

(37)

(38)

(39)

由此可將狀態(tài)微分方程轉(zhuǎn)化為+1組在LGL配點(diǎn)處的等式約束，但需要注意先將問(wèn)題P的定義域變換到區(qū)間[-1,+1]內(nèi)：

(40)

(41)

為了保證算法收斂時(shí)虛擬控制盡可能小，需要在目標(biāo)函數(shù)中施加懲罰項(xiàng)：

(42)

式中：為給定的虛擬控制懲罰權(quán)重，且懲罰函數(shù)選擇1范數(shù)是為了讓虛擬控制中具有盡可能多的0元素，也可以選擇2范數(shù)作為懲罰函數(shù)，即保證全部元素盡可能地趨近0。

注意到問(wèn)題P的目標(biāo)函數(shù)中并不含有積分項(xiàng)，因此并不需要進(jìn)行處理。若含有積分項(xiàng)，可使用Gauss-Lobatto求積公式進(jìn)行近似。

1.3.5 記錄手術(shù)時(shí)長(zhǎng)、術(shù)中液體出入量和術(shù)中發(fā)生不良反應(yīng)如低血壓、低血氧、牽拉反射引起的牽拉痛、惡心、嘔吐的例數(shù)等。

(43)

式中：Δ為給定的信賴域半徑。

最后，算法收斂時(shí)，應(yīng)滿足前后兩次迭代的解差別不大，同時(shí)虛擬控制足夠小，則可以定義收斂標(biāo)準(zhǔn)如下：

(44)

式中：和為給定的收斂標(biāo)準(zhǔn)。

注意到過(guò)程約束、以及狀態(tài)量或控制量相關(guān)的約束不涉及微分或積分項(xiàng)，因此離散化后約束表達(dá)式不變，只需將連續(xù)變量改為相應(yīng)的離散變量即可。

至此，離散Legendre偽譜凸優(yōu)化問(wèn)題P可定義為：

(45)

s.t.(10),(11),(30),(41),(43)

(46)

2.3 序列凸優(yōu)化

綜上所述，Legendre偽譜凸優(yōu)化算法的具體步驟為：

3 仿真校驗(yàn)

表1 仿真參數(shù)設(shè)置

根據(jù)上述仿真參數(shù)，對(duì)比了本文LPCP方法、自適應(yīng)偽譜法、以及施加虛擬控制的SCP方法。其中LPCP設(shè)置56個(gè)離散節(jié)點(diǎn)，SCP設(shè)置201個(gè)離散節(jié)點(diǎn)，二者的初始猜測(cè)軌跡均為-35°常值傾側(cè)角積分所得軌跡，懲罰權(quán)重=1×10，信賴域半徑Δ=0.5，收斂標(biāo)志=0.05，=1×10，且LPCP和SCP均采用Yalmip工具箱和Mosek求解器進(jìn)行仿真。自適應(yīng)偽譜法采用GPOPS軟件求解，設(shè)置求解器為IPOPT，網(wǎng)格收斂精度為1×10。仿真結(jié)果如圖2～圖5所示。

圖2 經(jīng)度-緯度剖面和速度-高度剖面

圖3 航跡傾角和航跡偏角剖面

圖4 時(shí)間-傾側(cè)角、傾側(cè)角速度剖面

圖5 過(guò)程約束剖面

同時(shí)，GPOPS產(chǎn)生的最優(yōu)軌跡雖然和LPCP、SCP的最優(yōu)軌跡一致，但傾側(cè)角、航跡傾角和航跡偏角剖面和其余兩個(gè)方法有一定區(qū)別，同時(shí)傾側(cè)角速度剖面的振蕩也較大。從算法上來(lái)說(shuō)，GPOPS將最優(yōu)控制問(wèn)題轉(zhuǎn)化為非線性規(guī)劃問(wèn)題，而LPCP和SCP則是通過(guò)序列凸優(yōu)化逼近原始最優(yōu)控制問(wèn)題，且離散格式也不相同，因此兩類方法近似解存在一定差異。

為了進(jìn)一步說(shuō)明本文LPCP算法的優(yōu)勢(shì)，表2中給出了三種方法的結(jié)果對(duì)比。從表2中可知，LPCP的計(jì)算耗時(shí)最少，為 6.195 s，而SCP和GPOPS分別耗時(shí)10.92 s和48.66 s。對(duì)比LPCP和SCP，兩種方法的差別僅有離散格式不同，且LPCP的離散節(jié)點(diǎn)數(shù)(56個(gè))遠(yuǎn)少于SCP的離散節(jié)點(diǎn)數(shù)(201個(gè))，這表明Legendre離散格式相比于均勻離散格式能夠加速算法收斂，同時(shí)LPCP的目標(biāo)函數(shù)值(10.868 km)和SCP的目標(biāo)函數(shù)值(10.853 km)非常接近。對(duì)比LPCP和GPOPS，二者分別采用Legendre和Radau偽譜離散格式，盡管GPOPS的目標(biāo)函數(shù)值(11.17 km)優(yōu)于LPCP(10.868 km)，但差距不大，約300 m，而LPCP的計(jì)算耗時(shí)遠(yuǎn)遠(yuǎn)小于GPOPS，這表明凸優(yōu)化方法比非線性規(guī)劃方法的計(jì)算效率更高，但在尋優(yōu)能力上略微下降。

表2 仿真結(jié)果統(tǒng)計(jì)

4 結(jié) 論

本文提出Legendre偽譜凸優(yōu)化方法求解火星進(jìn)入末端高度最大化問(wèn)題，可以快速精確地求出滿足各項(xiàng)約束條件的最優(yōu)軌跡。研究表明：

1)以縱向航程角為自變量建立動(dòng)力學(xué)模型，不必優(yōu)化末端時(shí)間或已知末端高度，可以直接求解末端高度最大化問(wèn)題；

2)提出的Legendre偽譜凸優(yōu)化方法相比一般序列凸優(yōu)化方法，收斂速度更快，且不損失最優(yōu)性；

3)凸優(yōu)化方法(包括Legendre偽譜凸優(yōu)化方法和一般序列凸優(yōu)化方法)相比自適應(yīng)偽譜法，最優(yōu)解接近，但求解速度更快。