考慮J2項攝動的小推力燃料最優轉移軌道設計

潘 迅，泮斌峰，唐 碩

(1.西北工業大學 航天學院，西安710072；2.航天飛行動力學技術國家級重點實驗室(西北工業大學)，西安 710072)

考慮J2項攝動的小推力燃料最優轉移軌道設計

潘 迅1,2，泮斌峰1,2，唐 碩1

(1.西北工業大學 航天學院，西安710072；2.航天飛行動力學技術國家級重點實驗室(西北工業大學)，西安 710072)

為優化得到考慮地球扁率J2項攝動影響的小推力燃料最優轉移軌道，提出了一種3次同倫方法.構造較簡單的采用“線性引力”，且不考慮J2項攝動的大推力能量最優轉移軌道作為同倫初始問題.引入3個同倫參數，分別對動力學模型、推力大小和性能指標進行同倫，根據極小值原理推導得到同倫過程中的最優控制律，并通過跟蹤同倫參數的連續變化求解一系列的同倫迭代子問題，分別得到J2攝動模型下的大推力能量最優轉移軌道和小推力能量最優轉移軌道，并最終優化得到小推力燃料最優轉移軌道.以航天器與位于太陽同步軌道的碎片的交會任務為算例進行數值仿真，驗證所提出的3次同倫方法在求解J2項攝動影響下的小推力燃料最優轉移軌道優化問題中的有效性.結果表明,利用打靶法容易對同倫初始問題進行求解，在同倫過程中能連續穩定地跟蹤同倫參數，進而得到所需的燃料最優小推力轉移軌道，利用該方法能有效地解決J2項攝動導致的非線性強、推力小、轉移圈數多等原因所導致的一般數值優化算法不易收斂的難題.

軌道轉移；J2項攝動；小推力推進；同倫方法；最優控制

航天器軌道動力學與控制是航天技術工程中的重要組成部分.在對航天器的軌道機動進行優化設計時，通常將地球視為質點模型，航天器在中心引力場中進行機動變軌.然而在實際飛行中，航天器受到空間環境各種攝動因素的影響，其軌道要素發生變化.這些攝動力包括地球形狀非球形的附加引力、大氣阻力、日月引力、太陽光壓等攝動力，其中地球扁率的J2項攝動遠大于其他攝動因素，是引起航天器軌道變化的主要因素.

考慮J2項攝動的軌道優化問題，逐漸成為國內外學者的研究熱點.文獻[1]提出了基于開普勒二體運動推算地球引力模型J2攝動的間接補償的方法; 對于采用脈沖推進的航天器，文獻[2] 在二體模型的基礎上，對J2項攝動下的Lambert問題進行求解；文獻[3]闡述了J2攝動的原因和模型，并對J2攝動的補償算法進行研究；文獻[4]通過理論分析確定兩次脈沖假設下，考慮地球扁率的機動軌道在地心轉移角接近π時對邊界條件比較敏感的原因，并提出相應解決辦法；文獻[5]研究了在地球扁率影響下，固定時間兩異面橢圓軌道間的多沖量最優交會控制問題；文獻[6]針對J2項攝動影響下的多脈沖最優機動的衛星編隊構型的建立和重構，提出了一種高效的在線計算方法.

與脈沖推進相比，小推力發動機具有比沖大、控制精度高等優點，利用小推力推進進行軌道機動是未來發展趨勢.由于小推力發動機推力小，持續時間長，其軌道優化存在較大困難.對于小推力轉移軌道優化，該過程屬于轉移軌道設計過程中的局部優化，優化方法主要可以分為間接法、直接法和混合法，其中間接法通過推導最優性一階必要條件，得到局部最優解，其最優性優于混合法和直接法，文獻[7]對小推力轉移軌道的優化方法進行了詳細的比較說明.文獻[8]通過間接法對深空探測中的燃料最優小推力轉移軌道進行研究，利用粒子群算法選取初值，結合同倫算法和開關函數檢測法克服Bang-Bang控制的困難.文獻[9]對轉移軌道的推力大小同倫過程進行研究，提出2次同倫方法，克服同倫過程不連續的問題.對于考慮J2項的小推力轉移.文獻[10]對J2項攝動影響下，衛星編隊重構的燃料最優機動策略進行研究，通過Gauss偽譜法進行離散，并利用非線性優化軟件TOMLAB/SNOPT進行優化. 文獻[11]利用直接法研究了考慮J2項攝動的平面內最小時間的小推力轉移. 文獻[12]對于推力加速度和J2項攝動加速度為同一量級的情況下，將小推力和地球扁率視為攝動力，研究了小推力和J2項攝動對軌道要素的影響. 文獻[13]以時間最短為性能指標，推力幅值恒定，利用Gauss偽譜法對攝動情況下的有限推力轉移軌道和交會聯合優化問題進行了研究.

本文針對發動機推力遠小于J2項攝動力的燃料最優轉移軌道問題，基于間接法和同倫法對其進行研究，通過引入多個同倫參數，分別對動力學模型、推力大小和性能指標進行同倫.首先采用“線性引力”假設，選取性能指標為能量最優，求解較大推力的轉移軌道；對動力學模型進行同倫，得到考慮J2攝動的較大推力能量最優轉移軌道；然后對推力大小進行同倫，得到小推力的能量最優轉移軌道；最后對性能指標進行同倫，克服Bang-Bang控制不連續的問題，最終得到所需的燃料最優轉移軌道.

1 考慮J2項攝動的動力學模型

對于在近地軌道運行的航天器，其受到的地球扁率所引起的J2項攝動力比三體引力、太陽光壓等其他攝動力大數個量級，是主要的攝動力，在此僅考慮航天器受到地球中心引力和J2項攝動的影響.此時，航天器的動力學方程為

(1)

式中：aJ2=[axJ2,ayJ2,azJ2]T；J2=0.001 082 63為地球扁率攝動常數；RE=6 378 137 m為地球平均半徑；r為航天器的地心距.

對于轉移軌道設計問題，可將性能指標表示為

結合動力學模型(1)，根據龐德里亞金極小值原理，引入拉格朗日乘子，構造哈密爾頓函數，有

為使哈密爾頓函數取得極小值，對其求一階偏導，可得歐拉方程為：

(2)

(3)

(4)

式中B2=[0,0,0;0,0,0;0,0,1]T.

根據最優控制理論，可推導得到最優推力方向為

在對推力大小u的最優控制律進行推導時，需要考慮不同性能指標的影響：

1)當性能指標為能量最優時，其最小化性能指標可以表示成

將其代入到哈密爾頓函數，并對哈密爾頓函數取極值，可得推力大小的最優控制律為

(5)

其中S為開關函數，其表達式為

2)當性能指標為燃料最優時，其最小化性能指標可以表示成

同理可得其最優控制律為

(6)

推導得到推力大小和推力方向的最優控制律后，可通過對動力學模型和歐拉方程進行積分，得到狀態量和協態變量的瞬時值，再根據控制律得到推力大小和方向，從而將軌道優化問題轉移成兩點邊值問題.

2 多次同倫方法

在將轉移軌道優化問題轉換為兩點邊值問題之后，雖然理論上可通過打靶法對其進行求解，然而在實際求解過程中，存在下列3個難點：

1)與只考慮地球中心引力的情況進行對比，由于J2項攝動的引入，動力學模型的非線性明顯增強，增加了打靶收斂的難度；

2)對于小推力轉移軌道，其推力加速度一般為10-4m/s2量級，遠小于地球引力，在軌道轉移過程中需要作用很長時間才能改變軌道完成轉移，從而在打靶過程中容易發散；

3)從式(6)可以看出，對于燃料最優的轉移軌道，其最優控制為Bang-Bang控制，推力u在開關函數S=0時刻存在突變，在對轉移軌道進行積分過程中推力發生階躍，不利于打靶求解.

針對上述3個問題，本文提出一種求解方法，通過多次同倫，依次解決存在的問題，從而得到最終的轉移軌道.

2.1 動力學方程的同倫

針對包含J2項攝動的動力學模型的復雜非線性問題，從簡單動力學模型出發，先對不包含J2項，且采用“線性引力”假設的動力學模型進行求解，得到初始解，然后將動力學模型同倫到包含J2項攝動的真實引力的動力學模型.

“線性引力”近似是指將航天器受到的地球引力近似為大小不變、方向變化的力.引入第1個同倫參數ε1，在同倫過程中航天器的動力學模型為

(7)

式中：r1為近似的航天器地心距，在轉移過程中為常值，不隨航天器的位置變化而改變.式(7)中，第2式的右邊第1項表示近似的線性引力，第2項表示地球中心引力和J2項攝動力，第3項為發動機推力.

針對同倫過程中的動力學模型，通過構造新的哈密爾頓函數，推導得到同倫過程中新的歐拉方程為:

(8)

(9)

(10)

相比較于原來的歐拉方程式(2)～(4)，推導得到同倫過程中的歐拉方程式(8)～(10)中，只有位置協態變量λr的積分方程式(8)發生變化，而關于速度協態變量λv的式(9)和質量協態變量λm的式(10)保持不變.

對于動力學模型式(7)，在同倫初始時刻，即ε1=0時，航天器只受到線性引力，在對其進行優化時不需要考慮J2攝動，則此時式(8)可簡化為

與包含J2項攝動時的歐拉方程相比，雖然只有λr的積分方程的非線性明顯降低，但由于λr和λv相互耦合，導致λv的非線性也相應降低，而最優推力方向只與λv相關，推力又影響航天器狀態量的變化，因此，采用線性引力假設，大幅降低了整個動力學系統的非線性，從而改善了初始問題的求解難度.

2.2 推力幅值的同倫

對于推力值過小導致的打靶不收斂問題，先在大推力幅值的情況下進行打靶求解，然后通過對推力值的大小進行同倫，得到小推力下的轉移軌道.

對動力學模型中的推力幅值，引入第2個同倫參數ε2，則在同倫過程中，推力幅值為

Tmax=(1-ε2)T1+ε2T2.

式中：T1為初始時的較大推力，T2為所需計算的小推力幅值.

在對推力大小進行同倫時，推力從大到小變化.若航天器的轉移時間自由，則隨著推力減小，航天器改變軌道狀態的能力減弱，完成相同的軌道轉移需要消耗更多的時間，轉移軌道圈數也隨之增加.文獻[9]中對于時間最優的小推力轉移軌道進行研究，在同倫過程中發現在圈數增加的某些時刻，同倫軌跡會發生突變，一般的同倫方法不能順利進行，因此提出二次同倫方法，克服了該同倫過程中存在的問題.

雖然二次同倫方法可以解決轉移圈數發生變化時產生的問題，但是所需計算時間較長，因此，在本文中，對轉移時間給定，在同倫過程中不改變轉移時間和轉移圈數，則不需要進行二次同倫方法即可完成該同倫過程.

在確定轉移時間時，先根據Lambert算法確定轉移所需的速度增量ΔV，再按照下式計算所需消耗的燃料為

Δm=m01-exp(-ΔV/(Ispg0)).

(11)

(12)

按照此原則選取得到的轉移時間，航天器該時間內能完成軌道轉移，同時不會消耗過多時間.

2.3 性能指標的同倫

針對燃料最優的Bang-Bang控制問題，根據Bertrand等[14]在2002年提出的同倫方法，通過構造從能量最優到燃料最優平滑過渡的性能指標，從而降低求解難度.

通過引入第3個同倫參數ε3，構造新的性能指標為

當ε3:0→1過程中，性能指標從能量最優轉變為燃料最優.根據最優性一階必要條件，推力大小u的控制律為

式中：q=1-ε3，當ε3≠1時，u為連續控制力，且ε3趨近于1時，即可得到燃料最優的轉移軌道.

得到性能指標為能量最優的轉移軌道之后，對同倫參數ε3按照一定步長進行迭代計算，最終可得到燃料最優的轉移軌道.

在實際求解過程中的同倫初始問題，單獨采用上述3種同倫的任意一種進行打靶求解，都存在很大難度.同時設置3個同倫參數均為0，即采用“線性引力”假設，不考慮J2項攝動，發動機推力幅值較大，性能指標為能量最優的情況下，打靶法才能快速收斂，然后再對這3個參數進行同倫.多次同倫方法的同倫過程如圖1所示，分別對同倫參數ε1、ε2和ε3依次同倫，最終得到原問題的小推力燃料最優轉移軌道.

圖1 多次同倫過程流程

3 數值仿真

以第八屆全國空間軌道設計競賽的甲題為例，航天器對位于太陽同步軌道的空間碎片進行交會，其動力學模型需考慮地球扁率J2項的攝動影響.航天器初始質量為1 000 kg，發動機最大推力Tmax=0.5 N，發動機比沖Isp=1 000 s.以半長軸為7 250 km，偏心率為0，軌道傾角為98.5°的軌道為例，航天器在該軌道上受到的地球中心引力加速度為aE=7.583 4 m/s2，由J2項攝動引起的最大引力加速度為aJ2=1.855 2×10-2m/s2，發動機推力產生的加速度為aT=5×10-4m/s2.對比可知，航天器的發動機推力遠小于受到的地球引力，甚至比J2項攝動力還要小兩個量級，若直接求解該動力學模型下的小推力最優轉移軌道，打靶法難以收斂，因此采用本文提出的多次同倫方法對其進行求解.

本文對航天器從前一個碎片出發，到后一個碎片交會的轉移軌道進行優化設計.對于本算例中的交會問題，交會精度在計算過程中體現為打靶精度，在對狀態量進行歸一化后，本文中的打靶精度設為1×10-8，即交會精度小于0.063 8m和7.905 4×10-5m/s.本算例中的計算在CPU為i5-4590，主頻為3.30 GHz，內存為8 G的臺式計算機上進行，計算平臺為MATLAB2015a，積分程序轉換為mexw32文件后利用fsolve函數進行打靶求解.

經過初步選取，得到兩個碎片的軌道要素見表1，其中歷元時刻是指對應于該軌道要素的約化儒略日.以碎片2的當前時刻為交會時刻，先將碎片1積分到交會時刻，得到碎片1在交會時刻的軌道；計算碎片1與碎片2的Lambert轉移所需速度增量ΔV，根據式(11)計算得到Δm=5.794 35 kg，再根據式(12)確定轉移時間為tf=1.973 0 d；然后對碎片1積分，得到航天器出發時刻的所在位置，從而確定轉移軌道初始時刻的狀態量、終端時刻的狀態量和轉移時間.轉移軌道初始狀態量和終端狀態量見表2.

表1 初步選取的兩個碎片的軌道要素

表2歸一化后的轉移軌道的初始狀態和終端狀態

Tab.2 The non-dimensionalized initial state and final state of the satellite

初始狀態量[-0．190071，1．111730，0．044457，0．130397，0．071664，-0．927521]終端狀態量[0．136987，0．182873，-1．107989，0．202075，-0．906749，-0．128812]

對于推力幅值的同倫，設初始時的最大推力為T1=30 N.當ε1=ε2=ε3=0時，利用MATLAB的fsolve函數進行打靶，很容易得到收斂解.先對動力學方程進行同倫，ε1從0增加到1，得到原動力學模型下的大推力能量最優轉移軌道.在此同倫過程中，同倫初始時(即ε1=ε2=ε3=0)的最優推力方向隨時間的變化曲線如圖2所示，同倫結束時(即ε1=1，ε2=ε3=0)的最優推力方向如圖3所示.從圖中可以看出，圖2中的推力方向近似于周期變化，圖3中的推力方向變化更為劇烈，這是同倫過程中動力學模型的非線性增強所導致的，同時這也可以解釋采用“線性引力”假設更容易得到收斂解的原因.

然后對推力最大幅值進行同倫.在此同倫過程中，始終保持ε1=1,ε3=0，ε2從0增加到1.對該同倫初始時刻和結束時刻的轉移軌道的推力大小隨時間的變化曲線進行對比，如圖4所示.當ε2=0時，Tmax=T1=30 N，由式(5)可知，在轉移過程中發動機推力在0和Tmax之間變化，而優化得到的推力幅值始終不大于1 N；當ε2=1時，Tmax=T2=0.5 N，在轉移過程中，推力多次達到最大值.

圖2 ε1=ε2=ε3=0時的推力方向隨時間的變化曲線

圖3 ε1=1,ε2=ε3=0時的推力方向隨時間的變化曲線

最后對性能指標進行同倫，令ε3:0→1，從而使性能指標從能力最優轉變為燃料最優.此同倫過程中的位置協態變量初值λr0和速度協態變量初值λv0隨同倫參數的變化曲線如圖5所示.從圖5中可以看出，協態變量的變化比較緩慢均勻，不存在突變跳躍的情況，同倫過程比較順利.同倫過程中航天器的燃料消耗質量變化曲線如圖6所示，從能量最優時需要消耗7.454 3 kg減小到燃料最優時的7.260 2 kg.當ε3=1時，得到燃料最優轉移軌道，此時的推力為Bang-Bang控制，開關函數和推力變化曲線如圖7所示，發動機在轉移過程中開機次數達到33次.燃料最優的小推力轉移軌道如圖8所示，航天器需要運行約28圈才能完成軌道轉移，從而對碎片2進行交會.

圖4 推力同倫過程中推力大小隨時間變化曲線

圖5性能指標同倫過程中協態變量隨同倫參數的變化曲線

Fig.5 The homotopic path of initial costates in the homotopy of index performance

圖6 性能指標同倫過程中航天器燃料消耗質量變化

Fig.6 The fuel consumption in the homotopy of index performance

本算例計算過程中，共用時1 055 s，其中在打靶過程中對動力學模型的積分時間為1 023 s，在3次同倫過程中，同倫次數分別為246次、99次、72次.與文獻[15]進行對比，其利用間接法進行優化，利用粒子群法得到打靶初值，再通過同倫得到所需解，利用C++語言編譯.其中水星交會的算例的轉移圈數為4圈，對于燃料最優的性能指標，粒子群法無法得到合適的初始值；對于能量最優，粒子群法得到初始值所需時間為338 s，然后再通過同倫得到燃料最優解，而本文中初始問題的打靶時間少于3 s；對于文獻[15]中的近地軌道長時間交會算例，轉移時間為15 d，計算用時為5.45～11.02 h，該算例與本文的不同之處在于其交會時間較長，在交會過程中未考慮J2項攝動，而本文中的動力學模型的復雜程度要遠高于該算例.綜合上述對比分析，本文所提出的多次同倫方法具有更好的普適性及更快的求解速度.

圖7 燃料最優時的開關函數和推力大小隨時間變化曲線

Fig.7 The switching function and optimal thrust profiles of fuel optimal transfer trajectory

圖8 燃料最優的轉移軌道

4 結 論

1)對動力學模型及性能指標同倫的相關公式進行推導，得到同倫過程中的最優性條件及最優控制律.

2)同倫初始問題為采用“線性引力”假設的較大推力轉移軌道優化問題，由于其非線性較弱，容易收斂，利用打靶法易于得到解，且得到的最優推力方向近似于周期變化.

3)通過引入3個同倫參數，分別對動力學模型、推力大小和性能指標進行同倫，可以克服由于推力小、J2攝動引起的模型非線性強等原因所導致的求解困難，最終得到所需的燃料最優轉移軌道.

[1] 張俊. 基于開普勒二體運動修正地球扁率J2攝動項算法[J]. 航天控制, 2014, 32(6):22-25. DOI:10.3969/j.issn.1006-3242.2014.06.005.

ZHANG Jun. The correction algorithm of J2 perturbation of the Earth oblateness based on Kepler two-body motion[J]. Aerospace Control, 2014, 32(6):22-25. DOI:10.3969/j.issn.1006-3242.2014.06.005.

[2] 張錦繡, 曹喜濱, 張劉, 等. J2項攝動下的Lambert問題算法研究[J]. 哈爾濱工業大學學報, 2009, 41(5):6-9.

ZHANG Jinxiu, CAO Xibin, ZHANG Liu, et al. An algorithm for lambert’s solution based on J2 perturbation[J]. Journal of Harbin institute of technology, 2009, 41 (5): 6-9.

[3] 王斌. 基于Lambert問題的J2攝動空間交會策略[D]. 哈爾濱：哈爾濱工業大學，2016.

WANG Bin. J2 perturbation space rendezvous strategy based onlambert problem[D]. Harbin: Harbin Institute of Technology, 2016.

[4] 張洪波, 鄭偉, 湯國建. 采用打靶法設計考慮地球扁率的機動軌道[J]. 宇航學報, 2008, 29(4):1177-1181.DOI:10.3873/j.issn.1000-1328.2008.04.017.

ZHANG Hongbo, ZHENG Wei, TANG Guojian. Maneuver trajectory design with J2 correction based on shooting procedure [J]. Journal of Astronautics, 2008, 29(4):1177-1181. DOI:10.3873/j.issn.1000-1328.2008.04.017.

[5] 周軍, 常燕. 考慮地球扁率J2攝動影響的異面橢圓軌道多沖量最優交會[J]. 宇航學報, 2008, 29(2):472-475.DOI:10.3873/j.issn.1000-1328.2008.02.015.

ZHOU Jun, CHANG Yan. Optimal multiple-impulse rendezvous between non-coplanar elliptic orbits considering the J2 perturbation effects [J]. Journal of Astronautics, 2008, 29(2):472-475.

[6] ROSCOE C W T, WESTPHAL J J, GRIESBACH J D, et al. Formation establishment and reconfiguration using differential elements in J2-perturbed orbits[J]. Journal of Guidance, Control, and Dynamics, 2015, 38(9): 1725-1740. DOI: 10.2514/1.G000999.

[7] 李俊峰, 蔣方華. 連續小推力航天器的深空探測軌道優化方法綜述[J].力學與實踐，2011,33(3):1-6.

LI Junfeng, JIANG Fanghua. Survey of low-thrust trajectory optimization methods for deep space exploration[J]. Mechanics in Engineering, 2011, 33(3):1-6.

[8] JIANG Fanghua, BAOYIN Hexi, LI Junfeng. Practical techniques for low-thrust trajectory optimization with homotopic approach[J]. Journal of Guidance, Control, and Dynamics, 2012, 35(1): 245-258. DOI: 10.2514/1.52476.

[9] PAN Binfeng, LU Ping, PAN Xun, et al. Double-homotopy method for solving optimal control problems[J]. Journal of Guidance, Control, and Dynamics, 2016, 39(8): 1706-1720. DOI: 10.2514/1.G001553.

[10]WU Baolin, WANG Danwei, POH E K, et al. Nonlinear optimization of low-thrust trajectory for satellite formation: Legendre pseudospectral approach[J]. Journal of Guidance, Control, and Dynamics, 2009, 32(4): 1371-1381. DOI: 10.2514/1.37675.

[11]SCHEEL W A, CONWAY B A. Optimization of very-low-thrust, many-revolution spacecraft trajectories[J]. Journal of Guidance, Control, and Dynamics, 1994, 17(6): 1185-1192. DOI: 10.2514/3.21331.

[12]李亮, 和興鎖, 張娟, 等. 考慮地球扁率影響的任意橢圓軌道小推力最優交會控制[J]. 西北工業大學學報, 2005, 23(3):401-405. DOI:10.3969/j.issn.1000-2758.2005.03.028.

LI Liang, HE Xingsuo, ZHANG Juan, et al. Optimal low thrust power rendezvous between neighboring elliptic orbits considering the oblateness of Earth [J]. Journal of Northwestern Polytechnical University, 2005, 23(3):401-405. DOI:10.3969/j.issn.1000-2758.2005.03.028.

[13]韓威華, 甘慶波, 王曉光, 等. 攝動情況下有限推力軌道轉移與交會聯合優化[J]. 系統工程與電子技術, 2013, 35(7):1486-1491. DOI：10.3969/j.issn.1001-506X.2013.07.22.

HAN Weihua, GAN Qingbo, WANG Xiaoguang, et al. Comprehensive optimization for orbit transfer and rendezvous of finite thrust with perturbation [J]. Systems Engineering and Electronics, 2013, 35(7):1486-1491. DOI：10.3969/j.issn.1001-506X.2013.07.22.

[14]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. DOI: 10.1002/oca.709.

[15]李俊峰, 寶音賀西, 蔣方華, 等. 深空探測動力學與控制[M]. 北京: 清華大學出版社, 2014.

LI Junfeng, BAOYIN Hexi, JIANG Fanghua, et al. Dynam ics and control of interplanetary flight [M]. Beijing: Tsinghua University Press, 2014.

Fuel-optimallowthrusttrajectorydesignwithJ2perturbation

PAN Xun1,2, PAN Binfeng1,2,TANG Shuo1

(1.School of Astronautics, Northwestern Polytechnical University, Xi’an 710072，China； 2.State Key Laboratory ofAerospace Flight Dynamics(Northwestern Polytechnical University), Xi’an 710072，China)

A multiple homotopy method is proposed to optimize the very low thrust fuel-optimal transfer trajectory under the effects of J2 perturbation. A simple problem of high thrust, energy-optimal transfers in the linear gravity without J2 perturbation is constructed as the homotopy initial problem. Three homotopic parameters are embedded in the kinetic equations, thrust magnitude and performance index, respectively, and the optimal control laws in the homotopy process are deduced according to the minimum principle. By solving the subproblems with iterative homotopic parameters, the high thrust energy-optimal transfers, low thrust energy-optimal transfers and fuel-optimal transfers are solved in turn. A numerical example about rendezvous mission of satellite and debris on sun-synchronous orbits is given to substantiate the effectiveness of the method in fuel-optimal low thrust trajectory design in the gravity with J2 perturbation. Using the proposed method, the difficulties of the highly nonlinearity of the dynamic system caused by J2 perturbation, the fuel optimal problem’s discontinuous structure of Bang-Bang control and many revolution transfers can be solved.

orbit transfer; J2 perturbation; low thrust; homotopy method; optimal control

10.11918/j.issn.0367-6234.201609030

V448.2

A

0367-6234(2017)10-0015-07

2016-09-08

國家自然科學基金(11672234)

潘 迅(1990—)，男，博士研究生；

唐 碩(1963—)，男，教授，博士生導師

潘 迅，932741825@qq.com

(編輯張 紅)