參數(shù)自適應(yīng)凸優(yōu)化下的月面著陸最優(yōu)軌跡規(guī)劃

穆榮軍，鄧雁鵬，吳 鵬

（哈爾濱工業(yè)大學(xué)航天學(xué)院，哈爾濱 150001）

0 引言

為順利實(shí)施月球著陸探測(cè)器自主避障、保證精確安全著陸，通常將月球著陸過(guò)程分為3個(gè)階段，即動(dòng)力下降段、接近段與著陸段，如圖1所示［1］。這些階段通常還可能根據(jù)任務(wù)不同進(jìn)一步細(xì)分，例如中國(guó)的嫦娥系列探測(cè)器［1-4］。以嫦娥三號(hào)為例，月球著陸器在動(dòng)力下降段前處于遠(yuǎn)地點(diǎn)100 km、近月點(diǎn)15 km的轉(zhuǎn)移軌道上。自近月點(diǎn)附近開(kāi)始進(jìn)入動(dòng)力下降段。動(dòng)力下降段的主要目標(biāo)為盡可能降低著陸器的高度和速度。在動(dòng)力下降段，著陸器高度從15 km降低到2.4 km，速度將從1.7 km∕s 降低到0 左右［5］。該階段通常會(huì)消耗掉約80%的燃料，因此該段制導(dǎo)算法要求在著陸器具有初始狀態(tài)誤差和模型不確定性的情況下盡量降低燃料消耗，同時(shí)具有足夠的位置和速度精度以達(dá)到后續(xù)接近段的入口狀態(tài)要求。

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

嫦娥系列探測(cè)器為在動(dòng)力下降階段盡可能地節(jié)省燃料，采用了具有燃料次優(yōu)特性的線性正切制導(dǎo)律。但線性正切制導(dǎo)律不能控制著陸器的航程，因此并不能用于需要精確定點(diǎn)著陸的應(yīng)用場(chǎng)景。例如，嫦娥三號(hào)的預(yù)定著陸區(qū)是位于月球虹灣的一個(gè)長(zhǎng)356 km、寬91 km 的帶狀區(qū)域［6］。這顯然不能滿足未來(lái)月球探測(cè)的任務(wù)需求。

中國(guó)計(jì)劃在未來(lái)10年內(nèi)開(kāi)展載人登月活動(dòng)，同時(shí)建立有人駐守的月球基地［7］，這就需要著陸器能夠在基地附近定點(diǎn)精確著陸；此外，隨著對(duì)月球探索的進(jìn)一步深入，月球著陸器還需要能夠在具有復(fù)雜地形的高科學(xué)價(jià)值目標(biāo)附近精確著陸［8-9］。以上均需要下一代月球探測(cè)器具有精確定點(diǎn)著陸的能力。而動(dòng)力下降段作為月面著陸過(guò)程中飛行時(shí)間最長(zhǎng)，燃料消耗最多，航程最遠(yuǎn)的主要階段，其軌跡規(guī)劃的質(zhì)量與制導(dǎo)的精度將直接決定最終著陸精度與燃料消耗，因此許多機(jī)構(gòu)和個(gè)人對(duì)此進(jìn)行了深入的研究［10-11］。

現(xiàn)有的月面著陸制導(dǎo)算法分為閉環(huán)顯式制導(dǎo)算法和標(biāo)稱軌跡跟蹤制導(dǎo)算法兩類［12］，其中顯式制導(dǎo)律以A2PDG［13］、基于分?jǐn)?shù)多項(xiàng)式的動(dòng)力下降制導(dǎo)［14］以及ZEM∕ZEV 制導(dǎo)律［15-17］為代表，在具有良好的自主性和實(shí)時(shí)性的同時(shí)也具有一定的抗干擾能力，能夠保證末端制導(dǎo)精度。

相較于閉環(huán)顯式制導(dǎo)算法而言，基于軌跡優(yōu)化算法的標(biāo)稱軌跡制導(dǎo)方法首先需要利用最優(yōu)化原理，按照一定的目標(biāo)及約束條件設(shè)計(jì)標(biāo)稱著陸軌跡。該方法通常消耗燃料更少，具有更強(qiáng)的魯棒性。其中，凸優(yōu)化算法由于其具有多項(xiàng)式級(jí)的計(jì)算復(fù)雜度與理論上的全局最優(yōu)性，而得到了廣泛的研究與應(yīng)用［18-21］。

然而，常見(jiàn)的飛行器軌跡規(guī)劃動(dòng)力學(xué)模型約束和控制量約束都是非凸的，因此需要先將該問(wèn)題通過(guò)凸化轉(zhuǎn)變?yōu)橥箚?wèn)題再進(jìn)行求解。常見(jiàn)的凸化方法主要有序列凸化和無(wú)損凸化兩種［22］。其中無(wú)損凸化技術(shù)通過(guò)引入松弛變量的方法將非凸約束松弛為凸的形式，再證明松弛問(wèn)題的最優(yōu)解也是原始問(wèn)題的最優(yōu)解，實(shí)現(xiàn)原始非凸問(wèn)題向凸問(wèn)題的等價(jià)轉(zhuǎn)換。該方法相較于序列凸化，具有無(wú)需原始參考軌跡，同時(shí)不需要添加附加的信賴域約束，收斂特性良好的特點(diǎn)。Acikmese 等［21］針對(duì)火星著陸問(wèn)題提出了一種用于動(dòng)力下降制導(dǎo)問(wèn)題的凸規(guī)劃算法，該算法將火星著陸問(wèn)題通過(guò)無(wú)損凸化最終轉(zhuǎn)化為一個(gè)二階錐規(guī)劃（Second order cone programming，SOCP）問(wèn)題，進(jìn)而保證了燃料最優(yōu)。隨后這一方法得到廣泛的關(guān)注，并在不同方面得到進(jìn)一步地改進(jìn)。有的學(xué)者為提高該方法的收斂性，收斂速度，精度等，針對(duì)該方法本身做了改進(jìn)。

文獻(xiàn)［23］進(jìn)一步明確了凸優(yōu)化所需要的最佳飛行時(shí)間（time of flight，TOF）的確定方法；文獻(xiàn)［24］將凸優(yōu)化與偽譜法結(jié)合，通過(guò)合理的選擇節(jié)點(diǎn)集合，減少了節(jié)點(diǎn)個(gè)數(shù)，節(jié)約了CPU 時(shí)間，提高了運(yùn)算精度。文獻(xiàn)［25］進(jìn)一步提出將該方法運(yùn)用到了交會(huì)對(duì)接方面。

另一方面，很多學(xué)者針對(duì)該方法對(duì)于復(fù)雜障礙環(huán)境，復(fù)雜動(dòng)力學(xué)環(huán)境的適應(yīng)能力做出了相應(yīng)改進(jìn)。例如文獻(xiàn)［4］對(duì)月面著陸器動(dòng)力下降制導(dǎo)過(guò)程中時(shí)變的慣性加速度和重力加速度難以估計(jì)與補(bǔ)償?shù)葐?wèn)題提出一種基于序列凸優(yōu)化的在線制導(dǎo)算法，將時(shí)變加速度剖面同倫地添加到原非線性模型中序列迭代求解，使得該算法不依賴初始參考軌跡，增強(qiáng)了算法的自主性和多種應(yīng)用場(chǎng)景適應(yīng)能力。文獻(xiàn)［23］通過(guò)同倫方法與不動(dòng)點(diǎn)迭代思想處理時(shí)變的氣動(dòng)約束，在此基礎(chǔ)上設(shè)計(jì)了基于凸優(yōu)化的火箭垂直在線軌跡規(guī)劃方法。文獻(xiàn)［26］將凸優(yōu)化與曲率調(diào)整策略結(jié)合在了一起，在保證避障能力的同時(shí)降低了油耗。文獻(xiàn)［27］考慮了在具有不確定性情況下的最優(yōu)軌跡規(guī)劃。文獻(xiàn)［28］通過(guò)對(duì)該算法進(jìn)行改進(jìn)，將其應(yīng)用于具有不規(guī)則形狀與不規(guī)則引力場(chǎng)的小行星著陸。

在另一方面，不確定性廣泛存在于航天系統(tǒng)設(shè)計(jì)中，對(duì)系統(tǒng)響應(yīng)、最終控制精度等影響很大，在系統(tǒng)不確定性超過(guò)設(shè)計(jì)冗余時(shí)可能導(dǎo)致設(shè)計(jì)失敗。因此，有必要針對(duì)系統(tǒng)不確定性建模，以提升系統(tǒng)精度及魯棒性。

為解決上述問(wèn)題，許多學(xué)者提出了先進(jìn)控制理論［29］，如自適應(yīng)神經(jīng)網(wǎng)絡(luò)控制［29-32］、非參數(shù)自適應(yīng)控制［33］和滑模控制［34-35］等方法。此外，部分學(xué)者針對(duì)系統(tǒng)不確定性進(jìn)行在線辨識(shí)，同時(shí)通過(guò)在線辨識(shí)結(jié)果對(duì)模型進(jìn)行在線補(bǔ)償以降低系統(tǒng)的不確定性。文獻(xiàn)［36］針對(duì)小推力航天器軌道保持問(wèn)題，設(shè)計(jì)了一種基于改進(jìn)型條件積分滑模面和徑向基函數(shù)神經(jīng)網(wǎng)絡(luò)估計(jì)器，并通過(guò)實(shí)時(shí)補(bǔ)償實(shí)現(xiàn)高精度軌道保持。文獻(xiàn)［37］針對(duì)空間探測(cè)問(wèn)題，對(duì)不確定擾動(dòng)按從輕到重4 個(gè)等級(jí)進(jìn)行劃分，通過(guò)模糊神經(jīng)網(wǎng)絡(luò)進(jìn)行評(píng)估并采取響應(yīng)策略，從而提高任務(wù)的可靠性及靈活性。文獻(xiàn)［38］利用范數(shù)自適應(yīng)控制算法和自適應(yīng)估計(jì)算法有效估計(jì)充液航天器擾動(dòng)的未知上界以及液體晃動(dòng)的位移變量，實(shí)現(xiàn)了對(duì)充液航天器的自適應(yīng)魯棒控制。在月球探測(cè)方面，文獻(xiàn)［39］針對(duì)月球著陸過(guò)程中需要高精度懸停的需求，設(shè)計(jì)了以位置為量測(cè)量的觀測(cè)器，以實(shí)現(xiàn)對(duì)質(zhì)量和推力不確定性的快速估計(jì)和補(bǔ)償。文獻(xiàn)［40］將月球上升段制導(dǎo)參數(shù)的在線估計(jì)問(wèn)題轉(zhuǎn)換為離散線性系統(tǒng)的最優(yōu)估計(jì)問(wèn)題，解決了動(dòng)力顯示制導(dǎo)的參數(shù)自適應(yīng)問(wèn)題。

本文以嫦娥五號(hào)探測(cè)器所攜帶的著陸器與上升器組合體為例，提出了一種利用現(xiàn)有的發(fā)動(dòng)機(jī)實(shí)現(xiàn)精確月球著陸算法。在標(biāo)稱狀態(tài)下，通過(guò)合適的轉(zhuǎn)移軌道和推力剖面規(guī)劃即可以實(shí)現(xiàn)精確著陸。然而，考慮到在下降著陸過(guò)程中，由于發(fā)動(dòng)機(jī)工作條件、加工工藝、氧化劑和燃料的狀態(tài)存在一定的不確定性，因此，著陸器的比沖存在一定的不確定性；同時(shí)，由于工質(zhì)揮發(fā)等因素，著陸器的質(zhì)量也具有一定的不確定性［26，39］，這些不確定性會(huì)影響最優(yōu)軌跡規(guī)劃效果，使得所規(guī)劃軌跡喪失最優(yōu)性，增加燃料消耗。為消除系統(tǒng)的不確定性，盡可能地節(jié)約燃料，提升著陸精度，本文提出了一種基于參數(shù)自適應(yīng)凸優(yōu)化的軌跡規(guī)劃算法。該算法以凸優(yōu)化軌跡規(guī)劃為基礎(chǔ)，考慮到飛行過(guò)程中比沖和實(shí)時(shí)質(zhì)量的不確定性，為實(shí)現(xiàn)參數(shù)的快速高精度收斂，通過(guò)設(shè)計(jì)最優(yōu)觀測(cè)器，利用加速計(jì)量測(cè)數(shù)據(jù)對(duì)以上參數(shù)進(jìn)行在線實(shí)時(shí)估計(jì)。在估計(jì)結(jié)果達(dá)到收斂后再利用凸優(yōu)化重新進(jìn)行軌跡規(guī)劃，以得到符合實(shí)際參數(shù)的最優(yōu)軌跡。

本文編排如下：在第1 節(jié)首先介紹了月面著陸動(dòng)力學(xué)模型，引入了基于凸優(yōu)化的最優(yōu)軌跡規(guī)劃方法；在第2節(jié)本文首先考慮了不確定性模型，并基于該模型設(shè)計(jì)了最優(yōu)狀態(tài)觀測(cè)器，利用該觀測(cè)器實(shí)現(xiàn)對(duì)著陸器參數(shù)的在線實(shí)時(shí)估計(jì)；最后通過(guò)仿真分析驗(yàn)證了該觀測(cè)器的收斂速度及最終收斂精度，并通過(guò)蒙特卡洛打靶試驗(yàn)證明了該方法僅僅需要發(fā)動(dòng)機(jī)具有離散的推力域，同時(shí)能夠提升制導(dǎo)精度。

1 基于凸優(yōu)化算法的著陸動(dòng)力學(xué)模型

1.1 總體技術(shù)方案

為實(shí)現(xiàn)月面精確著陸，著陸器所攜帶的發(fā)動(dòng)機(jī)需要具有一定的推力調(diào)節(jié)能力。目前嫦娥五號(hào)所攜帶的發(fā)動(dòng)機(jī)為一臺(tái)變推力發(fā)動(dòng)機(jī)，其采用系統(tǒng)自鎖閥與發(fā)動(dòng)機(jī)氣動(dòng)閥控制。該發(fā)動(dòng)機(jī)可以輸出7 500 N的額定推力，也可以輸出1 500 N 至5 000 N 范圍內(nèi)的任意推力，但該發(fā)動(dòng)機(jī)的推力在5 000 N和7 500 N之間無(wú)法連續(xù)調(diào)節(jié)［41-43］，即該著陸器的推力域約束可以描述為

式中：Fc為主發(fā)動(dòng)機(jī)輸出推力，F(xiàn)1=5 000 N，F(xiàn)2=7 500 N，F(xiàn)3=1 500 N。

在此基礎(chǔ)上，本文進(jìn)一步探討一種攜帶僅具有離散推力域發(fā)動(dòng)機(jī)的探測(cè)器在月球表面著陸的問(wèn)題。這種發(fā)動(dòng)機(jī)的推力沒(méi)有任何連續(xù)調(diào)節(jié)能力，僅能按照F1或F2進(jìn)行輸出，即：

而在F1與F2之間，推力并不連續(xù)可調(diào)。

根據(jù)最優(yōu)化理論，采用凸優(yōu)化算法進(jìn)行軌跡規(guī)劃能夠保證在盡量節(jié)省燃料的前提下從給定的初始狀態(tài)轉(zhuǎn)移到預(yù)定的終端狀態(tài)，同時(shí)燃料消耗具有最優(yōu)性，能夠滿足動(dòng)力下降段節(jié)約燃料的要求。同時(shí)，由最優(yōu)化理論所產(chǎn)生的推力指令遵循最大-最小-最大模式［21］，即先輸出最大推力，然后在某時(shí)刻切換至小推力輸出，最后再切換至最大推力輸出，如下式所示：

式中：t1，t2為兩個(gè)推力突變時(shí)刻；tf為動(dòng)力下降段飛行時(shí)間。

根據(jù)式（3），最優(yōu)推力曲線呈現(xiàn)分段恒定的形式，可以斷定，式（2）所形成的推力大小在一定區(qū)間內(nèi)不可用的約束本質(zhì)上并不影響最優(yōu)解。因此本文采用凸優(yōu)化算法進(jìn)行最優(yōu)軌跡規(guī)劃，將式（2）的離散推力限制首先松弛為連續(xù)的推力區(qū)間：

在此基礎(chǔ)上，結(jié)合著陸器初末狀態(tài)約束、動(dòng)力學(xué)約束等進(jìn)一步將該問(wèn)題進(jìn)行凸化和離散化，最后轉(zhuǎn)化為SOCP 問(wèn)題后進(jìn)行求解。在解出最優(yōu)推力剖面后，通過(guò)主發(fā)動(dòng)機(jī)在最大與最小推力之間的切換，實(shí)現(xiàn)對(duì)該最優(yōu)推力剖面的跟蹤，進(jìn)而實(shí)現(xiàn)月面高精度定點(diǎn)著陸。

此外，考慮到在實(shí)際飛行過(guò)程中，著陸器的比沖、質(zhì)量可能與標(biāo)稱值不一致，從而對(duì)所規(guī)劃的軌跡的最優(yōu)性造成影響，使得燃料消耗增加，著陸精度降低。為解決上述問(wèn)題，本文根據(jù)主減速段標(biāo)稱軌跡來(lái)源進(jìn)一步將主減速段劃分為使用標(biāo)稱參數(shù)凸優(yōu)化以及在線參數(shù)自適應(yīng)凸優(yōu)化的兩段，如圖2所示。著陸器首先將根據(jù)標(biāo)稱參數(shù)離線產(chǎn)生一條標(biāo)稱軌跡，并在主發(fā)動(dòng)機(jī)點(diǎn)火后沿著該利用標(biāo)稱參數(shù)進(jìn)行最優(yōu)軌跡規(guī)劃形成的軌跡飛行；與此同時(shí)，在飛行過(guò)程中，著陸器將根據(jù)發(fā)動(dòng)機(jī)推力指令，加速度計(jì)的測(cè)量結(jié)果，利用最優(yōu)觀測(cè)器實(shí)時(shí)估計(jì)著陸器質(zhì)量，比沖等相關(guān)參數(shù)，在參數(shù)達(dá)到收斂后再重新進(jìn)行凸優(yōu)化，生成新的最優(yōu)軌跡，并沿著該軌跡飛行。此時(shí)由于最優(yōu)軌跡是根據(jù)實(shí)際參數(shù)生成的，因此能夠達(dá)到更高的著陸精度，同時(shí)消耗更少的燃料。

本文考慮到在月心慣性坐標(biāo)系下建模不需要考慮慣性力項(xiàng)，因此具有模型簡(jiǎn)單的優(yōu)勢(shì)，故在月心慣性坐標(biāo)系下建立著陸動(dòng)力學(xué)模型。其中月心慣性如圖3 所示，坐標(biāo)系的原點(diǎn)O位于月心，OzI與月球平赤道面垂直，指向月球北極；OxI在月球平赤道平面上指向J2000 時(shí)刻的春分點(diǎn)方向；OyI與另外兩軸構(gòu)成右手直角坐標(biāo)系。

圖3 月球慣性坐標(biāo)系Fig.3 Lunar inertial coordinate system

圖4 月球固聯(lián)坐標(biāo)系Fig.4 Lunar fixed coordinate system

此外，考慮到月面固連坐標(biāo)系較為直觀，物理意義明確，本文將在該坐標(biāo)系下設(shè)置仿真初始條件及展示仿真結(jié)果。該坐標(biāo)系的原點(diǎn)OL為著陸器的目標(biāo)著陸位置，OLxL軸的指向平行于著陸器標(biāo)稱初始位置到目標(biāo)著陸點(diǎn)的矢量在月面的投影，OLzL豎直向上，OLyL軸與以上兩軸組成右手坐標(biāo)系。

從月球固聯(lián)坐標(biāo)系到月球慣性坐標(biāo)系下的轉(zhuǎn)換關(guān)系為：

式中：β為月球固聯(lián)坐標(biāo)系與當(dāng)?shù)卣龞|方向的夾角，以東偏北為正。

1.2 著陸器動(dòng)力學(xué)模型

在月球慣性坐標(biāo)系下，著陸器的動(dòng)力學(xué)模型為

式中：r(t)，v(t)分別為著陸器t時(shí)刻的位置和速度；Fc(t)為著陸器推力；m(t)為著陸器質(zhì)量；Isp為推力器比沖；ge為標(biāo)稱地球重力加速度；g(t)為月球重力加速度；α=1∕(Ispge)為質(zhì)量流量系數(shù)。

該問(wèn)題經(jīng)由凸化和離散化最終可以化為以下SOCP問(wèn)題［4，19-21］：

通過(guò)解式（9）這一SOCP 問(wèn)題，即可得到最優(yōu)控制量以及相應(yīng)的最優(yōu)軌跡。可以看到式（9）中含有質(zhì)量流量系數(shù)α，根據(jù)定義，取決于著陸器比沖；同時(shí)其狀態(tài)量和約束均與初始質(zhì)量mwet有關(guān)，因此有必要對(duì)這些量進(jìn)行估計(jì)，降低系統(tǒng)的不確定性，提升優(yōu)化效果。

2 參數(shù)自適應(yīng)算法

在凸優(yōu)化過(guò)程中需要用到著陸器的比沖、初始質(zhì)量等參數(shù)。如果凸優(yōu)化使用的參數(shù)與實(shí)際不一致將影響軌跡的最優(yōu)性，從而進(jìn)一步影響燃料消耗與著陸精度。因此，本文設(shè)計(jì)了一個(gè)基于參數(shù)自適應(yīng)的凸優(yōu)化算法。通過(guò)構(gòu)造最優(yōu)觀測(cè)器在線實(shí)時(shí)估計(jì)比沖，質(zhì)量等參數(shù)。并在這些參數(shù)達(dá)到收斂后重新進(jìn)行凸優(yōu)化以保證軌跡的最優(yōu)性。

2.1 不確定性模型

考慮在著陸過(guò)程中，由于發(fā)動(dòng)機(jī)工作條件，加工工藝以及氧化劑及燃料狀態(tài)等存在一定的不確定性，因此著陸器的比沖存在一定的不確定性。同時(shí)由于著陸器攜帶的燃料揮發(fā)，物資重量變化等等影響，初始質(zhì)量也具有一定的不確定性，這樣飛行過(guò)程中的實(shí)時(shí)質(zhì)量亦具有不確定性。因此，本文考慮利用加速度計(jì)的在線測(cè)量值對(duì)比沖及質(zhì)量進(jìn)行估計(jì)。

假設(shè)比沖在飛行過(guò)程中保持不變，則質(zhì)量流量系數(shù)也為一恒定值。則著陸器質(zhì)量和質(zhì)量流量系數(shù)滿足以下微分方程：

以加速度計(jì)的測(cè)量結(jié)果作為量測(cè)值，則量測(cè)方程為

令θ=[mα]=[θ1θ2]，則式（12），（13）可記為

2.2 最優(yōu)觀測(cè)器設(shè)計(jì)

下面針對(duì)系統(tǒng)設(shè)計(jì)最優(yōu)觀測(cè)器。考慮觀測(cè)器模型為：

則H(t)可表示為

式中：P(t)為如下Riccati方程的解

其中，P0，Q，R為正定對(duì)稱矩陣。

將P，Q，R記為

可將式（16）與式（19）聯(lián)立展開(kāi)為

當(dāng)觀測(cè)器運(yùn)行時(shí)間t大于T0，且協(xié)方差矩陣P滿足ρ(P)≤ε時(shí)，認(rèn)為已經(jīng)收斂到θ。其中，ρ(P)表示矩陣P的譜半徑，ε為收斂后P的譜半徑上限。在實(shí)現(xiàn)收斂后，即可按照如式（9）所示的模型重新進(jìn)行凸優(yōu)化，獲得新的參考軌跡。

綜上所述，參數(shù)自適應(yīng)凸優(yōu)化算法流程如圖5所示。著陸器首先依據(jù)離線凸優(yōu)所產(chǎn)生的標(biāo)稱軌跡進(jìn)行標(biāo)稱軌跡制導(dǎo)，同時(shí)通過(guò)制導(dǎo)系統(tǒng)輸出的實(shí)時(shí)推力指令以及加速度計(jì)測(cè)量數(shù)據(jù)對(duì)實(shí)時(shí)質(zhì)量和比沖進(jìn)行在線參數(shù)估計(jì)。通過(guò)協(xié)方差矩陣P的譜半徑ρ(P)判斷所估計(jì)參數(shù)的收斂性。在參數(shù)估計(jì)達(dá)到收斂后，根據(jù)系統(tǒng)當(dāng)前狀態(tài)及系統(tǒng)實(shí)際參數(shù)，通過(guò)解有限維SOCP 問(wèn)題在線生成新的標(biāo)稱軌跡。在此后著陸器依據(jù)該軌跡進(jìn)行標(biāo)稱軌跡制導(dǎo)直到達(dá)到目標(biāo)狀態(tài)。

圖5 參數(shù)自適應(yīng)凸優(yōu)化算法流程Fig.5 Parameter adaptive convex optimization algorithm

3 仿真校驗(yàn)

為說(shuō)明本文所設(shè)計(jì)算法的有效性，本節(jié)將首先就最優(yōu)觀測(cè)器的收斂性及參數(shù)估計(jì)精度進(jìn)行仿真校驗(yàn)，其次將對(duì)比分析經(jīng)典的凸優(yōu)化算法與本文所提出的參數(shù)自適應(yīng)凸優(yōu)化算法的燃料消耗及著陸精度。

3.1 參數(shù)估計(jì)

本節(jié)將利用第2 節(jié)設(shè)計(jì)的最優(yōu)估計(jì)器對(duì)比沖、質(zhì)量的不確定性進(jìn)行估計(jì)，以驗(yàn)證本文所設(shè)計(jì)估計(jì)器的收斂速度及收斂精度。參考文獻(xiàn)［1-4，33］，仿真參數(shù)設(shè)置如表1所示。

表1 著陸器仿真參數(shù)設(shè)置Table 1 Lander simulation parameter setting

表2 參數(shù)估計(jì)誤差Table 2 Parameter estimation error

最優(yōu)觀測(cè)器對(duì)質(zhì)量和比沖的估計(jì)結(jié)果如圖6～圖7 所示。其中圖6 為著陸器質(zhì)量估計(jì)曲線，由于量測(cè)值與質(zhì)量直接相關(guān)，因此對(duì)質(zhì)量的估計(jì)收斂很快，在0.5 s 內(nèi)即收斂到了真實(shí)值，然后跟蹤真實(shí)值隨著飛行過(guò)程燃料的消耗而逐漸減少。相比之下，由于比沖（或質(zhì)量流量系數(shù)）與量測(cè)值的導(dǎo)數(shù)相關(guān)，屬于間接量測(cè)，所以比沖的收斂速度要慢于質(zhì)量，如圖7 所示。比沖大約用了10 s 才實(shí)現(xiàn)收斂。在本算例下，最終質(zhì)量估計(jì)誤差的均值為-0.005 8 kg，標(biāo)準(zhǔn)差為0.224 3 kg；最終比沖估計(jì)誤差均值為0.028 6 s，標(biāo)準(zhǔn)差為0.578 s。圖8 顯示了協(xié)方差矩陣P的譜半徑ρ(P)隨時(shí)間變化情況，由于P中各元素量綱不一致，圖中縱軸采用SI 單位制。可以看到隨著待估計(jì)參數(shù)的逐漸收斂，P的譜半徑在震蕩中也逐漸減小，收斂時(shí)間不超過(guò)10 s，最終收斂值小于150。因此，在本算例給定的條件下，可取T0為10 s，ε取150作為收斂性判據(jù)。

圖6 質(zhì)量估計(jì)曲線Fig.6 Mass estimation curve

圖7 比沖估計(jì)曲線Fig.7 Specific impulse estimation curve

圖8 譜半徑曲線Fig.8 Spectral radius curve

此外，比沖估計(jì)曲線在460 s后出現(xiàn)了±1.5 s 的波動(dòng)，對(duì)應(yīng)的譜半徑曲線也出現(xiàn)了波動(dòng)。該波動(dòng)是由于推力的突然變化引起的，對(duì)應(yīng)的推力曲線如圖9 所示。由于推力Fc的突然變化，根據(jù)式（17），Jacobian 矩陣AJ與CJ也會(huì)發(fā)生突變，進(jìn)而影響協(xié)方差矩陣P與反饋增益H(t)，并最終引起了估計(jì)值的突變。如果推力突變點(diǎn)出現(xiàn)在被估計(jì)參數(shù)收斂后，則對(duì)于著陸下降沒(méi)有影響；若出現(xiàn)在收斂前，則圖5中收斂性判據(jù)部分需要對(duì)觀測(cè)器運(yùn)行時(shí)間t重新計(jì)時(shí)，使得該突變的影響得以被排除，進(jìn)而使得觀測(cè)器有足夠的時(shí)間實(shí)現(xiàn)收斂。一般而言，凸優(yōu)化算法產(chǎn)生的推力曲線會(huì)產(chǎn)生兩個(gè)推力突變點(diǎn)，從而引起估計(jì)值的突變。同時(shí)由圖9 可以看出，采用凸優(yōu)化方法后形成的推力曲線符合如式（3）描述，呈現(xiàn)分段恒定的特征，進(jìn)而符合式（2）的推力域約束要求。

圖9 推力曲線Fig.9 Thrust curve

3.2 綜合仿真校驗(yàn)

為說(shuō)明基于參數(shù)自適應(yīng)凸優(yōu)化的最優(yōu)軌跡規(guī)劃方法的有效性，本節(jié)進(jìn)行了1 000 次蒙特卡洛打靶試驗(yàn)。分別采用經(jīng)典凸優(yōu)化方法，ZEM∕ZEV 制導(dǎo)方法，自適應(yīng)凸優(yōu)化方法進(jìn)行著陸。仿真參數(shù)如表3所示，表3中未列舉出的參數(shù)與表1保持一致。

表3 蒙特卡洛打靶仿真參數(shù)設(shè)置Table 3 Monte Carlo shooting parameter Setting

蒙特卡洛打靶結(jié)果如圖10～圖15所示，其中高度誤差即為Z軸方向位置誤差，垂直速度誤差即為Z方向速度誤差，相對(duì)頻率為各結(jié)果對(duì)于總仿真次數(shù)的占比，對(duì)應(yīng)的統(tǒng)計(jì)如表4 所示。由圖10～圖15可以看出，采用經(jīng)典的凸優(yōu)化算法的最終著陸點(diǎn)三軸位置誤差分別為-33.7 m，1.4 m，16.0 m，速度誤差分別為0.600 m∕s，0.016 m∕s，0.370 m∕s；采用ZEM∕ZEV 制導(dǎo)方案三軸位置誤差分別為3.98 m，0.29 m，1.69 m，速度誤差均值分別為0.050 m∕s，0.012 m∕s，0.030 m∕s；采用自適應(yīng)凸優(yōu)化的最終著陸點(diǎn)三軸位置誤差分別為4.00 m，0.32 m，1.73 m，速度誤差均值分別為0.060 m∕s，0.012 m∕s，0.032 m∕s。可以看到，相較于經(jīng)典凸優(yōu)化方法，自適應(yīng)凸優(yōu)化方法在很大程度上提高了速度精度與位置精度，其最終著陸精度與ZEM∕ZEV制導(dǎo)相當(dāng)。

表4 最終著陸精度統(tǒng)計(jì)Table 4 Final landing accuracy statistics

圖10 X方向位置誤差統(tǒng)計(jì)結(jié)果Fig.10 X-axis position error statistical result

圖11 Y方向位置誤差統(tǒng)計(jì)結(jié)果Fig.11 Y-axis position error statistical result

圖12 高度誤差統(tǒng)計(jì)結(jié)果Fig.12 Altitude error statistical result

圖14 Y方向速度誤差統(tǒng)計(jì)結(jié)果Fig.14 Y-axis speed error statistical result

圖15 垂直速度誤差統(tǒng)計(jì)結(jié)果Fig.15 Vertical speed error statistical result

值得注意的是，無(wú)論是否采用何種制導(dǎo)方案，Y方向上位置，速度的精度均高于其他兩軸。這是由于著陸過(guò)程中著陸器主要在航向和豎直方向運(yùn)動(dòng)，因此主要運(yùn)動(dòng)集中在縱平面內(nèi)，橫向運(yùn)動(dòng)較小，因此引起的該方向上的絕對(duì)誤差也相應(yīng)較小。

圖16展示了1 000 次蒙特卡洛打靶的燃料消耗情況的統(tǒng)計(jì)結(jié)果。運(yùn)用經(jīng)典凸優(yōu)化算法所消耗燃料的均值為1 428.9 kg，標(biāo)準(zhǔn)差為17.29 kg；采用本文所提出的參數(shù)自適應(yīng)凸優(yōu)化算法所消耗的均值為1 427.6 kg，標(biāo)準(zhǔn)差為16.38 kg，而采用ZEM∕ZEV 制導(dǎo)的燃料消耗均值為1 435.7 kg，標(biāo)準(zhǔn)差為16.50 kg。可以看出，這3 種算法所消耗燃料質(zhì)量的標(biāo)準(zhǔn)差相似，這是由于所消耗的燃料分布在極大程度上取決于著陸器實(shí)際的初始質(zhì)量與比沖，而這是由本仿真算例的初始條件給定的，因此其標(biāo)準(zhǔn)差幾乎一致。但就均值而言，采用參數(shù)自適應(yīng)凸優(yōu)化算法所消耗的燃料相較于經(jīng)典凸優(yōu)化算法減少了2.26 kg，相較ZEM∕ZEV 減少了8.1 kg。從仿真結(jié)果可以看出，參數(shù)自適應(yīng)算法相較于凸優(yōu)化算法在精度提高的同時(shí)消耗燃料更少，這是由于該算在進(jìn)行最優(yōu)軌跡規(guī)劃時(shí)采用的參數(shù)更接近于實(shí)際值，因此軌跡的最優(yōu)性更好，從而更省燃料。而ZEM∕ZEV 制導(dǎo)算法依賴反饋調(diào)節(jié)不斷接近制導(dǎo)目標(biāo)，并未考慮參數(shù)偏差，因此雖然制導(dǎo)精度很高，但是燃料消耗顯著多余參數(shù)自適應(yīng)算法。

圖16 燃料消耗統(tǒng)計(jì)結(jié)果Fig.16 Fuel consumption statistical result

圖17 展示了這3 種算法在打靶運(yùn)行第一次時(shí)的推力指令剖面，可以看到ZEM∕ZEV 這一定點(diǎn)著陸制導(dǎo)算法不但在某些時(shí)間段內(nèi)指令推力超出了發(fā)動(dòng)機(jī)的最大推力，同時(shí)其還要求發(fā)動(dòng)機(jī)在一定范圍內(nèi)連續(xù)可調(diào)。但無(wú)論是經(jīng)典的凸優(yōu)化算法還是自適應(yīng)凸優(yōu)化算法均符合推力分段恒定的特征，滿足發(fā)動(dòng)機(jī)工作要求。

圖17 典型推力指令剖面Fig.17 Typical thrust command profile

4 結(jié)論

本文提出了一種基于參數(shù)自適應(yīng)凸優(yōu)化的月面著陸最優(yōu)軌跡設(shè)計(jì)方法。本文首先在月球著陸動(dòng)力學(xué)模型的基礎(chǔ)上將月面著陸最優(yōu)軌跡設(shè)計(jì)問(wèn)題轉(zhuǎn)化為有限維SOCP 問(wèn)題。針對(duì)該SOCP 求解過(guò)程中比沖，初始質(zhì)量等參數(shù)具有一定不確定性的問(wèn)題，首先利用標(biāo)稱參數(shù)獲得最優(yōu)軌跡，并在沿著該軌跡進(jìn)行最優(yōu)軌跡跟蹤制導(dǎo)時(shí)，通過(guò)本文所設(shè)計(jì)的最優(yōu)觀測(cè)器對(duì)以上參數(shù)進(jìn)行實(shí)時(shí)估計(jì)。待參數(shù)估計(jì)滿足收斂性指標(biāo)要求后，再將參數(shù)的估計(jì)值帶入以上SOCP 問(wèn)題重新求解，以此獲得新的標(biāo)稱最優(yōu)軌跡。通過(guò)對(duì)參數(shù)的實(shí)時(shí)估計(jì)以及利用該參數(shù)的估計(jì)值對(duì)月面著陸最優(yōu)軌跡進(jìn)行實(shí)時(shí)在線設(shè)計(jì)，使得系統(tǒng)對(duì)于參數(shù)不確定性具有了一定的適應(yīng)能力。仿真結(jié)果表明，本文所設(shè)計(jì)的最優(yōu)觀測(cè)器能夠完成快速收斂，實(shí)現(xiàn)對(duì)實(shí)時(shí)質(zhì)量，比沖的在線觀測(cè)。其中，在本算例下，質(zhì)量收斂速度小于0.5 s，比沖收斂速度小于10 s。質(zhì)量的平均估計(jì)誤差小于0.23 kg，比沖估計(jì)誤差小于0.6 s，極大的減小了參數(shù)不確定度。1 000 次蒙特卡洛打靶仿真表明，參數(shù)自適應(yīng)凸優(yōu)化算法具有推力輸出剖面嚴(yán)格分段恒定，不需要發(fā)動(dòng)機(jī)推力連續(xù)可調(diào)，僅需要發(fā)動(dòng)機(jī)具有離散的推力域即可的優(yōu)勢(shì)；同時(shí)相較于其他算法，該算法能夠在具有參數(shù)不確定性的前提下顯著提升著陸精度。