火星大氣進入軌跡偽譜凸優化設計方法

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

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

0 引 言

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

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

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

以上軌跡優化方法均為確定性優化方法，而面對真實環境中的大量不確定性，確定性優化方法的局限性日益顯現。為解決不確定條件下的軌跡優化問題，NASA已經在航空航天領域開展了不確定性量化研究，該方法主要環節包括不確定性建模、不確定性分析和不確定性優化。不確定性建模時，根據概率信息的完備情況，可以考慮隨機不確定性和認知不確定性兩類，其中隨機不確定性的建模可以采用概率論方法，而認知不確定性由于無法準確給出概率密度模型或高階統計矩的概率信息，因而只能用上下界表示。不確定性分析的主要目的是估算系統輸出不確定性分布的期望或方差，常見的方法主要有蒙特卡洛法、敏感度矩陣法、協方差分析法和廣義多項式混沌法、隨機配點法、響應面法、拉丁超立方采樣法等。而不確定性優化的研究主要集中于魯棒優化和可靠性優化方法，其中魯棒優化要求在最差工況下仍能滿足任務約束并最小化目標函數，而可靠性優化則給出最小化目標函數時違反任務約束的概率。在大氣進入軌跡不確定性量化方面，Halder等，Prabhakar等較早開展了相關研究，國內崔平遠、李海陽、李爽等人的課題組也開展了相關工作。

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

1 問題描述

1.1 火星大氣進入模型

假設火星大氣相對火星表面靜止，且火星為均勻球體，則彈道升力式火星探測器在火星大氣中的無動力飛行過程可以用一組微分方程表示：

(1)

(2)

(3)

sincoscos)

(4)

(5)

(6)

(7)

式中：、、和分別表示火星探測器的升力系數、阻力系數、參考面積和質量，表示火星大氣密度，表達式為：

=e-

(8)

式中：為火星表面大氣密度;為飛行高度;為標準大氣模型的參考高度;e為自然對數的底數。

火星大氣進入問題的典型約束條件包括過程約束和邊界約束，過程約束為熱流密度、動壓和過載約束：

(9)

()=,()=,≤≤

(10)

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

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

(11)

=-

(12)

1.2 新的自變量

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

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

縱向航程角的初始值為0，末端值則可以根據初始和末端時刻的經緯度確定：

cos=sinsin+coscoscos(-)

(13)

同時，縱向航程角始終單調遞增，相對時間的微分表達式也很簡潔，因此本文選擇作為新的自變量建立火星進入動力學模型。相比基于時間的模型，本文建模方法可以將末端時間自由問題轉化為末端縱向航程角固定問題，避免對末端時間進行優化。同時，相比基于能量的模型，不需要事先已知末端高度和速度，也不需要忽略火星自轉項，模型精度得以保證。以縱向航程角為自變量的狀態微分方程為：

(14)

1.3 狀態量擴充

注意到式(14)中不包含飛行時間，因此本文將無量綱時間擴充為狀態變量，所以本文方法也可以通過定義時間最優目標函數=來求解時間最優問題。引入飛行時間后，狀態量擴展為：

=[,,,,,,]

(15)

此外，選擇傾側角的導數作為控制變量可以抑制傾側角的高頻振蕩，則狀態微分方程可改寫為控制仿射形式：

′=()+()+()

(16)

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

(17)

(18)

(19)

因此，火星進入末端高度最大化問題P可以表述為：

(20)

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

(21)

2 偽譜凸優化方法

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

2.1 一階線性化

凸優化問題是指目標函數和約束條件都為凸函數的最優化問題，而火星進入末端高度最大化問題的非凸性來源于狀態微分方程和過程約束。因此，本文采用一階泰勒展開方法將這兩項約束線性化。

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

(22)

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

(23)

(24)

(25)

(26)

(27)

其余元素均為0，此外：

=-,=

(28)

=-,=

(29)

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

(30)

(31)

(32)

(33)

至此，問題P可以近似為連續凸優化問題P：

(34)

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

(35)

2.2 偽譜離散化

一般的序列凸優化(Sequential convex programming, SCP)方法采用均勻節點離散問題，并通過梯形積分處理狀態微分方程和目標函數中的積分項。為保證離散精度，通常需要布置較多離散點，因此問題規模較大。而偽譜法將問題在一系列全局正交配點上離散，并通過Gauss積分處理狀態微分方程和目標函數的積分項，可以用較少的離散節點保證離散精度。本文采用Legendre偽譜離散格式將問題在LGL節點處離散。

(36)

式中：()為Lagrange插值基函數,且有：

(37)

(38)

(39)

由此可將狀態微分方程轉化為+1組在LGL配點處的等式約束，但需要注意先將問題P的定義域變換到區間[-1,+1]內：

(40)

(41)

為了保證算法收斂時虛擬控制盡可能小，需要在目標函數中施加懲罰項：

(42)

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

注意到問題P的目標函數中并不含有積分項，因此并不需要進行處理。若含有積分項，可使用Gauss-Lobatto求積公式進行近似。

1.3.5 記錄手術時長、術中液體出入量和術中發生不良反應如低血壓、低血氧、牽拉反射引起的牽拉痛、惡心、嘔吐的例數等。

(43)

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

最后，算法收斂時，應滿足前后兩次迭代的解差別不大，同時虛擬控制足夠小，則可以定義收斂標準如下：

(44)

式中：和為給定的收斂標準。

注意到過程約束、以及狀態量或控制量相關的約束不涉及微分或積分項，因此離散化后約束表達式不變，只需將連續變量改為相應的離散變量即可。

至此，離散Legendre偽譜凸優化問題P可定義為：

(45)

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

(46)

2.3 序列凸優化

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

3 仿真校驗

表1 仿真參數設置

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

圖2 經度-緯度剖面和速度-高度剖面

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

圖4 時間-傾側角、傾側角速度剖面

圖5 過程約束剖面

同時，GPOPS產生的最優軌跡雖然和LPCP、SCP的最優軌跡一致，但傾側角、航跡傾角和航跡偏角剖面和其余兩個方法有一定區別，同時傾側角速度剖面的振蕩也較大。從算法上來說，GPOPS將最優控制問題轉化為非線性規劃問題，而LPCP和SCP則是通過序列凸優化逼近原始最優控制問題，且離散格式也不相同，因此兩類方法近似解存在一定差異。

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

表2 仿真結果統計

4 結 論

本文提出Legendre偽譜凸優化方法求解火星進入末端高度最大化問題，可以快速精確地求出滿足各項約束條件的最優軌跡。研究表明：

1)以縱向航程角為自變量建立動力學模型，不必優化末端時間或已知末端高度，可以直接求解末端高度最大化問題；

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

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