混合推進最省燃料軌道設計方法

張 晨，趙育善

(北京航空航天大學宇航學院，北京100191)

0 引言

在深空探測任務中，探測器使用的發動機可分為化學發動機和電推進發動機兩大類。化學發動機產生的推力大但耗費燃料。電推進發動機因其高比沖的特性能夠極大地節省燃料消耗，但其相對較小的推力幅值導致較長的任務時間。Mingotti等［1］提出的混合推進方式能夠兼容以上兩種發動機各自的優點，這種方法假定衛星上面級提供一次脈沖后與衛星分離，之后衛星攜帶的電推進發動機開機并最終進入目標軌道。混合推進方式相比單純使用化學發動機，能夠大幅降低燃料消耗，還有操控靈活，對發射窗口不敏感等優點;且相比單純使用電推進發動機，混合推進方式提供了更短的任務時間。

近十多年來，國內外對電推進軌道優化技術進行了豐富的研究。任遠等使用混合法設計了地火電推進轉移軌道［2］，避免直接求解兩點邊值問題，從而提高了迭代的收斂性。Jiang等［3］研究了同倫方法以及開關檢測技術在最省燃料問題當中的應用，并獲得了很高的收斂效率和打靶精度。Guo等［4］將同倫方法和粒子群算法結合起來研究了引力輔助在電推進軌道優化方面的應用，獲得了更高的剩余燃料比。黃鎬等［5］將協態歸一化方法和曲線擬合技術相結合，研究了最優電推進逃逸軌道設計問題，降低了多圈問題協態初值猜測的敏感性。上述研究探索了二體模型下的電推進軌道設計方法，近些年來針對深空探測任務提出的新需求，基于三體以及四體模型的電推進軌道設計問題在國內外得到了越來越多的關注。Ozimek等［6］研究了使用變比沖發動機的三體電推進軌道優化問題，基于數值方面的考慮，其推力和比沖沒有設定上下限。Caillau等［7］基于龐特里亞金極大值原理和同倫方法，探討了三體模型下使用定比沖發動機的最省燃料問題，為降低求解復雜度，固定初始發射軌道為地球靜止軌道。Russell［8］使用帕累托最優探索了多目標優化下的三體電推進問題，舍棄了一個打靶約束用于提高收斂性。而到目前為止探討混合推進方式軌道設計方法的研究還較少。Mingotti等［1］將電推進技術和多體引力場下的彈道捕獲技術相結合，研究了混合推進方式捕獲軌道的設計方法，但其所采用的方法并未將兩種發動機的燃料消耗總和進行統一優化。

本文探討了地－月圓型限制性三體模型下混合推進方式的最省燃料軌道設計問題。將兩種發動機的燃料消耗總和作為目標函數統一進行優化，并推導了一階必要條件，此外還在Russell［8］工作的基礎上，進一步推導了最省能量問題的狀態轉移矩陣。上述方法被用于設計從近地圓軌道出發，到達地－月L1附近Halo軌道的混合推進轉移軌道。仿真結果表明，通過設定初始發射脈沖為優化變量，可以進一步節省燃料消耗，并且飛行時間和最終剩余燃料不同的組合方式增加了任務設計的靈活性，而上述改進提高了同倫方法的魯棒性和計算效率。

1 圓型限制性三體問題

限制性三體問題是在兩個大天體P1和P2共同旋轉坐標系中研究質量可忽略不計的小天體P3的運動;如果兩個大天體繞其公共質心做圓周運動則稱為圓型限制性三體問題(Circular Restricted Three-Body Problem，CRTBP)。在地－月CRTBP下，P1代表地球，P2為月球，定義其質量分別為m1和m2，設定兩個天體之間的距離、總質量以及旋轉角速率都為1，系統的質量參數定義為μ=m2/(m1+m2)。容易驗證地球的質量為1－μ并位于(－μ，0，0)，而月球的質量為μ并位于(1－μ，0，0)。衛星P3在地－月引力場和電推進發動機共同作用下的動力學方程表示如下［8］:

式中:Tmax是最大推力幅值，u∈［0，1］為發動機推力幅值比，α是推力方向的單位矢量，c2=Isp2g0代表電推進發動機燃氣速度，(Isp2是電推進發動機比沖，g0是海平面重力加速度)，定義g(r)和h(v)為旋轉坐標系下位置與速度的函數，其具體表達式為［8］

r1和r2分別代表衛星到地球和月球的距離，表示為

2 混合推進軌道優化問題

假設衛星初始位于近地圓軌道，衛星上面級提供一次脈沖后電推進發動機開機，則混合推進方式最省燃料問題的性能指標可表示為

式(2)等號右端第一項和第二項分別代表化學推進發動機和電推進發動機的燃料消耗質量(需要注意的是，式(2)使用歸一化的質量單位，即脈沖前系統質量為單位1)。Δv為脈沖所產生的速度改變量。c1為化學推進發動機燃氣速度。t0和tf分別代表任務初始和最終時刻。對于最省燃料問題，推力幅值在任何時刻只能為輸出最大推力或者關機兩種狀態(不考慮奇異情況)，這種“bang-bang”性質使得動力學方程不連續并導致打靶過程非常敏感［3］。Bertrand等［9］使用一種平滑技術即同倫方法來改善收斂性。通過引入同倫參數ε，定義新的性能指標為

當ε=1時式(3)代表最省能量問題，當ε=0時式(3)表示最省燃料問題。

如圖1所示，假設衛星初始位于地心距為r0的圓軌道上。衛星上面級提供一次脈沖后電推進發動機開機，則施加脈沖后衛星在地－月旋轉系下的狀態x0可由兩個參數κ和θ確定，如下所示

圖1 混合推進方式的轉移軌道示意圖Fig．1 Hybrid transfer trajectory illustration

式中:κ∈［1，1．4］表示初始脈沖比，即施加脈沖后，衛星速度介于停泊軌道速度與地球逃逸速度之間。而θ∈［0，2π］表示脈沖發射相位角。此外式(4)中脈沖Δv可表示為κ的函數，即

衛星上面級施加脈沖后電推進發動機開機，電推進轉移軌道需要滿足的初始狀態約束表示為

可以發現當給定參數κ和θ時，式(5)的約束自然滿足，即x(t0)=x0(κ，θ)。假設衛星在tf時刻到達地－月系L1附近Halo軌道上固定的一點，且剩余質量m(tf)自由，電推進轉移軌道需要滿足的終端狀態約束為

為方便討論，定義式(3)中化學發動機燃料消耗質量為Θ。對于邊界約束，引入拉格朗日乘子υ0和，有約束的最優控制問題就被轉換成了無約束的最優控制問題，新的目標函數J1表示如下

假設飛行時間tf固定，J1取得極值的必要條件是其一階變分等于0，即

通過式(9)可推導協態變量的動力學方程為［8］

梯度矩陣G=?g(r)/?r，H=?h(v)/?v。此外從式(9)還能得到如下關系

式(11)也稱為橫截條件，由于 λ0，λf，υ0，υf均為未知的拉格朗日乘子，式(11)只表明如果終端時刻的質量m(tf)自由，則其所對應的協態變量應等于0，即λm(tf)=0。此外，由式(9)還能得到由于κ和θ自由而產生的兩個約束條件

龐特里亞金極大值原理(Pontryagin Maximum Principle，PMP)揭示了控制變量應在整個時間歷程上都使得哈密頓函數取極小值，則最優推力方向為

式中:S代表開關函數

推力幅值比u介于［0，1］之間，為使哈密頓函數H取極小值，需再次使用PMP。觀察式(14)可知H是u的二次函數，兩個根分別為u1=0以及u2=(ε－S)/ε，此外該二次函數開口向上。根據u2所在區間不同，最優推力幅值比u*一共有三種情況。由于u2=(ε－S)/ε，可進一步得到開關函數S和最優推力幅值比u*的關系(見式(16)和圖2)。需說明的是，對最省燃料問題(即±ε=0)，u*的選取只和開關函數S的符號有關，即當S為正時u*=0，而當S為負時u*=1。

圖2 最優推力幅值比u*根據開關函數S位置的不同具有三種情況Fig．2 Optimal thrust magnitude ratio u*with respect to the time-varying S

為方便討論，將狀態和協態變量合并為正則變量即y=［xT，λT］T，并將狀態方程(1)和協態方程(10)合并在一起稱為正則方程，表示為

至此，混合推進軌道優化問題被轉換成兩點邊值問題。當給定t0時刻9個未知的打靶變量［λT0，κ，θ］T，通過積分隱含 α*和u*的式(17)就能夠得到終端時刻tf的狀態和協態變量。而所選擇的9個打靶變量需使得9個打靶等式約束得到滿足，與之對應的打靶函數Z如式(18)所示。

一旦求解了最省能量問題(即ε=1)，從1到0逐步降低同倫參數 ε，并將當前步驟所得到的打靶變量作為初值帶入到下一步中，最終能夠得到所

對應最省燃料問題的解。需注意的是，為避免無窮轉移時間問題，需要將飛行時間固定，而最省時間問題給出了飛行時間的下限(即tf≥tfmin)，所以最省時間問題應先于最省能量和最省燃料問題進行求解。對于最省時間問題除了增加一個打靶變量tf之外，還需再增加一個打靶約束，即哈密頓函數在tf時刻的值H(tf)等于0。由于篇幅所限，最省時間問題的一階必要條件在此不再進行推導。

3 解析雅克比矩陣

如前所述，最優控制問題被轉換成兩點邊值問題并通過打靶法進行求解，但與打靶法相關的牛頓法對于初始點的選取非常敏感。在本文中，為提高牛頓法的魯棒性和收斂速度，解析的雅克比矩陣被用于提供打靶過程的梯度信息。狀態轉移矩陣(State Transition Matrix，STM)反映了正則變量y在t0時刻的擾動對t時刻所產生的影響。定義STM為Φ，其性質如下

1)在一個積分步長內如果發動機始終處于開機或者關機狀態(即u*為標量)，則表示為［8］

2)當在一個積分步長內發動機推力幅值始終處于中等推力時，由于u*是m，λv以及 λm的函數，參見式(15)和式(16)，通過鏈式法則可得

式中:

需注意的是，上述狀態轉移矩陣只對動力學方程連續的正則變量進行映射，即發動機始終保持開機、中等推力或者關機。然而如果推力幅值u*在積分步長內發生切換，則需考慮切換點tj處的狀態轉移矩陣 ψ(tj)［8］:

至此，式(23)中的狀態轉移矩陣Φ(tf，ti)就可用于構造9個打靶等式約束相對于9個打靶變量的雅可比矩陣，表示如下(參見式(18))

4 仿真和結果

假設衛星初始位于近地點高度為400 km的圓軌道，衛星上面級施加一次脈沖Δv后，電推進發動機開始工作并不斷提升衛星能量，最終衛星在tf時刻與地－月L1點附近Halo軌道上的一個固定點重合。表1給出了本文所用到的參數，其中長度、時間、速度以及質量都進行了歸一化處理。

表1 物理常量及含義Table 1 Simulation parameters

首先求解最省時間問題以確定飛行時間下限。假設衛星上面級提供一次固定大小的脈沖，求解到達L1附近Halo軌道的電推進最省時間轉移軌道，其仿真結果如表2所示。當κ=1時代表衛星上面級不提供初始脈沖而只依靠電推進發動機提供動力，整個任務耗時22．6467天，最終剩余質量比為0．7783，其轉移軌道如圖3所示。隨著脈沖比κ從1逐漸增加到1．4，任務時間從22．6467天縮短到了3．3692天，最終剩余質量從0．7783降低到0．3200。當 κ等于1．4時，衛星被施加了一個近似等于地球逃逸速度的脈沖，其轉移軌道如圖4所示。

表2 發射脈沖固定的電推進最省時間問題收斂解Table 2 Simulation results of minimum-time problem with fixed initial impulse and electronic propulsion

圖3 初始發射脈沖為零(即κ=1，θ=－2)的最省時間電推進轉移軌道Fig．3 Minimum-time low-thrust trajectory with initial impulse equal to zero(i．e．，κ=1，θ=－2)

圖4 初始發射脈沖約為地球逃逸速度(即κ=1．4，θ=－2)的最省時間電推進轉移軌道Fig．4 Minimum-time low-thrust trajectory with high initial impulse(i．e．，κ=1．4，θ=－2)

在得到不同κ對應的最短轉移時間tfmin之后，可進一步求解混合推進方式的最省能量和最省燃料問題。需說明的是，打靶過程中κ和θ的選取非常敏感，可先將其固定，找到收斂解后再將其依次設定為優化變量。此外為避免“無窮轉移時間”問題，可通過引入標量ctf來固定飛行時間，即tf=ctftfmin。

表3 混合推進方式最省燃料問題部分收斂解Table 3 Minimum-fuel solutions with hybrid propulsion system

表3列出了混合推進方式的部分收斂解，其中tf為飛行時間，mf為最終剩余質量。通過設定κ和θ為優化變量，mf得到了不同程度的提高。下面選取一組收斂解對比轉移軌道和電推進發動機的控制曲線。

1)固定 κ=1．24、θ=－2、飛行時間比ctf=1．2，此時在CRTBP下最省燃料轉移軌道如圖5所示。其中粗實線代表電推進發動機推進段，細實線代表無動力滑行段，虛線代表Halo軌道。可以發現電推進發動機在近地點附近開機，并在遠地點處關機用于節省燃料消耗。

圖5 地－月旋轉系下固定發射脈沖的最省燃料轉移軌道(κ=1．24，θ=－2且c tf=1．2)Fig．5 Minimum-fuel low-thrust trajectory in CRTBP with fixed initial impulse(κ=1．24，θ=－2 and c tf=1．2)

圖6為推力幅值比和開關函數隨時間變化的曲線，推力幅值呈現“bang-bang”控制結構。通過計算可知，哈密頓函數在整個時間歷程上都近似為常值，其最大誤差在10－10以內，從而證明了解的精確性。

圖6 固定發射脈沖時推力幅值比u和開關函數S隨時間變化曲線(κ=1．24，θ=－2且c tf=1．2)Fig．6 When the initial impulse is fixed，the time history of throttle factor u and switching function S(κ=1．24，θ=－2 and c tf=1．2)

2)設定κ和θ為優化變量，飛行時間比ctf=1．2，混合推進方式的轉移軌道在地－月旋轉系下的仿真圖如圖7所示。通過和圖5對比發現，采用混合推進方式的衛星在初始時刻施加了一個相對較小的初始脈沖，但在近地空間內保持連續開機以提升軌道能量。此外還需要說明的是，接近目標周期軌道的無動力滑行段非常逼近從周期軌道上出發的一條穩定流形，即說明最優軌道能夠自動沿著穩定流型運動以節省燃料。之所以在Halo軌道附近仍需施加一段推力弧段，是因為終端約束為Halo軌道上的一個固定點，如果最終不施加推力弧段，則衛星沿著穩定流形需要無窮長的時間才能最終進入Halo軌道，而在固定轉移時間的前提下，衛星只能先沿著近似穩定流形的轉移軌道逼近以節省燃料，最后發動機開機并最終進入Halo軌道。

圖7 地－月旋轉系下初始發射脈沖自由的最省燃料轉移軌道(κ和θ自由且c tf=1．2)Fig．7 Minimum-fuel low-thrust trajectory in CRTBP with free initial impulse(freeκandθ，c tf=1．2)

電推進發動機的推力幅值比和開關函數隨時間變化的曲線如圖8所示，可以發現設定κ和θ為優化變量之后，衛星在任務時間的兩端保持開機，而只在仿真時間的中間部分關機。

圖8 初始發射脈沖自由時推力幅值比u和開關函數S隨時間變化曲線(κ和θ自由且c tf=1．2)Fig．8 When the initial impulse is free，the time history of throttle factor u and switching function S(freeκandθ，c tf=1．2)

5 結論

本文基于混合推進方式探討了地－月圓型限制性三體模型下最優轉移軌道的設計問題。根據變分法推導了最優解的一階必要條件。針對復雜動力學模型下多圈轉移軌道帶來的打靶迭代速度慢和分段連續的推力幅值導致打靶難以收斂的問題，在前人工作基礎上，給出了最省能量問題解析的狀態轉移矩陣。上述方法被用于求解從近地圓軌道出發，到達地－月L1附近Halo軌道的混合推進轉移軌道。數值仿真表明，上述方法可以將化學發動機和電推進發動機的燃料消耗總和進行優化。此外由于給出了燃料消耗以及飛行時間不同的組合方式，混合推進方式使任務設計具有更大的靈活性。

［1］Mingotti G，Topputo F，Bernelli Z F．Earth–Mars transfers with ballistic escape and low-thrust capture［J］．Celestial Mech．and Dyna．Astro．2011，110(2):169－188．

［2］ 任遠，崔平遠，欒恩杰．利用混合法進行地球－火星電推進軌道設計［J］．哈爾濱工業大學學報，2007，39(3):359－362．［Ren Yuan，Cui Ping-yuan，Luan En-jie．An Earth-Mars low-thrust trajectory design based on hybrid method［J］．Journal of Harbin Institude of Technology，2007，39(3):359－362．］

［3］Jiang F H，Baoyin H X，Li J F．Practical techniques for lowthrust trajectory optimization with homotopic approach［J］．J．of Guid．Con．and Dyna．2012，35(1):245–258．

［4］Guo T D，Jiang F H，Baoyin H X，et al．Fuel optimal low thrust rendezvous with outer planets via gravity assist［J］．Sci．China Phy．Mech．and Astro．2011，54(4):756－769．

［5］ 黃鎬，韓潮．基于間接法的最優電推進逃逸軌道設計［J］．中國空間科學技術，2013，33(2):25－31．［Huang Hao，Han Chao．Indirect optimization of low thrust escape trajectories［J］．Chinese Space Sci．and Tech．2013，33(2):25－31．］

［6］Ozimek M，Howell K．Low-thrust transfers in the Earth-Moon system，including applications to libration point orbits［J］．J．of Guid．Con．and Dyna．2010，33(2):533–549．

［7］Caillau J，Daoud B，Gergaud J．Minimum fuel control of the planar circular restricted three-body problem［J］．Celestial Mech．and Dyn．Astro．2012，114(1－2):137–150．

［8］Russell R．Primer vector theory applied to global low-thrust trade studies［J］．J．of Guid．Con．and Dyna．2007，30(2):460–472．

［9］Bertrand R，Epenoy R．New smoothing techniques for solving bang–bang optimal control problems numerical results and statistical interpretation［J］．Optimal Control Applications and Methods，2002，23(4):171–197．

［10］Lawden．Optimal trajectories for space navigation［M］．London UK:Butterworths，1963．