潘迅，泮斌峰

（西北工業(yè)大學(xué) 航天學(xué)院，西安 710012）

基于同倫方法三體問題小推力推進轉(zhuǎn)移軌道設(shè)計

潘迅，泮斌峰

（西北工業(yè)大學(xué) 航天學(xué)院，西安 710012）

提出一種基于同倫方法限制性三體問題小推力推進轉(zhuǎn)移軌道設(shè)計方法。首先根據(jù)最優(yōu)控制原理分析了航天器在軌道轉(zhuǎn)移中不同性能指標時的最優(yōu)控制律，然后引入同倫參數(shù)構(gòu)造新的性能指標，在基于遺傳算法和打靶法得到能量最優(yōu)的解基礎(chǔ)上，采用偽弧長法跟蹤同倫軌跡，進而得到燃料最優(yōu)的轉(zhuǎn)移軌道。最后對地月系下從GEO軌道到L1點Lyapunov軌道的轉(zhuǎn)移軌道進行優(yōu)化。仿真結(jié)果表明：利用遺傳算法能優(yōu)化得到較為合適的流形拼接點和協(xié)態(tài)變量初始值，利用打靶法能有效地優(yōu)化得到小推力燃料最優(yōu)轉(zhuǎn)移軌道。

同倫算法；偽弧長法；小推力；軌道優(yōu)化；限制性三體問題

0 引 言

在深空探測中，小推力發(fā)動機具有比沖大、控制精度高等優(yōu)點，采用小推力推進系統(tǒng)能減少燃料消耗，增加航天器有效載荷和航天任務(wù)的科學(xué)回報。相比于脈沖推進的動力學(xué)模型，小推力推進系統(tǒng)的動力學(xué)模型具有更強的非線性，小推力發(fā)動機的推力很小，所以改變航天器運動狀態(tài)需要較長的時間，從而增加了軌道設(shè)計和優(yōu)化的難度。相比于傳統(tǒng)的二體模型，圓限制性三體問題的模型更精確，但其動力學(xué)模型更復(fù)雜，在三體模型下更難進行轉(zhuǎn)移軌道優(yōu)化。在航天器轉(zhuǎn)移軌道設(shè)計過程中，可將軌道優(yōu)化看作是一類復(fù)雜的、非線性并受到嚴格約束限制的最優(yōu)控制問題，其求解的理論基礎(chǔ)為最優(yōu)控制理論。根據(jù)其對動力學(xué)模型處理方式的不同主要分為直接法和間接法。目前對小推力轉(zhuǎn)移軌道的研究較為豐富，但其中大部分是基于二體模型的軌道優(yōu)化[1-5]。在限制性三體問題中，Howell等[6]研究了小推力變比沖發(fā)動機推進的地球到平動點的轉(zhuǎn)移；Daero等[7]利用配點法結(jié)合小推力和脈沖推進進行平動點周期軌道的轉(zhuǎn)移；Mingotti等基于直接配點法進行了小推力轉(zhuǎn)移軌道研究，并取得了較為豐富的成果[8-10]。相比于直接法，間接法雖然具有對初始條件十分敏感、收斂域小等缺點，但其可通過對連續(xù)系統(tǒng)分析得到精確的最優(yōu)解。Jiang等[11]結(jié)合粒子群算法和同倫算法進行了小推力轉(zhuǎn)移研究。目前針對間接法求解三體模型下小推力轉(zhuǎn)移軌道的研究較少，有必要對其進行更深入的研究。

本文以地月系限制性三體問題模型為研究對象，分析了不同性能指標下的最優(yōu)推力控制律。針對燃料最優(yōu)的bang-bang控制問題，通過引入同倫參數(shù)構(gòu)造新的性能指標，并利用偽弧長法跟蹤得到最優(yōu)控制律，最終得到燃料最優(yōu)的轉(zhuǎn)移軌道。最后以GEO軌道到L1點Lyapunov軌道的小推力轉(zhuǎn)移為例驗證了算法的有效性。

1 圓限制性三體問題

圓限制性三體問題（Circular-Restricted Three Bodies Problem，CRTBP）是指兩個主天體繞著系統(tǒng)質(zhì)心做勻速圓周運動，而第三體（一般指航天器）質(zhì)量很小，遠小于兩個主天體的質(zhì)量，對主天體的運動的影響較小，可以忽略不計。月球軌道偏心率0.054 9，圓限制性三體問題的動力學(xué)模型能較好地滿足實際情況，因此本文在地月系圓限制性三體問題模型下進行轉(zhuǎn)移軌道設(shè)計研究。

考慮地月圓限制性三體問題平面運動，以旋轉(zhuǎn)坐標系為參考坐標系，原點與地月系質(zhì)心重合，x軸為地月連線并指向月球，y軸與其構(gòu)成右手坐標系，建立其動力學(xué)模型為

其中

μ= 0.012 15為地月系三體系統(tǒng)唯一參數(shù)，月球質(zhì)量與地月質(zhì)量之和的比值；為航天器在坐標系中的位置矢量，為速度矢量，Isp為發(fā)動機比沖，g0為地球海平面重力加速度，Tmax為發(fā)動機最大推力，u為發(fā)動機工作效率為單位推力在坐標軸上的分量。

2 轉(zhuǎn)移軌道優(yōu)化問題描述

在轉(zhuǎn)移軌道設(shè)計中，主要考慮的性能指標包括能量最優(yōu)和燃料最優(yōu)，其只與推力大小u有關(guān)，與推力方向α無關(guān)，因此可將性能指標寫為如下形式

根據(jù)最優(yōu)控制理論，通過引入?yún)f(xié)態(tài)變量

構(gòu)造哈密爾頓函數(shù)

將哈密爾頓函數(shù)對狀態(tài)變量求一階偏導(dǎo)，可得到歐拉方程為

對于控制量α，為使H最小，α應(yīng)與速度協(xié)態(tài)變量λv方向相反，即

對于初始時刻和終端時刻航天器位置都給定的轉(zhuǎn)移，其在tf時刻滿足狀態(tài)約束

對于終端時刻航天器質(zhì)量和轉(zhuǎn)移時間不固定的情況，其橫截條件為

根據(jù)龐德里亞金極小值原理，對控制量u，需滿足最優(yōu)性一階必要條件。下面考慮不同性能指標時候的u的控制律。

1）性能指標為燃料最優(yōu)

根據(jù)哈密爾頓函數(shù)對u求偏導(dǎo)，有

由式（16）可知，此時的轉(zhuǎn)移軌道最優(yōu)控制為bang-bang控制，且在S= 0時奇異。由于S連續(xù)變化，一般不會出現(xiàn)S≡ 0，因此本文中不考慮奇異段。

2）性能指標為能量最優(yōu)

對哈密爾頓函數(shù)求偏導(dǎo)，可得

則最優(yōu)控制律為

由式（19）可知，此時的推力大小u為連續(xù)控制。

對于一般的航天任務(wù)，其任務(wù)的開展依賴于有效載荷，而節(jié)省燃料不僅對增加有效載荷有著重要作用，而且對航天器壽命起關(guān)鍵作用。對于燃料最優(yōu)的軌道轉(zhuǎn)移，其推力控制為bang-bang控制，利用間接法進行優(yōu)化時存在很大的難度。針對該問題，Bertrand等[12]提出了一種平滑技術(shù)，即所謂的同倫方法，從而降低了求解難度。本文在其基礎(chǔ)上，利用同倫算法進行限制性三體問題下的最優(yōu)燃料轉(zhuǎn)移軌道設(shè)計。

3 同倫算法

引入同倫參數(shù)p，構(gòu)造新的性能指標

其中：ε= 1 -p，當p≠ 1時，u為連續(xù)控制，且p趨近于1時，即可得到燃料最優(yōu)的轉(zhuǎn)移軌道。

在得到能量最優(yōu)問題的解之后，同倫參數(shù)p需按照一定步長進行迭代，并將當前得到的解作為下一步迭代的初值進行計算。該方法看似簡單，但是仍存在一些困難：第一，同倫迭代過程中可能存在奇異點，采用同倫參數(shù)遞增方法在奇異點將導(dǎo)致計算失?。坏诙?，當p趨近于1時，雖然理論上推力u是連續(xù)的，但其接近bang-bang控制，在開關(guān)切換區(qū)域變化劇烈，p按給定步長進行減小時，很難保證其精度，因此本文中利用偽弧長方法進行迭代求解。

弧長法最初由Riks和Wempne提出，后經(jīng)Ramm和Crisfield發(fā)展，其基本思想是引入弧長參數(shù)，通過加入該參數(shù)在解曲線上增加一個約束方程，從而有效地解決了非線性分析中的穩(wěn)定性與收斂性問題。將終端等式約束和橫截條件表示成

對求（22）一階偏導(dǎo)，有

其中：ds為偽弧長。當步長足夠小時，根據(jù)幾何關(guān)系，有

結(jié)合式（23）和（24），可表示成

4 轉(zhuǎn)移軌道設(shè)計

為驗證同倫算法在燃料最優(yōu)轉(zhuǎn)移軌道優(yōu)化中的有效性，本節(jié)進行從GEO軌道到L1點Lyapunov軌道的轉(zhuǎn)移軌道設(shè)計，相關(guān)參數(shù)如表 1所示，其中RE為地球半徑，hGEO為GEO軌道高度，g0為海平面重力加速度，Ax為Lyapunov軌道幅值。

表1 轉(zhuǎn)移軌道及發(fā)動機初始參數(shù)Table 1 The parameters of transfers and engine

航天器從GEO軌道到Lyapunov軌道的轉(zhuǎn)移過程可以描述為：航天器從GEO軌道出發(fā)，通過小推力發(fā)動機使航天器運動狀態(tài)發(fā)生改變，并在拼接點進入Lyapunov軌道的穩(wěn)定流形，沿流形無動力滑行，最終進入目標Lyapunov軌道。初始時刻航天器在GEO軌道上的初始位置可用與x軸之間的夾角ω表示，有

軌道推進段轉(zhuǎn)移時間記為tf，拼接點位于Lyapunov軌道的穩(wěn)定流形上，可表示為其中τ表示流形初始點在Lyapunov軌道上的位置，θ表示拼接點在該流形上的位置。根據(jù)時間在Lyapunov軌道上取360個點，分別計算得到360條流形，以與x軸負半軸相交為流形終點，再將流形按積分時間均勻取2 000個點，即τ∈[1，360]，θ∈[1，2 000]，且τ，θ∈N，如圖 1所示，圖中DU為無量綱化長度，1 DU = 38.44萬km。因此，初始時刻航天器位置和終端約束可用[ω，τ，θ]表示。確定轉(zhuǎn)移段的初始點和終點后，根據(jù)極小值原理，轉(zhuǎn)移軌道的優(yōu)化可轉(zhuǎn)化為對進行求解。

圖1 拼接點在流形上的位置Fig.1 Position of insertion points in manifolds

對于轉(zhuǎn)移軌道的設(shè)計優(yōu)化，其過程為：

1）針對能量最優(yōu)的轉(zhuǎn)移軌道，其控制律為式（19），利用遺傳算法對變量定遺傳算法的指標為流形段拼接點與小推力轉(zhuǎn)移段終點的位置速度誤差值最小，得到較為合理的初值。

2）以得到的流形拼接點為小推力轉(zhuǎn)移段終點，并在保持轉(zhuǎn)移段初始點不變的情況下，用打靶法對進一步優(yōu)化，得到能量最優(yōu)的轉(zhuǎn)移軌道。

3）采用偽弧長法跟蹤同倫軌跡，最終得到燃料最優(yōu)的小推力轉(zhuǎn)移軌道。先采用遺傳算法是因為其具有較好的全局收斂性，但由于其不能滿足嚴格的約束條件，因此需要多次計算，選取合適的值，并用打靶法進一步優(yōu)化，然后才能作為初值進行同倫演化。根據(jù)多次計算，選取遺傳算法得到的起點和終點參數(shù)為ω= 4.712 4，τ= 44，θ= 1 538，進一步優(yōu)化過程中的協(xié)態(tài)變量的值如表 2所示。同倫過程中，變量隨參數(shù)p變化的關(guān)系如圖 2所示。航天器剩余質(zhì)量隨著同倫參數(shù)p的增加而增加，燃料消耗質(zhì)量比從0.056 074減小到0.055 226，減小了1.51%。

表2 能量最優(yōu)和燃料最優(yōu)的變量值Table 2 The solutions of energy optimal transfer and fuel optimal transfer

圖2 優(yōu)化變量隨同倫參數(shù)p的變化關(guān)系Fig.2 The relationships between the optimization variable and the homotopy parameterp

旋轉(zhuǎn)坐標系下從GEO軌道到L1點Lyapunov軌道的轉(zhuǎn)移軌道如圖 3所示，由從GEO出發(fā)的優(yōu)化段和流形段組成，其中優(yōu)化段為虛線表示性能指標為能量最優(yōu)的轉(zhuǎn)移軌道，實線表示燃料最優(yōu)的轉(zhuǎn)移軌道。燃料最優(yōu)時的開關(guān)函數(shù)S及推力控制u如圖 4所示，其中1 TU = 37.576萬秒，為無量綱時間單位。此時推力控制為典型的bang-bang控制，發(fā)動機開機3次，關(guān)機2次。

圖3 旋轉(zhuǎn)坐標系下GEO軌道到L1點Lyapunov軌道的轉(zhuǎn)移軌道Fig.3 The transfers from GEO to L1 Lyapunov orbit in rotating frame

圖4 燃料最優(yōu)時的開關(guān)函數(shù)和推力控制曲線Fig.4 The switching curve and thrust magnitude about the fuel optimal transfer

5 結(jié) 論

采用間接法設(shè)計小推力推進轉(zhuǎn)移軌道優(yōu)化問題時，當性能指標為燃料最優(yōu)時存在控制變量的bangbang控制形式導(dǎo)致優(yōu)化問題高度敏感和難以求解。本文采用同倫方法實現(xiàn)從能量最優(yōu)軌道到燃料最優(yōu)軌道的連續(xù)變換，并采用擬弧長法實現(xiàn)同倫軌跡的跟蹤。本文以地月平面圓形限制性三體問題模型為研究對象并進行仿真驗證，仿真結(jié)果表明：本文提出的算法能夠自動確定轉(zhuǎn)移過程中發(fā)動機開關(guān)機次數(shù)和切換時間，實現(xiàn)限制性三體問題下的小推力燃料最優(yōu)轉(zhuǎn)移軌道的優(yōu)化設(shè)計。

[1]John T B.Practical methods for optimal control and estimation using nonlinear programming[M].USA：The Society for Industrial and Applied Mathematics，2010.

[2]尚海濱，崔平遠，徐瑞，等.基于高斯偽光譜的星際小推力轉(zhuǎn)移軌道快速優(yōu)化[J].宇航學(xué)報，2010，4（31）：1005-1011.Shang H B，Cui P Y，Xu R，et al.Fast optimization of interplanetary low-thrust transfer trajectory based on Gauss pseudospectral algorithm[J].Journal of Astronautics，2010，4（31）：1005-1011.

[3]李俊峰，蔣方華.連續(xù)小推力航天器的深空探測軌道優(yōu)化方法綜述[J].力學(xué)與實踐，2011，33（3）：1-6.Li J F，Jiang F H.Survy of low-thrust trajectory optimization methods for deep space exploration[J].Advances in Mechanics，2011，33（3）：1-6.

[4]Leonardo M V.Finite thrust orbital transfers[J].Acta Astronautica，2014，100（1）：107-128.

[5]鄭博，張澤旭，周浩，等.一種小推力借力飛行轉(zhuǎn)移軌道初始設(shè)計方法[J].深空探測學(xué)報，2015，2（3）：256-261.Zheng B，Zhang Z X， Zhou H，et al.A novel initial design for lowthrust transfer trajectories using gravity-assist[J].Journal of Deep Space Exploration，2015，2（3）：256-261.

[6]Ozimek M T，Howell K C.Low-thrust transfers in the Earth-Moon system including applications to libration point orbits[J].Journal of Guidance，Control，and Dynamics，33（2）：533-549.

[7]Daero L，Eric A B，Amit K S.Optimal interior Earth-Moon Lagrange point transfer trajectories using mixed impulsive and continuous thrust[J].Aerospace Science and Technology，2014，39：281-292.

[8]Mingotti G，Topputo F，Bernelli-Zazzera F.Low-energy，low-thrust transfers to the Moon[J].Celestial Mechanics and Dynamical Astronomy，2009，105（1-3）：61-74.

[9]Mingotti G，Topputo F，Bernelli-Zazzera F.Optimal low-thrust invariant manifold trajectories via attainable sets[J].Journal of Guidance，Control，and Dynamics，2011，34（6）：1644-1656.

[10]Mingotti G，Topputo F，Bernelli-Zazzera F.Efficient invariantmanifold，low-thrust planar trajectories to the Moon[J].Communications in Nonlinear Science and Numerical Simulation，2012，17（2）：817-831.

[11]Jiang F H，Baoyin H X，Li J F.Practical techniques for low-thrust trajectory optimization with homotopic approach[J].Journal of Guidance，Control，and Dynamics，2012，35（1）：245-258.

[12]Betrand R，Epenoy R.New smoothing techniques for solving bangbang optimal control problems-numerical results and statistical interpretation[J].Optimal Control Applications and Methods，2002，23（4）：171-197.

Optimization of Low-Thrust Transfers Using Homotopic Method in the Restricted Three-Body Problem

PAN Xun，PAN Binfeng
（School of Astronautics，Northwestern Polytechnical University，Xi’an 710072，China）

A method for optimization of low-thrust transfers in the restricted three-body problem is proposed.First，the optimal control laws of different performance index in trajectory optimization are deduced based on the optimal control theory.Then，a new parameter is used to construct a modified performance index.The genetic algorithm and single shooting method are used to obtain the solution of energy optimal transfer，and arc-length method is adopted to track the homotopic pathtill the fuel optimal transfer is obtained.Finally，the numerical example about transfers from GEO to Lyapunov orbit of L1 point in the Earth-Moon system is studied.Simulation results show that the initial guess can be obtained by genetic algorithm，and single shooting method can be used to obtain the fuel-optimal transfer trajectory.

homotopy；arc-length method；low-thrust；trajectory optimization；restricted three-body problem

V41

A

2095-7777（2017）03-0270-06

[責(zé)任編輯：高莎，英文審校：任樹芳]

10.15982/j.issn.2095-7717.2017.03.011

潘迅，泮斌峰.基于同倫方法三體問題小推力推進轉(zhuǎn)移軌道設(shè)計[J].深空探測學(xué)報，2017，4（3）：270-275.

Reference format：Pan X，Pan B F.Optimization of low-thrust transfers using homotopic method in the restricted three-body problem[J].Journal of Deep Space Exploration，2017，4（3）：270-275.

2016-09-25

2016-12-01

潘迅（1990- ），男，博士，主要研究方向：空天飛行器動力學(xué)與制導(dǎo)控制。

通訊地址：陜西省西安市碑林區(qū)友誼西路127號西北工業(yè)大學(xué)航天學(xué)院（710072）

E-mail：xpan2012@gmail.com