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

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

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

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

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

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

針對擾動環(huán)境,常采用2種方式進(jìn)行克服：① 采 用預(yù)先規(guī)劃加反饋跟蹤控制的思想進(jìn)行;② 在 下降過程中不斷進(jìn)行重規(guī)劃。其中,前者在擾動作用下,規(guī)劃軌跡在實際中可能不可行。后者雖然從規(guī)劃層次上考慮了擾動,相比于前者具由更強(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),會出現(xiàn)重規(guī)劃頻率與計算能力的矛盾。增加內(nèi)環(huán)反饋控制后可降低重規(guī)劃頻率,但同時需要考慮規(guī)劃問題與內(nèi)環(huán)控制的耦合。本文針對終端自由擾動環(huán)境下火星著陸動力下降段軌跡規(guī)劃問題,結(jié)合Tube-MPC思想[16]設(shè)計反饋控制律,考慮規(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)一步提高抗擾性能。同時給出了重規(guī)劃問題可行的必要條件,為實際工程應(yīng)用提供參考。

1 問題描述

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

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

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

1.1 著陸器動力學(xué)

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

(1)

(2)

1.2 耗能最優(yōu)

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

maxm(tf)-m(0)

(3)

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

(4)

1.3 過程約束

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

(5)

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

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)

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

(9)

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

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

燃料最優(yōu)火星著陸動力下降段軌跡規(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]對優(yōu)化問題P1進(jìn)行凸化處理與變量替換,通過引入新控制變量Γ,對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)行變量替換后,動力學(xué)約束式(1)變?yōu)榫€性形式如式(15)所示。

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

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

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

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

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

(21)

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

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

(22)

s.t.

(23)

(24)

(25)

(26)

r(t)∈X

(27)

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

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

(28)

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

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

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

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

(29)

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

(30)

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

(31)

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

(32)

代入線性反饋控制律得：

(33)

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

(34)

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

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

C1eλ1t+C2eλ2t

(35)

有

(36)

第1項展開為

(37)

(38)

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

(39)

放縮得到

(40)

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

(41)

因此式(36)可放縮為

(42)

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

(43)

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

(44)

(45)

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

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

(46)

所以

(47)

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

(48)

有

(49)

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

(50)

(51)

(52)

(53)

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

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

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

(54)

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

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

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

(55)

(56)

移項并代入uemax=kp⊙rtube+kd⊙vtube得：

(57)

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

(58)

□

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

(59)

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

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

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

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

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

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

3 終端時刻自由

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

3.1 飛行時域映射

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

(60)

(61)

進(jìn)而問題轉(zhuǎn)化為P5(注意,此時求導(dǎo)運(yùn)算為對τ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)、動力學(xué)約束式(63)～式(65) 均為非線性形式。

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

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

(72)

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

3.3 離散化

在時域[0,1]上,使用歐拉法將動力學(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ù)；符號≥κ表示在二階錐κ中的廣義不等式。式(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ù)先輸入一個基準(zhǔn)。

4.1 求解步驟

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

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

步驟2在yn處對非凸約束進(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é)果可行性分析

首先驗證算法收斂結(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)、對偶約束式(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)的對偶問題為

(105)

展開表示為如下形式：

(106)

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

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

(107)

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

(108)

則有

(109)

由于

(110)

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

(111)

因此有

(112)

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

(113)

將yn=y*代入得：

(114)

又因為

(115)

所以

λigi(y*)=0

(116)

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

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

5 仿真驗證

面向擾動環(huán)境下火星動力下降段自主軌跡規(guī)劃問題,開展數(shù)值仿真研究,針對具體任務(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è)置控制時間間隔為0.01 s,重規(guī)劃觸發(fā)條件選擇周期性觸發(fā),觸發(fā)時間間隔為5 s,每次規(guī)劃離散區(qū)間個數(shù)N=50,求解結(jié)果采用三次樣條插值進(jìn)行進(jìn)一步離散。隨機(jī)擾動上界設(shè)置為η=0.05(根據(jù)文獻(xiàn)[22]中火星大氣參數(shù),當(dāng)著陸器阻力系數(shù)為0.9,橫截面積為2 m2時,以80 m/s的速度下降空氣阻力所產(chǎn)生的加速度不超過0.049 8 m/s2),控制器參數(shù)選取為kp=0.2,

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

預(yù)規(guī)劃問題P4初始終端時間基準(zhǔn)設(shè)置為以火星重力加速度完成下降任務(wù)的時間的2倍,初始軌跡基準(zhǔn)選擇為連接始末位置的連線軌跡,初始速度基準(zhǔn)為始末速度的連線速度,初始控制量基準(zhǔn)在每個時刻設(shè)置為0。重規(guī)劃問題初始時刻為根據(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迭代收斂,求解耗時0.935 6 s；11次重規(guī)劃平均耗時0.565 8 s,最高耗時0.940 8 s。規(guī)劃結(jié)果如圖1 所示,不同顏色表示不同時刻的規(guī)劃結(jié)果,符號“+”表示重規(guī)劃觸發(fā)時刻的位置,錐形區(qū)域為滑降約束。通過局部放大可以看出,由于下降過程中存在擾動,每次重規(guī)劃結(jié)果的軌跡都會發(fā)生改變,但所有規(guī)劃結(jié)果均可抵達(dá)目標(biāo)位置。位置分量(rX、rY、rZ)、速度分量(vX、vY、vZ)、推力分量(TcX、TcY、TcZ)隨時間的變化如圖2所示,最終實際軌跡與多次重規(guī)劃結(jié)果幾乎完全相同。速度大小v變化與推力大小Tc變化如圖3所示,可以看出推力大小滿足約束。

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

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

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

實際位置與上次規(guī)劃最優(yōu)位置誤差隨時間的變化如圖4(a)所示,速度誤差隨時間變化如圖4(b) 所示。可以看出位置誤差穩(wěn)定在上述計算的rtube=0.25 m以內(nèi),僅在重規(guī)劃初始時候略有超出,這是受重規(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ā)后,一般會引起誤差突變,在下次重規(guī)劃觸發(fā)前,依靠線性反饋控制使誤差逐漸收斂至有界,直至下一次重規(guī)劃問題觸發(fā)。

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

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

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

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

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

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

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

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

6 結(jié) 論

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

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

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

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

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