潘 迅，泮斌峰，唐 碩

(1.西北工業(yè)大學(xué) 航天學(xué)院，西安710072；2.航天飛行動(dòng)力學(xué)技術(shù)國(guó)家級(jí)重點(diǎn)實(shí)驗(yàn)室(西北工業(yè)大學(xué))，西安 710072)

考慮J2項(xiàng)攝動(dòng)的小推力燃料最優(yōu)轉(zhuǎn)移軌道設(shè)計(jì)

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

(1.西北工業(yè)大學(xué) 航天學(xué)院，西安710072；2.航天飛行動(dòng)力學(xué)技術(shù)國(guó)家級(jí)重點(diǎn)實(shí)驗(yàn)室(西北工業(yè)大學(xué))，西安 710072)

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

軌道轉(zhuǎn)移；J2項(xiàng)攝動(dòng)；小推力推進(jìn)；同倫方法；最優(yōu)控制

航天器軌道動(dòng)力學(xué)與控制是航天技術(shù)工程中的重要組成部分.在對(duì)航天器的軌道機(jī)動(dòng)進(jìn)行優(yōu)化設(shè)計(jì)時(shí)，通常將地球視為質(zhì)點(diǎn)模型，航天器在中心引力場(chǎng)中進(jìn)行機(jī)動(dòng)變軌.然而在實(shí)際飛行中，航天器受到空間環(huán)境各種攝動(dòng)因素的影響，其軌道要素發(fā)生變化.這些攝動(dòng)力包括地球形狀非球形的附加引力、大氣阻力、日月引力、太陽(yáng)光壓等攝動(dòng)力，其中地球扁率的J2項(xiàng)攝動(dòng)遠(yuǎn)大于其他攝動(dòng)因素，是引起航天器軌道變化的主要因素.

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

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

本文針對(duì)發(fā)動(dòng)機(jī)推力遠(yuǎn)小于J2項(xiàng)攝動(dòng)力的燃料最優(yōu)轉(zhuǎn)移軌道問(wèn)題，基于間接法和同倫法對(duì)其進(jìn)行研究，通過(guò)引入多個(gè)同倫參數(shù)，分別對(duì)動(dòng)力學(xué)模型、推力大小和性能指標(biāo)進(jìn)行同倫.首先采用“線性引力”假設(shè)，選取性能指標(biāo)為能量最優(yōu)，求解較大推力的轉(zhuǎn)移軌道；對(duì)動(dòng)力學(xué)模型進(jìn)行同倫，得到考慮J2攝動(dòng)的較大推力能量最優(yōu)轉(zhuǎn)移軌道；然后對(duì)推力大小進(jìn)行同倫，得到小推力的能量最優(yōu)轉(zhuǎn)移軌道；最后對(duì)性能指標(biāo)進(jìn)行同倫，克服Bang-Bang控制不連續(xù)的問(wèn)題，最終得到所需的燃料最優(yōu)轉(zhuǎn)移軌道.

1 考慮J2項(xiàng)攝動(dòng)的動(dòng)力學(xué)模型

對(duì)于在近地軌道運(yùn)行的航天器，其受到的地球扁率所引起的J2項(xiàng)攝動(dòng)力比三體引力、太陽(yáng)光壓等其他攝動(dòng)力大數(shù)個(gè)量級(jí)，是主要的攝動(dòng)力，在此僅考慮航天器受到地球中心引力和J2項(xiàng)攝動(dòng)的影響.此時(shí)，航天器的動(dòng)力學(xué)方程為

(1)

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

對(duì)于轉(zhuǎn)移軌道設(shè)計(jì)問(wèn)題，可將性能指標(biāo)表示為

結(jié)合動(dòng)力學(xué)模型(1)，根據(jù)龐德里亞金極小值原理，引入拉格朗日乘子，構(gòu)造哈密爾頓函數(shù)，有

為使哈密爾頓函數(shù)取得極小值，對(duì)其求一階偏導(dǎo)，可得歐拉方程為：

(2)

(3)

(4)

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

根據(jù)最優(yōu)控制理論，可推導(dǎo)得到最優(yōu)推力方向?yàn)?/p>

在對(duì)推力大小u的最優(yōu)控制律進(jìn)行推導(dǎo)時(shí)，需要考慮不同性能指標(biāo)的影響：

1)當(dāng)性能指標(biāo)為能量最優(yōu)時(shí)，其最小化性能指標(biāo)可以表示成

將其代入到哈密爾頓函數(shù)，并對(duì)哈密爾頓函數(shù)取極值，可得推力大小的最優(yōu)控制律為

(5)

其中S為開關(guān)函數(shù)，其表達(dá)式為

2)當(dāng)性能指標(biāo)為燃料最優(yōu)時(shí)，其最小化性能指標(biāo)可以表示成

同理可得其最優(yōu)控制律為

(6)

推導(dǎo)得到推力大小和推力方向的最優(yōu)控制律后，可通過(guò)對(duì)動(dòng)力學(xué)模型和歐拉方程進(jìn)行積分，得到狀態(tài)量和協(xié)態(tài)變量的瞬時(shí)值，再根據(jù)控制律得到推力大小和方向，從而將軌道優(yōu)化問(wèn)題轉(zhuǎn)移成兩點(diǎn)邊值問(wèn)題.

2 多次同倫方法

在將轉(zhuǎn)移軌道優(yōu)化問(wèn)題轉(zhuǎn)換為兩點(diǎn)邊值問(wèn)題之后，雖然理論上可通過(guò)打靶法對(duì)其進(jìn)行求解，然而在實(shí)際求解過(guò)程中，存在下列3個(gè)難點(diǎn)：

1)與只考慮地球中心引力的情況進(jìn)行對(duì)比，由于J2項(xiàng)攝動(dòng)的引入，動(dòng)力學(xué)模型的非線性明顯增強(qiáng)，增加了打靶收斂的難度；

2)對(duì)于小推力轉(zhuǎn)移軌道，其推力加速度一般為10-4m/s2量級(jí)，遠(yuǎn)小于地球引力，在軌道轉(zhuǎn)移過(guò)程中需要作用很長(zhǎng)時(shí)間才能改變軌道完成轉(zhuǎn)移，從而在打靶過(guò)程中容易發(fā)散；

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

針對(duì)上述3個(gè)問(wèn)題，本文提出一種求解方法，通過(guò)多次同倫，依次解決存在的問(wèn)題，從而得到最終的轉(zhuǎn)移軌道.

2.1 動(dòng)力學(xué)方程的同倫

針對(duì)包含J2項(xiàng)攝動(dòng)的動(dòng)力學(xué)模型的復(fù)雜非線性問(wèn)題，從簡(jiǎn)單動(dòng)力學(xué)模型出發(fā)，先對(duì)不包含J2項(xiàng)，且采用“線性引力”假設(shè)的動(dòng)力學(xué)模型進(jìn)行求解，得到初始解，然后將動(dòng)力學(xué)模型同倫到包含J2項(xiàng)攝動(dòng)的真實(shí)引力的動(dòng)力學(xué)模型.

“線性引力”近似是指將航天器受到的地球引力近似為大小不變、方向變化的力.引入第1個(gè)同倫參數(shù)ε1，在同倫過(guò)程中航天器的動(dòng)力學(xué)模型為

(7)

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

針對(duì)同倫過(guò)程中的動(dòng)力學(xué)模型，通過(guò)構(gòu)造新的哈密爾頓函數(shù)，推導(dǎo)得到同倫過(guò)程中新的歐拉方程為:

(8)

(9)

(10)

相比較于原來(lái)的歐拉方程式(2)～(4)，推導(dǎo)得到同倫過(guò)程中的歐拉方程式(8)～(10)中，只有位置協(xié)態(tài)變量λr的積分方程式(8)發(fā)生變化，而關(guān)于速度協(xié)態(tài)變量λv的式(9)和質(zhì)量協(xié)態(tài)變量λm的式(10)保持不變.

對(duì)于動(dòng)力學(xué)模型式(7)，在同倫初始時(shí)刻，即ε1=0時(shí)，航天器只受到線性引力，在對(duì)其進(jìn)行優(yōu)化時(shí)不需要考慮J2攝動(dòng)，則此時(shí)式(8)可簡(jiǎn)化為

與包含J2項(xiàng)攝動(dòng)時(shí)的歐拉方程相比，雖然只有λr的積分方程的非線性明顯降低，但由于λr和λv相互耦合，導(dǎo)致λv的非線性也相應(yīng)降低，而最優(yōu)推力方向只與λv相關(guān)，推力又影響航天器狀態(tài)量的變化，因此，采用線性引力假設(shè)，大幅降低了整個(gè)動(dòng)力學(xué)系統(tǒng)的非線性，從而改善了初始問(wèn)題的求解難度.

2.2 推力幅值的同倫

對(duì)于推力值過(guò)小導(dǎo)致的打靶不收斂問(wèn)題，先在大推力幅值的情況下進(jìn)行打靶求解，然后通過(guò)對(duì)推力值的大小進(jìn)行同倫，得到小推力下的轉(zhuǎn)移軌道.

對(duì)動(dòng)力學(xué)模型中的推力幅值，引入第2個(gè)同倫參數(shù)ε2，則在同倫過(guò)程中，推力幅值為

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

式中：T1為初始時(shí)的較大推力，T2為所需計(jì)算的小推力幅值.

在對(duì)推力大小進(jìn)行同倫時(shí)，推力從大到小變化.若航天器的轉(zhuǎn)移時(shí)間自由，則隨著推力減小，航天器改變軌道狀態(tài)的能力減弱，完成相同的軌道轉(zhuǎn)移需要消耗更多的時(shí)間，轉(zhuǎn)移軌道圈數(shù)也隨之增加.文獻(xiàn)[9]中對(duì)于時(shí)間最優(yōu)的小推力轉(zhuǎn)移軌道進(jìn)行研究，在同倫過(guò)程中發(fā)現(xiàn)在圈數(shù)增加的某些時(shí)刻，同倫軌跡會(huì)發(fā)生突變，一般的同倫方法不能順利進(jìn)行，因此提出二次同倫方法，克服了該同倫過(guò)程中存在的問(wèn)題.

雖然二次同倫方法可以解決轉(zhuǎn)移圈數(shù)發(fā)生變化時(shí)產(chǎn)生的問(wèn)題，但是所需計(jì)算時(shí)間較長(zhǎng)，因此，在本文中，對(duì)轉(zhuǎn)移時(shí)間給定，在同倫過(guò)程中不改變轉(zhuǎn)移時(shí)間和轉(zhuǎn)移圈數(shù)，則不需要進(jìn)行二次同倫方法即可完成該同倫過(guò)程.

在確定轉(zhuǎn)移時(shí)間時(shí)，先根據(jù)Lambert算法確定轉(zhuǎn)移所需的速度增量ΔV，再按照下式計(jì)算所需消耗的燃料為

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

(11)

(12)

按照此原則選取得到的轉(zhuǎn)移時(shí)間，航天器該時(shí)間內(nèi)能完成軌道轉(zhuǎn)移，同時(shí)不會(huì)消耗過(guò)多時(shí)間.

2.3 性能指標(biāo)的同倫

針對(duì)燃料最優(yōu)的Bang-Bang控制問(wèn)題，根據(jù)Bertrand等[14]在2002年提出的同倫方法，通過(guò)構(gòu)造從能量最優(yōu)到燃料最優(yōu)平滑過(guò)渡的性能指標(biāo)，從而降低求解難度.

通過(guò)引入第3個(gè)同倫參數(shù)ε3，構(gòu)造新的性能指標(biāo)為

當(dāng)ε3:0→1過(guò)程中，性能指標(biāo)從能量最優(yōu)轉(zhuǎn)變?yōu)槿剂献顑?yōu).根據(jù)最優(yōu)性一階必要條件，推力大小u的控制律為

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

得到性能指標(biāo)為能量最優(yōu)的轉(zhuǎn)移軌道之后，對(duì)同倫參數(shù)ε3按照一定步長(zhǎng)進(jìn)行迭代計(jì)算，最終可得到燃料最優(yōu)的轉(zhuǎn)移軌道.

在實(shí)際求解過(guò)程中的同倫初始問(wèn)題，單獨(dú)采用上述3種同倫的任意一種進(jìn)行打靶求解，都存在很大難度.同時(shí)設(shè)置3個(gè)同倫參數(shù)均為0，即采用“線性引力”假設(shè)，不考慮J2項(xiàng)攝動(dòng)，發(fā)動(dòng)機(jī)推力幅值較大，性能指標(biāo)為能量最優(yōu)的情況下，打靶法才能快速收斂，然后再對(duì)這3個(gè)參數(shù)進(jìn)行同倫.多次同倫方法的同倫過(guò)程如圖1所示，分別對(duì)同倫參數(shù)ε1、ε2和ε3依次同倫，最終得到原問(wèn)題的小推力燃料最優(yōu)轉(zhuǎn)移軌道.

圖1 多次同倫過(guò)程流程

3 數(shù)值仿真

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

本文對(duì)航天器從前一個(gè)碎片出發(fā)，到后一個(gè)碎片交會(huì)的轉(zhuǎn)移軌道進(jìn)行優(yōu)化設(shè)計(jì).對(duì)于本算例中的交會(huì)問(wèn)題，交會(huì)精度在計(jì)算過(guò)程中體現(xiàn)為打靶精度，在對(duì)狀態(tài)量進(jìn)行歸一化后，本文中的打靶精度設(shè)為1×10-8，即交會(huì)精度小于0.063 8m和7.905 4×10-5m/s.本算例中的計(jì)算在CPU為i5-4590，主頻為3.30 GHz，內(nèi)存為8 G的臺(tái)式計(jì)算機(jī)上進(jìn)行，計(jì)算平臺(tái)為MATLAB2015a，積分程序轉(zhuǎn)換為mexw32文件后利用fsolve函數(shù)進(jìn)行打靶求解.

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

表1 初步選取的兩個(gè)碎片的軌道要素

表2歸一化后的轉(zhuǎn)移軌道的初始狀態(tài)和終端狀態(tài)

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

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

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

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

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

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

最后對(duì)性能指標(biāo)進(jìn)行同倫，令ε3:0→1，從而使性能指標(biāo)從能力最優(yōu)轉(zhuǎn)變?yōu)槿剂献顑?yōu).此同倫過(guò)程中的位置協(xié)態(tài)變量初值λr0和速度協(xié)態(tài)變量初值λv0隨同倫參數(shù)的變化曲線如圖5所示.從圖5中可以看出，協(xié)態(tài)變量的變化比較緩慢均勻，不存在突變跳躍的情況，同倫過(guò)程比較順利.同倫過(guò)程中航天器的燃料消耗質(zhì)量變化曲線如圖6所示，從能量最優(yōu)時(shí)需要消耗7.454 3 kg減小到燃料最優(yōu)時(shí)的7.260 2 kg.當(dāng)ε3=1時(shí)，得到燃料最優(yōu)轉(zhuǎn)移軌道，此時(shí)的推力為Bang-Bang控制，開關(guān)函數(shù)和推力變化曲線如圖7所示，發(fā)動(dòng)機(jī)在轉(zhuǎn)移過(guò)程中開機(jī)次數(shù)達(dá)到33次.燃料最優(yōu)的小推力轉(zhuǎn)移軌道如圖8所示，航天器需要運(yùn)行約28圈才能完成軌道轉(zhuǎn)移，從而對(duì)碎片2進(jìn)行交會(huì).

圖4 推力同倫過(guò)程中推力大小隨時(shí)間變化曲線

圖5性能指標(biāo)同倫過(guò)程中協(xié)態(tài)變量隨同倫參數(shù)的變化曲線

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

圖6 性能指標(biāo)同倫過(guò)程中航天器燃料消耗質(zhì)量變化

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

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

圖7 燃料最優(yōu)時(shí)的開關(guān)函數(shù)和推力大小隨時(shí)間變化曲線

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

圖8 燃料最優(yōu)的轉(zhuǎn)移軌道

4 結(jié) 論

1)對(duì)動(dòng)力學(xué)模型及性能指標(biāo)同倫的相關(guān)公式進(jìn)行推導(dǎo)，得到同倫過(guò)程中的最優(yōu)性條件及最優(yōu)控制律.

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

3)通過(guò)引入3個(gè)同倫參數(shù)，分別對(duì)動(dòng)力學(xué)模型、推力大小和性能指標(biāo)進(jìn)行同倫，可以克服由于推力小、J2攝動(dòng)引起的模型非線性強(qiáng)等原因所導(dǎo)致的求解困難，最終得到所需的燃料最優(yōu)轉(zhuǎn)移軌道.

[1] 張俊. 基于開普勒二體運(yùn)動(dòng)修正地球扁率J2攝動(dòng)項(xiàng)算法[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項(xiàng)攝動(dòng)下的Lambert問(wèn)題算法研究[J]. 哈爾濱工業(yè)大學(xué)學(xué)報(bào), 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問(wèn)題的J2攝動(dòng)空間交會(huì)策略[D]. 哈爾濱：哈爾濱工業(yè)大學(xué)，2016.

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

[4] 張洪波, 鄭偉, 湯國(guó)建. 采用打靶法設(shè)計(jì)考慮地球扁率的機(jī)動(dòng)軌道[J]. 宇航學(xué)報(bào), 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攝動(dòng)影響的異面橢圓軌道多沖量最優(yōu)交會(huì)[J]. 宇航學(xué)報(bào), 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] 李俊峰, 蔣方華. 連續(xù)小推力航天器的深空探測(cè)軌道優(yōu)化方法綜述[J].力學(xué)與實(shí)踐，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]李亮, 和興鎖, 張娟, 等. 考慮地球扁率影響的任意橢圓軌道小推力最優(yōu)交會(huì)控制[J]. 西北工業(yè)大學(xué)學(xué)報(bào), 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]韓威華, 甘慶波, 王曉光, 等. 攝動(dòng)情況下有限推力軌道轉(zhuǎn)移與交會(huì)聯(lián)合優(yōu)化[J]. 系統(tǒng)工程與電子技術(shù), 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]李俊峰, 寶音賀西, 蔣方華, 等. 深空探測(cè)動(dòng)力學(xué)與控制[M]. 北京: 清華大學(xué)出版社, 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

國(guó)家自然科學(xué)基金(11672234)

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

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

潘 迅，932741825@qq.com

(編輯張 紅)