太陽(yáng)系多目標(biāo)探測(cè)軌道優(yōu)化設(shè)計(jì)*

楊大林 徐 波 高有濤 張 磊

1.南京航空航天大學(xué)航天學(xué)院，南京 210016 2.南京大學(xué)天文與空間科學(xué)學(xué)院，南京 210093

行星際探測(cè)已成為航天領(lǐng)域的一大熱點(diǎn)，備受航天大國(guó)的關(guān)注，但行星際探測(cè)具有能耗高，持續(xù)時(shí)間長(zhǎng)，任務(wù)約束多，發(fā)射窗口小等特點(diǎn)，為此任務(wù)設(shè)計(jì)過(guò)程中要求完成科學(xué)探測(cè)目標(biāo)的同時(shí)，盡量降低任務(wù)的成本和周期。因此攜帶盡可能多的有效載荷，以最小的能量消耗，在最短的飛行時(shí)間內(nèi)，到訪多個(gè)目標(biāo)成為太陽(yáng)系行星際探測(cè)發(fā)展的趨勢(shì)。

迄今為止，行星際探測(cè)軌道設(shè)計(jì)領(lǐng)域的主要研究方向有多天體交會(huì)的行星借力飛行技術(shù)、行星際低能轉(zhuǎn)移軌道優(yōu)化設(shè)計(jì)技術(shù)、小推力軌道設(shè)計(jì)技術(shù)等。行星借力飛行技術(shù)[1]即飛越大行星時(shí)借助大行星引力場(chǎng)效應(yīng)，獲取所需的速度增量，而飛越大行星時(shí)也為探測(cè)器探測(cè)該行星提供了機(jī)會(huì)，探測(cè)器進(jìn)入借力行星影響球到離開(kāi)該行星影響球所持續(xù)的時(shí)間大約為20多個(gè)小時(shí)[2]，雖然相對(duì)于整個(gè)任務(wù)時(shí)間而言非常小，可以忽略不計(jì)，但這段時(shí)間足夠獲取該行星的基本資料，并可以在此期間向該行星釋放低質(zhì)量且不回收的科學(xué)探測(cè)設(shè)備。

國(guó)內(nèi)外的許多學(xué)者針對(duì)行星際軌道優(yōu)化設(shè)計(jì)問(wèn)題進(jìn)行了研究，提出了多種行之有效的優(yōu)化方法[3-7]，這些方法可以分為解析法和數(shù)值解法。大多數(shù)解析法以變分法得到的必要條件為基礎(chǔ)，將軌跡優(yōu)化問(wèn)題轉(zhuǎn)換成兩點(diǎn)邊界值問(wèn)題，求解這類問(wèn)題富有挑戰(zhàn)并且需要數(shù)值解法，同時(shí)該問(wèn)題對(duì)終端條件敏感，求解過(guò)程中容易使目標(biāo)軌道陷于局部極小，這些特點(diǎn)使得此類方法的實(shí)際應(yīng)用受到了限制[3-5]。數(shù)值法則是將軌跡優(yōu)化問(wèn)題簡(jiǎn)化為參數(shù)優(yōu)化問(wèn)題，再通過(guò)求解非線性規(guī)劃問(wèn)題來(lái)解決[6-7]，這里提到的“數(shù)值算法”是指不明確使用變分法中的解析必要條件，而是使用打靶法等數(shù)值解法求解兩點(diǎn)邊界值問(wèn)題，求解過(guò)程中需要給該算法提供一個(gè)初始“估計(jì)值”。在軌跡優(yōu)化問(wèn)題中，初值通常不易確定，并且求解過(guò)程中容易陷入給定初值附近的局部極小點(diǎn)。Bett[8]全面分析了上述2種方法的特點(diǎn)以及2種方法的聯(lián)系。

進(jìn)化算法是航天器軌跡優(yōu)化的一種有效方法，是一種模擬自然界中生物進(jìn)化過(guò)程與機(jī)制的數(shù)值優(yōu)化方法。與解析法和傳統(tǒng)的數(shù)值解法相比，進(jìn)化算法主要有以下優(yōu)勢(shì)：進(jìn)化算法相對(duì)簡(jiǎn)單，不需要提供初始“估計(jì)值”(隨機(jī)產(chǎn)生一群初始解)，而且比傳統(tǒng)的優(yōu)化方法更容易找到全局最小點(diǎn)。Cage將遺傳算法(GA)應(yīng)用到行星際軌跡優(yōu)化設(shè)計(jì)中[9]，設(shè)計(jì)了地球—火星和地球—木星的轉(zhuǎn)移軌道。Bessette和Spencer將粒子群優(yōu)化(PSO)技術(shù)應(yīng)用到時(shí)間受限的轉(zhuǎn)移軌道和行星際脈沖轉(zhuǎn)移軌道的優(yōu)化設(shè)計(jì)中[10-11]。但是經(jīng)典的進(jìn)化算法也具有明顯的缺點(diǎn)，即由于隨機(jī)性，并不能保證算法一定收斂到全局最優(yōu)?？紤]到太陽(yáng)系多目標(biāo)探測(cè)軌道優(yōu)化問(wèn)題的復(fù)雜性和多峰性，本文融合改進(jìn)的網(wǎng)格搜索算法和差分進(jìn)化算法，提出一種組合優(yōu)化算法來(lái)求解太陽(yáng)系多目標(biāo)探測(cè)軌道優(yōu)化設(shè)計(jì)問(wèn)題。以2018～2020年太陽(yáng)系多目標(biāo)探測(cè)為例，通過(guò)對(duì)最終環(huán)繞土星，中途探測(cè)太陽(yáng)系內(nèi)多顆大行星的深空軌道的優(yōu)化設(shè)計(jì)，表明了該算法可以明顯地減少計(jì)算時(shí)間，并有效收斂到全局范圍內(nèi)的最優(yōu)解。

1 太陽(yáng)系多目標(biāo)探測(cè)動(dòng)力學(xué)模型

太陽(yáng)系多目標(biāo)探測(cè)行星際飛越軌道可以用包含非線性約束的有限維全局優(yōu)化問(wèn)題來(lái)描述。在J2000日心黃道坐標(biāo)系中，只考慮太陽(yáng)的中心引力，不考慮其他天體的引力和各天體的非球形引力攝動(dòng)，探測(cè)器的動(dòng)力學(xué)方程如下：

(1)

式中：μs為太陽(yáng)引力常數(shù)；rs為探測(cè)器在J2000日心黃道坐標(biāo)系中的位置矢量；rs為位置標(biāo)量。探測(cè)器與行星交會(huì)時(shí)，在交會(huì)行星的影響球內(nèi)只考慮該行星的引力效應(yīng)，并且考慮動(dòng)力飛越的情形，即在探測(cè)器與行星交會(huì)時(shí)，可以對(duì)探測(cè)器施加速度脈沖(參考1.1節(jié))。此時(shí)，探測(cè)器的動(dòng)力學(xué)方程如下：

(2)

其中：μp為行星引力常數(shù)；rp為探測(cè)器在以行星為中心的坐標(biāo)系中的位置矢量；rp為位置標(biāo)量。

1.1 動(dòng)力飛越軌道

(3)

類似地，在離開(kāi)交點(diǎn)處

(4)

(5)

(6)

由于進(jìn)入段和離開(kāi)段的近地點(diǎn)半徑相同，則

rper=a1(1-e1)=a2(1-e2)

(7)

同時(shí)轉(zhuǎn)角之間還存在如下關(guān)系

(8)

僅利用式(7)和(8)無(wú)法有效的求解e1和e2，因此需要對(duì)上述表達(dá)式作適當(dāng)轉(zhuǎn)換

(9)

對(duì)方程(9)采用牛頓迭代法可以求得e2，再根據(jù)式(7)即可求得e1和rper。根據(jù)活力公式可以求得進(jìn)入和離開(kāi)雙曲線軌道的近星點(diǎn)速度大小。

(10)

由式(11)可得近地點(diǎn)所需施加的速度脈沖大小為

(11)

1.2 問(wèn)題描述

(12)

因此，為使整個(gè)過(guò)程中消耗的燃料最小，本文引入下述性能指標(biāo)：

(13)

為了求解的方便，可以將上述問(wèn)題轉(zhuǎn)化為一般的非線性規(guī)劃問(wèn)題[12]：

優(yōu)化指標(biāo)函數(shù)：

(14)

服從的動(dòng)力學(xué)特性：

(15)

需滿足的邊界條件：

τ∈I

(16)

其中τ=[t0,T1,...,TN+1]，I∈RN+2。

需滿足的段匹配約束：

?i=1,...,N+1

(17)

需滿足的約束條件：

(18)

圖1 太陽(yáng)系多目標(biāo)探測(cè)行星際飛越軌道求解流程圖

2 軌道優(yōu)化設(shè)計(jì)算法

優(yōu)化算法的目的是使與時(shí)變動(dòng)力學(xué)系統(tǒng)相關(guān)的目標(biāo)方程最小化，問(wèn)題可以簡(jiǎn)化成非線性規(guī)劃問(wèn)題。通過(guò)對(duì)比現(xiàn)有的數(shù)值優(yōu)化方法(如列舉法、隨機(jī)試驗(yàn)法、復(fù)合形法、可行方向法、進(jìn)化算法)，列舉法是唯一可以保證搜索到全局最優(yōu)的算法(前提是步長(zhǎng)要選擇的合適)，但列舉法所需的時(shí)間隨著搜索空間的維數(shù)成指數(shù)增長(zhǎng)，針對(duì)太陽(yáng)系多目標(biāo)探測(cè)軌跡優(yōu)化問(wèn)題的復(fù)雜性，為了減少計(jì)算時(shí)間，引入下述改進(jìn)的網(wǎng)格搜索算法。

2.1 改進(jìn)后的網(wǎng)格搜索算法

將改進(jìn)網(wǎng)格搜索算法的每一步分別與太陽(yáng)系多目標(biāo)探測(cè)的每一段轉(zhuǎn)移軌道相對(duì)應(yīng)。通過(guò)不斷調(diào)整時(shí)間序列，可以得到任務(wù)的起始時(shí)間和每一階段理想的轉(zhuǎn)移時(shí)間，進(jìn)而得到滿足整個(gè)任務(wù)要求的完整軌道。應(yīng)用改進(jìn)的網(wǎng)格搜索算法，定義太陽(yáng)系多目標(biāo)探測(cè)每階段的優(yōu)化模型如下：

1)地球到第1個(gè)飛越天體，時(shí)間序列：τ1=[t0,T1](其中t0為地球發(fā)射時(shí)間，T1為從地球到第1個(gè)飛越天體的轉(zhuǎn)移時(shí)間)

(19)

(20)

(21)

圖2 改進(jìn)的網(wǎng)格搜索算法流程圖

2.2 差分進(jìn)化算法

差分進(jìn)化算法是一種群體進(jìn)化算法，具有記憶個(gè)體最優(yōu)解和種群內(nèi)信息共享的特點(diǎn)，該算法每次隨機(jī)選擇2個(gè)或多個(gè)群體元素，比較它們性能指標(biāo)的差異，從而確定解在結(jié)果空間中的走向，其本質(zhì)是一種基于實(shí)數(shù)編碼的具有保優(yōu)思想的貪婪遺傳算法。

由于太陽(yáng)系多目標(biāo)探測(cè)軌跡優(yōu)化問(wèn)題的多峰性，單純的使用差分進(jìn)化算法并不能保證得到結(jié)果的最優(yōu)性。因此，本文基于改進(jìn)后的網(wǎng)格搜索算法和差分進(jìn)化算法，提出了一種組合優(yōu)化算法。該算法首先利用改進(jìn)的網(wǎng)格搜索算法作全局搜索，保留較理想的結(jié)果，得到的理想結(jié)果分布在整個(gè)解空間中，然后用差分進(jìn)化算法作局部尋優(yōu)，以改進(jìn)的網(wǎng)格搜索算法得到的結(jié)果收縮搜索空間，將原來(lái)的搜索空間分割成n(改進(jìn)后的網(wǎng)格搜索算法得到的結(jié)果數(shù))份，使用差分進(jìn)化算法逐一搜索各個(gè)子空間，再通過(guò)全局范圍內(nèi)比較得到的各個(gè)極小值，就可以在全局范圍內(nèi)找到一些比較理想的結(jié)果。

3 仿真算例

以2018～2020年太陽(yáng)系多目標(biāo)探測(cè)為例。任務(wù)過(guò)程中先后探測(cè)金星和木星，最后被土星引力場(chǎng)捕獲，形成近地點(diǎn)半徑rp=108950km，偏心率e=0.98的大橢圓軌道。任務(wù)過(guò)程中借助金星，地球，木星的引力，既可以達(dá)到探測(cè)器變軌的目的，節(jié)省發(fā)射能量，又能實(shí)現(xiàn)對(duì)多個(gè)天體的探測(cè)。

該任務(wù)是上述非線性規(guī)劃問(wèn)題的一個(gè)實(shí)例，一共可以分為5個(gè)階段，因此時(shí)間序列的形式為：τ=[t0,T1,T2,T3,T4,T5]，包括地球的出發(fā)時(shí)刻(單位為MJD2000)，每階段轉(zhuǎn)移時(shí)間(單位為d)。優(yōu)化的目標(biāo)函數(shù)為：

(22)

其中，Δv0為地球發(fā)射時(shí)的剩余速度，Δvi(i=1,2,3,4)為探測(cè)器飛行途中飛越大行星所需的速度脈沖，Δv5為木星捕獲所需的速度脈沖。時(shí)間序列的取值范圍如表1所示。

表1 時(shí)間序列參數(shù)取值范圍

發(fā)射能量約束：

Δv0≤4km/s

近地點(diǎn)半徑約束：

rp1>6351.8km

rp2>6351.8km

rp3>6778.1km

rp4>671492km

為了求解方便可以在性能指標(biāo)中引入懲罰項(xiàng)，將有約束的軌跡優(yōu)化問(wèn)題轉(zhuǎn)換成無(wú)約束的軌跡優(yōu)化問(wèn)題。(其中ωi為約束的權(quán)重大小)

1)發(fā)射能量約束

(23)

2)近地點(diǎn)半徑約束

(24)

在仿真過(guò)程中發(fā)現(xiàn)使用上述介紹的改進(jìn)的網(wǎng)格搜索算法，遍歷整個(gè)搜索空間只需很少的時(shí)間，例如同樣使用matlab作為工具，以相同的步長(zhǎng)，列舉法需要大約一周的時(shí)間，而改進(jìn)的網(wǎng)格搜索算法只需要2h，并且得到的理想結(jié)果與列舉法相似(步長(zhǎng)一致)。改進(jìn)的網(wǎng)格搜索算法得到的結(jié)果投影到一個(gè)平面上的結(jié)果如圖3所示，x軸為解空間的分布，y軸為解空間的每一點(diǎn)對(duì)應(yīng)的太陽(yáng)系多目標(biāo)轉(zhuǎn)移軌道所消耗的能量。

圖3 高維解在一個(gè)切平面上的投影

表2 2018～2021年太陽(yáng)系多目標(biāo)探測(cè)行星際飛越軌道

從表2中可以看出前3組數(shù)據(jù)所對(duì)應(yīng)的太陽(yáng)系多目標(biāo)探測(cè)行星際飛越軌道雖然消耗的燃料很少，但轉(zhuǎn)移時(shí)間卻達(dá)到了19年左右，這是工程任務(wù)所無(wú)法接受的。后面4組數(shù)據(jù)中盡管燃料消耗相對(duì)多一點(diǎn)，但任務(wù)時(shí)間均縮短一半以上。探索土星的國(guó)際卡西尼太空任務(wù)經(jīng)過(guò)大約7年的飛行時(shí)間最終抵達(dá)土星[13]，表2中后4組數(shù)據(jù)所對(duì)應(yīng)的太陽(yáng)系多目標(biāo)探測(cè)行星際飛越軌道的轉(zhuǎn)移時(shí)間大都在7年左右，并且所需的燃料也在實(shí)際工程約束范圍之內(nèi)，因此，后4組軌道均是可實(shí)際應(yīng)用的。

表3給出了表2中第7組數(shù)據(jù)所對(duì)應(yīng)的行星際飛越軌道的詳細(xì)參數(shù)，并且表3中的速度單位都是(km/s)，從表3中可以看出，上述結(jié)果滿足發(fā)射能量約束和近地點(diǎn)半徑約束，這說(shuō)明以第7組數(shù)據(jù)對(duì)應(yīng)的軌道參數(shù)設(shè)計(jì)的太陽(yáng)系多目標(biāo)轉(zhuǎn)移軌道有效。使用上述參數(shù)進(jìn)行仿真，最終得到太陽(yáng)系多目標(biāo)探測(cè)行星際飛越軌道的三維圖像如圖4所示。圖5為表2中第一組數(shù)據(jù)所對(duì)應(yīng)的太陽(yáng)系多目標(biāo)探測(cè)行星際飛越軌道的三維圖像。圖4和圖5再次說(shuō)明了本文所提出算法的有效性。

表3 發(fā)射窗口為2018年9月18日的太陽(yáng)系多目標(biāo)探測(cè)行星際飛越軌道的詳細(xì)參數(shù)

圖4 發(fā)射窗口為2018/09/12的太陽(yáng)系多目標(biāo)探測(cè)行星際飛越軌道的三維圖像

圖5 發(fā)射窗口為2018/07/25的太陽(yáng)系多目標(biāo)探測(cè)行星際飛越軌道的三維圖像

4 結(jié)束語(yǔ)

行星際探測(cè)能耗高，持續(xù)時(shí)間長(zhǎng)，通過(guò)一次任務(wù)探測(cè)太陽(yáng)系多顆行星，不僅可以節(jié)省發(fā)射成本，同時(shí)結(jié)合新技術(shù)驗(yàn)證，可以推動(dòng)航天科技的發(fā)展。根據(jù)我國(guó)的實(shí)際情況，考慮到我國(guó)可能于2020年左右進(jìn)行的太陽(yáng)系多目標(biāo)深空探測(cè)任務(wù)，本文對(duì)其軌道優(yōu)化設(shè)計(jì)進(jìn)行了預(yù)先研究。首先根據(jù)太陽(yáng)系多目標(biāo)探測(cè)行星際飛越的基本理論，建立了太陽(yáng)系多目標(biāo)探測(cè)的軌道優(yōu)化模型，然后將軌跡優(yōu)化問(wèn)題轉(zhuǎn)換成非線性規(guī)劃問(wèn)題,并且針對(duì)太陽(yáng)系多目標(biāo)探測(cè)軌跡優(yōu)化設(shè)計(jì)問(wèn)題的復(fù)雜性和多峰性，提出了一種融合了改進(jìn)的網(wǎng)格搜索算法和差分進(jìn)化算法的組合優(yōu)化算法。該算法將整個(gè)階段分成多個(gè)子設(shè)計(jì)段，從而可以對(duì)每個(gè)子設(shè)計(jì)段進(jìn)行單獨(dú)設(shè)計(jì)，減少了全局優(yōu)化所需要的計(jì)算量并可以有效的收斂到全局范圍內(nèi)的最優(yōu)值。仿真結(jié)果所給出的多組2018～2020年最終環(huán)繞土星，中途探測(cè)金星，木星的轉(zhuǎn)移軌道參數(shù)可以為我國(guó)未來(lái)的太陽(yáng)系多目標(biāo)探測(cè)任務(wù)設(shè)計(jì)提供有價(jià)值的參考。

參 考 文 獻(xiàn)

[1] Broucke R A. The Celestial Mechanics of Gravity Assist [C].AIAA/AAS Astrodynamics Conference, Minnesota, Aug 15-17, 1998.

[2] Louis A Damario, Larrye Bright. Galileo Trajectory Design [J]. Space Science Reviews, 1992, 60(4): 23-78.

[3] Lawden D F. Optimal Trajectories for Space Navigation [M]. Butterworths, London, 1963.

[4] Lion P M, Handelsman M. The Primer Vector on Fixed-time Impulsive Trajectories [J]. AIAA Journal, 1968, 1(6): 127-132.

[5] Petropoulos A E, Longuski J M. Shape-based Algorithm for Automated Design of Low-Thrust Gravity-Assist Trajectories [J]. Journal of Spacecraft and Rockets, 2004, 5(41): 787-796.

[6] 張旭輝, 劉竹生.火星探測(cè)無(wú)動(dòng)力借力飛行軌道研究[J]. 宇航學(xué)報(bào), 2008, 29(6): 1739-1746.(Zhang Xu-hui, Liu Zhu-sheng. Unpowered Swingby Flight Orbit Design for Mars Exploration[J].Journal of Astronautics, 2008, 29(6): 1739-1746.)

[7] 喬棟, 崔平遠(yuǎn), 徐瑞.星際探測(cè)借力飛行軌道的混合設(shè)計(jì)方法研究[J].宇航學(xué)報(bào), 2010, 31(3): 655-661. (Qiao Dong, Cui Ping-yuan, Xu Rui. Research on Hybrid Design Method of Gravity Assists Trajectory for Interplanetary Mission [J]. Journal of Astronautics, 2010, 31(3):655-661.)

[8] Betts J T. Survey of Numerical Methods for Trajectory Optimization [J]. Journal of Guidance, Control and Dynamics, 1998, 2(21): 193-207.

[9] Cage P G. Interplanetary Trajectory Optimization Using a Genetic Algorithm [C]. AIAA/AAS Astrodynamics Conference, Scottsdale, Aug 1-3, 1994.

[10] Bessette C R, Spencer D B. Identifying Optimal Interplanetary Trajectories Through a Genetic Approach[C].Astrodynamics Specialist Conference and Exhibit, Keystone, 2006.

[11] Bessette C R, Spencer D B. Optimal Space Trajectory Design: A Heuristic-Based Approach [J]. Advances in the Astronautical Sciences, 2006, 6(124): 1611-1628.

[12] V M Izzo, D Myatt, D R Nasuto, S J Bishop. Search Space Pruning and Global Optimisation of Multiple Gravity Assist Spacecraft Trajectories[J]. Journal of Global Optimization, 2007, 2(38): 283-296.

[13] Peralta F, Flanagan S. Cassini Interplanetary Trajectory Mission [J]. Control Engineering Practice, 1995, 3(11): 1603-1610.