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

潘迅，泮斌峰

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

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

潘迅，泮斌峰

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

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

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

0 引 言

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

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

1 圓限制性三體問題

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

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

其中

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 同倫算法

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3）采用偽弧長法跟蹤同倫軌跡，最終得到燃料最優(yōu)的小推力轉(zhuǎn)移軌道。先采用遺傳算法是因?yàn)槠渚哂休^好的全局收斂性，但由于其不能滿足嚴(yán)格的約束條件，因此需要多次計(jì)算，選取合適的值，并用打靶法進(jìn)一步優(yōu)化，然后才能作為初值進(jìn)行同倫演化。根據(jù)多次計(jì)算，選取遺傳算法得到的起點(diǎn)和終點(diǎn)參數(shù)為ω= 4.712 4，τ= 44，θ= 1 538，進(jìn)一步優(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)坐標(biāo)系下從GEO軌道到L1點(diǎn)Lyapunov軌道的轉(zhuǎn)移軌道如圖 3所示，由從GEO出發(fā)的優(yōu)化段和流形段組成，其中優(yōu)化段為虛線表示性能指標(biāo)為能量最優(yōu)的轉(zhuǎn)移軌道，實(shí)線表示燃料最優(yōu)的轉(zhuǎn)移軌道。燃料最優(yōu)時(shí)的開關(guān)函數(shù)S及推力控制u如圖 4所示，其中1 TU = 37.576萬秒，為無量綱時(shí)間單位。此時(shí)推力控制為典型的bang-bang控制，發(fā)動(dòng)機(jī)開機(jī)3次，關(guān)機(jī)2次。

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

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

5 結(jié) 論

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

[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]尚海濱，崔平遠(yuǎn)，徐瑞，等.基于高斯偽光譜的星際小推力轉(zhuǎn)移軌道快速優(yōu)化[J].宇航學(xué)報(bào)，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é)與實(shí)踐，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ì)方法[J].深空探測學(xué)報(bào)，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

潘迅，泮斌峰.基于同倫方法三體問題小推力推進(jìn)轉(zhuǎn)移軌道設(shè)計(jì)[J].深空探測學(xué)報(bào)，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- ），男，博士，主要研究方向：空天飛行器動(dòng)力學(xué)與制導(dǎo)控制。

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

E-mail：xpan2012@gmail.com