行星動力下降多約束最優(yōu)反饋制導(dǎo)方法

張曉文，王天舒，王大軼，李 驥

(1. 清華大學(xué)航天航空學(xué)院，北京 100084；2. 北京控制工程研究所，北京 100190；3. 北京空間飛行器總體設(shè)計部，北京 100094)

0 引 言

未來行星著陸探測任務(wù)普遍提出定點著陸的工程目標(biāo),一般要求落點精度為100 m。實現(xiàn)此目標(biāo)主要依靠動力制導(dǎo)。以往月球和火星著陸器的動力制導(dǎo)方法研究和工程應(yīng)用多采用多項式制導(dǎo),例如阿波羅制導(dǎo)和E制導(dǎo)。阿波羅制導(dǎo)的優(yōu)點是終端三維位置、速度和加速度全部受控,缺點是固定終端時刻后未對推進(jìn)劑消耗作優(yōu)化。Lu[1-2]為改進(jìn)阿波羅制導(dǎo),先是參考隱式制導(dǎo)在原制導(dǎo)律中增加了一個控制增益作為可調(diào)參數(shù),之后應(yīng)用分?jǐn)?shù)多項式理論又增加了分?jǐn)?shù)指數(shù)參數(shù),并將新方法命名為分?jǐn)?shù)多項式動力下降制導(dǎo)。阿波羅制導(dǎo)中發(fā)動機(jī)推力加速度為時間的二次多項式形式。經(jīng)典的E制導(dǎo)及與其等同的基本型ZEM/ZEV制導(dǎo)[3]的推力加速度的制導(dǎo)指令為時間的一次多項式形式。E制導(dǎo)已被證明是常值引力場假設(shè)下的固定終端時刻的加速度模平方積分最優(yōu)制導(dǎo)[4],即能量最優(yōu)制導(dǎo),雖然不能完全等同于推進(jìn)劑消耗最優(yōu),但是已經(jīng)對推進(jìn)劑消耗做了顯著優(yōu)化。E制導(dǎo)的一項缺點是終端三維加速度不完全受控。E制導(dǎo)和阿波羅制導(dǎo)都未考慮推力范圍限制,因此實際應(yīng)用中指令推力可能觸發(fā)推力飽和。

E制導(dǎo)或基本型ZEM/ZEV制導(dǎo)盡管是最優(yōu)反饋制導(dǎo),但是面對實際工程應(yīng)用時還是存在一些不足,因此對其改進(jìn)的方法不斷被提出。這些改進(jìn)包括完善最優(yōu)制導(dǎo)時長確定方法,提高制導(dǎo)算法的魯棒性和碰撞規(guī)避能力等。Guo等[5]使用數(shù)值優(yōu)化程序GPOPS選擇中間航路點作為ZEM/ZEV制導(dǎo)的中間目標(biāo),實現(xiàn)了火星動力下降約束下滑角的碰撞規(guī)避。Zhou等[6]針對火星動力下降段制導(dǎo)規(guī)避碰撞地平問題,在ZEM/ZEV制導(dǎo)原性能指標(biāo)中添加了加權(quán)的高度項,然后推導(dǎo)出有地平碰撞規(guī)避能力的改進(jìn)制導(dǎo)律。Wibben等[7]為了提高ZEM/ZEV制導(dǎo)在火星動力下降段的干擾適應(yīng)性,加入最優(yōu)滑動制導(dǎo)并證明了新方法的全局穩(wěn)定性。Zhang等[8]利用平均加速度修正基本型ZEM/ZEV制導(dǎo)加速度的來解決月面著陸地平碰撞規(guī)避問題,還應(yīng)用模型預(yù)測靜態(tài)規(guī)劃(MPSP)方法迭代計算滿足終端狀態(tài)約束的解。Zhang等[9]對Zhou等[6]的方法繼續(xù)改進(jìn),將性能指標(biāo)中高度項的常值加權(quán)系數(shù)改為與高度平方相關(guān)的時變系數(shù),進(jìn)一步提升碰撞規(guī)避能力。Ahn等[10]研究了軌道攔截或交會機(jī)動問題的由最小能量和最小時間組成的性能指標(biāo),利用最優(yōu)控制理論和橫截條件,推導(dǎo)出確定ZEM/ZEV制導(dǎo)最優(yōu)制導(dǎo)時長的代數(shù)方程。Wang等[11]提出一種適用于火星精確著陸的兩段ZEM/ZEV反饋制導(dǎo)方法,通過虛擬終端速度和在線參數(shù)搜索來提高制導(dǎo)的障礙規(guī)避能力。高峰等[12]利用脈寬脈頻(PWPF)調(diào)制器將ZEM/ZEV制導(dǎo)的連續(xù)指令轉(zhuǎn)換為常值推力器的開關(guān)指令。

碰撞規(guī)避是行星著陸制導(dǎo)必須具備的能力。除了上面提到方法外,還有一種解決思路是控制著陸軌跡的曲率。Cui等[13]先是提出了一種基于多項式的曲率制導(dǎo)策略,通過調(diào)整曲率使著陸器沿凸形軌跡著陸,降低障礙物碰撞的風(fēng)險。之后進(jìn)一步改進(jìn),將曲率約束凸化為二階錐約束,并利用非線性規(guī)劃工具求解凸規(guī)劃問題來優(yōu)化推進(jìn)劑消耗[14]。

近年來數(shù)值優(yōu)化方法和人工智能技術(shù)被越來越多的應(yīng)用于行星著陸制導(dǎo)問題的理論研究中。林曉輝等[15]應(yīng)用凸優(yōu)化理論將月球動力下降段軌道優(yōu)化問題轉(zhuǎn)化為二階錐問題,并采用內(nèi)點法求解了最優(yōu)標(biāo)稱軌跡。鄧雁鵬等[16]應(yīng)用序列凸優(yōu)化方法解決月面動力下降過程中時變加速度剖面實時估計和補(bǔ)償問題。Furfaro等[17]將深度強(qiáng)化學(xué)習(xí)方法與ZEM/ZEV制導(dǎo)方法相結(jié)合,通過動態(tài)調(diào)節(jié)2個制導(dǎo)增益和制導(dǎo)時長來滿足行星著陸問題的路徑等約束并優(yōu)化推進(jìn)劑消耗。數(shù)值優(yōu)化制導(dǎo)計算量大,依賴高效穩(wěn)定的求解器。人工智能制導(dǎo)目前處于理論研究階段,工程應(yīng)用同樣面臨計算量大和算法穩(wěn)定性方面的挑戰(zhàn)。

總體而言,傳統(tǒng)最優(yōu)反饋制導(dǎo)方法有著計算量小、穩(wěn)定性好、節(jié)省推進(jìn)劑的優(yōu)勢,但是僅滿足初始和終端位置速度約束。實際工程應(yīng)用中還需要考慮推力范圍受限、連續(xù)變推力、碰撞規(guī)避和終端加速度等多方面約束。單臺變推力發(fā)動機(jī)推力的大小和方向在物理上只能連續(xù)變化,因此要求制導(dǎo)的指令推力連續(xù)變化。著陸終端希望著陸器為豎直姿態(tài),并且合加速度接近為零,這就對制導(dǎo)的終端加速度提出了約束?，F(xiàn)有最優(yōu)反饋制導(dǎo)方法研究欠缺對上述約束的綜合考慮。

本文考慮推力范圍受限、連續(xù)變推力、碰撞規(guī)避和終端加速度共4方面約束,對動力下降應(yīng)用的傳統(tǒng)最優(yōu)反饋制導(dǎo)方法進(jìn)行改進(jìn)和完善。內(nèi)容安排為,首先建立行星定點著陸動力下降制導(dǎo)問題模型,然后介紹和分析E制導(dǎo)和阿波羅制導(dǎo)方法,接著詳細(xì)闡述提出的多約束最優(yōu)反饋制導(dǎo)方法,最后通過數(shù)學(xué)仿真驗證所提方法有效性。

1 問題模型

1.1 制導(dǎo)坐標(biāo)系定義

制導(dǎo)坐標(biāo)系原點O定義在目標(biāo)著陸點,從原點指向天向定義為Z軸,垂直Z軸從原點指向著陸器初始位置方向為X軸,Y軸按照右手螺旋法則確定。在該坐標(biāo)系下,XOZ平面即為理想著陸軌跡面,或又稱為縱向平面,Y軸即為軌跡面負(fù)法線方向。X軸位移體現(xiàn)著陸器航程變化,因此可以將X向稱之為航向,Z軸位移體現(xiàn)著陸器高度變化,Y軸位移體現(xiàn)著陸器垂直軌跡面位置變化,因此可以將Y向稱之為橫向。

1.2 行星動力下降動力學(xué)方程

月球沒有大氣?；鹦谴髿庀”?在著陸器動力下降段氣動力顯著小于制導(dǎo)推力。因此以往分析行星動力下降制導(dǎo)問題時忽略大氣影響,僅考慮中心天體引力和發(fā)動機(jī)推力。制導(dǎo)坐標(biāo)系下的行星動力下降動力學(xué)方程如下:

(1)

(2)

a=F/m

(3)

(4)

1.3 定點著陸制導(dǎo)問題

定點著陸制導(dǎo)問題即控制著陸器以指定終端速度到達(dá)指定終端位置。研究定點著陸制導(dǎo)問題時一般將終端位置和速度取為零。定點著陸制導(dǎo)問題屬于軌跡優(yōu)化問題,軌跡優(yōu)化問題又屬于最優(yōu)控制問題。設(shè)初始時刻為0,下面依次列出定點著陸制導(dǎo)問題的初始條件、終端條件、性能指標(biāo)和推力約束表述。

初始條件:

r(0)=r0

(5)

v(0)=v0

(6)

m(0)=m0

(7)

終端條件:

r(tf)=rf

(8)

v(tf)=vf

(9)

性能指標(biāo):

(10)

推力約束:

Fmin≤F≤Fmax

(11)

2 多項式制導(dǎo)

多項式制導(dǎo)即制導(dǎo)指令加速度為時間的多項式形式。常見的多項式制導(dǎo)有E制導(dǎo)和阿波羅制導(dǎo)。下面對E制導(dǎo)和阿波羅制導(dǎo)展開原理介紹和性質(zhì)分析。在做性質(zhì)分析時假設(shè)引力加速度為常值。

2.1 E制導(dǎo)

在E制導(dǎo)方法中,制導(dǎo)推力加速度三維矢量取為時間的一次多項式,包含2個待確定系數(shù)矢量,即指令推力加速度為

a=c0+c1t,t∈[0,tf]

(12)

將式(12)代入動力學(xué)方程(2)中,然后連續(xù)積分兩次,并代入位置速度的初值和終值確定制導(dǎo)律系數(shù),結(jié)果如下:

(13)

(14)

將式(13)代入式(12)中后,得到E制導(dǎo)的初始時刻推力加速度計算公式:

(15)

定義零控位移偏差rzem和零控速度偏差vzev如下:

(16)

vzev?vf-v0-tfg

(17)

代入式(15)中,得到零控位移偏差rzem和零控速度偏差vzev形式的制導(dǎo)指令推力加速度計算公式,即基本型ZEM/ZEV制導(dǎo)律公式:

(18)

上式說明E制導(dǎo)與基本型ZEM/ZEV制導(dǎo)完全相同。將零控位移偏差rzem和零控速度偏差vzev前的常值控制增益6和-2改為可調(diào)控制增益kr和kv,則得到了一般型的ZEM/ZEV制導(dǎo)律公式[18]:

(19)

利用龐特里亞金極小值原理可以證明,E制導(dǎo)是常值引力場假設(shè)下的加速度模平方積分最優(yōu)制導(dǎo),即能量最優(yōu)制導(dǎo),使得如下性能指標(biāo)最小:

(20)

E制導(dǎo)性質(zhì)1.E制導(dǎo)的單個方向推力加速度為單調(diào)變化。

(21)

E制導(dǎo)性質(zhì)3.已知初始高度大于零,即rz0>0,在t∈[0,tf),滿足rz>0的條件為

(22)

2.2 阿波羅制導(dǎo)

在阿波羅制導(dǎo)方法中,制導(dǎo)加速度三維矢量取為時間的二次多項式,包含3個待確定系數(shù)矢量,即指令推力加速度為

a=c0+c1t+c2t2,t∈[0,tf]

(23)

將式(23)代入動力學(xué)方程(2)中,然后連續(xù)積分兩次,并代入位置、速度和加速度的初值和終值確定制導(dǎo)律系數(shù),結(jié)果如下:

(24)

(25)

(26)

將式(24)代入式(23)中后,得到E制導(dǎo)的初始時刻推力加速度計算公式:

(27)

由阿波羅制導(dǎo)的求解過程可知,求解常值制導(dǎo)參數(shù)時用到了終端位置、速度和加速度條件,因此與E制導(dǎo)相比,阿波羅制導(dǎo)除了可以滿足終端位置和速度約束外,還可以滿足終端加速度約束。

3 多約束最優(yōu)反饋制導(dǎo)

3.1 制導(dǎo)時長尋優(yōu)算法

多約束最優(yōu)反饋制導(dǎo)的基礎(chǔ)是E制導(dǎo),也即基本型ZEM/ZEV制導(dǎo)。E制導(dǎo)給出了終端時刻固定的指令加速度計算公式,沒有給出終端時刻的確定方法。如果制導(dǎo)加速度允許范圍不受約束,則E制導(dǎo)的終端時刻可由E制導(dǎo)性質(zhì)2,即式(21)確定,保證了制導(dǎo)加速度模平方的積分最小。工程應(yīng)用中制導(dǎo)加速度允許范圍一定受限。一般來講,制導(dǎo)時長越短越節(jié)省推進(jìn)劑[6],但是推力變化范圍也會隨之變大,最終指令推力超過了發(fā)動機(jī)實際允許范圍,即指令推力飽和?？紤]推力范圍約束后,實際輸出推力矢量滿足如下公式:

(28)

E制導(dǎo)有解析的軌跡預(yù)報公式,理論上只要終端指令推力不飽和,則規(guī)劃軌跡滿足終端狀態(tài)約束。因此,將終端指令推力不飽和作為制導(dǎo)算法基本原則。在此基礎(chǔ)上盡量縮短制導(dǎo)時長作為優(yōu)化推進(jìn)劑消耗的方向。最優(yōu)制導(dǎo)時長參數(shù)難以通過解析計算公式直接獲得,因此這里采用數(shù)值方法來尋找。優(yōu)化制導(dǎo)時長首先需要給出其一個合理的初值tf0。該初值應(yīng)盡量接近最優(yōu)值又要便于計算。這里提出將使得初始指令推力正好等于允許最大推力的制導(dǎo)時長作為優(yōu)化初值tf0。設(shè)rf=0和vf=0,由式(15)推導(dǎo)得到E制導(dǎo)的推力加速度模平方的初值為

(29)

由最大允許推力Fmax和初始質(zhì)量m0得到最大允許初始加速度為a0max=Fmax/m0。將最大允許初始加速度a0max代入式(29)得到制導(dǎo)時長tf的4次方程:

(30)

利用解析公式求解該方程,將其中最大實數(shù)根作為制導(dǎo)時長初值tf0。將制導(dǎo)時長初值tf0代入E制導(dǎo)公式,進(jìn)行首次數(shù)值積分預(yù)報,得到終端時刻標(biāo)稱推力指令。如果終端指令推力飽和,則需要增加制導(dǎo)時長。這里采用預(yù)測校正方法來計算增加量。記數(shù)值積分預(yù)報給出的指令推力在制導(dǎo)末端持續(xù)飽和的時長為tfs,校正系數(shù)記為ka,則給出制導(dǎo)時長第1個校正公式如下:

tf,i+1=tf,i+katfs

(31)

式中:下標(biāo)i為迭代計算拍數(shù)。將校正后制導(dǎo)時長再次代入制導(dǎo)方程,重新進(jìn)行數(shù)值積分預(yù)報。當(dāng)終端指令推力不再飽和時認(rèn)為制導(dǎo)時長優(yōu)化完成。

如果首次數(shù)值積分預(yù)報的終端指令推力未飽和,則說明制導(dǎo)時長還可以縮短。記制導(dǎo)時長縮短系數(shù)為ks,則制導(dǎo)時長縮短公式如下:

tf,i+1=kstf,i

(32)

應(yīng)保證制導(dǎo)時長縮短足夠充分,使得終端指令推力一定飽和。

接下來再次進(jìn)行數(shù)值積分預(yù)報。仍然采用預(yù)測校正方法來計算增加量。記校正系數(shù)為kb,則給出制導(dǎo)時長第2個校正公式如下:

tf,i+1=tf,i+kbtfs

(33)

校正系數(shù)ka和kb取值偏小的話,則單次修正量過小,需要的迭代計算次數(shù)多。校正系數(shù)取值偏大的話,則單次修正量過大,計算得到制導(dǎo)時長比其最優(yōu)值偏大較多。

定義推力俯仰角為推力方向和天向的夾角。為避免著陸器大幅度姿態(tài)機(jī)動帶來多方面不利影響需要約束推力俯仰角。通常不希望推力俯仰角超過90°,對應(yīng)的等效條件是垂向推力加速度非負(fù)。由垂向推力加速度非負(fù)得到的制導(dǎo)時長的不等式約束為:

(34)

顯然該式給出了一個制導(dǎo)時長的最小值約束。

3.2 碰撞規(guī)避算法

碰撞地面風(fēng)險規(guī)避要求飛行高度不低于地面。由E制導(dǎo)性質(zhì)3得到制導(dǎo)時長需滿足不等式(22),因此定義碰撞規(guī)避時長如下

(35)

式中:vzmin為足夠小的速度閾值,避免除零。

為規(guī)避碰撞地平風(fēng)險,計算垂向制導(dǎo)加速度時,制導(dǎo)時長計算公式如下:

(36)

依據(jù)上式得到新的制導(dǎo)加速度,進(jìn)行數(shù)值積分預(yù)報。如果推力沒有飽和,那么用上式方式來制導(dǎo)可以實現(xiàn)臨界碰撞規(guī)避。但是實際推力很可能初始就飽和,且上式取值臨界,此時碰撞規(guī)避就很可能已經(jīng)來不及。這種情況下需要進(jìn)一步優(yōu)化推力加速度的分配。由于碰撞規(guī)避的優(yōu)先級高,因此推力飽和時,優(yōu)先滿足垂向的推力加速度需求。給出碰撞規(guī)避控制的推力優(yōu)化分配律如下:

(37)

(38)

式中:m為著陸器當(dāng)前質(zhì)量;amax為當(dāng)前最大允許推力加速度;ax、ay和az為加速度優(yōu)化分配前指令加速度三軸分量,且az≥0;F為優(yōu)化分配后的指令推力矢量。

為了進(jìn)一步提高著陸安全性,增加碰撞規(guī)避的裕度,可以再增加下滑角約束,即著陸器與目標(biāo)著陸點連線,與地平的夾角需大于閾值。在滿足地平碰撞規(guī)避的基礎(chǔ)上,下滑角約束通過坐標(biāo)系旋轉(zhuǎn)的方法即可方便的得到滿足[11]。記下滑角閾值為α,則由制導(dǎo)坐標(biāo)系到下滑角坐標(biāo)系的旋轉(zhuǎn)陣為

(39)

將初始位置矢量r0和速度矢量v0,終端目標(biāo)位置矢量rf和速度矢量vf,引力加速度矢量g,都左乘坐標(biāo)系旋轉(zhuǎn)陣Cgl,得到下滑角坐標(biāo)系對應(yīng)的分量。將經(jīng)坐標(biāo)系旋轉(zhuǎn)變換后的變量作為新的初始約束、終端約束和引力加速度矢量代入制導(dǎo)計算公式便可得到滿足下滑角約束的解。滿足下滑角約束后既提高了碰撞規(guī)避能力又增加了著陸點可見性。

3.3 勻速轉(zhuǎn)動制導(dǎo)

著陸終端一般希望著陸器為豎直姿態(tài)并且合加速度近似為零,即推力加速度應(yīng)為負(fù)的引力加速度。這就對制導(dǎo)的終端加速度提出了約束。但是如前所述,E制導(dǎo)的終端三維加速度不完全受控,即不能滿足該約束。為了滿足終端加速度約束,在E制導(dǎo)結(jié)束之后,引入勻速轉(zhuǎn)動制導(dǎo),迅速調(diào)整推力加速度到目標(biāo)值。

勻速轉(zhuǎn)動制導(dǎo)的功能是依據(jù)最大角速度和最大推力變化率約束實現(xiàn)最短時間內(nèi)調(diào)整推力加速度矢量的方向和大小。勻速轉(zhuǎn)動制導(dǎo)的指令推力在由初始和終端加速度確定的平面內(nèi)勻速轉(zhuǎn)動。勻速轉(zhuǎn)動制導(dǎo)的輸入為初始制導(dǎo)加速度a0和終端目標(biāo)加速度af,輸出為指令推力F、制導(dǎo)時長tf、位置增量Δr和速度增量Δv。初始制導(dǎo)加速度a0和終端目標(biāo)加速度af的夾角為:

(40)

(41)

定義勻速轉(zhuǎn)動制導(dǎo)的平面坐標(biāo)系為:原點為推力作用點,X軸為初始制導(dǎo)加速度方向,X軸向終端目標(biāo)加速度方向旋轉(zhuǎn)90°得到Y(jié)軸。

勻速轉(zhuǎn)動制導(dǎo)的原理是首先計算平面坐標(biāo)系下推力加速度矢量,然后將其轉(zhuǎn)換為三維坐標(biāo)系下矢量。最終的勻速轉(zhuǎn)動制導(dǎo)的指令推力計算公式為:

(42)

F0=m0||a0||

(43)

Ff=mf||af||

(44)

θ=θft/tf

(45)

式中:M為將平面坐標(biāo)系矢量轉(zhuǎn)為三維坐標(biāo)系矢量的轉(zhuǎn)換陣;mf為著陸器在制導(dǎo)終端的質(zhì)量,可以用初始質(zhì)量m0近似代替。

對合加速度做一次和二次積分分別得到位置增量Δr和速度增量Δv的解析表達(dá)式:

(46)

(47)

圖1 組合制導(dǎo)參數(shù)預(yù)測校正規(guī)劃方法示意Fig.1 Diagram of predictor-corrector planning method for integrated guidance parameters

最優(yōu)反饋制導(dǎo)和勻速轉(zhuǎn)動制導(dǎo)交接處的位置rm和速度vm的校正公式如下:

rm=rf-Δr

(48)

vm=vf-Δv

(49)

4 仿真校驗

以火星探測任務(wù)的動力下降段為例,通過數(shù)學(xué)仿真方法來驗證和分析新提出的制導(dǎo)方法。仿真用例選取時參考了文獻(xiàn)[5,7,9,11]中的共同用例。著陸器的初始質(zhì)量取為1 905 kg,發(fā)動機(jī)比沖取為2 205 Ns/kg,制導(dǎo)坐標(biāo)系下的初始位置為[2 000, 0, 1 500]Tm,初始速度的多組工況在后面分別給出;目標(biāo)終端位置為[0, 0, 0]Tm,目標(biāo)終端速度為[0, 0, 0]Tm/s,目標(biāo)終端加速度為[0, 0, 0]Tm/s2。發(fā)動機(jī)推力最大值為13.258 kN,最小值為4.971 kN?；鹦且铀俣热閇0, 0, -3.711 4]Tm/s2。校正系數(shù)ka和kb的取值分別為0.7和0.3。制導(dǎo)時長縮短系數(shù)ks取為0.5。

4.1 有效性驗證

設(shè)計制導(dǎo)有效性驗證仿真工況對3.1節(jié)中方法做驗證,設(shè)定了4組初始速度取值,詳見表1。

表1 有效性驗證仿真工況的初始速度Table 1 Initial velocities of effectiveness verification simulation cases

有效性驗證仿真工況的軌跡見圖2。從圖中可知,4條軌跡的初始位置相同,由于初始速度不同導(dǎo)致中途軌跡不同,不過最終都匯集到了坐標(biāo)系原點,即目標(biāo)終端位置。比較工況2到3的制導(dǎo)飛行軌跡,當(dāng)航向初始速度為負(fù)值時,制導(dǎo)飛行軌跡偏左上,航向初始速度為正值時,制導(dǎo)飛行軌跡偏右下。比較工況3和4的制導(dǎo)飛行軌跡,當(dāng)垂向和航向初始速度相同,有橫向初始速度時,制導(dǎo)飛行軌跡將偏右下。有效性驗證仿真工況4的橫向位置速度曲線見圖3。從中可知,橫向位置偏差先增大然后逐漸減小到零,橫向速度偏差先由正值變?yōu)樨?fù)值然后逐漸減小到零。有效性驗證仿真工況的推力比見圖4。從中可知,初始航向和橫向速度越大則制導(dǎo)時長越長,初始推力飽和時長也越長。工況1、3和4的推力由飽和段、減小段和增加段組成,即推力大致按照大小大的趨勢變化,近似于最優(yōu)控制理論給出的變推力制導(dǎo)的最大最小最大推力剖面[19]。工況2雖然沒有推力飽和段,但是推力符合大、小、大的變化趨勢。有效性驗證仿真工況的推力俯仰角見圖5。從圖中可發(fā)現(xiàn)推力俯仰角的2個變化規(guī)律。第1個規(guī)律是推力俯仰角先大后小然后再增大,原因是著陸器為了實現(xiàn)水平方向的位移,需要在水平方向上先加速然后再減速。第2個規(guī)律是初始航向和橫向速度越大則初始推力俯仰角越大,原因是制導(dǎo)需要在水平方向分配更多的推力來抵消初始的反向水平速度。

圖2 有效性驗證仿真工況的軌跡Fig.2 Trajectories of effectiveness verification simulation cases

圖3 有效性驗證仿真工況4的橫向位置速度Fig.3 Cross position and velocity of effectiveness verification simulation case 4

圖4 有效性驗證仿真工況的推力比Fig.4 Thrust throttles of effectiveness verification simulation cases

圖5 有效性驗證仿真工況的推力俯仰角Fig.5 Thrust pitch angles of effectiveness verification simulation cases

制導(dǎo)有效性驗證仿真工況的制導(dǎo)時長和推進(jìn)劑消耗仿真結(jié)果見表2。從表中可知,初始航向和橫向速度越大則制導(dǎo)時長和推進(jìn)劑消耗越長和越多。工況2由于初始航向速度與目標(biāo)水平移動方向相同,因而制導(dǎo)時長最短并且推進(jìn)劑消耗最少。表中給出了預(yù)測校正法和遍歷法兩種方法的推進(jìn)劑消耗來對比。遍歷法的搜索步長取為1 s?？梢钥闯?預(yù)測校正法的推進(jìn)劑消耗十分接近遍歷法,最大僅相差1.1%。表中還給出了兩種預(yù)測校正的實際迭代次數(shù)。工況1實際執(zhí)行式(31)所示的預(yù)測校正0次,執(zhí)行式(33)所示的預(yù)測校正2次。工況2實際執(zhí)行式(31)所示的預(yù)測校正1次,執(zhí)行式(33)所示的預(yù)測校正0次。所有4個工況僅需執(zhí)行1到2次預(yù)測校正便可獲得滿意的結(jié)果,說明所提出的預(yù)測校正方法不僅有效而且高效。

表2 有效性驗證仿真工況的制導(dǎo)時長和推進(jìn)劑消耗Table 2 Time-to-go and fuel consumption of effectiveness verification simulation cases

4.2 碰撞規(guī)避驗證

文獻(xiàn)[5,9,11]中為了測試驗證制導(dǎo)算法碰撞規(guī)避性能,設(shè)計了惡劣的碰撞規(guī)避仿真工況,初始速度取為[100, 0, -75]Tm/s。在此初始速度條件下,分別按照無碰撞規(guī)避、有地平規(guī)避、有8°下滑角約束和有9°下滑角約束進(jìn)行制導(dǎo)方法的測試仿真,對比驗證3.2節(jié)中的碰撞規(guī)避方法。碰撞規(guī)避仿真工況的軌跡見圖6。從中可知,無碰撞規(guī)避的制導(dǎo)飛行軌跡有部分位于了地平線下,該結(jié)果意味著著陸器將與地面碰撞?？紤]地平規(guī)避的制導(dǎo)飛行軌跡在著陸后期緊緊貼著地平,實現(xiàn)了臨界的碰撞規(guī)避,該結(jié)果具有理論分析意義,但是不能用于實際飛行,因為其安全裕度不夠高,實際地面必然是具有一定幅度的起伏。有8°和9°下滑角約束的著陸軌跡可以明顯看到著陸后期保持了一定的飛行高度,與地平線形成了固定夾角,顯然該飛行軌跡的碰撞規(guī)避裕度更高,對一定幅度的地形突起也可以實現(xiàn)碰撞規(guī)避,因而更具有工程應(yīng)用價值。工況2到4的高度是單調(diào)變化的,另外著陸后期的軌跡線為直線,原因是著陸后期在下滑角坐標(biāo)系下垂向速度已被控制到零。文獻(xiàn)[5]利用航路點方法作碰撞規(guī)避的后期航跡也近似直線,然而航路點前后制導(dǎo)推力不連續(xù)。文獻(xiàn)[9,11]中碰撞規(guī)避段的高度是先探底然后上升之后再下降,也就是高度非單調(diào)下降,并且文獻(xiàn)[11]中也有制導(dǎo)推力躍變點。碰撞規(guī)避仿真工況的垂向速度見圖7,從中可知,工況1的垂向速度先為負(fù)值后為正值,工況2和3的垂向速度全程為非正值,工況4的垂向速度先為負(fù)值中間為正值后為負(fù)值。工況4的垂向速度出現(xiàn)正值的原因是在碰撞規(guī)避控制段的推力始終飽和。碰撞規(guī)避仿真工況的推力比曲線見圖8,從中可知,全部工況的推力由飽和段、減小段和增加段組成,即推力大致按照大小大的趨勢變化。比較工況2到4可知,碰撞規(guī)避的約束下滑角越大,制導(dǎo)時長越長。另外,考慮碰撞規(guī)避后制導(dǎo)初期的推力飽和段時長增大,并且碰撞規(guī)避裕度越大,制導(dǎo)初期推力飽和段時長增大的越多。碰撞規(guī)避仿真工況的推力俯仰角見圖9,從中可知,碰撞規(guī)避要求越嚴(yán)格則推力的初始俯仰角越小,如果初始推力飽和則初始俯仰角等于負(fù)的下滑角。該現(xiàn)象的原因是,為了實現(xiàn)碰撞規(guī)避就需要減小下降速度,因此制導(dǎo)算法在垂向分配的推力更大,因此推力俯仰角就會相對變小。工況4的推力俯仰角在30 s時發(fā)生3.1°的躍變,原因是碰撞規(guī)避控制結(jié)束時刻已到,但是推力還處于飽和狀態(tài),因此導(dǎo)致制導(dǎo)加速度變化不連續(xù)。工況4說明如果希望制導(dǎo)推力連續(xù),則在碰撞規(guī)避控制結(jié)束時制導(dǎo)推力應(yīng)處于非飽和狀態(tài)。工況3在保證指令推力連續(xù)的前提下滿足了8°下滑角約束,大于文獻(xiàn)[5,11]中相同條件下實現(xiàn)的4°。

圖6 碰撞規(guī)避仿真工況的軌跡Fig.6 Trajectories of collision avoidance simulation cases

圖7 碰撞規(guī)避仿真工況的垂向速度Fig.7 Vertical velocities of collision avoidance simulation cases

圖8 碰撞規(guī)避仿真工況的推力比Fig.8 Thrust throttles of collision avoidance simulation cases

圖9 碰撞規(guī)避仿真工況的推力俯仰角Fig.9 Thrust pitch angles of collision avoidance simulation cases

碰撞規(guī)避仿真工況的制導(dǎo)時長和推進(jìn)劑消耗見表3。從中可知,有地平規(guī)避與無碰撞規(guī)避相比制導(dǎo)時長減少3.2 s推進(jìn)劑消耗增加6.5 kg,有8°下滑角約束與無碰撞規(guī)避相比制導(dǎo)時長增加0.7 s推進(jìn)劑消耗增加18.3 kg,有9°下滑角約束與無碰撞規(guī)避相比制導(dǎo)時長增加13.9 s推進(jìn)劑消耗增加75.6 kg。本節(jié)仿真結(jié)果說明新制導(dǎo)方法的碰撞規(guī)避策略有效,并且滿足碰撞規(guī)避要求后付出的推進(jìn)劑代價相對小。

表3 碰撞規(guī)避仿真工況的制導(dǎo)時長和推進(jìn)劑消耗Table 3 Time-to-go and fuel consumptions of collision avoidance simulation cases

4.3 有終端加速度約束

設(shè)計有終端加速度約束仿真工況來驗證3.3節(jié)中方法,初始速度取為[0, 0, -75]Tm/s,目標(biāo)終端加速度為[0, 0, 0]Tm/s2。姿態(tài)機(jī)動最大角速度和最大推力變化率取值見表4。預(yù)測校正迭代計算的終止條件為終端位置偏差小于1 m,速度偏差小于0.1 m/s。

表4 有終端加速度約束仿真工況的最大角速度和最大推力變化率Table 4 Maximum angular rate and maximum thrust change rate of terminal acceleration constrained simulation cases

有終端加速度約束仿真工況的軌跡見圖10,從中可知,著陸前中期3條軌跡基本重合,著陸后期3條軌跡略有差異,終端時刻3條軌跡都滿足了終端位置約束。有終端加速度約束仿真工況的推力比見圖11,以工況1為例,前40.2 s為最優(yōu)反饋制導(dǎo),之后為勻速轉(zhuǎn)動制導(dǎo),推力比從0.95降至0.47。有終端加速度約束仿真工況的俯仰角見圖12,以工況1為例,在40.2 s處可以清楚的看到俯仰角發(fā)生轉(zhuǎn)折,這是因為制導(dǎo)律切換為了勻速轉(zhuǎn)動制導(dǎo),俯仰角迅速從54.7°轉(zhuǎn)動到0°。推力比和俯仰角的仿真結(jié)果表明制導(dǎo)終端狀態(tài)滿足了終端加速度約束。

圖10 有終端加速度約束仿真工況的軌跡Fig.10 Trajectories of terminal acceleration constrained simulation cases

圖11 有終端加速度約束仿真工況的推力比Fig.11 Thrust throttles of terminal acceleration constrained simulation cases

圖12 有終端加速度約束仿真工況的俯仰角Fig.12 Thrust pitch angles of terminal acceleration constrained simulation cases

有終端加速度約束仿真工況的數(shù)值結(jié)果見表5。從表中可知,雖然勻速轉(zhuǎn)動制導(dǎo)時長隨著最大角速度和最大推力變化率增大而減小,但是總制導(dǎo)時長并沒有出現(xiàn)等幅度變化,而是變化很小。勻速轉(zhuǎn)動制導(dǎo)時長越短則推進(jìn)劑消耗越少。工況1的推進(jìn)劑消耗最多,與4.1節(jié)中無終端加速度約束工況相比僅增加了6.3 kg,相對增量為2.7%。工況1到3的預(yù)測校正迭代次數(shù)依次減少,說明勻速轉(zhuǎn)動制導(dǎo)允許的最大角速度和最大推力變化率越大,則預(yù)測校正需要的迭代次數(shù)越少。以上結(jié)果表明最優(yōu)反饋制導(dǎo)和勻速轉(zhuǎn)動制導(dǎo)的組合使用方法有效,使得著陸終端既滿足了位置速度約束,又滿足了加速度約束。

表5 有終端加速度約束仿真工況的制導(dǎo)時長和推進(jìn)劑消耗Table 5 Time-to-go and fuel consumption of terminal acceleration constrained simulation cases

5 結(jié) 論

針對行星定點著陸動力下降段制導(dǎo)問題,從工程實際出發(fā),考慮推力范圍受限、連續(xù)變推力、碰撞規(guī)避和終端加速度共4方面約束,對傳統(tǒng)的最優(yōu)反饋制導(dǎo)方法進(jìn)行改進(jìn)和完善。首先,提出推力范圍受限且連續(xù)變推力約束下的最優(yōu)反饋制導(dǎo)的制導(dǎo)時長高效尋優(yōu)方法,經(jīng)1到2次預(yù)測校正即可獲得推進(jìn)劑消耗量接近遍歷法的結(jié)果,推力剖面為大小大形式。其次,提出基于碰撞規(guī)避時長和推力優(yōu)化分配律的碰撞規(guī)避方法,結(jié)合坐標(biāo)系旋轉(zhuǎn)可滿足下滑角約束,在連續(xù)變推力前提下顯著提高碰撞規(guī)避能力和著陸點可見性,著陸全程高度單調(diào)下降,碰撞規(guī)避控制結(jié)束后的縱向平面軌跡為直線。最后,提出勻速轉(zhuǎn)動制導(dǎo)和基于預(yù)測校正的組合制導(dǎo)參數(shù)規(guī)劃方法,連續(xù)銜接最優(yōu)反饋制導(dǎo)和勻速轉(zhuǎn)動制導(dǎo),充分發(fā)揮推力調(diào)節(jié)能力,使得終端狀態(tài)擴(kuò)充滿足加速度約束。新制導(dǎo)方法改進(jìn)原理清晰,多約束適應(yīng)性強(qiáng),雖然使用了數(shù)值積分預(yù)報,但是總體還屬于解析制導(dǎo)方法,計算量顯著低于凸優(yōu)化等數(shù)值優(yōu)化方法和人工智能方法,適合工程應(yīng)用于行星定點著陸任務(wù)中。