基于形狀的行星際小推力轉移軌道初始設計方法

尚海濱,崔平遠,喬 棟

(北京理工大學深空探測技術研究所,北京100081)

0 引 言

在星際探測任務中,采用小推力發動機可有效提高探測器的有效載荷,增加探測任務科學回報[1]。星際小推力轉移軌道的設計通常分為初始設計和精確設計兩步,初始設計的目的是為精確設計提供可行的設計初值[2]。與脈沖轉移軌道不同,小推力轉移弧段不能用傳統的圓錐曲線來刻畫,許多脈沖軌道設計中的理論與方法不再適用,這給行星際小推力轉移軌道的設計帶來了困難,尋求一種快速有效的初始設計方法成為目前研究的熱點[3]。

基于形狀的星際小推力軌道設計理論是目前較為有效的一類方法。基于形狀逼近思想,Petropoulos提出了一種基于正弦指數曲線逼近的設計方法[4-5],由于正弦指數曲線不能滿足始末端速度邊界條件,該方法只能用于設計發射及射入能量不為零的軌道。針對此,Wall又提出了一種基于六次逆多項式曲線逼近的設計方法[6],該函數曲線是通過對最優轉移軌道數據擬合得到的,由于具有七個待定參數,能夠同時滿足始末端位置速度邊界條件以及飛行時間約束,因此可對任務類型轉移軌道進行設計。然而,由于沒有考慮轉移軌道的實際動力學約束要求,導致該方法經常出現遺漏可行軌道的情況。

本文針對以上問題,對逆多項式逼近軌道的動力學可行性進行了深入分析,推導了關鍵系數的可行范圍;然后,基于改進的逆多項式逼近理論,建立了一種簡單有效的行星際小推力轉移軌道初始設計模型,同時,為了獲得全局最優轉移機會,采用了微分進化算法對搜索參數進行尋優計算。最后,以地球-火星和地球-金星小推力轉移為例,對所提算法進行了仿真驗證。

1 逆多項式逼近策略

在日心引力場中,采用極坐標系描述的探測器運動方程為:

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

由于推力加速度項的存在,方程(1)不存在閉環解析解。為便于求解,Petropoulos和Wall先后提出兩種函數來描述探測器的運動軌跡:

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

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

由文獻[6]可知,在給定始末端邊界條件的情況下,要獲得解析的轉移軌道,只存在a3一個未知量,該參數可根據時間約束條件通過迭代求解。然而,文獻[6]所給方法在迭代過程中經常會出現˙θ在某一相角范圍內有˙θ＜0,導致所求得的轉移時間t出現虛值,無法得到滿足條件的a3,從而可能丟失可行的轉移軌道。

2 動力學可行性分析

對于可行的轉移軌道,˙θ在轉移過程中的變化趨勢應該保持不變。如果始終有˙θ≥0,則對應順行軌道;如果始終有˙θ＜0,則對應逆行軌道。本文針對順行軌道展開討論,根據式(4),由于 μr4＞0,若要使˙θ≥0,必須滿足如下條件:

為便于描述,假設探測器初始相角 θ1=0,末端相角為0≤θ2≤2π,則探測器轉過的相角為 θf=2nrevπ+θ2(nrev為探測器在轉移過程中繞太陽轉過的完整圈數)。分別令

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

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

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

根據式(5),f是相角θ和逆多項式系數a3的函數,可以改寫成:

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

可見,系數 a3的可行范圍主要是由函數g(θ)的取值確定的。考察g(θ)發現,g(θ)具有如下性質:

進一步分析可知函數g(θ)實根個數n隨轉移相角θf的變化特點為:

根據以上討論,針對不同的轉移相角,為了保證逆多項式函數(3)描述的軌道滿足動力學要求,系數a3實際的可行范圍可選定為:

圖1 函數 g(θ)與θf的關系Fig.1 Relation between g(θ)and θf

其中:θ′為g(θ)在區間[0,θf]的實根,xmin和 xmax分別為

其中:x min可以表示為

其中:xmin可以表示為

根據確定的可行范圍,圖2給出了逆多項式曲線描述軌道的飛行時間隨a3的變化曲線。

由圖2可以看出,轉移時間 ΔT隨系數a3具有良好的單調遞減關系。因此,對于給定的發射時間和到達時間,采用逆多項式逼近軌道時只存在無解和單解兩種情況,為保證計算的高效性,避免搜索過程中的冗余計算,可首先通過式(20)判斷是否有解:

圖2 飛行時間隨a3的變化曲線Fig.2 Relation betweenΔT and a3

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

根據確定的可行范圍,如果E＞0,則不存在可行軌道,如果 E≤0,則存在可行軌道,采用Brent算法可以簡單的尋找到滿足時間飛行時間的系數 a3,并且在迭代過程中不會出現˙θ＜0的情況。

3 軌道設計方法

在上述理論研究基礎上,本文發展了一種基于逆多項式逼近的行星際軌道初始化設計方法。

3.1 軌道初始設計模型

基于逆多項式逼近理論,燃料最省星際小推力轉移軌道設計問題可以描述如下:尋求最優的發射時間 t1、到達目標星體時間t2和轉移圈數 nrev,使性能指標(21)達到最小。

這里,尋優參數可以表示為:

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

此外,由于逆多項式理論是在加速度大小可調且無界假設下推導的。為了滿足實際的推進器模型,轉移過程中加速度大小必須滿足如下約束條件:

其中:Fmax為推進器能夠提供的最大加速度。

方程(24)描述的約束可看作是尋優問題的路徑約束,為無限維的連續約束,為便于計算,可簡化為如下一維形式:

綜上所述,基于逆多項式逼近理論,星際小推力轉移軌道設計問題可以轉化為一個帶有整數變量nrev的混合非線性規劃問題。

3.2 微分進化算法

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

值得注意的是,由于微分進化沒有約束處理機制,為此,這里采用模擬退火罰函數策略將混合非線性規劃問題轉化為無約束優化問題,構造新的性能指標為:

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

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

其中:ε為提前定義的收斂參數,本文中選取ε=10-5,在每一步迭代過程中,優化程序都將記錄所有群體中的最優性能指標,當迭代次數G＞15時,優化程序根據式(26)計算性能指標的改善情況,如果滿足式(28),則終止優化程序,否則繼續迭代。

4 數值仿真與分析

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

4.1 固定時間地球火星交會任務

在算例1中,假定探測器2015年1月1日從地球發射,飛行500天后與火星成功交會,探測器轉移過程中繞太陽圈數nrev=1。針對以上問題,采用逆多項式方法求解的結果如表1所示。

表1 地球-火星轉移軌道設計結果Table 1 Design result of Earth-Mars transfer trajectory

針對這一問題,采用文獻[6]所給方法進行求解時,無法得到可行的解。這是因為,由于文獻[6]中沒有給出搜索參數的可行范圍,只簡單的令其初值為零進行計算,導致搜索過程中參數a3跳出其可行范圍,得到的飛行時間為虛值,不滿足要求。由此可見,采用改進后的擬多項式逼近理論,可以有效避免傳統算法的失效情況,提高逆多項式逼近理論的魯棒性。圖3給出了采用改進逆多項式逼近理論求得的轉移軌道示意圖。

圖3 地球-火星的轉移軌道Fig.3 Earth-Mars transfer trajectory

4.2 非固定時間地球金星交會任務

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

針對該問題,微分進化算法迭代81步達到了收斂,即滿足式(27)給出的終止條件。圖4給出了地球-金星轉移軌道的性能指標隨發射時間和到達時間的分布情況。

表2 地球-金星探測任務約束Table2 Mission constraints of Earth-Venus transfer

表3 微分進化算法參數設置Table3 Parameters set of DE algorithm

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

為驗證設計結果的有效性,作者針對該問題,采用高斯偽光譜配點法進行了精確優化[8]。優化過程中,探測任務的約束條件同表2,探測器采用平面運動方程,假設推力大小和方向都可調,性能指標取為式(21)。高斯偽光譜配點算法的離散點數取為50,則優化問題有302個尋優參數,250個約束條件。兩種求解方式的計算結果如表4所示。

表4 計算結果比較Table 4 Comparison between two algorithms

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

圖5 地球-金星轉移軌道示意圖Fig.5 Earth-Venus transfer trajectory

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

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

5 結 論

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

[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.