擾動(dòng)環(huán)境下火星精確著陸自主軌跡規(guī)劃方法

鄧云山,夏元清,*,孫中奇,沈剛輝

1. 北京理工大學(xué) 自動(dòng)化學(xué)院,北京 100081

2. 復(fù)雜系統(tǒng)智能控制與決策國家重點(diǎn)實(shí)驗(yàn)室,北京 100081

未來的火星探測(cè)任務(wù)要求探測(cè)器具有精確著陸的能力[1-2]。探測(cè)器需要經(jīng)歷大氣進(jìn)入(Entry)、下降(Descent)以及最終著陸(Landing)過程(也稱為EDL過程)[3],其中下降段包含傘降段與動(dòng)力下降段(Powered-Descent Guidance, PDG),而動(dòng)力下降段的軌跡規(guī)劃是動(dòng)力下降段制導(dǎo)控制問題中的重要研究內(nèi)容。因此提高軌跡規(guī)劃速度,使著陸器在線實(shí)時(shí)生成著陸軌跡對(duì)提高著陸器的生存能力和著陸精度具有重要意義[4]。

著陸器拋棄降落傘后,探測(cè)器進(jìn)入動(dòng)力下降段,需要利用火箭反推技術(shù),進(jìn)行姿態(tài)調(diào)整及障礙規(guī)避,最終懸停于預(yù)期著陸點(diǎn)正上方,釋放火星車[3]。整個(gè)過程需要滿足最大、最小推力約束,最小下降角等多個(gè)約束條件。文獻(xiàn)[5]將能量最優(yōu)動(dòng)力下降段制導(dǎo)問題轉(zhuǎn)化為二階錐優(yōu)化問題。在給定精度下,此類型優(yōu)化問題可使用內(nèi)點(diǎn)法在多項(xiàng)式時(shí)間內(nèi)求解[6-8]。文獻(xiàn)[9]在文獻(xiàn)[5]的基礎(chǔ)上,針對(duì)著陸點(diǎn)不可達(dá)問題,提出了最小著陸誤差的軌跡規(guī)劃算法,首先求解最小著陸誤差問題得到距離預(yù)期著陸點(diǎn)最小距離的可降落區(qū)域,然后在該降落區(qū)域內(nèi)求解能量最優(yōu)規(guī)劃問題。文獻(xiàn)[4] 將能量最優(yōu)下降問題轉(zhuǎn)化為有限個(gè)推力為常值的分段函數(shù),給出了每段內(nèi)微分方程的解析表達(dá)式,并利用序列二次規(guī)劃(Sequential Quadratic Program,SQP)算法對(duì)待求參數(shù)進(jìn)行迭代求解。這種方法在一定程度上提高了計(jì)算效率。文獻(xiàn)[10]基于凸優(yōu)化方法研究了給定著陸精度下的約束可控集與約束可達(dá)集,通過對(duì)數(shù)據(jù)集的量化,評(píng)估給定著陸器著陸的可行性。文獻(xiàn)[11]結(jié)合凸優(yōu)化方法與偽譜法,對(duì)比了Radau偽譜法、Lobatto偽譜法與標(biāo)準(zhǔn)凸優(yōu)化方法的求解結(jié)果,結(jié)果表明,結(jié)合偽譜法的求解精度較高,但求解時(shí)間增加,在少量譜點(diǎn)情況下具有工程實(shí)用性。文獻(xiàn)[12] 研究了不規(guī)則小行星動(dòng)力著陸軌跡規(guī)劃問題,考慮引力模型,對(duì)固定時(shí)間著陸問題進(jìn)行了無損凸化。

上述方法只考慮了固定終端問題,實(shí)際問題中,終端時(shí)間也為未知變量。文獻(xiàn)[5,9]采用外嵌黃金分割法(Golden Search Procedure, GSP)進(jìn)行終端時(shí)間的線性搜索,這種方法可以解決終端時(shí)間自由問題,但是以犧牲計(jì)算時(shí)間為代價(jià)。文獻(xiàn)[13]考慮阻力模型,針對(duì)終端時(shí)間自由問題,對(duì)非線性動(dòng)力學(xué)直接進(jìn)行線性化,初步解決了終端時(shí)間自由軌跡規(guī)劃問題。文獻(xiàn)[14]針對(duì)直接線性化處理中,終端時(shí)間初值的估計(jì)不準(zhǔn)確時(shí),問題可能在迭代前幾次不可行進(jìn)而造成求解失敗的情況,選擇對(duì)高度變量進(jìn)行求導(dǎo)并設(shè)計(jì)高度離散方法。所設(shè)計(jì)的高度離散方法可盡可能使對(duì)應(yīng)離散點(diǎn)時(shí)間間隔相同,但無法保證統(tǒng)一的時(shí)間步長,同時(shí)無法對(duì)擾動(dòng)進(jìn)行建模。

針對(duì)擾動(dòng)環(huán)境,常采用2種方式進(jìn)行克服：① 采 用預(yù)先規(guī)劃加反饋跟蹤控制的思想進(jìn)行;② 在 下降過程中不斷進(jìn)行重規(guī)劃。其中,前者在擾動(dòng)作用下,規(guī)劃軌跡在實(shí)際中可能不可行。后者雖然從規(guī)劃層次上考慮了擾動(dòng),相比于前者具由更強(qiáng)的魯棒性,但重規(guī)劃問題的可行性無法從理論上進(jìn)行論證,可能出現(xiàn)重規(guī)劃不可行的情況。值得說明的是,文獻(xiàn)[15]結(jié)合本質(zhì)魯棒MPC(Model Predict Control)思想,提出了保證可行性的重規(guī)劃問題,巧妙解決了重規(guī)劃不可行問題。但其未考慮2次重規(guī)劃間隔間的抗擾,抗擾能力與2次重規(guī)劃間隔相關(guān),重規(guī)劃頻率越高,抗擾能力越強(qiáng),但依靠高頻率優(yōu)化進(jìn)行制導(dǎo),會(huì)出現(xiàn)重規(guī)劃頻率與計(jì)算能力的矛盾。增加內(nèi)環(huán)反饋控制后可降低重規(guī)劃頻率,但同時(shí)需要考慮規(guī)劃問題與內(nèi)環(huán)控制的耦合。本文針對(duì)終端自由擾動(dòng)環(huán)境下火星著陸動(dòng)力下降段軌跡規(guī)劃問題,結(jié)合Tube-MPC思想[16]設(shè)計(jì)反饋控制律,考慮規(guī)劃問題與內(nèi)環(huán)控制的耦合,提出了魯棒意義下的自主軌跡規(guī)劃框架并分析了重規(guī)劃問題的可行性。外環(huán)采用序列凸優(yōu)化(Sequential Convex Programming, SCP)方法進(jìn)行求解,并分析了求解結(jié)果的最優(yōu)性；內(nèi)環(huán)采用反饋控制律直接控制,進(jìn)一步提高抗擾性能。同時(shí)給出了重規(guī)劃問題可行的必要條件,為實(shí)際工程應(yīng)用提供參考。

1 問題描述

火星著陸動(dòng)力下降段任務(wù)高度在距離火星5 km 高度內(nèi)[5],重力加速度可視為常數(shù)。由于火星大氣稀薄,在動(dòng)力下降段最優(yōu)軌跡規(guī)劃中通常將氣動(dòng)力等其他力的作用視作擾動(dòng)。

動(dòng)力下降段開始前,通過傘降段充分減速,此時(shí)著陸器速度量級(jí)為100 m/s。動(dòng)力下降段開始時(shí),著陸器與降落傘分離,點(diǎn)燃火箭反推裝置,從距離火星表面約1.5 km開始下降到期望著陸點(diǎn)正上方20 m處進(jìn)行懸停。

動(dòng)力下降段自主軌跡規(guī)劃在初始時(shí)刻進(jìn)行一次規(guī)劃,隨著著陸器的下降進(jìn)行實(shí)時(shí)重規(guī)劃。整個(gè)過程火箭反推推力需要位于允許最小推力與允許最大推力之間,著陸軌跡需要避免下降過程中與火星表面發(fā)生碰撞,并盡可能使整個(gè)過程能耗最小。

1.1 著陸器動(dòng)力學(xué)

以期望著陸點(diǎn)為原心,X軸沿火星經(jīng)線指向東方,Z軸豎直向上,Y軸滿足右手定則,建立火星表面固連坐標(biāo)系。著陸器動(dòng)力學(xué)為

(1)

(2)

1.2 耗能最優(yōu)

能耗最優(yōu)的需求可描述為動(dòng)力下降任務(wù)開始(0時(shí)刻)到結(jié)束(tf時(shí)刻),所消耗燃料的質(zhì)量最小,也即最大化任務(wù)結(jié)束時(shí)刻著陸器質(zhì)量與任務(wù)初始時(shí)刻著陸器質(zhì)量的差：

maxm(tf)-m(0)

(3)

同時(shí),由于著陸器質(zhì)量滿足式(2),因此上述指標(biāo)等價(jià)于最小化任務(wù)過程中推力大小的積分：

(4)

1.3 過程約束

下降過程中,火箭發(fā)動(dòng)機(jī)推力大小處于最小允許推力ρmin與最大允許推力ρmax之間,即

(5)

同時(shí)為避免與火星表面的山峰狀障礙物相撞,下降軌跡需要滿足滑降約束：

r(t)∈X

(6)

X={r∈R3:

(7)

式中:rX、rY、rZ分別為著陸器位置r沿坐標(biāo)軸方向的投影分量;rfX、rfY、rfZ分別為預(yù)期懸停位置rf沿坐標(biāo)軸方向的投影分量;γ為一給定角度常數(shù)。

滑降約束幾何含義為著陸器下降軌跡位于以預(yù)期懸停位置為頂點(diǎn),開口方向平行于Z軸向上的圓錐型區(qū)域內(nèi),屬于二階錐約束。

1.4 起始、終端約束

任務(wù)起始于著陸器滿載燃料的質(zhì)量mwet,初始位置r0與速度v0：

(8)

要求終止時(shí)刻著陸器的質(zhì)量大于燃料耗盡時(shí)著陸器的質(zhì)量mdry,且終止位置、速度達(dá)到預(yù)期懸停位置rf與速度vf：

(9)

1.5 無擾燃料最優(yōu)軌跡預(yù)規(guī)劃問題

著陸過程可描述為以能量最優(yōu)為指標(biāo)的最優(yōu)控制問題,需要滿足著陸器動(dòng)力學(xué)約束、狀態(tài)約束、控制約束。

燃料最優(yōu)火星著陸動(dòng)力下降段軌跡規(guī)劃問題以下降過程中消耗燃料最少為性能指標(biāo),表示為[9]

P1(非凸燃料最優(yōu)問題)：

(10)

s.t. 式(1),式(2),式(5),式(6),式(8),式(9)

其中:式(1)～式(2)為非線性等式約束,式(5)為非凸不等式約束。

文獻(xiàn)[5]與文獻(xiàn)[9]對(duì)優(yōu)化問題P1進(jìn)行凸化處理與變量替換,通過引入新控制變量Γ,對(duì)P1進(jìn)行無損凸化,并給出了無損凸化的充分條件。并通過如下變量替換,得到問題P2。

(11)

(12)

z(t)=lnm(t)

(13)

P2(松弛燃料最優(yōu)問題)：

(14)

s.t.

(15)

(16)

(17)

0≤ρmine-z(t)≤σ(t)≤ρmaxe-z(t)

(18)

r(t)∈X

(19)

(20)

引入新變量并進(jìn)行變量替換后,動(dòng)力學(xué)約束式(1)變?yōu)榫€性形式如式(15)所示。

2 擾動(dòng)環(huán)境下軌跡重規(guī)劃

實(shí)際下降過程中,由于受以氣動(dòng)擾動(dòng)為主的擾動(dòng)影響,預(yù)先規(guī)劃的軌跡在實(shí)際下降過程中可能不可行,因此需要實(shí)時(shí)進(jìn)行重規(guī)劃,克服擾動(dòng)的影響。對(duì)于重規(guī)劃問題,本節(jié)將分析其可行性,并構(gòu)造出重規(guī)劃問題的可行解。

首先,給出擾動(dòng)環(huán)境下的動(dòng)力學(xué)模型,針對(duì)線性反饋控制律,參考Tube魯棒思想[16],提出擾動(dòng)環(huán)境下的軌跡重規(guī)劃優(yōu)化問題。其次,證明重規(guī)劃問題的可行性。最后,提出擾動(dòng)環(huán)境下火星著陸軌跡規(guī)劃架構(gòu)。

2.1 擾動(dòng)環(huán)境下的軌跡重規(guī)劃優(yōu)化模型

由于火星大氣啟動(dòng)參數(shù)難以在著陸前進(jìn)行精確測(cè)量,因此在已處理動(dòng)力學(xué)模型式(15)的基礎(chǔ)上建立擾動(dòng)動(dòng)力學(xué)模型：

(21)

式中:dp(t)為火星大氣等客觀因素引起的擾動(dòng),擾動(dòng)具有上界|dp(t)|≤η, |·|表示向量按位取絕對(duì)值。記m⊙n為2個(gè)向量按位相乘。tk時(shí)刻擾動(dòng)環(huán)境下重規(guī)劃問題如下。

P3(k)(擾動(dòng)環(huán)境下重規(guī)劃問題)：

(22)

s.t.

(23)

(24)

(25)

(26)

r(t)∈X

(27)

考慮單一方向上的特征根分布情況,特征方程在復(fù)平面左半平面內(nèi)有2個(gè)共軛復(fù)根λ1=λRe+λImi,λ2=λRe-λImi(λRe<0,λIm≥0)時(shí),rtube=η(2λIm-λRe)/(λImλ1λ2),vtube=η/λIm；當(dāng)特征方程在復(fù)平面左半平面內(nèi)有2個(gè)不同實(shí)根λ1≠λ2時(shí),有rtube=η/(λ1λ2),vtube=η/|λ1-λ2|;當(dāng)特征方程在復(fù)平面左半平面內(nèi)有兩相同實(shí)根λ1=λ2=λ時(shí),rtube=2η/λ2,vtube=-η/eλ,e為自然常數(shù)。

注意,重規(guī)劃問題P3(k)的初始狀態(tài)被作為一個(gè)決策變量,實(shí)際控制量采用線性反饋形式進(jìn)行控制,控制律如下：

(28)

事實(shí)上,預(yù)先規(guī)劃的本質(zhì)是開環(huán)的,無法從規(guī)劃的角度上考慮實(shí)時(shí)擾動(dòng)；重規(guī)劃對(duì)下降過程中的軌跡進(jìn)行修正,可以達(dá)到更優(yōu)的性能。重規(guī)劃中,若不對(duì)優(yōu)化問題進(jìn)行處理,則無法保證問題的可行性。

2.2 重規(guī)劃可行性分析

與直接以P2問題進(jìn)行重規(guī)劃相比,P3(k)可以在理論上保證重規(guī)劃問題的可行性,并可以根據(jù)上次規(guī)劃結(jié)果構(gòu)造出下次重規(guī)劃問題的可行解。本節(jié)將在連續(xù)時(shí)間情況下,分析重規(guī)劃問題的可行性,首先給出引理1。

1) 單一方向上特征根λ1、λ2為不相等共軛復(fù)根,即λ1=λRe+λImi,λ2=λRe-λImi(λRe<0,λIm>0)時(shí):

(29)

2) 單一方向上特征根為不相等實(shí)根,即λ1≠λ2時(shí):

(30)

3) 單一方向上特征根為相等實(shí)根,即λ1=λ2=λ時(shí):

(31)

證明：考慮單個(gè)方向的誤差狀態(tài),由于

(32)

代入線性反饋控制律得：

(33)

即誤差狀態(tài)方程為

(34)

下面就特征根進(jìn)行分類討論：

1) 當(dāng)特征根λ1、λ2為不相等共軛復(fù)根時(shí)：

C1eλ1t+C2eλ2t

(35)

有

(36)

第1項(xiàng)展開為

(37)

(38)

對(duì)(35)式求導(dǎo)有

(39)

放縮得到

(40)

2) 當(dāng)特征根λ1、λ2為兩不等的實(shí)根時(shí),由于：

(41)

因此式(36)可放縮為

(42)

式(40)變?yōu)?/p>

(43)

3) 當(dāng)特征根λ1=λ2=λ<0為相等的實(shí)根時(shí)：

(44)

(45)

求導(dǎo)可得(λt-1)eλt單調(diào)遞增,因此有

-1≤(λt-1)eλt≤0

(46)

所以

(47)

對(duì)式(44)求導(dǎo)得

(48)

有

(49)

求導(dǎo)可得teλt在t=-1/λ時(shí)取最大值,結(jié)合teλt≥0可得

(50)

(51)

(52)

(53)

下面對(duì)重規(guī)劃問題P3(k)進(jìn)行可行性分析。

不難驗(yàn)證部分解滿足P3(k+1)問題的式(6),式(20),式(23)～(26)約束。

考察部分解的初始狀態(tài),由引理1可知：

(54)

由于重歸化問題可行性需求,需要對(duì)預(yù)規(guī)劃問題P2的約束式(18)進(jìn)行修正,用P3(k)中的約束式(26)替代,得到修正后的預(yù)規(guī)劃問題P4。

通過數(shù)學(xué)歸納法可得,對(duì)于連續(xù)時(shí)間系統(tǒng),當(dāng)修正后的預(yù)規(guī)劃問題P4可行時(shí),后續(xù)重規(guī)劃問題可行,且最優(yōu)性更強(qiáng)。下面給出預(yù)規(guī)劃問題P4的必要性條件,為參數(shù)篩選提供參考。

定理2當(dāng)預(yù)規(guī)劃問題P4可行時(shí),控制參數(shù)滿足如下不等式：

(55)

(56)

移項(xiàng)并代入uemax=kp⊙rtube+kd⊙vtube得：

(57)

依據(jù)式(13)進(jìn)行變量反代換得：

(58)

□

其逆否命題可作為控制參數(shù)初步篩選的依據(jù)：滿足式(59)的參數(shù)一定使預(yù)規(guī)劃問題P4不可行。

(59)

注意,滿足不等式(55)的參數(shù)不一定使問題P4可行。

2.3 擾動(dòng)環(huán)境燃料最優(yōu)軌跡規(guī)劃算法框架

擾動(dòng)環(huán)境燃料最優(yōu)軌跡規(guī)劃框架如算法1所示。

算法1 擾動(dòng)環(huán)境燃料最優(yōu)軌跡規(guī)劃方法1. 初始化:k=1, t=t02. 在t0時(shí)刻檢測(cè)著陸器實(shí)際狀態(tài),構(gòu)建并求解修正軌跡預(yù)規(guī)劃問題P4,獲得最優(yōu)控制序列;3. 根據(jù)線性反饋控制律式(28)計(jì)算實(shí)際控制信號(hào),對(duì)著陸器進(jìn)行實(shí)時(shí)控制;4. 檢測(cè)當(dāng)前時(shí)刻t是否滿足重規(guī)劃觸發(fā)條件,若滿足,進(jìn)入步驟5;若不滿足,返回步驟3;5. 令tk=t,構(gòu)建并求解軌跡重規(guī)劃問題P3(k),獲得最優(yōu)控制序列,k=k+1;6. 檢測(cè)著陸任務(wù)是否完成:若未完成,返回步驟3;若完成,則退出算法。

重規(guī)劃觸發(fā)條件可以根據(jù)任務(wù)實(shí)際需求進(jìn)行設(shè)置,可以設(shè)置為一定時(shí)間間隔Δt后自動(dòng)開始重規(guī)劃,也可設(shè)置為誤差相關(guān)的事件觸發(fā)條件。

注意,若直接在P2中代入當(dāng)前狀態(tài)信息進(jìn)行重規(guī)劃,難以構(gòu)造可行解,并可能出現(xiàn)優(yōu)化不收斂的情況,這極大降低了重規(guī)劃的可靠性。但對(duì)于問題P3(k),可以根據(jù)上次規(guī)劃結(jié)果直接構(gòu)造重規(guī)劃問題可行解,當(dāng)求解器因?yàn)槠渌蚯蠼馐r(shí),可依據(jù)上次規(guī)劃結(jié)果構(gòu)造可行解,直至某次重規(guī)劃成功。同時(shí),可設(shè)計(jì)事件觸發(fā)重規(guī)劃條件,進(jìn)一步降低不必要的計(jì)算量。

3 終端時(shí)刻自由

對(duì)于終端時(shí)刻固定問題而言,P3(k)與P4可經(jīng)過線性化、離散化轉(zhuǎn)化為二階錐優(yōu)化問題,現(xiàn)有求解器可在多項(xiàng)式時(shí)間內(nèi)進(jìn)行求解。對(duì)于終端時(shí)刻自由問題,P3(k)與P4中的動(dòng)力學(xué)約束為非線性約束,本節(jié)將以P4為例,對(duì)飛行時(shí)域進(jìn)行映射,進(jìn)一步將問題凸化為有限維二階錐優(yōu)化問題。

3.1 飛行時(shí)域映射

P4中,由于終端時(shí)刻tf是未知量,因此將飛行時(shí)域[0,tf]映射到區(qū)間[0,1],有：

(60)

(61)

進(jìn)而問題轉(zhuǎn)化為P5(注意,此時(shí)求導(dǎo)運(yùn)算為對(duì)τt求導(dǎo))。

P5：

(62)

s.t.

(63)

(64)

(65)

(66)

(67)

(68)

r(τt)∈X

(69)

z(0)=lnmwet,r(0)=r0,v(0)=v0

(70)

z(1)≥lnmdry,r(1)=rf,v(1)=vf

(71)

其中,目標(biāo)函數(shù)式(62)、動(dòng)力學(xué)約束式(63)～式(65) 均為非線性形式。

3.2 目標(biāo)函數(shù)凸化

非線性目標(biāo)函數(shù)可轉(zhuǎn)化為線性目標(biāo)函數(shù)外加一個(gè)非線性不等式約束的形式：

(72)

式中：δ(·)為中間變量。

3.3 離散化

在時(shí)域[0,1]上,使用歐拉法將動(dòng)力學(xué)離散為N段,得到P6。

P6：

(73)

s.t.

(74)

(75)

zi-zi-1=-ασitf

(76)

(77)

(78)

(79)

σitf≤δi

(80)

ri=0,ri=1,…,ri=N∈X

(81)

zi=0=lnmwet,ri=0=r0,vi=0=v0

(82)

zi=N≥lnmdry,ri=N=rf,vi=N=vf

(83)

P6包含非線性等式約束式(74)～式(76)以及非凸不等式約束式(78)～式(80),其余約束為二階錐形式(包含半平面)。因此可寫為標(biāo)準(zhǔn)形式P6*。

P6*：

(84)

s.t.

Ay-b≥κ0

(85)

Cy-d=0

(86)

gi(y)≤0

(87)

hj(y)=0

(88)

式中：y為優(yōu)化變量；c為式(73)轉(zhuǎn)化為標(biāo)準(zhǔn)形式后的加權(quán)系數(shù)；符號(hào)≥κ表示在二階錐κ中的廣義不等式。式(85)表示旋轉(zhuǎn)二次錐約束[17],相關(guān)參數(shù)A、b、κ由約束式(77)和式(81)整理得到；C、d由線性等式約束式(82)～式(83)整理得到；式(87) 為非線性不等式約束的一般形式,函數(shù)gi(y)由約束式(78)～式(80)整理得到；式(88)為非線性等式約束,函數(shù)hj(y)由約束式(74)～式(76) 整理得到。

3.4 線性近似

s.t.

(89)

(90)

(91)

(92)

(93)

(94)

ri=0,ri=1,…,ri=N∈X

zi=0=lnmwet,ri=0=r0,vi=0=v0

zi=N≥lnmdry,ri=N=rf,vi=N=vf

(95)

其中:約束式(95)為信賴域約束,是二階錐約束,確保凸化的精確性,是算法魯棒性的保證[18]。

4 序列凸優(yōu)化方法

序列凸優(yōu)化將通過將原始非凸問題轉(zhuǎn)化為一系列凸優(yōu)化子問題進(jìn)行迭代求解。每次需求解的凸優(yōu)化子問題以上次迭代結(jié)果為基準(zhǔn),第1次迭代需要預(yù)先輸入一個(gè)基準(zhǔn)。

4.1 求解步驟

序列凸優(yōu)化方法求解分為如下5個(gè)步驟：

步驟1初始化設(shè)計(jì)變量y0,及迭代次數(shù)n=0；

步驟2在yn處對(duì)非凸約束進(jìn)行凸化,得到子問題PP6(yn)；

步驟3求解子問題,得到解記為yn+1；

步驟4檢查是否滿足收斂準(zhǔn)則,若是,則算法結(jié)束,輸出解yn+1,若否,則繼續(xù)下一步；

步驟5n=n+1,返回步驟2。

s.t.

Ay-b≥κ0

Cy-d=0

(96)

(97)

4.2 結(jié)果可行性分析

首先驗(yàn)證算法收斂結(jié)果是否滿足P6問題的所有約束,即解是否可行。給出如下定理。

定理3若算法得到的解數(shù)列{yn}收斂于解y*,則y*為P6*可行解。

證明：解數(shù)列{yn}收斂于y*,因此

y*=PP*(y*)

將y*代入PP*(y*)約束得：

(98)

即y*為滿足P6*約束,因此y*為P6*可行解。□

4.3 結(jié)果最優(yōu)性分析

定理4若PP*(yn)嚴(yán)格可行,解數(shù)列{yn}收斂于解y*,則y*為P6*的KKT(Karush-Kuhn-Tucker)解。

證明：P6*的拉格朗日函數(shù)為

(99)

式中：λ0、λi、υ0、υj為拉格朗日乘子。需要證明y*滿足P6*的KKT條件[19],即原始約束式(100)、對(duì)偶約束式(101)、互補(bǔ)松弛式(102)以及拉格朗日函數(shù)導(dǎo)數(shù)為零等式(式(103))。

Ay-b≥κ0,Cy=d,gi(y)≤0,hj(y)=0

(100)

(101)

(102)

(103)

由定理3可知,y*滿足原始約束式(100)。

PP*(yn)的拉格朗日函數(shù)為

(104)

PP*(yn)的對(duì)偶問題為

(105)

展開表示為如下形式：

(106)

由于y*為PP*(yn)最優(yōu)解,因此y*滿足對(duì)偶約束式(101)及拉格朗日函數(shù)導(dǎo)數(shù)為零等式(式(103))。

由于y*滿足PP*(yn)約束,因此以下表達(dá)式成立：

(107)

由于PP*(yn)嚴(yán)格可行,y*滿足PP*(yn)問題的SLATER條件,因此對(duì)偶間隙為0,即原問題最優(yōu)值P*與對(duì)偶問題最優(yōu)值D*的差為0,將式(107) 與式(103)代入得：

(108)

則有

(109)

由于

(110)

根據(jù)對(duì)偶錐定義[19]得：

(111)

因此有

(112)

由于對(duì)偶間隙為0且拉格朗日函數(shù)導(dǎo)數(shù)為0,因此以下表達(dá)式成立：

(113)

將yn=y*代入得：

(114)

又因?yàn)?/p>

(115)

所以

λigi(y*)=0

(116)

結(jié)合式(112),y*滿足互補(bǔ)松弛條件式(102)。

綜上,y*為P6*的KKT解。但實(shí)際求解中,由于精度限制,算法收斂結(jié)果接近P6*的KKT解。□

5 仿真驗(yàn)證

面向擾動(dòng)環(huán)境下火星動(dòng)力下降段自主軌跡規(guī)劃問題,開展數(shù)值仿真研究,針對(duì)具體任務(wù),進(jìn)行仿真分析。硬件環(huán)境為Intel Core i7-6700 3.4 GHz PC,采用Yalmip進(jìn)行問題描述與轉(zhuǎn)化,采用Mosek求解器[20-21]求解,編程環(huán)境為MATLAB 2019b。

5.1 任務(wù)參數(shù)與算法參數(shù)

任務(wù)參數(shù)如表1所示。火星重力加速度g=[0,0,-3.114]TN/kg,初始位置r0=[0,2 000,1 500]Tm,初始速度v0=[0,50,-75]Tm/s,終端位置rf=[0,0,10]Tm,終端速度vf=[0,0,0]Tm/s。

表1 任務(wù)參數(shù)

設(shè)置控制時(shí)間間隔為0.01 s,重規(guī)劃觸發(fā)條件選擇周期性觸發(fā),觸發(fā)時(shí)間間隔為5 s,每次規(guī)劃離散區(qū)間個(gè)數(shù)N=50,求解結(jié)果采用三次樣條插值進(jìn)行進(jìn)一步離散。隨機(jī)擾動(dòng)上界設(shè)置為η=0.05(根據(jù)文獻(xiàn)[22]中火星大氣參數(shù),當(dāng)著陸器阻力系數(shù)為0.9,橫截面積為2 m2時(shí),以80 m/s的速度下降空氣阻力所產(chǎn)生的加速度不超過0.049 8 m/s2),控制器參數(shù)選取為kp=0.2,

kd=2.1,根據(jù)引理1可以計(jì)算出rtube=0.25 m,vtube=0.026 3 m/s。

預(yù)規(guī)劃問題P4初始終端時(shí)間基準(zhǔn)設(shè)置為以火星重力加速度完成下降任務(wù)的時(shí)間的2倍,初始軌跡基準(zhǔn)選擇為連接始末位置的連線軌跡,初始速度基準(zhǔn)為始末速度的連線速度,初始控制量基準(zhǔn)在每個(gè)時(shí)刻設(shè)置為0。重規(guī)劃問題初始時(shí)刻為根據(jù)上次規(guī)劃結(jié)果構(gòu)造的可行解。

5.2 自主規(guī)劃結(jié)果

著陸器從預(yù)先位置觸發(fā),經(jīng)過1次預(yù)規(guī)劃,11次重規(guī)劃后,抵達(dá)終點(diǎn)。預(yù)規(guī)劃經(jīng)過5次SCP迭代收斂,求解耗時(shí)0.935 6 s；11次重規(guī)劃平均耗時(shí)0.565 8 s,最高耗時(shí)0.940 8 s。規(guī)劃結(jié)果如圖1 所示,不同顏色表示不同時(shí)刻的規(guī)劃結(jié)果,符號(hào)“+”表示重規(guī)劃觸發(fā)時(shí)刻的位置,錐形區(qū)域?yàn)榛导s束。通過局部放大可以看出,由于下降過程中存在擾動(dòng),每次重規(guī)劃結(jié)果的軌跡都會(huì)發(fā)生改變,但所有規(guī)劃結(jié)果均可抵達(dá)目標(biāo)位置。位置分量(rX、rY、rZ)、速度分量(vX、vY、vZ)、推力分量(TcX、TcY、TcZ)隨時(shí)間的變化如圖2所示,最終實(shí)際軌跡與多次重規(guī)劃結(jié)果幾乎完全相同。速度大小v變化與推力大小Tc變化如圖3所示,可以看出推力大小滿足約束。

圖1 下降規(guī)劃結(jié)果

圖2 狀態(tài)及控制量變化(重規(guī)劃與實(shí)際對(duì)比)

圖3 速度大小與推力大小(重規(guī)劃與實(shí)際對(duì)比)

實(shí)際位置與上次規(guī)劃最優(yōu)位置誤差隨時(shí)間的變化如圖4(a)所示,速度誤差隨時(shí)間變化如圖4(b) 所示。可以看出位置誤差穩(wěn)定在上述計(jì)算的rtube=0.25 m以內(nèi),僅在重規(guī)劃初始時(shí)候略有超出,這是受重規(guī)劃問題最優(yōu)解離散精度的影響。同理,速度誤差穩(wěn)定在vtube=0.026 3 m/s以內(nèi),符合引理1結(jié)論。

圖4 狀態(tài)誤差變化

由于重規(guī)劃問題與當(dāng)前著陸器狀態(tài)相關(guān),因此重規(guī)劃觸發(fā)后,一般會(huì)引起誤差突變,在下次重規(guī)劃觸發(fā)前,依靠線性反饋控制使誤差逐漸收斂至有界,直至下一次重規(guī)劃問題觸發(fā)。

自主規(guī)劃的最優(yōu)值與最優(yōu)終端時(shí)間情況如圖5 所示。最優(yōu)值呈現(xiàn)線性下降趨勢(shì),這是由于重規(guī)劃中的代價(jià)函數(shù)沒有考慮已經(jīng)完成的軌跡代價(jià)。最后終端時(shí)間大體相同,但相比于預(yù)規(guī)劃,略有減少,這是因?yàn)樵谥匾?guī)劃中修正了前次規(guī)劃中次優(yōu)的部分軌跡。

圖5 自主規(guī)劃的最優(yōu)值與最優(yōu)終端時(shí)間

5.3 規(guī)劃結(jié)果對(duì)比

為進(jìn)一步對(duì)比有無重規(guī)劃的差異,在相同設(shè)置下進(jìn)行了僅有預(yù)規(guī)劃與反饋控制的仿真。僅有預(yù)規(guī)劃時(shí),下降時(shí)間為60.83 s,消耗燃料質(zhì)量為292.035 5 kg；有重規(guī)劃時(shí),下降時(shí)間縮短至56.7 s,消耗燃料質(zhì)量為272.607 0 kg。

實(shí)際上,在擾動(dòng)作用下,預(yù)規(guī)劃出的軌跡可能在出發(fā)不久就不是最優(yōu)軌跡了,通過不斷地重規(guī)劃可以克服擾動(dòng)影響并在一定程度上提高實(shí)際下降結(jié)果的最優(yōu)性,減少下降時(shí)間。

將提出方法與文獻(xiàn)[15]方法(PFGMPG)進(jìn)行對(duì)比。重規(guī)劃觸發(fā)間隔為5 s,控制間隔為0.01 s,狀態(tài)懲罰項(xiàng)系數(shù)設(shè)置為1。進(jìn)行2組PFGMPG仿真,PFGMPG 1組中含有反饋控制律；PFGMPG 2組中不含反饋控制律,由規(guī)劃結(jié)果直接控制。對(duì)比如圖6所示。

從圖6(a)中看出,含有反饋控制律的PFGMPG 1組在實(shí)際下降時(shí)推力違反了約束,這是

由于規(guī)劃時(shí)忽略了控制環(huán)節(jié)的耦合作用,沒有給內(nèi)環(huán)控制器預(yù)留足夠的抗擾控制量。從圖6(b)中可以看出,在2次規(guī)劃間隔內(nèi),提出方法位置誤差在反饋控制律的作用下呈現(xiàn)逐漸減少趨勢(shì),而不加反饋控制律的PFGMPG 2組由于沒有反饋,位置誤差隨著擾動(dòng)影響逐漸增加；從圖6(c)中可以看出不加反饋控制律的PFGMPG 2組的速度誤差在2次規(guī)劃間隔中可能超出約束,這是2次規(guī)劃間隔中開環(huán)控制的直接結(jié)果。

6 結(jié) 論

針對(duì)擾動(dòng)環(huán)境下火星著陸動(dòng)力下降段自主軌跡規(guī)劃問題,本文考慮終端時(shí)刻自由情況,建立預(yù)規(guī)劃問題與重規(guī)劃問題模型,使用SCP方法實(shí)現(xiàn)軌跡規(guī)劃。給出了預(yù)規(guī)劃問題可行的必要性條件,為控制器參數(shù)的選取提供參考。得到以下結(jié)論：

1) 在連續(xù)時(shí)間系統(tǒng)下,當(dāng)預(yù)規(guī)劃問題可行時(shí),重規(guī)劃問題可行,且可以構(gòu)造出重規(guī)劃問題可行解。

2) 誤差線性反饋控制律參數(shù)在滿足一定條件時(shí)可將實(shí)際軌跡約束在最優(yōu)軌跡附近。

3) 在本文提出的重規(guī)劃框架下,每次規(guī)劃的最優(yōu)值呈線性下降趨勢(shì),實(shí)際下降軌跡的時(shí)間略長于預(yù)規(guī)劃出的下降時(shí)間結(jié)果。

4) 仿真結(jié)果表明,實(shí)際下降軌跡既可以保證在2次規(guī)劃間隔內(nèi)誤差有界,又可保證推力滿足推力大小約束。