尚海濱,崔平遠(yuǎn),喬 棟

(北京理工大學(xué)深空探測(cè)技術(shù)研究所,北京100081)

0 引 言

在星際探測(cè)任務(wù)中,采用小推力發(fā)動(dòng)機(jī)可有效提高探測(cè)器的有效載荷,增加探測(cè)任務(wù)科學(xué)回報(bào)[1]。星際小推力轉(zhuǎn)移軌道的設(shè)計(jì)通常分為初始設(shè)計(jì)和精確設(shè)計(jì)兩步,初始設(shè)計(jì)的目的是為精確設(shè)計(jì)提供可行的設(shè)計(jì)初值[2]。與脈沖轉(zhuǎn)移軌道不同,小推力轉(zhuǎn)移弧段不能用傳統(tǒng)的圓錐曲線(xiàn)來(lái)刻畫(huà),許多脈沖軌道設(shè)計(jì)中的理論與方法不再適用,這給行星際小推力轉(zhuǎn)移軌道的設(shè)計(jì)帶來(lái)了困難,尋求一種快速有效的初始設(shè)計(jì)方法成為目前研究的熱點(diǎn)[3]。

基于形狀的星際小推力軌道設(shè)計(jì)理論是目前較為有效的一類(lèi)方法?；谛螤畋平枷?Petropoulos提出了一種基于正弦指數(shù)曲線(xiàn)逼近的設(shè)計(jì)方法[4-5],由于正弦指數(shù)曲線(xiàn)不能滿(mǎn)足始末端速度邊界條件,該方法只能用于設(shè)計(jì)發(fā)射及射入能量不為零的軌道。針對(duì)此,Wall又提出了一種基于六次逆多項(xiàng)式曲線(xiàn)逼近的設(shè)計(jì)方法[6],該函數(shù)曲線(xiàn)是通過(guò)對(duì)最優(yōu)轉(zhuǎn)移軌道數(shù)據(jù)擬合得到的,由于具有七個(gè)待定參數(shù),能夠同時(shí)滿(mǎn)足始末端位置速度邊界條件以及飛行時(shí)間約束,因此可對(duì)任務(wù)類(lèi)型轉(zhuǎn)移軌道進(jìn)行設(shè)計(jì)。然而,由于沒(méi)有考慮轉(zhuǎn)移軌道的實(shí)際動(dòng)力學(xué)約束要求,導(dǎo)致該方法經(jīng)常出現(xiàn)遺漏可行軌道的情況。

本文針對(duì)以上問(wèn)題,對(duì)逆多項(xiàng)式逼近軌道的動(dòng)力學(xué)可行性進(jìn)行了深入分析,推導(dǎo)了關(guān)鍵系數(shù)的可行范圍;然后,基于改進(jìn)的逆多項(xiàng)式逼近理論,建立了一種簡(jiǎn)單有效的行星際小推力轉(zhuǎn)移軌道初始設(shè)計(jì)模型,同時(shí),為了獲得全局最優(yōu)轉(zhuǎn)移機(jī)會(huì),采用了微分進(jìn)化算法對(duì)搜索參數(shù)進(jìn)行尋優(yōu)計(jì)算。最后,以地球-火星和地球-金星小推力轉(zhuǎn)移為例,對(duì)所提算法進(jìn)行了仿真驗(yàn)證。

1 逆多項(xiàng)式逼近策略

在日心引力場(chǎng)中,采用極坐標(biāo)系描述的探測(cè)器運(yùn)動(dòng)方程為:

其中:r為探測(cè)器軌道的矢徑大小,θ為在日心慣性系中的相角,μ為太陽(yáng)的引力常數(shù),F為推力加速度大小,α為推力方向角。

由于推力加速度項(xiàng)的存在,方程(1)不存在閉環(huán)解析解。為便于求解,Petropoulos和Wall先后提出兩種函數(shù)來(lái)描述探測(cè)器的運(yùn)動(dòng)軌跡:

理論推導(dǎo)表明,方程(3)提出的逆多項(xiàng)式曲線(xiàn)可以用于包括交會(huì)軌道在內(nèi)的任務(wù)類(lèi)型的軌道轉(zhuǎn)移問(wèn)題,相比方程(2)具有更廣的應(yīng)用范圍。另外,由于方程(3)是基于最優(yōu)轉(zhuǎn)移軌道數(shù)據(jù)擬合得到的函數(shù),因此具有更優(yōu)的軌道轉(zhuǎn)移特性。Bradley基于方程(3),通過(guò)使探測(cè)器軌道分別滿(mǎn)足始末端的位置速度邊界條件,得到了解析的逆多項(xiàng)式曲線(xiàn)的參數(shù)a0～a2及a4～ a6。

進(jìn)一步,假定轉(zhuǎn)移過(guò)程中推力矢量始終沿飛行路徑角方向或其反方向,即令 α=γ+mπ(m為0或1),可以得到相角 θ變化率為:

由文獻(xiàn)[6]可知,在給定始末端邊界條件的情況下,要獲得解析的轉(zhuǎn)移軌道,只存在a3一個(gè)未知量,該參數(shù)可根據(jù)時(shí)間約束條件通過(guò)迭代求解。然而,文獻(xiàn)[6]所給方法在迭代過(guò)程中經(jīng)常會(huì)出現(xiàn)˙θ在某一相角范圍內(nèi)有˙θ＜0,導(dǎo)致所求得的轉(zhuǎn)移時(shí)間t出現(xiàn)虛值,無(wú)法得到滿(mǎn)足條件的a3,從而可能丟失可行的轉(zhuǎn)移軌道。

2 動(dòng)力學(xué)可行性分析

對(duì)于可行的轉(zhuǎn)移軌道,˙θ在轉(zhuǎn)移過(guò)程中的變化趨勢(shì)應(yīng)該保持不變。如果始終有˙θ≥0,則對(duì)應(yīng)順行軌道;如果始終有˙θ＜0,則對(duì)應(yīng)逆行軌道。本文針對(duì)順行軌道展開(kāi)討論,根據(jù)式(4),由于 μr4＞0,若要使˙θ≥0,必須滿(mǎn)足如下條件:

為便于描述,假設(shè)探測(cè)器初始相角 θ1=0,末端相角為0≤θ2≤2π,則探測(cè)器轉(zhuǎn)過(guò)的相角為 θf(wàn)=2nrevπ+θ2(nrev為探測(cè)器在轉(zhuǎn)移過(guò)程中繞太陽(yáng)轉(zhuǎn)過(guò)的完整圈數(shù))。分別令

則 a′4、a′5和 a′6分別可以表示為 :

其中:輔助變量 w′1、w′2和 w′3可以表示為

其中:r2為探測(cè)器軌道末端矢徑大小,γ2為軌道末端飛行路徑角。

根據(jù)式(5),f是相角θ和逆多項(xiàng)式系數(shù)a3的函數(shù),可以改寫(xiě)成:

其中:g(θ)和 H(θ)可以分別表示為

可見(jiàn),系數(shù) a3的可行范圍主要是由函數(shù)g(θ)的取值確定的。考察g(θ)發(fā)現(xiàn),g(θ)具有如下性質(zhì):

進(jìn)一步分析可知函數(shù)g(θ)實(shí)根個(gè)數(shù)n隨轉(zhuǎn)移相角θf(wàn)的變化特點(diǎn)為:

根據(jù)以上討論,針對(duì)不同的轉(zhuǎn)移相角,為了保證逆多項(xiàng)式函數(shù)(3)描述的軌道滿(mǎn)足動(dòng)力學(xué)要求,系數(shù)a3實(shí)際的可行范圍可選定為:

圖1 函數(shù) g(θ)與θf(wàn)的關(guān)系Fig.1 Relation between g(θ)and θf(wàn)

其中:θ′為g(θ)在區(qū)間[0,θf(wàn)]的實(shí)根,xmin和 xmax分別為

其中:x min可以表示為

其中:xmin可以表示為

根據(jù)確定的可行范圍,圖2給出了逆多項(xiàng)式曲線(xiàn)描述軌道的飛行時(shí)間隨a3的變化曲線(xiàn)。

由圖2可以看出,轉(zhuǎn)移時(shí)間 ΔT隨系數(shù)a3具有良好的單調(diào)遞減關(guān)系。因此,對(duì)于給定的發(fā)射時(shí)間和到達(dá)時(shí)間,采用逆多項(xiàng)式逼近軌道時(shí)只存在無(wú)解和單解兩種情況,為保證計(jì)算的高效性,避免搜索過(guò)程中的冗余計(jì)算,可首先通過(guò)式(20)判斷是否有解:

圖2 飛行時(shí)間隨a3的變化曲線(xiàn)Fig.2 Relation betweenΔT and a3

其中:amin3和amax3分別為系數(shù)a3的取值上下限,δ為避免ΔT奇異的常值參數(shù),這里取為δ=10-4。

根據(jù)確定的可行范圍,如果E＞0,則不存在可行軌道,如果 E≤0,則存在可行軌道,采用Brent算法可以簡(jiǎn)單的尋找到滿(mǎn)足時(shí)間飛行時(shí)間的系數(shù) a3,并且在迭代過(guò)程中不會(huì)出現(xiàn)˙θ＜0的情況。

3 軌道設(shè)計(jì)方法

在上述理論研究基礎(chǔ)上,本文發(fā)展了一種基于逆多項(xiàng)式逼近的行星際軌道初始化設(shè)計(jì)方法。

3.1 軌道初始設(shè)計(jì)模型

基于逆多項(xiàng)式逼近理論,燃料最省星際小推力轉(zhuǎn)移軌道設(shè)計(jì)問(wèn)題可以描述如下:尋求最優(yōu)的發(fā)射時(shí)間 t1、到達(dá)目標(biāo)星體時(shí)間t2和轉(zhuǎn)移圈數(shù) nrev,使性能指標(biāo)(21)達(dá)到最小。

這里,尋優(yōu)參數(shù)可以表示為:

F(θ)為轉(zhuǎn)移過(guò)程中的加速度,可以表示為:

此外,由于逆多項(xiàng)式理論是在加速度大小可調(diào)且無(wú)界假設(shè)下推導(dǎo)的。為了滿(mǎn)足實(shí)際的推進(jìn)器模型,轉(zhuǎn)移過(guò)程中加速度大小必須滿(mǎn)足如下約束條件:

其中:Fmax為推進(jìn)器能夠提供的最大加速度。

方程(24)描述的約束可看作是尋優(yōu)問(wèn)題的路徑約束,為無(wú)限維的連續(xù)約束,為便于計(jì)算,可簡(jiǎn)化為如下一維形式:

綜上所述,基于逆多項(xiàng)式逼近理論,星際小推力轉(zhuǎn)移軌道設(shè)計(jì)問(wèn)題可以轉(zhuǎn)化為一個(gè)帶有整數(shù)變量nrev的混合非線(xiàn)性規(guī)劃問(wèn)題。

3.2 微分進(jìn)化算法

本文采用微分進(jìn)化(DE)算法對(duì)歸結(jié)的混合非線(xiàn)性規(guī)劃問(wèn)題進(jìn)行求解。微分進(jìn)化是由Storn和Price提出的一種新穎的基于達(dá)爾文進(jìn)化理論的隨機(jī)搜索算法[7]。與遺傳算法相比,微分進(jìn)化是基于雙精度實(shí)數(shù)編碼,并且可以選擇采用多種不同變異和選擇策略,因此,具有更強(qiáng)的自適應(yīng)能力,DE算法具體步驟參考文獻(xiàn)[7]。

值得注意的是,由于微分進(jìn)化沒(méi)有約束處理機(jī)制,為此,這里采用模擬退火罰函數(shù)策略將混合非線(xiàn)性規(guī)劃問(wèn)題轉(zhuǎn)化為無(wú)約束優(yōu)化問(wèn)題,構(gòu)造新的性能指標(biāo)為:

其中:T為退火因子,隨著迭代進(jìn)行由初始溫度 T0變化到凝結(jié)溫度 Tf,其遞推公式如下:

采用微分進(jìn)化求解小推力轉(zhuǎn)移軌道問(wèn)題的另一個(gè)關(guān)鍵問(wèn)題就是迭代的終止條件,如果微分進(jìn)化已經(jīng)尋找到最優(yōu)解或者已經(jīng)對(duì)當(dāng)前解不能有所改善,則都應(yīng)該對(duì)優(yōu)化程序進(jìn)行終止。微分進(jìn)化算法作為一種隨機(jī)非確定算法,每一次運(yùn)行結(jié)果并不相同,因此,必須針對(duì)以上兩種情況給出合理的終止條件。本文中作者采用如下終止策略:

其中:ε為提前定義的收斂參數(shù),本文中選取ε=10-5,在每一步迭代過(guò)程中,優(yōu)化程序都將記錄所有群體中的最優(yōu)性能指標(biāo),當(dāng)?shù)螖?shù)G＞15時(shí),優(yōu)化程序根據(jù)式(26)計(jì)算性能指標(biāo)的改善情況,如果滿(mǎn)足式(28),則終止優(yōu)化程序,否則繼續(xù)迭代。

4 數(shù)值仿真與分析

為了充分驗(yàn)證作者對(duì)逆多項(xiàng)式理論改進(jìn)以及所提轉(zhuǎn)移軌道初始設(shè)計(jì)方法的有效性,這里將給出兩個(gè)仿真算例。算例1將以固定時(shí)間地球-火星的轉(zhuǎn)移為例,驗(yàn)證作者所提改進(jìn)擬多項(xiàng)式理論的正確性,算例2將以非固定時(shí)間地球-水星的轉(zhuǎn)移為例,驗(yàn)證所提軌道初始設(shè)計(jì)方法的有效性。

4.1 固定時(shí)間地球火星交會(huì)任務(wù)

在算例1中,假定探測(cè)器2015年1月1日從地球發(fā)射,飛行500天后與火星成功交會(huì),探測(cè)器轉(zhuǎn)移過(guò)程中繞太陽(yáng)圈數(shù)nrev=1。針對(duì)以上問(wèn)題,采用逆多項(xiàng)式方法求解的結(jié)果如表1所示。

表1 地球-火星轉(zhuǎn)移軌道設(shè)計(jì)結(jié)果Table 1 Design result of Earth-Mars transfer trajectory

針對(duì)這一問(wèn)題,采用文獻(xiàn)[6]所給方法進(jìn)行求解時(shí),無(wú)法得到可行的解。這是因?yàn)?由于文獻(xiàn)[6]中沒(méi)有給出搜索參數(shù)的可行范圍,只簡(jiǎn)單的令其初值為零進(jìn)行計(jì)算,導(dǎo)致搜索過(guò)程中參數(shù)a3跳出其可行范圍,得到的飛行時(shí)間為虛值,不滿(mǎn)足要求。由此可見(jiàn),采用改進(jìn)后的擬多項(xiàng)式逼近理論,可以有效避免傳統(tǒng)算法的失效情況,提高逆多項(xiàng)式逼近理論的魯棒性。圖3給出了采用改進(jìn)逆多項(xiàng)式逼近理論求得的轉(zhuǎn)移軌道示意圖。

圖3 地球-火星的轉(zhuǎn)移軌道Fig.3 Earth-Mars transfer trajectory

4.2 非固定時(shí)間地球金星交會(huì)任務(wù)

在算例2中,將以非固定時(shí)間地球-金星的燃料最省交會(huì)任務(wù)為例對(duì)所給的軌道初始設(shè)計(jì)方法進(jìn)行驗(yàn)證。金星探測(cè)任務(wù)的約束條件如表2所示。微分進(jìn)化算法參數(shù)設(shè)置如表3所示。

針對(duì)該問(wèn)題,微分進(jìn)化算法迭代81步達(dá)到了收斂,即滿(mǎn)足式(27)給出的終止條件。圖4給出了地球-金星轉(zhuǎn)移軌道的性能指標(biāo)隨發(fā)射時(shí)間和到達(dá)時(shí)間的分布情況。

表2 地球-金星探測(cè)任務(wù)約束Table2 Mission constraints of Earth-Venus transfer

表3 微分進(jìn)化算法參數(shù)設(shè)置Table3 Parameters set of DE algorithm

圖4 地球-金星轉(zhuǎn)移軌道解空間分布Fig.4 Solution space for Earth-Venus transfer

為驗(yàn)證設(shè)計(jì)結(jié)果的有效性,作者針對(duì)該問(wèn)題,采用高斯偽光譜配點(diǎn)法進(jìn)行了精確優(yōu)化[8]。優(yōu)化過(guò)程中,探測(cè)任務(wù)的約束條件同表2,探測(cè)器采用平面運(yùn)動(dòng)方程,假設(shè)推力大小和方向都可調(diào),性能指標(biāo)取為式(21)。高斯偽光譜配點(diǎn)算法的離散點(diǎn)數(shù)取為50,則優(yōu)化問(wèn)題有302個(gè)尋優(yōu)參數(shù),250個(gè)約束條件。兩種求解方式的計(jì)算結(jié)果如表4所示。

表4 計(jì)算結(jié)果比較Table 4 Comparison between two algorithms

由表4可以看出:(1)采用作者所提的基于微分進(jìn)化的初始設(shè)計(jì)方法,設(shè)計(jì)結(jié)果與高斯偽光譜方法的結(jié)果較為接近,其中從地球發(fā)射時(shí)間與最優(yōu)解僅相差6天,飛行時(shí)間相差8.733天,所需要的總的速度脈沖相差約1.24%,同時(shí),初始設(shè)計(jì)方法求得軌道的最大推力約束為0.25N,完全滿(mǎn)足任務(wù)約束,這充分表明作者所提的設(shè)計(jì)方法是有效的,可以獲得滿(mǎn)足任務(wù)約束的次優(yōu)轉(zhuǎn)移軌道;(2)從另一個(gè)角度,既然初始設(shè)計(jì)方法的設(shè)計(jì)結(jié)果與高斯偽光譜方法得到的最優(yōu)解差別很小,也表明了初始設(shè)計(jì)方法獲得的次優(yōu)解可以為精確的轉(zhuǎn)移軌道設(shè)計(jì)提供可行有效的設(shè)計(jì)初值。由此可以看出,作者所給的行星際小推力轉(zhuǎn)移軌道初始設(shè)計(jì)方法是行之有效的,采用其獲得的初值猜測(cè)可以大大降低行星際小推力轉(zhuǎn)移軌道的設(shè)計(jì)難度。為了充分比較作者所提算法與高斯偽光譜所獲的最優(yōu)解,圖5和圖6分別給出了兩種算法獲得的轉(zhuǎn)移軌道示意圖和推力大小時(shí)間歷程。

圖5 地球-金星轉(zhuǎn)移軌道示意圖Fig.5 Earth-Venus transfer trajectory

由圖5和圖6可以看出,兩種算法獲得的轉(zhuǎn)移軌道和推力大小時(shí)間歷程盡管存在一定差別,但總體趨勢(shì)是一致的,這進(jìn)一步驗(yàn)證了作者所給設(shè)計(jì)方法的正確性和有效性。造成兩種算法所獲軌道差別的主要原因是:在逆多項(xiàng)式逼近算法中,假定推力方向始終沿速度方向或其反方向,而在高斯偽光譜方法中,認(rèn)為其推力方向是自由可調(diào)的。

圖6 推力大小時(shí)間歷程Fig.6 Time history of thrust magnitude

5 結(jié) 論

針對(duì)行星際小推力轉(zhuǎn)移軌道初始設(shè)計(jì)問(wèn)題,提出了一種基于改進(jìn)逆多項(xiàng)式逼近的設(shè)計(jì)方法。該方法具有以下三個(gè)優(yōu)點(diǎn):(1)推導(dǎo)了關(guān)鍵系數(shù)的動(dòng)力學(xué)可行范圍,相比傳統(tǒng)擬多項(xiàng)逼近算法具有更高的魯棒性;(2)所提設(shè)計(jì)方法的模型簡(jiǎn)單,設(shè)計(jì)參數(shù)少;(3)采用結(jié)合微分進(jìn)化與模擬退火的尋優(yōu)算法,可以快速有效獲得逆多項(xiàng)式逼近下的全局最優(yōu)解。作者所提初始設(shè)計(jì)方法可為精確設(shè)計(jì)提供行之有效的初值,從而有效提高行星際小推力轉(zhuǎn)移軌道的設(shè)計(jì)效率,降低其設(shè)計(jì)難度。本文的研究可以為行星際小推力探測(cè)任務(wù)的設(shè)計(jì)與規(guī)劃提供有價(jià)值的參考。

[1] Bertrand R,Bernussou J,Geffroy S.Electric transfer optimization for mars sample return mission[J].Acta Astronautica,2001,48(5-12):651-660.

[2] Strange N J,Longuski JM.Graphical method for gravity-assist trajectory design[J].Journal of Spacecraft and Rockets,2002,39(1):9-16.

[3] Pascale D P,Vasile M.Preliminary design of low-thrust multiplegravity-assist trajectories[J].Journal of Spacecraft and Rockets,2006,43(5):1065-1076.

[4] Petropoulos A E,Longuski JM.Shape-based analytic representations of low-Thrust trajectories for gravity-assist applications[C].AAS/AIAA Astrodynamics Specialist Conference,16-19August 1999:1-19.

[5] Petropoulos A E,Longuski JM.Shape-based algorithm for automated design of low-thrust,gravity-assist trajectories[J].Journal of Spacecraft and Rockets,2004,41(5):87-796.

[6] Wall B J,Conway B A.Shape-based approach to low-thrust rendezvous trajectory design[J].Journal of Guidance,Control,and Dynamics,2009,32(1):95-101.

[7] Storn R,Price K.Differential evolution:Asimple and efficient adaptive schemefor global optimization over continuous spaces[R].Technical Report TR-95-012,Berkeley,USA,1995.

[8] Benson D.A gauss pseudospectral transcription for optimal control[D].Department of Aeronautics and Astronautics,MIT,2004.