鄧雁鵬,穆榮軍,彭 娜,吳 鵬

(1. 哈爾濱工業(yè)大學(xué)航天學(xué)院,哈爾濱 150001；2. 上海航天電子技術(shù)研究所,上海 201109)

0 引 言

為順利實(shí)施月球著陸探測器自主避障、保證精確安全著陸,通常將月球著陸過程細(xì)分為4個(gè)階段,如圖1所示。月球著陸器在動(dòng)力下降段之前處于遠(yuǎn)地點(diǎn)100 km、近地點(diǎn)15 km的轉(zhuǎn)移軌道上。自近地點(diǎn)附近開始動(dòng)力下降,然后開始進(jìn)行姿態(tài)調(diào)整,經(jīng)過10～20 s的姿態(tài)調(diào)整后開始避障。動(dòng)力下降段的主要目標(biāo)為盡可能降低著陸器的高度和速度,通常歷時(shí)約500 s,航程500 km左右；著陸器高度從15 km降低到3 km,速度將從1.7 km/s降低到0左右。該階段通常會(huì)消耗掉著陸器約80%的燃料,因此該段制導(dǎo)算法要求,在著陸器具有初始狀態(tài)誤差和模型不確定性的情況下盡量減少燃料消耗并在終端達(dá)到后續(xù)姿態(tài)調(diào)整、避障著陸的入口狀態(tài)要求。由于著陸器在經(jīng)過環(huán)月軌道離軌后存在霍曼轉(zhuǎn)移軌道偏差和點(diǎn)火時(shí)間偏差,進(jìn)而造成初始位置偏差。因此,本文充分考慮初始位置偏差的影響,提出了一種基于最優(yōu)軌跡序列凸優(yōu)化的在線制導(dǎo)方法,在節(jié)省燃料的基礎(chǔ)上實(shí)現(xiàn)高精度動(dòng)力下降制導(dǎo)。

圖1 典型月面著陸過程Fig.1 Typical lunar landing process

月面著陸制導(dǎo)算法主要分為三種類型：重力轉(zhuǎn)彎制導(dǎo)、顯式制導(dǎo)和標(biāo)稱軌跡制導(dǎo)。其中重力轉(zhuǎn)彎制導(dǎo)算法原理簡單,但其不具備燃料經(jīng)濟(jì)性,同時(shí)對末端著陸位置沒有約束,故無法實(shí)現(xiàn)精確著陸。

顯式制導(dǎo)算法能夠?qū)χ懳恢脤?shí)現(xiàn)精確控制,同時(shí)由于其具有自主性,實(shí)時(shí)性與燃料消耗次優(yōu)性,因此近些年得到了廣泛的研究與應(yīng)用。文獻(xiàn)[1]依據(jù)動(dòng)力學(xué)模型將月面軟著陸下降制導(dǎo)分為縱向平面運(yùn)動(dòng)和橫向運(yùn)動(dòng)并分別設(shè)計(jì)制導(dǎo)律,其縱向平面運(yùn)動(dòng)采用改進(jìn)顯式制導(dǎo)算法,橫向運(yùn)動(dòng)采用ZEM/ZEV制導(dǎo)律并對脈寬脈頻進(jìn)行脈沖調(diào)制。文獻(xiàn)[2]針對動(dòng)力下降段的制導(dǎo)問題,在顯式制導(dǎo)算法的基礎(chǔ)上,提出了自適應(yīng)動(dòng)力顯式制導(dǎo)方法。通過對動(dòng)力下降段目標(biāo)的自適應(yīng)調(diào)整使其滿足接近段入口條件,通過在線實(shí)時(shí)估計(jì)推重比與比沖等參數(shù)實(shí)現(xiàn)制導(dǎo)參數(shù)的自適應(yīng)調(diào)整,由此實(shí)現(xiàn)月面高精度高可靠性著陸。

顯式制導(dǎo)雖算法簡潔、應(yīng)用廣泛,但通常未考慮月面著陸制導(dǎo)過程中所關(guān)注的燃料最優(yōu)性問題。相對顯式制導(dǎo)算法而言,基于最優(yōu)化原理的數(shù)值優(yōu)化算法按照一定的目標(biāo)及約束條件設(shè)計(jì)標(biāo)稱軌道。在飛行過程中或根據(jù)偏差進(jìn)行實(shí)時(shí)最優(yōu)軌跡規(guī)劃,或通過實(shí)時(shí)比較著陸器的位置和速度測量信息,導(dǎo)引著陸器飛向著陸目標(biāo)點(diǎn)。該方法通常消耗更少的燃料,具有更強(qiáng)的魯棒性。但其通常具有計(jì)算量較大,實(shí)時(shí)性難以保證等缺點(diǎn)。

傳統(tǒng)軌跡優(yōu)化方法分為直接法和間接法。間接法利用極小值原理,通過引入拉格朗日乘子法進(jìn)行求解。但其初值難以獲得,該解法對于初值又比較敏感,限制了該方法的應(yīng)用。直接法是將著陸軌跡優(yōu)化問題轉(zhuǎn)化為多種路徑約束和邊界約束下的靜態(tài)參數(shù)優(yōu)化問題,進(jìn)而采用數(shù)值優(yōu)化方法求解。數(shù)值優(yōu)化算法主要分為二次規(guī)劃算法和凸優(yōu)化兩類。第一類是將原問題通過偽譜法進(jìn)行離散化,將其轉(zhuǎn)化為非線性規(guī)劃問題,利用非線性規(guī)劃方法進(jìn)行求解。例如文獻(xiàn)[4]利用高斯偽譜法實(shí)現(xiàn)了著陸器在二維情況下的軟著陸軌道規(guī)劃；在此基礎(chǔ)上文獻(xiàn)[5]利用高斯偽譜法+序列二次規(guī)劃算法設(shè)計(jì)最優(yōu)著陸軌跡。文獻(xiàn)[6]提出了一種具有路標(biāo)點(diǎn)約束的多階段非線性優(yōu)化制導(dǎo)方法。這一算法在最大程度減少燃料消耗的同時(shí),可確保著陸器在著陸過程中通過兩個(gè)規(guī)劃的路標(biāo)點(diǎn),使視覺探測器能夠在通過這些路標(biāo)點(diǎn)時(shí)進(jìn)行著陸點(diǎn)檢測,進(jìn)一步確保任務(wù)的安全性。

偽譜法+序列二次規(guī)劃算法隨著問題規(guī)模增大,計(jì)算量將顯著增長,規(guī)劃效率難以滿足實(shí)時(shí)應(yīng)用需求。相比于偽譜法+序列二次規(guī)劃算法,基于凸優(yōu)化的軌跡規(guī)劃算法通常具有更小的計(jì)算量,同時(shí)凸問題能夠保證局部最優(yōu)解即為全局最優(yōu)解。因此,第二類算法將非線性軌跡規(guī)劃問題轉(zhuǎn)換為凸問題后,再進(jìn)行離散化,進(jìn)而對該問題進(jìn)行求解。例如文獻(xiàn)[8]針對小型天體弱引力場所導(dǎo)致的動(dòng)力學(xué)非線性問題,在凸優(yōu)化過程中采用局部線性化的方式對原問題進(jìn)行凸化處理。除采用局部線性化外,在著陸軌跡規(guī)劃中,還常采用變量松弛的方式將原非凸問題轉(zhuǎn)化為凸問題求解。文獻(xiàn)[9-10]分別針對月面與行星著陸問題,通過松弛變量將原非凸問題轉(zhuǎn)化為二階錐規(guī)劃(Second-order cone programming, SOCP)問題,引入內(nèi)點(diǎn)法求解,并通過數(shù)值模擬證明了所提出方法的性能。針對非線性氣動(dòng)阻力不可忽略的可重復(fù)使用運(yùn)載器著陸過程中的燃料優(yōu)化問題,文獻(xiàn)[11]提出了一種快速彈道優(yōu)化方法。該算法通過同倫迭代策略來凸化火箭的非線性系統(tǒng)動(dòng)力學(xué),通過逐次迭代逐步逼近時(shí)變阻力剖面。該算法不需要參考軌跡或初始猜測,極大地增強(qiáng)了算法的自主性。除采用松弛算法對原問題進(jìn)行無損凸化外,文獻(xiàn)[12]針對火箭高空再入定點(diǎn)回收問題,提出了一種考慮氣動(dòng)力和推力控制的多約束軌跡優(yōu)化方法,使用無損凸化對升力約束和推力約束進(jìn)行松弛,最終將問題轉(zhuǎn)化為SOCP問題進(jìn)行序列迭代求解。但其將重力加速度建模為常矢量,且未考慮慣性加速度項(xiàng)的影響?；诖?針對月面著陸動(dòng)力下降段最優(yōu)軌跡優(yōu)化問題,有必要建立強(qiáng)非線性時(shí)變項(xiàng)影響下的精細(xì)化動(dòng)力學(xué)模型并對該方法做出對應(yīng)調(diào)整,實(shí)現(xiàn)動(dòng)力下降段最優(yōu)軌跡的在線快速迭代求解。

以上文獻(xiàn)對于時(shí)變加速度剖面通常采用兩種處理方式。部分文獻(xiàn)采用在月心慣性坐標(biāo)系下建模的方式消除動(dòng)力學(xué)模型中的慣性力項(xiàng),如文獻(xiàn)[1,13]。在此坐標(biāo)系下的著陸動(dòng)力學(xué)模型較為簡單,但其三個(gè)方向上的坐標(biāo)不能直接反映著陸器高度、航程等信息,也難以添加相應(yīng)的姿態(tài)約束以滿足導(dǎo)航與避障傳感器的需求。同時(shí),如圖2所示,著陸器飛行過程中重力加速度的大小雖基本不變,但由于月球曲率的存在,其方向在不斷改變；如不對這一部分的時(shí)變加速度進(jìn)行精細(xì)化建模,將會(huì)影響最終制導(dǎo)精度及燃料消耗。此外,部分文獻(xiàn)選擇在月面固連坐標(biāo)系下建模,將該坐標(biāo)系下的重力加速度作為定常矢量處理,同時(shí)忽略在此坐標(biāo)系下難以建模和補(bǔ)償?shù)膽T性加速度項(xiàng)的影響。在此坐標(biāo)系下的動(dòng)力學(xué)模型中,著陸器三軸坐標(biāo)分別代表了高度與水平兩個(gè)方向的航程,三軸姿態(tài)角可直觀反映著陸器相對于月面的姿態(tài),具體物理意義明確,從而易于根據(jù)導(dǎo)航與避障傳感器需求添加相應(yīng)的姿態(tài)約束。但由于建模時(shí)通常將重力加速度考慮為常矢量,同時(shí)月面曲率與自轉(zhuǎn)引起的慣性力項(xiàng)的非線性特性難以求解,一般僅能應(yīng)用于慣性力較小可忽略、位置變動(dòng)較小的末制導(dǎo)階段。

圖2 著陸器動(dòng)力學(xué)模型Fig.2 Dynamics model of the lunar lander

同時(shí),經(jīng)典凸優(yōu)化算法對于非線性項(xiàng)一般通過序列凸優(yōu)化算法進(jìn)行處理,即通過一系列凸近似手段將原復(fù)雜非凸問題轉(zhuǎn)化為凸問題,從而快速獲得原優(yōu)化問題的近似最優(yōu)解。盡管該方法計(jì)算效率較高、實(shí)用性較強(qiáng),但仍需要依賴初始參考軌跡,利用Taylor展開使模型實(shí)現(xiàn)線性化,將其轉(zhuǎn)化為凸問題后開始迭代求解。與此同時(shí),該算法需要對該問題施加額外的信賴域約束,以保證線性化的有效性。此外,該算法的收斂性與參考軌跡的質(zhì)量具有一定的相關(guān)性,這類似于非線性規(guī)劃中的初值敏感性問題,目前很少有文獻(xiàn)針對這一內(nèi)容有所討論。序列凸優(yōu)化依賴于初始參考軌跡的特性弱化了局部最優(yōu)即為全局最優(yōu)這一凸優(yōu)化的主要優(yōu)點(diǎn)。因此,為提高算法的求解精度及自主性,本文提出改進(jìn)序列凸優(yōu)化算法,將時(shí)變加速度剖面同倫地添加到原非線性模型中序列迭代求解,使得該算法不依賴初始參考軌跡,增強(qiáng)了算法的自主性和多種應(yīng)用場景適應(yīng)能力。

本文結(jié)合凸優(yōu)化計(jì)算量小、實(shí)時(shí)性高的優(yōu)點(diǎn),提出了一種基于序列凸優(yōu)化的月球探測器動(dòng)力下降段在線制導(dǎo)方法。綜合考慮月球探測器動(dòng)力下降段在線軌跡規(guī)劃的快速性、制導(dǎo)精度與燃料消耗等需求,針對月球著陸過程中引力加速度與慣性加速度時(shí)變特性及模型復(fù)雜、難以求解等問題,首先在月面固聯(lián)坐標(biāo)系下建立了考慮月面曲率和月球自轉(zhuǎn)的精細(xì)化著陸動(dòng)力學(xué)模型,在此基礎(chǔ)上提出考慮時(shí)變加速度剖面的序列凸優(yōu)化算法。區(qū)別于經(jīng)典序列凸優(yōu)化算法通過對動(dòng)力學(xué)模型的序列線性化將原問題轉(zhuǎn)化為凸問題并進(jìn)行求解,本文以燃料消耗為性能指標(biāo),綜合考慮推力大小等約束條件,針對月面著陸飛行過程中,慣性加速度剖面逐漸趨于0,重力加速度逐漸趨于定常的特點(diǎn)。在提出迭代算法的基礎(chǔ)上,首先求取考慮重力加速度為常數(shù),慣性加速度為0時(shí)原問題的解。然后保持問題為凸問題,將時(shí)變加速度同倫地添加到該問題中進(jìn)行序列迭代求解。由于該方法在不需要求解Jacobi矩陣,減少了計(jì)算量；同時(shí)非線性項(xiàng)的凸化不是基于線性化,因此不需要參考軌跡或初始猜測,增強(qiáng)了算法的自主性,有利于在線應(yīng)用。最終,通過與凸優(yōu)化、多項(xiàng)式制導(dǎo)算法進(jìn)行仿真對比,并針對多種初始偏差開展打靶分析,驗(yàn)證了算法性能。

1 精細(xì)化月球著陸動(dòng)力學(xué)模型

1.1 月球著陸器動(dòng)力下降段動(dòng)力學(xué)模型

本文考慮到在月面固連參考坐標(biāo)系下三個(gè)坐標(biāo)軸方向上的坐標(biāo)可直接反映著陸器高度航程等信息,也易于添加相應(yīng)的姿態(tài)約束以滿足導(dǎo)航與避障傳感器的需求的優(yōu)點(diǎn)。在月面固連參考系下建立精細(xì)化月球著陸動(dòng)力學(xué)模型,并通過序列凸優(yōu)化算法實(shí)時(shí)的估計(jì)時(shí)變加速度剖面。此時(shí),著陸器的動(dòng)力學(xué)模型為：

(1)

式中：(),()分別為著陸器的位置和速度；()為著陸器推力；()為著陸器質(zhì)量；為推力器比沖；為標(biāo)稱地球重力加速度；(),()分別為慣性加速度與重力加速度,它們均隨時(shí)間變化。因此,采用考慮時(shí)變加速度剖面的著陸器動(dòng)力學(xué)模型,通過序列凸優(yōu)化算法實(shí)現(xiàn)對(),()的估計(jì)與補(bǔ)償,進(jìn)而保證算法的最優(yōu)性。

由式(1),有

(2)

兩邊積分可得著陸器剩余質(zhì)量：

(3)

為使得最終質(zhì)量()達(dá)到最大,需要以下目標(biāo)函數(shù)達(dá)到最小,

(4)

在著陸過程中,為避免著陸器撞擊到月面,需要引入與著陸安全性相關(guān)的約束。首先引入高度約束：

()≥

(5)

式中：為著陸器被允許到達(dá)的最低高度,低于此高度則有撞擊月面的風(fēng)險(xiǎn)。

基于安全考慮,在末制導(dǎo)階段通常還需要引入視線角約束。令

()=()-

(6)

則要求

(7)

式(7)要求著陸器在著陸軌跡上任意位置處相對于目標(biāo)點(diǎn)的連線與水平面的夾角要大于。由于本文主要針對著陸器動(dòng)力下降段應(yīng)用場景,該段末了點(diǎn)的高度足夠高,通常遠(yuǎn)高于著陸場地形起伏的程度；同時(shí)考慮到著陸器在動(dòng)力下降段航程較長,相對航程而言下降高度較小,視線角無法起到約束效果。因此,在動(dòng)力下降段可忽略式(7),此時(shí)并不會(huì)影響著陸安全性。

除著陸安全性帶來的約束之外,還需考慮推力約束與剩余質(zhì)量約束,其中推力約束如下：

(8)

式中：,為推力器能提供的最小與最大推力。

剩余質(zhì)量約束如下：

()≥

(9)

最后還有著陸器初末狀態(tài)約束：

(10)

以上動(dòng)力學(xué)模型和約束共同構(gòu)成了如下非凸優(yōu)化問題：

(11)

1.2 無損凸化

在式(1)表示的動(dòng)力學(xué)模型中,由于()出現(xiàn)在分母上,因此優(yōu)化問題式(11)為非凸模型。為使用凸優(yōu)化方法求解,需要將模型進(jìn)行凸化,轉(zhuǎn)化為一個(gè)凸問題,具體而言為一個(gè)SOCP問題。參考文獻(xiàn)[16]凸化方法,首先引入對數(shù)質(zhì)量

()=ln()

(12)

同時(shí)引入新的控制變量

(13)

在此基礎(chǔ)上,動(dòng)力學(xué)模型變?yōu)?/p>

(14)

目標(biāo)函數(shù)變?yōu)椋?/p>

(15)

推力幅值約束為

e-()≤()≤e-()

(16)

剩余質(zhì)量約束變?yōu)椋?/p>

()≥

(17)

還需引入控制量約束

(18)

當(dāng)取得最優(yōu)解時(shí)有

(19)

為了避免指數(shù)形非凸約束,將式(16)進(jìn)行Taylor展開有

(20)

式中:()滿足

()=ln(-)

(21)

式中：為著陸器初始質(zhì)量。

對數(shù)質(zhì)量()還需要滿足以下約束

ln(-)≤()≤ln(-)

(22)

故原問題化為以下SOCP問題

(23)

2 問題的離散化

對式(23),還需進(jìn)行離散化再進(jìn)行求解。定義狀態(tài)變量如下

(24)

動(dòng)力學(xué)模型式(14)變?yōu)?/p>

(25)

式中：

離散化的方法為：將時(shí)間區(qū)間[0,]離散化為個(gè)區(qū)間,在每個(gè)區(qū)間內(nèi),假設(shè)第個(gè)時(shí)間區(qū)間內(nèi)的控制量為()滿足

(26)

式中:(·)為特定的基函數(shù)。在本文中,取基函數(shù)(·)為常數(shù)函數(shù)。離散化的動(dòng)力學(xué)模型為

+1=++(+I)

(27)

式中：,,,即為(),(),(),I()。

=Δ+,=Δ,=Δ

式中:Δ為離散化的時(shí)間步長Δ=。

離散化后,問題變?yōu)?/p>

(28)

式中：=(),=(),=(),0=(),3=()。式(28)中的第一式為著陸性能指標(biāo),描述為各節(jié)點(diǎn)加速度幅值的和；第二式為動(dòng)力學(xué)約束；第三、四式為推力大小約束；第五、六式為質(zhì)量約束；第七式為高度約束。

3 考慮時(shí)變加速度剖面的制導(dǎo)模型

慣性加速度的來源為離心力與科里奧利力,其滿足：

(29)

式中：()為著陸器與月球的相對速度;()為著陸器相對于月球的角速度;為月球的自轉(zhuǎn)角速度;為的叉乘矩陣形式：

在本問題中,前兩個(gè)量(),()均為時(shí)變量,因此慣性加速度()也為時(shí)變量。

除慣性加速度()之外,月球引力加速度在月面固連坐標(biāo)系下也有微小的變化。如考慮月球引力模型到項(xiàng),則引力加速度隨著著陸器緯度與高度的變化會(huì)有微小的變化。在本文中,取月球引力常數(shù)=4.902801056×10m/s,歸一化項(xiàng)引力攝動(dòng)系數(shù)為=-9.08×10。

根據(jù)上述模型,對時(shí)變加速度剖面也采用離散化處理,具體計(jì)算過程如圖3所示。在每一個(gè)離散點(diǎn)處,均采用上一制導(dǎo)周期的計(jì)算結(jié)果作為標(biāo)稱的()與()初值,計(jì)算時(shí)變加速度剖面,并將該加速度剖面用于本制導(dǎo)周期凸優(yōu)化中；在進(jìn)行一次凸優(yōu)化之后,將本次凸優(yōu)化結(jié)果中的(),()與此前(),()進(jìn)行比較。如果二者差值在一定范圍內(nèi),則認(rèn)為收斂,輸出凸優(yōu)化結(jié)果；否則根據(jù)此次結(jié)果重新計(jì)算時(shí)變加速度剖面(),()并再次進(jìn)行凸優(yōu)化。如此反復(fù)進(jìn)行序列凸優(yōu)化,直至收斂,完成本次迭代周期的計(jì)算。

圖3 基于序列凸優(yōu)化的在線制導(dǎo)算法流程圖Fig.3 Flow chart of the online guidance algorithm based on sequence convex optimization method

4 仿真校驗(yàn)

以嫦娥系列著陸器典型參數(shù)為例,假設(shè)著陸器初始質(zhì)量為3000 kg,發(fā)動(dòng)機(jī)最大推力為7500 N,最小推力900 N,真空比沖309 s；初始速度=1700 m/s,==0,末速度=50 m/s,=0,=-50 m/s；初始高度15 km,最終高度3 km。本文將所提出的序列凸優(yōu)化算法與常規(guī)凸優(yōu)化、多項(xiàng)式制導(dǎo)相比,以表明其在燃料消耗等方面的優(yōu)勢；同時(shí)將本文方法與LGR偽譜法進(jìn)行對比,以明確其在計(jì)算效率上的優(yōu)越性。本文數(shù)值仿真計(jì)算機(jī)采用Intel Core i7-10700U 2.90 GHz CPU,CVX工具箱和GPOPS-II求解器進(jìn)行解算。

4.1 仿真結(jié)果

各方法制導(dǎo)結(jié)果由表1與圖4～圖9所示。由序列凸優(yōu)化方法,LGR偽譜法,凸優(yōu)化方法與多項(xiàng)式制導(dǎo)均實(shí)現(xiàn)了制導(dǎo)目標(biāo)。參照圖4～圖9,LGR偽譜法與序列凸優(yōu)化算法計(jì)算結(jié)果一致、同樣取得了最優(yōu)解。序列凸優(yōu)化,其中序列凸優(yōu)化與多項(xiàng)式制導(dǎo)方法用時(shí)578 s,凸優(yōu)化用時(shí)643 s,節(jié)約65 s。這是由于凸優(yōu)化模型未考慮離心慣性力的影響,因而離心慣性力在著陸器飛行過程中作為干擾項(xiàng)存在,在每個(gè)制導(dǎo)周期重新規(guī)劃最優(yōu)著陸曲線時(shí),為使凸優(yōu)化問題有界,對著陸時(shí)間調(diào)整導(dǎo)致的,由于用時(shí)較短,序列凸優(yōu)化方法與多項(xiàng)式制導(dǎo)方法的燃料消耗也顯著小于凸優(yōu)化方法。

表1 各制導(dǎo)方法結(jié)果對比Table 1 Comparison of results of each guidance method

圖4 質(zhì)量變化曲線Fig.4 Mass variation

圖5 推力變化曲線Fig.5 Thrust variation

圖6 軟著陸過程速度曲線Fig.6 Velocity during soft landing

圖7 軟著陸過程高度曲線Fig.7 Altitude during soft landing

圖8 軟著陸過程航程曲線Fig.8 Flying range during soft landing

圖9 軟著陸過程俯仰角曲線Fig.9 Pitch angle during soft landing

如圖4所示,采用序列凸優(yōu)化算法最終剩余質(zhì)量1666.4 kg,采用多項(xiàng)式制導(dǎo)方法最終剩余質(zhì)量1664.2 kg,序列凸優(yōu)化方法節(jié)約燃料2.2 kg。

如圖5所示,采用序列凸優(yōu)化算法構(gòu)造的制導(dǎo)指令推力值近似符合傳統(tǒng)變分法推導(dǎo)的max-min模式,推力值近似在最大與最小之間切換,切換次數(shù)為2。而采用凸優(yōu)化方法則完全符合max-min模式,但由于每一個(gè)制導(dǎo)周期重新規(guī)劃的制導(dǎo)指令與上一周期的制導(dǎo)指令有較大偏差,造成推力指令在最大與最小之間反復(fù)切換,這對執(zhí)行機(jī)構(gòu)執(zhí)行該制導(dǎo)指令帶來一定困難；對于多項(xiàng)式制導(dǎo)方法,推力值為一條平滑曲線,其后段與序列凸優(yōu)化方法幾乎重合,這由于在動(dòng)力下降段末期,時(shí)變的慣性加速度帶來的影響幾乎可以忽略,而多項(xiàng)式方法在此時(shí)也具有最優(yōu)性,因此二者重合。

三軸速度如圖6所示。四種算法在東向速度曲線和北向速度曲線上表現(xiàn)相似,而在天向速度方面,凸優(yōu)化算法與其余算法相差較大,這是因?yàn)橥箖?yōu)化算法沒有考慮離心加速度的影響；同時(shí)序列凸優(yōu)化算法形成的速度剖面具有分段光滑特性,這是由于max-min的控制模式導(dǎo)致的。圖7為高度曲線,由圖可以看出采用凸優(yōu)化方法的飛行高度顯著大于其余方法,這是由于凸優(yōu)化算法未考慮離心加速度導(dǎo)致的；而序列凸優(yōu)化算法的飛行高度僅略高于多項(xiàng)式算法。圖8為著陸器航程曲線,可以看到在本文所示的典型參數(shù)下,四種方法航程一致,均為557 km,結(jié)合圖7可以看出著陸器典型下降高度為12 km,對應(yīng)中的為1.2°,如此小的起不到約束效果。因此,本算例進(jìn)一步說明了在動(dòng)力下降段,可以忽略約束式(7),此時(shí)并不會(huì)影響著陸安全性。

俯仰角曲線如圖9所示,由圖可知,凸優(yōu)化算法的俯仰角在制導(dǎo)各個(gè)周期反復(fù)的大幅度切換,這對控制系統(tǒng)實(shí)現(xiàn)帶來了一定的問題；而對于序列凸優(yōu)化方法,在飛行過程中俯仰角僅有兩次幅度較小的突變,控制系統(tǒng)較容易實(shí)現(xiàn)；而多項(xiàng)式方法則在全程呈連續(xù)變化。

4.2 求解效率分析

按第4節(jié)仿真初始條件,對改進(jìn)序列凸優(yōu)化算法和LGR偽譜法的求解效率進(jìn)行對比驗(yàn)證。

由圖4～圖9及表1可看出序列凸優(yōu)化算法與LGR偽譜法所得結(jié)果一致,均達(dá)到了最優(yōu)解,實(shí)現(xiàn)了燃料消耗最優(yōu)。下面就序列凸優(yōu)化算法的收斂性及計(jì)算效率進(jìn)行分析,并與LGR偽譜法進(jìn)行比較。

圖10為每次序列迭代求解消耗時(shí)間。由圖可知,平均消耗時(shí)間為0.2343 s,其中耗時(shí)最長的一次計(jì)算耗時(shí)0.2678 s,出現(xiàn)在第二步,這是由時(shí)變加速度剖面在第二步被突然地加入其中所致。

圖10 序列計(jì)算耗時(shí)Fig.10 Time consumption of sequential calculation

圖11為各次序列計(jì)算所得高度曲線?？梢钥闯?前兩次計(jì)算曲線相差較大,而后曲線很快地收斂,第3次與第8次所得曲線已幾乎重合。為進(jìn)一步衡量算法的收斂性,以得到收斂速度,計(jì)算前后兩次計(jì)算所得高度剖面()的差,將其平均值作為衡量指標(biāo),即：

圖11 每次計(jì)算所得高度曲線Fig.11 Curves of calculated altitude

(30)

式中:3,()為第次序列計(jì)算所得著陸器高度剖面;為著陸器飛行時(shí)間。類似地,還可以定義：

(31)

式中:c,()為第次序列計(jì)算所得著陸器推力指令。若式(30)與(31)趨于0,則說明軌跡及推力指令不再隨序列凸優(yōu)化發(fā)生變化,算法收斂。圖12與圖13分別為高度剖面以及推力指令的收斂過程,由于Δ與Δ變化范圍較大,縱軸采用對數(shù)坐標(biāo)。由圖可以看出,前9次序列計(jì)算軌跡和指令均迅速收斂。其中Δ由111.3 m收斂到6.9×10m, Δ由2151 N收斂到8.2×10N,而后在這附近小范圍波動(dòng)。最大波動(dòng)范圍Δ<3.3×10m, Δ< 2.5×10N。就導(dǎo)航系統(tǒng)精度及推力系統(tǒng)精度而言,軌跡規(guī)劃只需要收斂至1.0×10m以內(nèi),推力指令收斂到0.1 N以內(nèi)即可。由圖可知,計(jì)算至第6次即可滿足精度要求,此時(shí)序列凸優(yōu)化算法總用時(shí)1.446 s,而偽譜法用時(shí)27.45 s,計(jì)算效率提高94.73%。以上為對初始軌跡規(guī)劃所用時(shí)間的分析。此外,除在初始軌跡規(guī)劃外的每個(gè)制導(dǎo)周期重新進(jìn)行軌跡規(guī)劃時(shí),可利用上一制導(dǎo)周期的結(jié)果作為序列計(jì)算初始值,此時(shí)僅需要計(jì)算1至2次即可符合收斂性要求。對歷次制導(dǎo)周期的收斂次數(shù)以及收斂時(shí)間進(jìn)行統(tǒng)計(jì),可知約78.3%的制導(dǎo)周期僅需要計(jì)算一次,21.7%僅需要計(jì)算兩次即可滿足收斂性要求,平均耗時(shí)0.2851 s,最長耗時(shí)0.5801 s。

圖12 高度剖面收斂過程Fig.12 Altitude profile convergence process

圖13 推力指令收斂過程Fig.13 Thrust command convergence process

5 精度分析

為分析本方法精度,考慮軌道轉(zhuǎn)移階段的制導(dǎo)及定軌精度、點(diǎn)火時(shí)間偏差等因素主要影響動(dòng)力下降段初始位置偏差。因此,本文確定一定的位置偏差作為初始條件進(jìn)行1000次蒙特卡洛打靶仿真驗(yàn)證。各工況初始偏差設(shè)置如表2所示。

表2 著陸器初始位置偏差Table 2 Lander initial position deviation

初始位置在徑向上的分布滿足：

(32)

各工況下著陸器初始的位置誤差分布如圖14所示。

圖14 著陸器初始位置誤差分布Fig.14 Initial position error distribution of the lander

仿真結(jié)果如圖15所示,關(guān)于著陸器精度具體的統(tǒng)計(jì)結(jié)果見表3。其中高度、垂直速度和水平速度這三個(gè)指標(biāo)決定了著陸器是否能夠?qū)崿F(xiàn)安全的月面軟著陸；落點(diǎn)分布則展示了著陸器著陸精度。圖15(a)為各個(gè)工況下最終高度散布情況,由圖可以看出,最終的高度誤差分布基本符合正態(tài)分布；且初始誤差越大,最終的高度誤差也越大,但即使對于工況3這一最惡劣的情況,最終的高度誤差也不超過±3.147 m；圖15(b)與(c)分別表示各工況下的垂直速度分布和水平速度分布,可以看出垂直速度誤差和水平速度誤差分別在0.036 m/s與0.003 m/s以內(nèi)。以上結(jié)果表明,即使當(dāng)初始誤差為±2500 m時(shí),序列凸優(yōu)化制導(dǎo)算法仍然能夠使著陸器在動(dòng)力下降段完成時(shí)達(dá)到預(yù)定的高度和速度,從而順利銜接姿態(tài)調(diào)整段,實(shí)現(xiàn)安全著陸。圖15(d)為最終落點(diǎn)示意圖,由圖可知在三種工況下著陸器動(dòng)力下降段落點(diǎn)誤差分別為0.8670 m,2.5336 m,3.6655 m。上述仿真結(jié)果表明,當(dāng)初始位置誤差在±2500 m的范圍內(nèi)時(shí),著陸器能夠達(dá)到預(yù)定的位置,實(shí)現(xiàn)高精度月面動(dòng)力下降。

表3 著陸器末狀態(tài)精度統(tǒng)計(jì)Table 3 Precision statistics of the lander in its final state

圖15 最終誤差分布及落點(diǎn)情況Fig.15 Distribution of the final errors and landing points

6 結(jié) 論

本文針對傳統(tǒng)優(yōu)化算法應(yīng)用于月面著陸制導(dǎo)時(shí),時(shí)變慣性加速度和重力加速度難以估計(jì)補(bǔ)償問題,提出了一種基于改進(jìn)序列凸優(yōu)化的在線制導(dǎo)方法。在考慮月面曲率及月球自轉(zhuǎn)的前提下建立了月面著陸器動(dòng)力學(xué)模型；首先對該模型和約束條件進(jìn)行凸化,得到一個(gè)SOCP問題；并通過逐次進(jìn)行序列凸優(yōu)化對動(dòng)力學(xué)模型中難以精確估計(jì)的時(shí)變加速度項(xiàng)進(jìn)行迭代優(yōu)化,直至收斂。在此過程中,時(shí)變加速度被同倫地添加到該問題中進(jìn)行序列迭代求解,使得存在時(shí)變加速度擾動(dòng)時(shí),著陸器仍能在節(jié)約燃料的同時(shí)實(shí)現(xiàn)高精度著陸。仿真結(jié)果表明,相比于經(jīng)典的顯式制導(dǎo)律,序列凸優(yōu)化算法所消耗的燃料更少,相比于經(jīng)典的偽譜法,序列凸優(yōu)化算法的計(jì)算效率更高。針對不同初始位置偏差的蒙特卡洛打靶結(jié)果表明,著陸器初始位置誤差即使高達(dá)±2500 m,著陸器的末狀態(tài)速度精度依然在±0.037 m/s以內(nèi),位置精度在±3.6655 m內(nèi)。綜上所述,基于序列凸優(yōu)化的在線制導(dǎo)算法能夠高精度地完成動(dòng)力下降制導(dǎo),使著陸器能順利地銜接姿態(tài)調(diào)整段,實(shí)現(xiàn)高精度安全著陸。