吳愛國，董 西

(哈爾濱工業(yè)大學(xué)(深圳)機(jī)電工程與自動(dòng)化學(xué)院· 深圳·518055)

0 引 言

深空探測(cè)是指人類對(duì)月球、火星及其他天體或空間環(huán)境所進(jìn)行的探測(cè)活動(dòng)，是航天領(lǐng)域重要的發(fā)展趨勢(shì)之一[1]。深空探測(cè)活動(dòng)不僅可以幫助人類了解宇宙和太陽系的發(fā)展和演變，而且在進(jìn)行空間資源開發(fā)和技術(shù)創(chuàng)新方面有著重要的科學(xué)和經(jīng)濟(jì)價(jià)值。深空軌道的設(shè)計(jì)和優(yōu)化方法不僅決定著航天器能否成功轉(zhuǎn)移，也會(huì)影響轉(zhuǎn)移過程中航天器的轉(zhuǎn)移時(shí)間、能量消耗等，是行星探索活動(dòng)中的一項(xiàng)關(guān)鍵技術(shù)。

改進(jìn)的逆五階和六階多項(xiàng)式法[5]改善了逆多項(xiàng)式形狀法無法實(shí)現(xiàn)多圈轉(zhuǎn)移的缺點(diǎn)。文獻(xiàn)[5]在二體模型中將長(zhǎng)半軸表示為逆五階或逆六階多項(xiàng)式的形式，通過在初始和末端位置引入一個(gè)單調(diào)系數(shù)，滿足使得軌跡半徑始終增大的要求，以滿足軌跡進(jìn)行安全轉(zhuǎn)移，適用于小推力共面軌道的多圈轉(zhuǎn)移問題。但針對(duì)交會(huì)時(shí)間短的情況而言，可能不存在解。同時(shí)，隨著轉(zhuǎn)移圈數(shù)的增多，速度增量也會(huì)越來越大，存在不收斂的問題。Petropoulos等人提出的指數(shù)正弦形狀法[6]是一種含有4個(gè)參數(shù)因子的初始設(shè)計(jì)方法，但是后驗(yàn)決定的推力并不能由設(shè)計(jì)任務(wù)精確描述。文獻(xiàn)[7]將形狀法應(yīng)用于速度向量而不是徑向向量上，對(duì)速度向量進(jìn)行有限Fourier分解，然后對(duì)速度向量進(jìn)行積分并得到徑向向量的表達(dá)式，再結(jié)合遺傳算法搜索到合適的初始猜測(cè)值。文獻(xiàn)[8]根據(jù)火星探測(cè)器的初始與末端的位置、速度矢量和轉(zhuǎn)移時(shí)間，利用基于逆六階多項(xiàng)式的二維形狀法建立了火星探測(cè)器的異面轉(zhuǎn)移軌跡。在小推力交會(huì)問題中，也可以將末端速度分量作為優(yōu)化參數(shù)來求解所要求的推力加速度和推力指向角β。文獻(xiàn)[9]對(duì)指數(shù)正弦形狀法的收斂性和偽二分法的可行性進(jìn)行了分析，針對(duì)小推力重力協(xié)助問題，提出了一種新的算法并進(jìn)行了最優(yōu)分析。文獻(xiàn)[10]分析了利用指數(shù)正弦形狀法求解小推力固定時(shí)間軌道轉(zhuǎn)移情況下的Lambert問題的動(dòng)力學(xué)可行性，證明了推力方向角的正切值在某一范圍內(nèi)且推力大小無限制的情況下是可行的。

相比形狀法，基于Fourier級(jí)數(shù)的設(shè)計(jì)方法可以對(duì)大量的自由參數(shù)進(jìn)行求解，通過找到相應(yīng)的Fourier展開式系數(shù)來找到滿足限制條件的解。其缺點(diǎn)是，當(dāng)沒有解時(shí)，有限制條件的Fourier變換問題需要更多的時(shí)間來確定沒有解，使得計(jì)算時(shí)間大大增加。文獻(xiàn)[11]在球坐標(biāo)系下建立了基于Fourier級(jí)數(shù)的近似模型進(jìn)行初始設(shè)計(jì)，再利用遺傳算法通過種群的迭代來得到轉(zhuǎn)移軌道的初始猜測(cè)值，使得即使在Fourier展開式系數(shù)很少的情況下，也可以得到次優(yōu)的火星轉(zhuǎn)移的三維初始軌跡，節(jié)省了計(jì)算時(shí)間。高斯偽譜法[12]和各種優(yōu)化控制方法[13-14]也被應(yīng)用于求解航天器的轉(zhuǎn)移軌跡。

本文提出了一種基于有限Fourier級(jí)數(shù)的形狀法來進(jìn)行航天器轉(zhuǎn)移軌跡的初始和優(yōu)化設(shè)計(jì)，將轉(zhuǎn)移軌跡中徑向大小和轉(zhuǎn)移角之間的關(guān)系用有限Fourier級(jí)數(shù)展開，根據(jù)已知的初始和末端參數(shù)進(jìn)行兩點(diǎn)邊值問題的求解。在初始設(shè)計(jì)中，選擇最小的有限Fourier級(jí)數(shù)項(xiàng)數(shù)求解轉(zhuǎn)移軌跡的初始猜測(cè)值。在優(yōu)化設(shè)計(jì)中，將Fourier級(jí)數(shù)的未知項(xiàng)數(shù)、轉(zhuǎn)移時(shí)間、最大推力加速度等作為約束條件，將最小速度增量作為目標(biāo)函數(shù)，利用所得的初始猜測(cè)值進(jìn)行小推力多約束多圈轉(zhuǎn)移問題的求解。該方法的創(chuàng)新性表現(xiàn)為如下幾個(gè)方面：

(1)利用有限Fourier級(jí)數(shù)求解了航天器小推力多圈轉(zhuǎn)移軌跡的轉(zhuǎn)移徑向大小和轉(zhuǎn)移角之間的關(guān)系，有別于將轉(zhuǎn)移軌跡徑向和轉(zhuǎn)移時(shí)間、轉(zhuǎn)移角度和轉(zhuǎn)移時(shí)間之間的關(guān)系進(jìn)行分開研究的方法；

(2)提出的基于有限Fourier級(jí)數(shù)的設(shè)計(jì)方法比較了在不同參數(shù)取值下的所得出的結(jié)果，利用數(shù)值結(jié)果形象直觀地說明了轉(zhuǎn)移軌道初始設(shè)計(jì)的好壞對(duì)優(yōu)化設(shè)計(jì)的影響；

(3)將提出的方法與逆五階多項(xiàng)式形狀法進(jìn)行了比較，在利用更小的最大推力加速度的條件下，本文所提出的方法可以減少航天器在轉(zhuǎn)移過程中75%的速度增量。

1 數(shù)學(xué)模型的建立

航天器在極坐標(biāo)下的運(yùn)動(dòng)模型為

(1)

式中，r為徑向大小，θ為極角，μ為地球引力參數(shù)，γ為飛行路徑角，Ta為推力加速度大小，α為推力方向角。徑向r對(duì)轉(zhuǎn)移時(shí)間t求一階導(dǎo)數(shù)可得

(2)

(3)

將飛行路徑角的正切值tanγ對(duì)極角θ進(jìn)行一階求導(dǎo)和二階求導(dǎo)，分別可得到

(4)

(5)

(6)

由式(1)和式(6)，可得

(7)

將式(7)代入式(1)，可得

(8)

假定推力方向是與速度方向一致的正切推力，即滿足α=γ+nπ，可得角速度滿足式(9)

(9)

(10)

代入航天器的運(yùn)動(dòng)模型式(1)，可得正切推力加速度為

(11)

火星探測(cè)器到達(dá)最終位置所需要的飛行時(shí)間tf和速度增量ΔV可以表示為式(12)和式(13)

(12)

(13)

2 設(shè)計(jì)方法

2.1 轉(zhuǎn)移軌跡的初始設(shè)計(jì)

本節(jié)采用三階有限Fourier級(jí)數(shù)來模擬正切小推力作用下所求得的轉(zhuǎn)移軌跡中航天器的徑向大小和轉(zhuǎn)移角之間的關(guān)系。假設(shè)在不存在時(shí)間約束的情況下，通過初始和末端的徑向大小、轉(zhuǎn)移角度、切向角六個(gè)參數(shù)來設(shè)計(jì)初始軌跡形狀。由于初始項(xiàng)未知數(shù)的存在，三階有限Fourier級(jí)數(shù)中含有7個(gè)未知數(shù)，則存在一個(gè)不確定的參數(shù)取值，本文選擇并比較了四種不同參數(shù)進(jìn)行初始設(shè)計(jì)時(shí)帶來的不同的結(jié)果。

三階有限Fourier級(jí)數(shù)的表示形式如下

(14)

其中，a0,ai,bi,i=1，2，3是需要根據(jù)初始和末端條件求得的7個(gè)未知數(shù)。將初始位置和末端位置的轉(zhuǎn)移角度的取值θ=θ0(θ0=0)和θ=θf分別代入式(14)及其一階、二階導(dǎo)數(shù)后，得到下式

(15)

(16)

由式(3)和式(9)，將徑向大小r對(duì)轉(zhuǎn)移角θ求一階和二階導(dǎo)數(shù)后可得

(17)

當(dāng)θ=0時(shí)，即當(dāng)航天器處于初始位置時(shí)，可得

(18)

當(dāng)θ=θf時(shí)，即當(dāng)航天器處于末端位置時(shí)，可得

(19)

為簡(jiǎn)化運(yùn)算的表達(dá)式，令

(20)

由式(15)～式(20)整理得到僅含a0,ai,bi,i=1，2，3的六個(gè)等式，利用解方程組的方法可以得到

(21)

剩下的3個(gè)待求未知數(shù)可以表示為

b1+2b2+3b3=C

-b1+2b2-3b3=D

(22)

可以得知，b2有唯一解b2=(C+D)/4。本文所提出的方法(簡(jiǎn)寫為FFS-3初始設(shè)計(jì))在進(jìn)行航天器小推力多圈轉(zhuǎn)移軌跡的設(shè)計(jì)時(shí)，應(yīng)滿足等式b1+3b3=(C-D)/4。為了計(jì)算的簡(jiǎn)便，選取了與式(20)相關(guān)的4組參數(shù)值，在第4節(jié)中將利用這4組滿足條件的參數(shù)值進(jìn)行仿真。

2.2 轉(zhuǎn)移軌跡的優(yōu)化設(shè)計(jì)

在航天器小推力多圈轉(zhuǎn)移過程中，轉(zhuǎn)移軌跡呈螺旋狀且具有周期性。根據(jù)法國數(shù)學(xué)家Fourier所提出的理論可知，任何周期函數(shù)都可以用正弦函數(shù)和余弦函數(shù)構(gòu)成的無窮級(jí)數(shù)來表示。根據(jù)初始設(shè)計(jì)所得出的初始迭代值，選擇恰當(dāng)?shù)某跏疾聹y(cè)值進(jìn)行航天器小推力多圈轉(zhuǎn)移軌跡的優(yōu)化設(shè)計(jì)。將有限Fourier級(jí)數(shù)的未知項(xiàng)數(shù)的個(gè)數(shù)、轉(zhuǎn)移時(shí)間和最大正切小推力等作為約束條件，以速度增量最小作為優(yōu)化設(shè)計(jì)的目標(biāo)函數(shù)，即將轉(zhuǎn)移過程中速度增量的最小值作為優(yōu)化目標(biāo)，得出多約束條件下的小推力多圈轉(zhuǎn)移軌跡。基于有限Fourier級(jí)數(shù)的形狀法可將航天器在轉(zhuǎn)移過程中的徑向大小r和轉(zhuǎn)移角θ之間的關(guān)系表示為

(23)

對(duì)式(23)中的徑向大小r對(duì)轉(zhuǎn)移角度θ分別求一階導(dǎo)數(shù)和二階導(dǎo)數(shù)，并得到結(jié)果如下

(24)

將初始位置和末端位置的轉(zhuǎn)移角度的取值θ=θ0(θ0=0)和θ=θf分別代入式(23)和式(24)，可得

(25)

(26)

根據(jù)式(25)和式(26)確定有限Fourier級(jí)數(shù)的6個(gè)系數(shù)，即還存在2nr-5個(gè)系數(shù)沒有被確定。利用所提出的優(yōu)化方法(簡(jiǎn)寫為FFS-nr優(yōu)化設(shè)計(jì))，聯(lián)立式(14)和式(23)可得關(guān)系式

將轉(zhuǎn)移角度θ等分為nr份，將M作為三角函數(shù)中角度的系數(shù)，將FFS-nr表達(dá)式中的近似系數(shù)用矩陣來表示

X=F-1Y

其中

X=[a0,a1,a2…anr,b1,b2…bnr]T

F=

利用猜測(cè)矩陣可以輔助求解余下的2nr-5個(gè)未知的系數(shù)。

3 仿真試驗(yàn)及結(jié)果分析

3.1 仿真結(jié)果

選取初始發(fā)射點(diǎn)的軌道根數(shù)為a1=3DU,e1=0.6,ω1=30°，末端終點(diǎn)的軌道根數(shù)為a2=6DU,e2=0.6,ω2=30°，并給定初始位置為330°，終端位置為120°，最大轉(zhuǎn)移圈數(shù)為Nrev=5。通過MATLAB進(jìn)行2.1節(jié)中初始設(shè)計(jì)的仿真，得出4種參數(shù)的取值情況所計(jì)算出的不同參數(shù)下的轉(zhuǎn)移速度增量ΔV、轉(zhuǎn)移時(shí)間tf和最大正切小推力加速度Ta(max)的值如表1(FFS-3初始設(shè)計(jì)仿真結(jié)果)所示。假設(shè)2.2節(jié)優(yōu)化設(shè)計(jì)中有限Fourier級(jí)數(shù)的項(xiàng)數(shù)為nr=10，以轉(zhuǎn)移時(shí)間大于200TU小于1200TU、最大推力加速度小于0.02(DU/TU2)為約束條件，以最小速度增量ΔV作為目標(biāo)函數(shù)進(jìn)行了優(yōu)化，目標(biāo)函數(shù)fmin為

fmin=minΔV

約束條件為

利用MATLAB中的Fmincon函數(shù)得到的仿真結(jié)果如表1(FFS-10優(yōu)化設(shè)計(jì)仿真結(jié)果)所示。

選取相同的參數(shù)、利用逆五階多項(xiàng)式形狀法進(jìn)行仿真，并與表1中b1=C/2,b3=-D/6的結(jié)果進(jìn)行比較，結(jié)果見表2。

表1 基于Fourier級(jí)數(shù)的小推力多圈轉(zhuǎn)移軌跡的設(shè)計(jì)(Nrev=5,nr=10)

表2 本文提出的設(shè)計(jì)方法與逆五階多項(xiàng)式形狀法的比較結(jié)果

三種設(shè)計(jì)方法所得到的正切推力加速度變化如圖1所示。

(a)逆五階多項(xiàng)式形狀法

3.2 結(jié)果分析

從表1中的數(shù)據(jù)比較可以發(fā)現(xiàn)，不同的參數(shù)取值會(huì)導(dǎo)致不同的初始設(shè)計(jì)結(jié)果，而由初始設(shè)計(jì)得到的不同軌跡表達(dá)式的初始猜測(cè)值將會(huì)影響所求目標(biāo)軌跡的優(yōu)化設(shè)計(jì)。在進(jìn)行航天器轉(zhuǎn)移軌跡的設(shè)計(jì)時(shí)，相對(duì)精確的初始猜測(cè)值不僅可以減少轉(zhuǎn)移時(shí)間，還能利用較小的推力值實(shí)現(xiàn)轉(zhuǎn)移目標(biāo)，并利用較小的速度增量達(dá)到轉(zhuǎn)移的目的。由表2可知，相比逆五階多項(xiàng)式設(shè)計(jì)方法，本文提出的基于有限Fourier項(xiàng)數(shù)的航天器小推力轉(zhuǎn)移軌跡形狀法的速度增量更小。當(dāng)航天器小推力多圈轉(zhuǎn)移的轉(zhuǎn)移圈數(shù)為5圈、有限Fourier級(jí)數(shù)的未知項(xiàng)數(shù)為10時(shí)，可以減少近75%的轉(zhuǎn)移速度增量。

4 結(jié) 論

本文提出了一種基于有限Fourier級(jí)數(shù)的形狀法，描述了航天器轉(zhuǎn)移軌跡的徑向大小和轉(zhuǎn)移角之間的關(guān)系，并進(jìn)行了航天器小推力多圈轉(zhuǎn)移軌跡的設(shè)計(jì)。通過選擇不同的參數(shù)值，完成了小推力航天器多圈轉(zhuǎn)移軌跡的初始設(shè)計(jì)和優(yōu)化設(shè)計(jì)，從數(shù)值上比較了由轉(zhuǎn)移軌跡初始設(shè)計(jì)得出的初始猜測(cè)值對(duì)優(yōu)化設(shè)計(jì)結(jié)果的影響，得出了精確的初始猜測(cè)值可以實(shí)現(xiàn)更佳轉(zhuǎn)移目標(biāo)的結(jié)論。但是，本文未考慮具有時(shí)間限制的轉(zhuǎn)移軌跡的設(shè)計(jì)問題，也沒有考慮非共面、非同軸、具有不同離心率等情況下的轉(zhuǎn)移軌跡初始設(shè)計(jì)和優(yōu)化問題，因此后續(xù)將繼續(xù)基于本文提出的有限Fourier級(jí)數(shù)形狀法進(jìn)行進(jìn)一步研究。