行星著陸動(dòng)力下降軌跡優(yōu)化的中心差分凸化方法

李明翔，泮斌峰

（1.西北工業(yè)大學(xué) 航天學(xué)院，西安 710072；2.航天飛行動(dòng)力學(xué)技術(shù)國家級(jí)重點(diǎn)實(shí)驗(yàn)室，西安 710072）

引 言

行星探測(cè)是人類探索太空最重要的途徑之一，月球和火星探測(cè)是當(dāng)前行星探測(cè)的熱點(diǎn)。當(dāng)目標(biāo)行星大氣比較稀薄或者沒有大氣時(shí)，探測(cè)器將無法使用降落傘降至安全速度，基于反推火箭的動(dòng)力下降成為實(shí)現(xiàn)行星表面軟著陸的一種重要手段。由于行星動(dòng)力下降過程存在諸多不確定性因素，其在線軌跡優(yōu)化是實(shí)現(xiàn)探測(cè)器安全、精確著陸所需要解決的重要問題。

近年來，序列凸優(yōu)化方法成為求解上述在線軌跡優(yōu)化的有效手段[1-2]，其本質(zhì)是將無限維空間中的最優(yōu)控制問題通過離散轉(zhuǎn)化為有限維空間中的非線性規(guī)劃問題，并利用序列凸優(yōu)化理論通過凸化處理后迭代求解。Szmuk等[3]通過引入虛擬控制量并且在動(dòng)力學(xué)方程中定義了新的時(shí)間狀態(tài)量，提出了一種終端時(shí)間自由的序列凸優(yōu)化方法計(jì)算框架，能夠較好地得到最優(yōu)燃料軌跡；其不足之處是對(duì)初始猜想敏感且迭代次數(shù)多。Blackmore等[4]通過將燃料最優(yōu)問題凸化為二階錐規(guī)劃（Second-Order-Cone-Programming，SOCP）問題，并考慮了坡度約束，能夠較快求解出燃料最優(yōu)的下降軌跡。上述方法對(duì)于終端時(shí)間固定的情況能夠得到較好的解，但由于其終端時(shí)間采用線性搜索方法計(jì)算，效率較低。此外，Liu等[5]考慮了氣動(dòng)力作用下的一級(jí)火箭回收問題，通過無損凸化技術(shù)和動(dòng)力學(xué)方程歸一化提高了算法收斂性和快速性。該方法在凸化過程中重定義了控制量，公式推導(dǎo)復(fù)雜且引入了額外的約束和計(jì)算量，只能處理固定終端時(shí)間問題。

有別于上述研究方法，本文探索使用預(yù)標(biāo)記中心差分替代解析表達(dá)式計(jì)算動(dòng)力學(xué)矩陣，并對(duì)動(dòng)力學(xué)方程采用直接線性化替代變量代換后再線性化的凸化方法，最后研究行星動(dòng)力下降燃料最優(yōu)控制問題中近似最優(yōu)終端時(shí)間的確定方法，以提高序列凸優(yōu)化在線軌跡優(yōu)化方法的工程實(shí)用性。

1 行星動(dòng)力下降軌跡優(yōu)化問題建模

以著陸點(diǎn)為原點(diǎn)建立空間直角坐標(biāo)系，著陸點(diǎn)一般不直接選在地表，而是距離地面稍有一定距離，方便提供安全緩沖余量。x軸方向預(yù)先約定，z軸垂直于地表向上，y軸滿足右手螺旋關(guān)系。行星動(dòng)力下降可以表述為經(jīng)典的Bolza型最優(yōu)控制問題

其中終端性能指標(biāo)為

式（1）中約束分別為初值約束、動(dòng)力學(xué)約束、控制量幅值約束和防撞約束，其中反推發(fā)動(dòng)機(jī)為了避免多次開關(guān)機(jī)帶來的可靠性問題，在實(shí)際工作中一般不將其完全關(guān)閉，故需限定最小推力。行星動(dòng)力下降動(dòng)力學(xué)方程為

其中：g=[00–g（rz）]T是行星的重力加速度矢量。

考慮到動(dòng)力下降段高度在2～5 km之間，高度相比行星半徑尺度較小，重力加速度變化也較小，可以直接通過海拔高度計(jì)算。本文以火星引力場(chǎng)為例，其表面重力加速度為gM=0.377、gE=3.7 m/s2、RM=3 397 km，距地表任意高度的重力加速度為

此外，d=[dxdydz]T是3個(gè)方向的氣動(dòng)加速度，分為可建模的氣動(dòng)（如飛行器的氣動(dòng)阻力）和無法建模的隨機(jī)風(fēng)場(chǎng)擾動(dòng)。由于火星大氣較為復(fù)雜，不適合在在線軌跡優(yōu)化中應(yīng)用[1]，因此本文采用文獻(xiàn)[3～4]的處理方式，在動(dòng)力下降段中以忽略氣動(dòng)力作用下簡化的動(dòng)力學(xué)進(jìn)行分析，氣動(dòng)擾動(dòng)可以通過后期設(shè)計(jì)控制器加以克服[6]。

2 序列凸優(yōu)化問題建立

2.1 控制量凸化

采用直接法求解優(yōu)化問題式（1）時(shí)，需要先將其離散到有限維空間中。當(dāng)終端時(shí)間tf確定時(shí)，將飛行時(shí)間[0，tf]等距離散為N個(gè)時(shí)刻，即N–1個(gè)子區(qū)間，并在這N個(gè)時(shí)刻對(duì)狀態(tài)量和控制量進(jìn)行離散，方便施加過程約束。

對(duì)于控制力幅值約束

因最小推力幅值約束是二階錐補(bǔ)集，為將非凸約束式（5）凸化，在各個(gè)離散點(diǎn)處引入推力松弛因子Fη

此時(shí)，控制幅值約束可以表示為

至此非凸約束式（5）被等價(jià)轉(zhuǎn)換為凸約束式（6）和式（7），積分型性能指標(biāo)也可以表達(dá)為關(guān)于Fη的線性組合

其中：Δt=（tf–t0）/（N–1），即使用矩形公式將積分化為求和。

2.2 中心差分法改進(jìn)動(dòng)力學(xué)凸化

為將非線性動(dòng)力學(xué)微分方程約束處理為凸約束，需進(jìn)行凸化處理。文獻(xiàn)[5]和文獻(xiàn)[7]中處理上述動(dòng)力學(xué)方程通常使用2種凸化技巧：將推力矢量分解為推力大小和方向余弦、對(duì)質(zhì)量取對(duì)數(shù)重定義新的狀態(tài)量。前者存在的問題是，方向余弦約束為非凸等式約束，雖然文獻(xiàn)[3]中通過極小值原理證明其可無損轉(zhuǎn)化為凸錐約束，但實(shí)際數(shù)值求解中經(jīng)常出現(xiàn)不穩(wěn)定因素，造成收斂速度較慢，且控制變量從原來的3維上升到4維，問題復(fù)雜度上升。后者對(duì)質(zhì)量取對(duì)數(shù)后重定義新的狀態(tài)量，再對(duì)推力約束進(jìn)行一階泰勒展開線性化和直接對(duì)動(dòng)力學(xué)泰勒展開線性化，本質(zhì)上都是線性化，并未降低原問題的非線性，僅將非線性從動(dòng)力學(xué)方程轉(zhuǎn)移到過程約束中。取對(duì)數(shù)后推力不等式約束中存在指數(shù)函數(shù)，相比原來僅含加、乘法的公式增加了計(jì)算量，且人為增加了推導(dǎo)公式和編程的工作量。

對(duì)于動(dòng)力學(xué)矩陣A、B、C的分量使用差分直接計(jì)算，第i個(gè)微分方程，對(duì)第n個(gè)離散點(diǎn)處的狀態(tài)量（控制量）的第j個(gè)分量的偏導(dǎo)的差分計(jì)算公式為

為了消除不必要的計(jì)算量，在設(shè)計(jì)程序時(shí)，提前對(duì)每個(gè)動(dòng)力學(xué)微分方程顯含的狀態(tài)量和控制量索引進(jìn)行標(biāo)注，僅計(jì)算這些位置的偏導(dǎo)數(shù)，其余元素直接賦值為0。將式（1）中的微分方程約束轉(zhuǎn)換為線性代數(shù)方程約束，為表示方便，簡記為

為兼顧計(jì)算精度和計(jì)算效率，使用梯形法施加動(dòng)力學(xué)約束

啟動(dòng)第一次迭代時(shí)，動(dòng)力學(xué)矩陣A，B，C的計(jì)算需要提供k=0時(shí)每個(gè)離散點(diǎn)n處的狀態(tài)量x（0）和控制量u（0），即初始猜想。對(duì)于狀態(tài)量，如果同時(shí)存在初值約束和終端約束，則狀態(tài)量的初始猜想直接設(shè)計(jì)為兩個(gè)狀態(tài)的線性插值

如果僅存在初值約束（如質(zhì)量m），則狀態(tài)量初始猜想設(shè)計(jì)為常值向量

對(duì)于控制量，初始猜測(cè)全部初始化為0，即

式（12）～（14）簡潔的設(shè)計(jì)過程直接消除了利用數(shù)值積分設(shè)計(jì)的軌跡引入的人為因素和計(jì)算量，雖然初始猜想并不符合動(dòng)力學(xué)方程約束，但在每次迭代求解后得到的軌跡將完全滿足動(dòng)力學(xué)約束，并在不斷迭代中使軌跡收斂到最優(yōu)解。

3 數(shù)值仿真

3.1 固定終端時(shí)間算例仿真

為驗(yàn)證前文所設(shè)計(jì)算法的有效性，給定飛行狀態(tài)和航天器參數(shù)如表1所示，求解行星動(dòng)力下降段燃料最優(yōu)的軌跡優(yōu)化問題。

表1 仿真設(shè)定參數(shù)Table 1 Simulation parameters

圖1 性能指標(biāo)及其絕對(duì)偏差和相對(duì)偏差隨迭代次數(shù)變化圖Fig.1 Performance index and its absolute deviation，relative deviation versus iterations

本文采用C++調(diào)用SOCP求解器ECOS[9]進(jìn)行求解。計(jì)算環(huán)境建立在CPU為i5-8300H（基頻2.2 GHz，睿頻3.8 GHz）Windows系統(tǒng)的PC上，耗時(shí)統(tǒng)計(jì)從進(jìn)入main函數(shù)開始到滿足序列凸優(yōu)化迭代條件er＜ 1%退出為止（不含優(yōu)化結(jié)果輸出時(shí)間），求解500次，均為迭代3步收斂，耗時(shí)統(tǒng)計(jì)如圖2所示，平均耗時(shí)0.130 1 s，優(yōu)化所得控制量與狀態(tài)量如圖3～6所示。

圖2 求解表1算例500次耗時(shí)統(tǒng)計(jì)Fig.2 Computing time statistic of solving the case in Table 1 500 times

優(yōu)化得到的推力曲線如圖3所示，可以看到推力幅值被精確地約束在[100，1 000]N的區(qū)間內(nèi)，表現(xiàn)為bang-bang控制。速度曲線如圖4所示，可看到終端速度被精確約束在0 m/s，圖5中質(zhì)量曲線m（t）的導(dǎo)數(shù)是控制力幅值函數(shù)||F（t）||2，可以看到m（t）隨||F（t）||2也相應(yīng)分為3段。

圖3 發(fā)動(dòng)機(jī)3個(gè)方向推力曲線Fig.3 Thrust component in three directions

圖4 飛行器速度曲線Fig.4 Velocity profiles of optimal trajectory

圖5 飛行器質(zhì)量曲線Fig.5 Mass profile of optimal trajectory.

圖1證明了用中心差分法建立序列凸優(yōu)化模型的有效性，良好收斂性和對(duì)初始猜測(cè)弱敏感性，事實(shí)上由圖1還可以看出如果以相對(duì)偏差作為性能指標(biāo)，可認(rèn)為第3步已經(jīng)收斂，軌跡隨迭代次數(shù)變化如圖7所示，3～5次迭代產(chǎn)生的軌跡基本重合，算例對(duì)應(yīng)的三維軌跡和各離散點(diǎn)處推力矢量如圖8所示。

圖6 飛行器位置相對(duì)于著陸點(diǎn)變化曲線Fig.6 Position profiles of optimal trajectory

圖7 tf=60 s軌跡隨迭代次數(shù)變化Fig.7 Trajectory versus iterations for case of tf=60 s

圖8 tf=60 s優(yōu)化軌跡和推力矢量Fig.8 Optimal trajectory and thrust for case oftf=60 s

3.2 最優(yōu)終端時(shí)間計(jì)算

圖9 tf∈[50，190]s最優(yōu)軌跡簇Fig.9 Optimal trajectories attf∈[50，190]s

圖10 不同tf下最優(yōu)軌跡對(duì)應(yīng)的剩余燃料質(zhì)量mfFig.10 mfof optimal trajectories at differenttf

由圖11可見，降低離散點(diǎn)個(gè)數(shù)對(duì)最優(yōu)終端時(shí)間的確定幾乎不產(chǎn)生影響，但可以顯著降低計(jì)算量，如圖12所示。在實(shí)際工程應(yīng)用中，可以通過降低離散點(diǎn)個(gè)數(shù)的方式或者使用低精度的動(dòng)力學(xué)約束進(jìn)行最優(yōu)控制模型的求解，以快速獲?。╰f，mf）樣本點(diǎn)。樣本點(diǎn)個(gè)數(shù)不能少于3個(gè)，使用最小二乘逼近算法得到a、b、c 3個(gè)參數(shù)，利用式（15）計(jì)算出最優(yōu)終端時(shí)間后，再增大離散點(diǎn)個(gè)數(shù)進(jìn)行精規(guī)劃。

圖11 不同離散點(diǎn)個(gè)數(shù)對(duì)應(yīng)的mf-tf曲線Fig.11 mf-tfcurves under different number of discrete points

圖12 不同離散點(diǎn)個(gè)數(shù)與動(dòng)力學(xué)約束下的計(jì)算耗時(shí)Fig.12 Computing time under different dynamics constraint versus discrete points

為驗(yàn)證上述方法的有效性，用計(jì)算機(jī)隨機(jī)產(chǎn)生100組初始狀態(tài)和終端狀態(tài)，計(jì)算出每組條件下的mf-tf曲線，仿真參數(shù)如表2所示。仿真結(jié)果如圖13所示，在100個(gè)飛行條件下，除去2個(gè)苛刻的飛行情況無解外，其余98個(gè)條件下都找到了最優(yōu)燃料控制的規(guī)劃軌跡，且有較高的規(guī)劃成功率。

表2 隨機(jī)仿真參數(shù)設(shè)置Table 2 Simulation parameters

圖13 隨機(jī)初始狀態(tài)下最優(yōu)軌跡簇的mf-tf曲線Fig.13 mf-tfcurves of optimal trajectories under random initial state

3.3 算法對(duì)比與誤差分析

記本文中不對(duì)動(dòng)力學(xué)改造直接采用中心差分進(jìn)行線性化建立的序列凸優(yōu)化算法為P0，另外兩種凸化方法分別為P1[10]和P2[11]

圖14 3種凸化方法最優(yōu)軌跡對(duì)比Fig.14 Comparison of optimal trajectories obtained by three convexification methods

圖15 3種凸化方法求解得到的推力幅值曲線對(duì)比Fig.15 Comparison of thrust magnitude obtained by three convexification methods

表3 3種凸化方法的計(jì)算結(jié)果Table 3 Results of three convexification methods

圖16 P0、P1、P2求得規(guī)劃軌跡與積分軌跡的誤差隨時(shí)間變化曲線Fig.16 Error versus time between the integral trajectories and the optimal trajectories obtained by methods P0，P1，P2.

為將非線性動(dòng)力學(xué)模型轉(zhuǎn)化為凸問題進(jìn)行求解，三者都使用了等距離散的分段線性動(dòng)力學(xué)進(jìn)行近似逼近，線性化誤差和離散誤差不可避免。使用控制量積分得到的軌跡和規(guī)劃所得軌跡進(jìn)行對(duì)比，積分步長為0.01 s。發(fā)現(xiàn)3種凸化方式規(guī)劃解和積分解誤差隨時(shí)間變化都在同一數(shù)量級(jí)內(nèi)。

在表4設(shè)置的參數(shù)下，產(chǎn)生500組隨機(jī)動(dòng)力下降任務(wù)，分別使用3種方法進(jìn)行求解，觀察終端誤差散布。仿真結(jié)果如圖17～20所示。

表4 終端誤差散布仿真參數(shù)Table 4 Simulation parameters

圖17 3種方法積分所得終端位置誤差散布圖Fig.17 Terminal position error of integral trajectories

圖18 P0和P2積分所得終端位置誤差散布圖Fig.18 Terminal position error of integral trajectories obtained by convexification methods P0 and P2

圖19 3種方法積分所得終端速度誤差散布圖Fig.19 Terminal velocity error of integral trajectories obtained by three methods

圖20 P0和P2積分所得終端位置誤差散布圖Fig.20 Terminal velocity error of integral trajectories obtained by convexification methods P0 and P2

根據(jù)對(duì)比，可以發(fā)現(xiàn)P1方法的數(shù)值穩(wěn)定性較差，部分算例達(dá)到最大迭代次數(shù)20次后仍然沒有收斂，導(dǎo)致終端誤差較大，P0和P2方法終端誤差數(shù)量級(jí)接近。P0相比于P2的位置終端誤差散布更加集中，且速度終端誤差在z方向上更小。

另外，離散步長對(duì)終端誤差也都有一定影響。以P0為例，計(jì)算表1任務(wù)在不同離散點(diǎn)個(gè)數(shù)下得到的規(guī)劃軌跡，并通過積分得到終端誤差如圖21和圖22所示?？梢钥吹?，終端誤差隨離散點(diǎn)個(gè)數(shù)的增加而不斷減少。

圖21 積分得到終端位置誤差隨離散點(diǎn)個(gè)數(shù)變化圖Fig.21 Terminal position error of integral trajectory versus the number of discrete points

圖22 積分得到終端速度誤差隨離散點(diǎn)個(gè)數(shù)變化圖Fig.22 Terminal velocity error of integral trajectory versus the number of discrete points

固定離散點(diǎn)個(gè)數(shù)，針對(duì)表5設(shè)置的隨機(jī)參數(shù)，采用P0在不同終端時(shí)間下各打靶200次，如圖23和圖24所示，積分所得終端誤差散布也處于可接受范圍內(nèi)。

圖23 不同固定終端時(shí)間下積分所得終端位移誤差三維散布圖及xy平面圖Fig.23 3D and projection on xy plane of terminal position error dispersion of integral trajectories at different fixed terminal times

圖24 不同固定終端時(shí)間下積分所得終端速度誤差三維散布圖及xy平面圖Fig.24 3D and projection on XY plane of terminal velocity error dispersion of integral trajectories at different fixed terminal times

表5 終端誤差散布仿真參數(shù)Table 5 Simulation parameters

對(duì)于收斂性，200個(gè)隨機(jī)初始狀態(tài)中，P0、P1都求出198個(gè)解（第43和181個(gè)失?。?，P2求出199個(gè)解（第43個(gè)失?。呤諗啃詿o顯著差別。兩任務(wù)的初始速度大，高度低，規(guī)劃也失敗屬于合理現(xiàn)象。雖然第181個(gè)任務(wù)的解被P2找到，但這條軌跡十分危險(xiǎn)，在18.18 s處離地面較近且仍然有較大的水平速度，考慮存在離散誤差和環(huán)境擾動(dòng)的情況下可能會(huì)導(dǎo)致任務(wù)失敗，兩任務(wù)參數(shù)如下。

圖25 P2解得任務(wù)181的規(guī)劃軌跡Fig.25 The optimal trajectory of case 181 by using P2 convexification method

圖26 P2解得任務(wù)181的速度軌線Fig.26 Velocity profiles of case 181 by using convexification method P2

圖27 200組隨機(jī)任務(wù)打靶耗時(shí)統(tǒng)計(jì)圖Fig.27 Computing time statistics of 200 random cases

圖28 200組隨機(jī)任務(wù)中指標(biāo)J的統(tǒng)計(jì)圖Fig.28 The index J statistics of 200 random cases

以上求解過程都未引入信賴域約束，若增加狀態(tài)量信賴域約束，在表1仿真參數(shù)下采用P0進(jìn)行求解，每種信賴域各計(jì)算100次，圖29顯示對(duì)解最優(yōu)性并未有明顯改善。圖30為不同信賴域下求解耗時(shí)統(tǒng)計(jì)。前3條軌線和無信賴域約束求解結(jié)果基本重合，第4條軌跡消耗燃料相比前3條節(jié)省了0.032 7 kg，僅占總質(zhì)量的0.065 4%，但耗時(shí)相對(duì)較長。在本問題中基本可視信賴域?yàn)槿哂嗉s束。

圖29 不同信賴域約束下優(yōu)化得到的軌線Fig.29 Optimal trajectories under different trust regions

圖30 不同信賴域下求解耗時(shí)統(tǒng)計(jì)Fig.30 Computing time statistics under different trust regions

表6 不同信賴域下的求解結(jié)果Table 6 Results under different trust regions.

4 結(jié) 論

本文使用中心差分逐次線性化動(dòng)力學(xué)，提高了序列凸優(yōu)化算法的規(guī)范性和工程性，避免因使用推力方向余弦無損凸化和質(zhì)量方程取對(duì)數(shù)凸化而引入額外的優(yōu)化變量，降低了公式推導(dǎo)工作量和編程復(fù)雜性。通過對(duì)比證明，相對(duì)于傳統(tǒng)動(dòng)力學(xué)凸化技巧，直接線性化的求解結(jié)果并未對(duì)算法收斂性、解的最優(yōu)性和終端誤差造成顯著影響。同時(shí)，采用性能指標(biāo)相對(duì)偏差作為迭代收斂條件，使收斂條件和具體模型解耦，具有更好的普適性。

針對(duì)動(dòng)力下降段最優(yōu)控制問題，給出最優(yōu)軌跡簇剩余燃料和終端時(shí)間的擬合函數(shù)，可解析求出近似最優(yōu)終端時(shí)間，能夠有效提高在線軌跡方法的效率。

本文使用的動(dòng)力學(xué)模型較為簡單，未考慮更多的過程約束和環(huán)境條件，如發(fā)動(dòng)機(jī)推力方向約束、落地傾角約束和氣動(dòng)力作用，但可為復(fù)雜約束下的軌跡規(guī)劃算法提供參考初始軌跡，也可為軌跡跟蹤控制算法提供標(biāo)稱軌跡。

由于存在線性化誤差和離散誤差二者的復(fù)合作用，主要采用大批次數(shù)值仿真計(jì)算來評(píng)估算法收斂能力，對(duì)序列凸優(yōu)化算法收斂性的嚴(yán)格數(shù)學(xué)證明尚比較困難，需要進(jìn)一步進(jìn)行深入研究。