本發(fā)明涉及一種小天體附著探測(cè)下降軌跡優(yōu)化方法，屬于航天航空領(lǐng)域。

背景技術(shù)：

小天體探測(cè)作為人類(lèi)了解太陽(yáng)系形成與演化、生命起源與進(jìn)化以及防御外來(lái)天體撞擊的重要途徑，將是未來(lái)深空探測(cè)活動(dòng)的主要內(nèi)容之一，附著探測(cè)是未來(lái)一段時(shí)間內(nèi)人類(lèi)探索小天體的主要方式。下降階段是著陸器附著小天體或完成采樣返回探測(cè)的關(guān)鍵階段，對(duì)能否安全、準(zhǔn)確到達(dá)預(yù)設(shè)的具有科學(xué)探測(cè)價(jià)值的目標(biāo)區(qū)域起著決定性的作用，這對(duì)下降階段的軌跡設(shè)計(jì)、導(dǎo)航和制導(dǎo)控制都提出了很高的要求。小天體著陸器的下降過(guò)程可以轉(zhuǎn)化為一個(gè)軌跡優(yōu)化問(wèn)題以及對(duì)標(biāo)稱軌跡的跟蹤控制問(wèn)題。設(shè)定的標(biāo)稱軌跡需要能夠安全、準(zhǔn)確地到達(dá)指定著陸點(diǎn)，滿足始末狀態(tài)約束、路徑約束、控制約束等多重約束，同時(shí)使某項(xiàng)重要的性能指標(biāo)最優(yōu)化，如燃耗、時(shí)間等。軌跡優(yōu)化方法主要包括基于參數(shù)化方法的直接法和基于龐特里亞金極小值原理的間接法。直接法無(wú)需推導(dǎo)優(yōu)化問(wèn)題的橫截條件，避免了求解兩點(diǎn)邊值問(wèn)題的協(xié)態(tài)初值敏感困難，因而得到了廣泛的應(yīng)用。此外，對(duì)于小天體附近的軌跡優(yōu)化問(wèn)題，還需要建立精確描述小天體附近引力加速度的引力場(chǎng)模型。

在已發(fā)展的小天體附著探測(cè)下降軌跡優(yōu)化方法中，在先技術(shù)[1](參見(jiàn)Lantoine G,Braun R.Optimal trajectories for soft landing on asteroids.Space Systems Design Lab,Georgia Institute of Technology,Atlanta,GA,AE8900 MS Special Problems Report,Dec.2006.)采用多面體模型求解目標(biāo)小天體附近的引力加速度，以能量作為優(yōu)化指標(biāo)，運(yùn)用直接法進(jìn)行大步長(zhǎng)優(yōu)化，從而估計(jì)得到優(yōu)化問(wèn)題的協(xié)態(tài)初值，然后基于龐特里亞金原理進(jìn)行間接法軌跡優(yōu)化計(jì)算。多面體模型求引力加速度效率低，優(yōu)化算法繁瑣，耗時(shí)長(zhǎng)。

在先技術(shù)[2](參見(jiàn)Ren,Y.and Shan,J.,“Reliability-Based Soft Landing Trajectory Optimization near Asteroid with Uncertain Gravitational Field,”Journal of Guidance,Control,and Dynamics,Vol.38,No.9,2015,pp.1810-1820.)，首先構(gòu)建收斂區(qū)域更明確的能量最優(yōu)問(wèn)題，然后通過(guò)調(diào)整同倫系數(shù)，序列求解，最終將能量最優(yōu)問(wèn)題轉(zhuǎn)化為燃料最優(yōu)問(wèn)題，采用同倫算法仍然需解兩點(diǎn)邊值問(wèn)題，優(yōu)化過(guò)程耗時(shí)較長(zhǎng)。

技術(shù)實(shí)現(xiàn)要素：

本發(fā)明目的是為了解決現(xiàn)有小天體附著探測(cè)下降軌跡優(yōu)化方法因采用多面體引力場(chǎng)模型求引力加速度效率差；以及采用間接法進(jìn)行軌跡優(yōu)化解算，因協(xié)態(tài)初值難以估計(jì)求解困難的問(wèn)題，提供一種小天體附著探測(cè)下降軌跡優(yōu)化方法。

一種小天體附著探測(cè)下降軌跡優(yōu)化方法，采用內(nèi)球諧引力場(chǎng)模型估計(jì)目標(biāo)小天體附近的引力加速度，采用凸優(yōu)化算法解軌跡優(yōu)化方程。

所述凸優(yōu)化算法包括松弛、線性化、離散化和內(nèi)點(diǎn)法。

一種小天體附著探測(cè)下降軌跡優(yōu)化方法，包括以下步驟：

步驟一、由最小二乘法估計(jì)內(nèi)球諧引力場(chǎng)模型的球諧系數(shù)：

在小天體外部構(gòu)造內(nèi)布里淵球，所構(gòu)造內(nèi)布里淵球應(yīng)該與小天體表面的目標(biāo)著陸點(diǎn)相切，球心選取應(yīng)保證著陸器下降軌跡包含在內(nèi)布里淵球內(nèi)部。在內(nèi)布里淵球內(nèi)部任取N_data個(gè)點(diǎn)，由多面體模型計(jì)算N_data個(gè)點(diǎn)的引力加速度，并由所述引力加速度構(gòu)造3N_data×1維矩陣由最小二乘法求取內(nèi)球諧引力模型球諧系數(shù)。

其中，n和m分別表示勒讓德多項(xiàng)式的階次和冪次，為一個(gè)(n²+2n)×1維的矩陣，包含了所要求取的各階內(nèi)球諧系數(shù)

為引力加速度關(guān)于內(nèi)球諧系數(shù)的導(dǎo)數(shù)構(gòu)成的矩陣，維度為(n²+2n)×3N_data。W為一個(gè)3N_data維單位陣；

步驟二、構(gòu)造小天體附著探測(cè)下降軌跡優(yōu)化方程：

在小天體固連坐標(biāo)系下，著陸器滿足如下的動(dòng)力學(xué)方程

其中，r＝[x,y,z]^T和分別為固連坐標(biāo)系下著陸器的位置矢量和速度矢量；ω＝[0,0,ω]^T為小天體自旋角速度矢量；T＝[T_x,T_y,T_z]^T為著陸器推力矢量；m_e為著陸器的質(zhì)量；I_sp為發(fā)動(dòng)機(jī)比沖；g_e＝9.807為地球引力加速度常數(shù)；為引力加速度矢量，由以下公式計(jì)算得到：

其中，G為萬(wàn)有引力常數(shù)，M為目標(biāo)小天體質(zhì)量，R為步驟一中內(nèi)布里淵球半徑，δ_0,m為Kronecker delta(克勞內(nèi)克)函數(shù)(當(dāng)m＝0時(shí)，δ_0,m＝1)。為步驟一所求的內(nèi)球諧系數(shù)；為步驟一中內(nèi)球諧引力場(chǎng)模型基函數(shù)，滿足如下遞推關(guān)系。

著陸器滿足如下的邊界條件

其中，r₀，v₀和m₀分別為起始位置處著陸器的位置、速度和質(zhì)量。t_f為終端時(shí)間，r_f為目標(biāo)著陸點(diǎn)，v_f為終端速度約束，對(duì)于軟著陸問(wèn)題，v_f＝0。

著陸器推力矢量滿足如下約束條件

0≤||T||≤T_max (7)

燃料最優(yōu)性能指標(biāo)表示如下

公式(3)、(6)、(7)、(8)構(gòu)成了下降軌跡優(yōu)化方程；

步驟三、向步驟二所得的下降軌跡優(yōu)化方程中引入松弛變量Γ，對(duì)所述下降軌跡優(yōu)化方程進(jìn)行松弛：

引入松弛變量Γ替代下降軌跡優(yōu)化方程中的||T||，則松弛后的軌跡優(yōu)化方程為：

步驟四、對(duì)步驟三所得方程線性化處理：

定義新變量如下

將以上變量代入步驟三的優(yōu)化方程中，得到線性化的優(yōu)化方程為

其中，t表示時(shí)間變量。

步驟五、對(duì)步驟四所得的優(yōu)化方程進(jìn)行離散化處理：

將時(shí)間區(qū)間[0,t_f]等分成N份，得到對(duì)步驟四的優(yōu)化方程進(jìn)行離散化處理，經(jīng)過(guò)離散化以后，軌跡優(yōu)化方程轉(zhuǎn)換為參數(shù)優(yōu)化方程，參數(shù)優(yōu)化方程的表達(dá)式如下：

其中，M＝[I₆,0_6×1]，t_k表示第k個(gè)時(shí)間節(jié)點(diǎn)，u_k和σ_k分別表示t_k時(shí)刻u和σ的取值。

步驟六、采用內(nèi)點(diǎn)法對(duì)步驟五中的參數(shù)優(yōu)化方程進(jìn)行迭代求解，具體求解過(guò)程如下：

1)令引力加速度為一常值▽V₀，采用內(nèi)點(diǎn)法求解步驟五中的參數(shù)優(yōu)化方程，得到一條軌跡；

2)將步驟1)所得到的軌跡作為參考軌跡，計(jì)算參考軌跡中每個(gè)節(jié)點(diǎn)(即k＝0,…,N)處的引力加速度，并帶入步驟五中的參數(shù)優(yōu)化方程中，采用內(nèi)點(diǎn)法進(jìn)行求解，得到一條新的軌跡，將得到的新軌跡作為下一次迭代的參考軌跡；

3)當(dāng)?shù)玫降能壽E收斂，則得到最優(yōu)解，即得到探測(cè)最優(yōu)下降軌跡。

步驟五所述離散化處理方法采用顯式四階龍格庫(kù)塔積分公式法；

有益效果

本發(fā)明所給出的一種小天體附著探測(cè)下降軌跡優(yōu)化方法，采用內(nèi)球諧引力場(chǎng)模型估計(jì)小天體附近引力加速度，具有計(jì)算效率高的優(yōu)點(diǎn)。采用凸優(yōu)化算法求解燃耗最優(yōu)控制問(wèn)題，避免了間接法的推導(dǎo)復(fù)雜，以及協(xié)態(tài)變量無(wú)物理意義不易猜測(cè)的問(wèn)題，推導(dǎo)步驟較為簡(jiǎn)便，節(jié)約了計(jì)算時(shí)間，是一種快速、高精度的估計(jì)優(yōu)化方法，且所得結(jié)果符合初始和末端狀態(tài)約束、動(dòng)力學(xué)約束和控制約束。

附圖說(shuō)明

圖1是本發(fā)明方法的流程圖；

圖2是仿真結(jié)果示意圖，共經(jīng)過(guò)四次迭代得到最有軌跡；其中(a)是小天體固連坐標(biāo)系下著陸器沿x軸方向的軌跡，(b)是沿y軸方向軌跡示意圖，(c)是沿z軸方向軌跡示意圖，(d)是推力大小示意圖。

具體實(shí)施方式

下面結(jié)合附圖與實(shí)施例對(duì)本發(fā)明作進(jìn)一步說(shuō)明。

實(shí)施例1

一種小天體附著探測(cè)下降軌跡優(yōu)化方法，包括以下步驟：

步驟一、由最小二乘法估計(jì)內(nèi)球諧引力場(chǎng)模型的球諧系數(shù)：

在小天體外部構(gòu)造內(nèi)布里淵球，所構(gòu)造內(nèi)布里淵球應(yīng)該與小天體表面的目標(biāo)著陸點(diǎn)相切，球心選取應(yīng)保證著陸器下降軌跡包含在內(nèi)布里淵球內(nèi)部。在內(nèi)布里淵球內(nèi)部任取N_data個(gè)點(diǎn)，由多面體模型計(jì)算N_data個(gè)點(diǎn)的引力加速度，并由所述引力加速度構(gòu)造3N_data×1維矩陣由最小二乘法求取內(nèi)球諧引力模型球諧系數(shù)。

其中，n和m分別表示勒讓德多項(xiàng)式的階次和冪次，為一個(gè)(n²+2n)×1維的矩陣，包含了所要求取的各階內(nèi)球諧系數(shù)

為引力加速度關(guān)于內(nèi)球諧系數(shù)的導(dǎo)數(shù)構(gòu)成的矩陣，維度為(n²+2n)×3N_data。W為一個(gè)3N_data維單位陣；

步驟二、構(gòu)造小天體附著探測(cè)下降軌跡優(yōu)化方程：

在小天體固連坐標(biāo)系下，著陸器滿足如下的動(dòng)力學(xué)方程

其中，r＝[x,y,z]^T和分別為固連坐標(biāo)系下著陸器的位置矢量和速度矢量；ω＝[0,0,ω]^T為小天體自旋角速度矢量；T＝[T_x,T_y,T_z]^T為著陸器推力矢量；m_e為著陸器的質(zhì)量；I_sp為發(fā)動(dòng)機(jī)比沖；g_e＝9.807為地球引力加速度常數(shù)；為引力加速度矢量，由以下公式計(jì)算得到：

其中，G為萬(wàn)有引力常數(shù)，M為目標(biāo)小天體質(zhì)量，R為步驟一中內(nèi)布里淵球半徑，δ_0,m為(克勞內(nèi)克)Kronecker delta函數(shù)(當(dāng)m＝0時(shí)，δ_0,m＝1)。為步驟一所求的內(nèi)球諧系數(shù)；為步驟一中內(nèi)球諧引力場(chǎng)模型基函數(shù)，滿足如下遞推關(guān)系。

著陸器滿足如下的邊界條件

其中，r₀，v₀和m₀分別為起始位置處著陸器的位置、速度和質(zhì)量。t_f為終端時(shí)間，r_f為目標(biāo)著陸點(diǎn)，v_f為終端速度約束，對(duì)于軟著陸問(wèn)題，v_f＝0。

著陸器推力矢量滿足如下約束條件

0≤||T||≤T_max (19)

燃料最優(yōu)性能指標(biāo)表示如下

公式(15)、(18)、(19)、(20)構(gòu)成了下降軌跡優(yōu)化方程；

步驟三、向步驟二所得的下降軌跡優(yōu)化方程中引入松弛變量Γ，對(duì)所述下降軌跡優(yōu)化方程進(jìn)行松弛：

引入松弛變量Γ替代下降軌跡優(yōu)化方程中的||T||，則松弛后的軌跡優(yōu)化方程為：

步驟四、對(duì)步驟三所得方程線性化處理：

定義新變量如下

將以上變量代入步驟三的優(yōu)化方程中，得到線性化的優(yōu)化方程為

其中，t表示時(shí)間變量。

步驟五、對(duì)步驟四所得的優(yōu)化方程進(jìn)行離散化處理：

將時(shí)間區(qū)間[0,t_f]等分成N份，可以得到由顯式四階龍格庫(kù)塔積分公式，對(duì)步驟四的優(yōu)化方程進(jìn)行離散化

其中，r_k，v_k和p_k分別表示r，v和p在時(shí)間節(jié)點(diǎn)t_k處的值，U_k和▽V_k分別表示狀態(tài)變量、控制輸入和引力加速度在時(shí)間節(jié)點(diǎn)t_k處的值。

u_k和σ_k分別表示t_k時(shí)刻u和σ的取值。經(jīng)過(guò)離散化以后，軌跡優(yōu)化方程轉(zhuǎn)換為一個(gè)參數(shù)優(yōu)化方程，其表達(dá)式如下：

其中，M＝[I₆,0_6×1]。

步驟六、采用內(nèi)點(diǎn)法對(duì)步驟五中的參數(shù)優(yōu)化方程進(jìn)行迭代求解，具體求解過(guò)程如下：

1)令引力加速度為一常值▽V₀，采用內(nèi)點(diǎn)法求解步驟五中的參數(shù)優(yōu)化方程，得到一條軌跡；

2)將步驟1)所得到的軌跡作為參考軌跡，計(jì)算參考軌跡中每個(gè)節(jié)點(diǎn)(即k＝0,…,N)處的引力加速度，并帶入步驟五中的參數(shù)優(yōu)化方程中，采用內(nèi)點(diǎn)法進(jìn)行求解，得到一條新的軌跡，將得到的新軌跡作為下一次迭代的參考軌跡；

3)當(dāng)?shù)玫降能壽E收斂，則得到最優(yōu)解，即得到探測(cè)最優(yōu)下降軌跡。

圖2是在小天體216Kleopatra上進(jìn)行的仿真示意圖。著陸器的初始位置為[-150108,6010,-1034]m，目標(biāo)位置為[-112600,6340,-12520]m。采用本發(fā)明所提出的小天體附著探測(cè)下降軌跡優(yōu)化方法，共經(jīng)過(guò)四次迭代，得到最終的燃耗最優(yōu)軌跡。由圖2可知，優(yōu)化結(jié)果滿足各項(xiàng)約束條件，著陸器以零速度軟著陸于預(yù)設(shè)著陸點(diǎn)，整個(gè)計(jì)算過(guò)程耗時(shí)11.4秒，同時(shí)優(yōu)化得到的控制力始終在設(shè)定的發(fā)動(dòng)機(jī)推力上限之內(nèi)。