姚 瑋，羅建軍，謝劍鋒，袁建平

(1.西北工業(yè)大學(xué)航天學(xué)院，西安 710072；2. 航天飛行動(dòng)力學(xué)技術(shù)重點(diǎn)實(shí)驗(yàn)室，西安 710072；3. 北京航天飛行控制中心，北京 100094)

0 引 言

隨著電推進(jìn)技術(shù)的應(yīng)用以及全電航天器和太陽(yáng)帆航天器技術(shù)的發(fā)展，連續(xù)推力軌道在航天任務(wù)中的應(yīng)用越來(lái)越廣泛[1-8]。連續(xù)推力軌道機(jī)動(dòng)的核心問(wèn)題是軌道設(shè)計(jì)。對(duì)于連續(xù)推力軌道設(shè)計(jì)，目前已有的正方法難以滿(mǎn)足軌道狀態(tài)約束，而反方法由于轉(zhuǎn)移軌道方程局限性太強(qiáng)，結(jié)果難以令人滿(mǎn)意。反方法又叫形狀法，學(xué)者們針對(duì)形狀法做了大量研究。Petropoulos和Longuski[9]提出了一種利用指數(shù)正弦曲線(xiàn)的形狀方法來(lái)處理平面內(nèi)連續(xù)推力重力輔助軌跡。該方法只能滿(mǎn)足位置一致，而無(wú)法實(shí)現(xiàn)速度一致。Izzo[10]利用基于指數(shù)正弦曲線(xiàn)的形狀方法解決多圈蘭伯特問(wèn)題，以地火轉(zhuǎn)移問(wèn)題為例進(jìn)行了研究。對(duì)于連續(xù)推力軌道攔截問(wèn)題，Petropoulos[11]研究了很多其他的可行方法，包括對(duì)數(shù)螺旋線(xiàn)，笛卡爾卵形線(xiàn)和卡西尼卵形線(xiàn)。崔平遠(yuǎn)等[12]提出了一種深空探測(cè)任務(wù)發(fā)射機(jī)會(huì)的快速搜索算法。尚海濱等[13]針對(duì)行星際連續(xù)推力Lambert軌道問(wèn)題，提出了一種用線(xiàn)性化策略簡(jiǎn)化問(wèn)題的標(biāo)稱(chēng)軌道設(shè)計(jì)方法。Wall和Conway[14]采用5階和6階逆多項(xiàng)式形狀方法(IP)為各種航天器軌跡問(wèn)題提供近似最優(yōu)解。然而，逆多項(xiàng)式方法的最大缺陷是6階多項(xiàng)式，最多只有5個(gè)極值。因此，當(dāng)轉(zhuǎn)移軌道周期大于3時(shí)，逆多項(xiàng)式法不能對(duì)橢圓軌道機(jī)動(dòng)問(wèn)題中轉(zhuǎn)移軌道遠(yuǎn)地點(diǎn)和近地點(diǎn)的振蕩特性進(jìn)行建模。尚海濱等對(duì)Wall的6階多項(xiàng)式逼近方法開(kāi)展進(jìn)一步研究，證明其結(jié)果可以為軌道優(yōu)化過(guò)程提供有效的優(yōu)化初值[15]，并把階數(shù)提高到了15階和20階[16]。Wall等[17]又研究了一種余弦逆多項(xiàng)式方法，該方法設(shè)計(jì)下的軌道矢徑可以有更多的極值。另外，Abdelkhalik和Taheri[18-19]還提出了利用有限傅里葉級(jí)數(shù)來(lái)解決未知機(jī)動(dòng)的軌道重構(gòu)問(wèn)題的新方法。近幾年，Xie等[20]提出了一種解決共面橢圓軌道間連續(xù)推力機(jī)動(dòng)問(wèn)題的簡(jiǎn)單形狀近似方法，該方法利用多項(xiàng)式與橢圓軌道形狀方程結(jié)合的復(fù)合函數(shù)完成軌道設(shè)計(jì)。本文研究的貝塞爾曲線(xiàn)是一種可以通過(guò)控制點(diǎn)靈活設(shè)計(jì)的曲線(xiàn)。學(xué)者們通過(guò)大量研究，應(yīng)用貝塞爾曲線(xiàn)解決了機(jī)械手抓捕區(qū)域問(wèn)題[21]，無(wú)人機(jī)路徑規(guī)劃問(wèn)題[22]，軌跡生成問(wèn)題[23]，彈道設(shè)計(jì)問(wèn)題[24]，智能車(chē)軌跡跟蹤[25]，車(chē)道檢測(cè)[26]，工業(yè)無(wú)人運(yùn)輸軌跡規(guī)劃[27]和多導(dǎo)彈協(xié)同航跡規(guī)劃[28]等問(wèn)題中的研究。然而，目前還沒(méi)有研究將貝塞爾曲線(xiàn)用于航天器機(jī)動(dòng)軌道設(shè)計(jì)和描述中。

本文首先給出了貝塞爾曲線(xiàn)的基本方程，并利用復(fù)合函數(shù)對(duì)機(jī)動(dòng)軌道進(jìn)行描述，給出了相應(yīng)的約束條件和控制點(diǎn)設(shè)計(jì)，作為設(shè)計(jì)變量；然后，給出了平面切向機(jī)動(dòng)的燃耗計(jì)算，作為優(yōu)化指標(biāo)；接著，針對(duì)推力峰值過(guò)大的問(wèn)題，提出了多段式貝塞爾曲線(xiàn)法，并通過(guò)自由控制點(diǎn)的梯度下降修正解決了平面交會(huì)問(wèn)題；最后相應(yīng)的案例仿真表明新方法在連續(xù)推力機(jī)動(dòng)軌道設(shè)計(jì)中優(yōu)勢(shì)明顯。

1 貝塞爾曲線(xiàn)

軌道設(shè)計(jì)方法的核心是設(shè)計(jì)一個(gè)合理的函數(shù)作為機(jī)動(dòng)軌道的軌道方程。而貝塞爾曲線(xiàn)法的軌道方程中包含有貝塞爾曲線(xiàn)方程，一個(gè)包含n+1個(gè)控制點(diǎn)的n階貝塞爾曲線(xiàn)方程可表示為n階Bernstein多項(xiàng)式的線(xiàn)性組合[22]，如式(1)所示。

(1)

式中：Pi=[ri，θi]為貝塞爾曲線(xiàn)的控制點(diǎn)參數(shù)，(1-s)n-isi為多項(xiàng)式部分，s為貝塞爾曲線(xiàn)的比例參數(shù)，當(dāng)s從0到1改變時(shí)，貝塞爾曲線(xiàn)上的點(diǎn)對(duì)應(yīng)的從起點(diǎn)移動(dòng)至終點(diǎn)。

貝塞爾曲線(xiàn)表達(dá)式對(duì)s的各階導(dǎo)數(shù)可表述為式(2)～(4)，其中C和A分別代表組合數(shù)和排列數(shù)。

(2)

(3)

(4)

由式(2)～(4)可得貝塞爾曲線(xiàn)矢徑rB對(duì)相位方向θB的各階導(dǎo)數(shù)如式(5)～(7)所示。

(5)

(6)

(7)

2 貝塞爾曲線(xiàn)法

本文通過(guò)將初始軌道方程、目標(biāo)軌道方程與貝塞爾曲線(xiàn)結(jié)合構(gòu)成的復(fù)合函數(shù)對(duì)機(jī)動(dòng)軌道的表達(dá)式進(jìn)行設(shè)計(jì)，并給出具體的轉(zhuǎn)移任務(wù)起止點(diǎn)約束條件，然后完成控制點(diǎn)的設(shè)計(jì)以及燃耗計(jì)算。

2.1 復(fù)合函數(shù)設(shè)計(jì)

本文設(shè)計(jì)中復(fù)合函數(shù)設(shè)計(jì)為初始軌道方程H1、目標(biāo)軌道方程H2和兩個(gè)貝塞爾曲線(xiàn)B1，B2加權(quán)相加形式[20]，即轉(zhuǎn)移軌道方程如式(8)所示：

r=rB1rH1+rB2rH2

(8)

機(jī)動(dòng)軌道設(shè)計(jì)的物理特性實(shí)際上包含兩部分：首先是需要滿(mǎn)足軌道高度的逐漸提升或降低；其次是機(jī)動(dòng)軌道必須滿(mǎn)足轉(zhuǎn)移過(guò)程中橢圓機(jī)動(dòng)軌道的物理特性，即近地點(diǎn)和遠(yuǎn)地點(diǎn)的交替出現(xiàn)。而復(fù)合函數(shù)的設(shè)計(jì)中，貝塞爾曲線(xiàn)B1，B2的設(shè)計(jì)使得轉(zhuǎn)移軌道可以靈活設(shè)計(jì)與優(yōu)化，并為軌道高度提供了一個(gè)全局變化趨勢(shì)；而H1，H2作為開(kāi)普勒橢圓軌道方程則為轉(zhuǎn)移軌道引入了橢圓軌道極坐標(biāo)系下固有的震蕩特性。當(dāng)貝塞爾曲線(xiàn)參數(shù)s從0到1緩慢變化時(shí)，航天器相應(yīng)的就從轉(zhuǎn)移軌道的起點(diǎn)逐漸移動(dòng)到終點(diǎn)，在軌道高度逐漸提升或降低的過(guò)程中伴隨有橢圓軌道的周期震蕩，從而更優(yōu)地完成軌道轉(zhuǎn)移任務(wù)。通過(guò)上述的復(fù)合函數(shù)矢徑表達(dá)式可得r對(duì)θ的一階導(dǎo)數(shù)和二階導(dǎo)數(shù)如式(9)和式(10)所示：

(9)

(10)

2.2 約束條件

軌道轉(zhuǎn)移任務(wù)要求轉(zhuǎn)移軌道的位置參數(shù)r、位置對(duì)相位角一階導(dǎo)數(shù)?r/?θ和位置對(duì)相位角二階導(dǎo)數(shù)?2r/?θ2與初始軌道H1在初始位置θ0，目標(biāo)軌道H2在終止位置θf(wàn)匹配，則

(11)

將式(11)代入式(8)～ (10)，軌道轉(zhuǎn)移任務(wù)的約束條件可改寫(xiě)為式(12)的形式。

(12)

2.3 控制點(diǎn)設(shè)計(jì)

在起止點(diǎn)處，貝塞爾曲線(xiàn)參數(shù)的n階導(dǎo)數(shù)只受相鄰的n個(gè)控制點(diǎn)影響，而軌道轉(zhuǎn)移任務(wù)要求1階、2階導(dǎo)數(shù)匹配，即至少需要與起止點(diǎn)相鄰的兩個(gè)控制點(diǎn)被約束才能完成轉(zhuǎn)移任務(wù)，因此貝塞爾曲線(xiàn)一共需要至少6個(gè)控制點(diǎn)。為了使設(shè)計(jì)結(jié)果更為靈活精確，本文額外引入一個(gè)自由控制點(diǎn)來(lái)提高貝塞爾曲線(xiàn)階數(shù)，嘗試獲得更優(yōu)的軌道設(shè)計(jì)結(jié)果。

首先給出6階貝塞爾曲線(xiàn)B1方程為

(13)

6階貝塞爾曲線(xiàn)r和θ關(guān)于s的各階導(dǎo)數(shù)為

(14)

(15)

(16)

由式(12)可知，6階貝塞爾曲線(xiàn)的第一個(gè)控制點(diǎn)和最后一個(gè)控制點(diǎn)參數(shù)如式(17)～(18)所示。

P11=[rB1|s=0,θB1|s=0]=[1,0]

(17)

P17=[rB1|s=1,θB1|s=1]=[0,θf(wàn)]

(18)

由式(5)～(7)和式(14)～(15)可得，在s=0和s=1處，一階導(dǎo)數(shù)如式(19)、(21)所示，二階導(dǎo)數(shù)如式(20)、(22)所示。

(19)

(20)

(21)

(22)

引入關(guān)于控制點(diǎn)P12，P13，P15，P16的四個(gè)中間變量k12，k13，k15，k16，分別為

k12=θ12-θ11

(23)

k16=θ17-θ16

(24)

k13=θ11-2θ12+θ13

(25)

k15=θ15-2θ16+θ17

(26)

考慮如式(12)所示的約束條件，P12，P13，P15，P16的控制點(diǎn)參數(shù)可表示為

P12=[r12,θ12]=

(27)

P13=[r13,θ13]=

(28)

P15=[r15,θ15]=

(29)

P16=[r16,θ16]=

(30)

P14作為一個(gè)自由移動(dòng)控制點(diǎn)，其參數(shù)為

P14=[r14,θ14]

(31)

同理，第二條6階貝塞爾曲線(xiàn)B2的7個(gè)控制點(diǎn)P21，P22，P23，P24，P25，P26，P27的參數(shù)分別為

P21=[r21,θ21]=[0,θ0]

(32)

P22=[r22,θ22]=[0,θ21+k22]

(33)

P23=[r23,θ23]=[0,k23+2θ22-θ21]

(34)

P24=[r24,θ24]

(35)

P25=[r25,θ25]=[1,k25+2θ26-θ27]

(36)

P26=[r26,θ26]=[1,θ27-k26]

(37)

P27=[r27,θ27]=[1,θf(wàn)]

(38)

至此，貝塞爾曲線(xiàn)法的控制點(diǎn)設(shè)計(jì)就已完成。這樣，控制點(diǎn)的參數(shù)k12,k13,r14,θ14,k15,k16和k22,k23,r24,θ24,k25,k26將會(huì)作為后續(xù)軌道優(yōu)化的優(yōu)化變量。

2.4 相位角歸一化

得到控制點(diǎn)的12個(gè)參數(shù)之后，最終的轉(zhuǎn)移軌道可以通過(guò)式(8)獲得。需要注意的是，當(dāng)計(jì)算復(fù)合函數(shù)時(shí)，θB1必須與θB2一致才能做簡(jiǎn)單的疊加處理。然而，前面的兩條貝塞爾曲線(xiàn)是獨(dú)立設(shè)計(jì)的，r和θ都是由s求得，對(duì)于兩條貝塞爾曲線(xiàn)同一個(gè)s值將會(huì)求出不同的θ。這里就需要做相位歸一化處理，即貝塞爾曲線(xiàn)B2的s應(yīng)由θB1反解獲得，進(jìn)而求得第二條貝塞爾曲線(xiàn)的參數(shù)rB2，這樣就可以實(shí)現(xiàn)θ=θB1=θB2的統(tǒng)一。其中，由θB1反解s的方程可通過(guò)式(13)獲得，具體表達(dá)式如式(39)所示。

c1s6+c2s5+c3s4+c4s3+c5s2+c6s1+c7s0=θB1

(39)

其中，系數(shù)矩陣可表示為

(40)

2.5 速度增量計(jì)算

前文給出了平面機(jī)動(dòng)軌道設(shè)計(jì)優(yōu)化問(wèn)題中具體的優(yōu)化變量，而本節(jié)將給出優(yōu)化問(wèn)題的指標(biāo)函數(shù)，即累積速度增量ΔV。其具體推導(dǎo)和計(jì)算流程如下:

在小推力機(jī)動(dòng)軌道形狀近似的方法中，方程通常定義在極坐標(biāo)系下。航天器在機(jī)動(dòng)軌道起止點(diǎn)的速度和位置必須分別與初始軌道和目標(biāo)軌道一致。

軌道形狀方程為

(41)

由式(41)可求得rH關(guān)于相位角θ的各階導(dǎo)數(shù)如下

(42)

(43)

(44)

航天器的速度可分解成兩個(gè)分量：切向速度Vθ和徑向速度Vr，如式(45)、(46)所示。

(45)

(46)

飛行航跡角γ和速度大小V由式(47)～(48)給出:

(47)

(48)

(49)

(50)

(51)

(52)

極坐標(biāo)系下軌道動(dòng)力學(xué)方程為

(53)

(54)

(55)

將式(55)代入式(53)并假設(shè)δ=γ，即推力方向與速度方向保持一致，δ為發(fā)動(dòng)機(jī)安裝角，則θ對(duì)t的一階導(dǎo)數(shù)為

(56)

從式(47)～(49)和(56)可以看出，轉(zhuǎn)移任務(wù)要求航天器在機(jī)動(dòng)軌道起止點(diǎn)的位置r，位置對(duì)相位角的一階偏導(dǎo)?r/?θ和位置對(duì)相位角的二階偏導(dǎo)?2r/?θ2分別與初始軌道和目標(biāo)軌道相匹配，這就是機(jī)動(dòng)任務(wù)的約束條件。

相位角θ對(duì)時(shí)間t的二階導(dǎo)數(shù)可由式(56)推導(dǎo)獲得，表達(dá)式如下:

(57)

將式(56)～(57)代入式(47)和式(54)可得推力表達(dá)式[20]為

(58)

則相應(yīng)地做積分可知，平面軌道轉(zhuǎn)移時(shí)間Δt和累積速度增量ΔV分別為

(59)

(60)

以上得到的累積速度增量就是優(yōu)化問(wèn)題的指標(biāo)函數(shù)Q=ΔV，再結(jié)合前面機(jī)動(dòng)軌道設(shè)計(jì)給出的具體優(yōu)化變量，接下來(lái)就可以利用優(yōu)化算法通過(guò)迭代尋找最優(yōu)機(jī)動(dòng)軌道。

綜上，平面轉(zhuǎn)移軌道設(shè)計(jì)與優(yōu)化的基本流程如圖1所示，首先通過(guò)復(fù)合函數(shù)設(shè)計(jì)給出軌道形狀方程的通式；接著，分別利用復(fù)合函數(shù)和控制點(diǎn)設(shè)計(jì)給出具體軌道方程和優(yōu)化變量；最后，以累積速度增量為指標(biāo)，通過(guò)優(yōu)化給出最優(yōu)結(jié)果。

圖1 平面轉(zhuǎn)移軌道設(shè)計(jì)與優(yōu)化流程圖Fig.1 The flow chart of coplanar transfer problem

3 多段貝塞爾曲線(xiàn)法

為了降低推力峰值，本文進(jìn)一步將整個(gè)過(guò)程用多段貝塞爾曲線(xiàn)進(jìn)行描述，分段點(diǎn)參數(shù)通過(guò)優(yōu)化獲得，完成綜合最優(yōu)機(jī)動(dòng)軌道設(shè)計(jì)。

多段貝塞爾曲線(xiàn)法將在每條貝塞爾曲線(xiàn)中引入一個(gè)自由設(shè)計(jì)的分段點(diǎn)將曲線(xiàn)拆分為兩條，且在分段點(diǎn)處二階連續(xù)。機(jī)動(dòng)軌道優(yōu)化中指標(biāo)函數(shù)計(jì)算與優(yōu)化算法與前文相同，但在優(yōu)化過(guò)程中額外添加Tmax<0.018的約束(該值為文獻(xiàn)[20]設(shè)計(jì)結(jié)果中的最大推力峰值)。如表1所示，多段貝塞爾曲線(xiàn)法中包含四條貝塞爾曲線(xiàn)。其中兩條曲線(xiàn)B1和B2分別由Ba，Bb和Bc，Bd組成，這樣多段貝塞爾曲線(xiàn)法中的優(yōu)化參數(shù)共計(jì)有32個(gè)。

4 貝塞爾曲線(xiàn)法交會(huì)軌道設(shè)計(jì)

前文利用貝塞爾曲線(xiàn)復(fù)合函數(shù)實(shí)現(xiàn)了機(jī)動(dòng)軌道的設(shè)計(jì)與優(yōu)化，而這里的機(jī)動(dòng)軌道實(shí)際上是轉(zhuǎn)移軌道，只要求航天器速度矢量和位置矢量滿(mǎn)足要求，然而對(duì)機(jī)動(dòng)時(shí)間沒(méi)有任何約束。當(dāng)考慮固定機(jī)動(dòng)時(shí)間為tfixed，軌道轉(zhuǎn)移問(wèn)題就變成了軌道交會(huì)問(wèn)題，為了滿(mǎn)足這個(gè)額外的機(jī)動(dòng)時(shí)間約束，本節(jié)就將在圖1的基礎(chǔ)上，通過(guò)在軌道設(shè)計(jì)和軌道優(yōu)化過(guò)程之間引入梯度下降修正算法來(lái)實(shí)施修正，滿(mǎn)足時(shí)間匹配約束，從而實(shí)現(xiàn)交會(huì)軌道的設(shè)計(jì)與優(yōu)化。軌道機(jī)動(dòng)時(shí)間與軌道周期直接相關(guān)，周期受軌道高度影響，而在本文中軌道高度又受自由控制點(diǎn)的矢徑值影響，因此需要對(duì)自由控制點(diǎn)的矢徑值進(jìn)行修正以滿(mǎn)足交會(huì)任務(wù)。

表1 多段貝塞爾曲線(xiàn)法控制點(diǎn)屬性表Table 1 List of the control points for multi-segment method

圖2 平面交會(huì)軌道設(shè)計(jì)與優(yōu)化流程圖Fig.2 The flow chart of coplanar rendezvous problem

如圖2所示，交會(huì)軌道優(yōu)化設(shè)計(jì)的具體流程與轉(zhuǎn)移軌道設(shè)計(jì)十分相似，唯一的區(qū)別是額外的梯度下降修正算法部分。在轉(zhuǎn)移任務(wù)中，機(jī)動(dòng)軌道設(shè)計(jì)部分給出的是可行轉(zhuǎn)移軌道，機(jī)動(dòng)軌道優(yōu)化部分則是在這些可行軌道中優(yōu)化出最優(yōu)轉(zhuǎn)移軌道；而對(duì)于交會(huì)任務(wù)，首先要把機(jī)動(dòng)軌道設(shè)計(jì)部分給出的可行轉(zhuǎn)移軌道，通過(guò)梯度下降修正算法變成可行交會(huì)軌道，即實(shí)現(xiàn)機(jī)動(dòng)時(shí)間匹配，然后再通過(guò)機(jī)動(dòng)軌道優(yōu)化部分優(yōu)化出燃料最優(yōu)的交會(huì)軌道。

當(dāng)通過(guò)梯度下降對(duì)貝塞爾曲線(xiàn)法轉(zhuǎn)移軌道設(shè)計(jì)的設(shè)計(jì)結(jié)果中兩條貝塞爾曲線(xiàn)的自由控制點(diǎn)的矢徑值r14和r24進(jìn)行修正時(shí)，如圖3所示，兩個(gè)自由控制點(diǎn)的位置從P14和P24移動(dòng)到P′14和P′24，對(duì)應(yīng)的兩條貝塞爾曲線(xiàn)也相應(yīng)改變。通過(guò)這樣的修正處理使得機(jī)動(dòng)時(shí)間達(dá)到任務(wù)需要的值，從而實(shí)現(xiàn)交會(huì)任務(wù)，具體設(shè)計(jì)過(guò)程如下(1)～(4)所示。

(1)在轉(zhuǎn)移任務(wù)每個(gè)個(gè)體確定之后，根據(jù)式(61)～(63)給出梯度下降算法的第一代個(gè)體，其中下標(biāo)r和t分別代表交會(huì)任務(wù)和轉(zhuǎn)移任務(wù)。

Δt(r)(1)=Δt(t)

(61)

r14(r)(1)=r14(t)

(62)

r24(r)(1)=r24(t)

(63)

(2)接著，通過(guò)式(64)～(65)，利用較小的數(shù)值偏差σ生成第二代個(gè)體，通常σ取1即可。

r14(r)(2)=r14(r)(1)+σ

(64)

r24(r)(2)=r24(r)(1)+σ

(65)

(3)根據(jù)r14(r)(2)和r24(r)(2)以及轉(zhuǎn)移任務(wù)優(yōu)化結(jié)果的其他參數(shù)計(jì)算可得Δt(r)(2)，進(jìn)而利用如式(66)～(67)進(jìn)行迭代生成之后的每一代個(gè)體，直到滿(mǎn)足如式(68)所示的終止條件后停止，并將最后一代的ΔV(r)作為優(yōu)化算法的指標(biāo)函數(shù)。

r14(r)(n)=r14(r)(n-1)+

(66)

r24(r)(n)=r24(r)(n-1)+

(67)

|Δtfixed-Δt(r)(n)|<10-5

(68)

(4)最后，通過(guò)SA-DE-RM[29]算法的參數(shù)優(yōu)化即可獲得滿(mǎn)足時(shí)間約束的最優(yōu)交會(huì)軌道。

5 數(shù)值仿真

5.1 貝塞爾曲線(xiàn)法轉(zhuǎn)移軌道設(shè)計(jì)仿真

優(yōu)化過(guò)程中，每一代個(gè)體數(shù)設(shè)為10，迭代代數(shù)為100。貝塞爾曲線(xiàn)法的目標(biāo)函數(shù)為minQ=ΔV，這是為了尋找燃耗最優(yōu)的轉(zhuǎn)移軌道。參數(shù)N代表轉(zhuǎn)移所需軌道周期數(shù)，T代表每一時(shí)刻所需的推力，是一個(gè)沒(méi)有任何約束的參數(shù)，完全由設(shè)計(jì)結(jié)果確定。貝塞爾曲線(xiàn)的前三個(gè)控制點(diǎn)的矢徑參數(shù)r1，r2，r3均為1，而后三個(gè)控制點(diǎn)的矢徑參數(shù)r4，r5，r6均為0，這是因?yàn)檫@6個(gè)控制點(diǎn)矢徑參數(shù)均受式(12)所示的約束條件限制。設(shè)計(jì)結(jié)果中兩條貝塞爾曲線(xiàn)的具體參數(shù)和變化趨勢(shì)如表3和圖4所示，其中k12，k13，r14，θ14，k15，k16和k22，k23，r24，θ24，k25，k26均為SA-DE-RM[29]算法的優(yōu)化變量。在這些變量中，r14，θ14，r24，θ24是不受任何約束的自由變量。從圖4可以看出，控制點(diǎn)P24的矢徑參數(shù)r24為負(fù)數(shù)，這是因?yàn)樗且粋€(gè)完全自由的優(yōu)化變量，不受任何約束限制。圖5和圖6分別給出了貝塞爾曲線(xiàn)法設(shè)計(jì)結(jié)果與文獻(xiàn)[20]方法在極坐標(biāo)系和直角坐標(biāo)系下的軌道對(duì)比圖，盡管兩種方法都獲得了一個(gè)可行的轉(zhuǎn)移軌道，但設(shè)計(jì)結(jié)果差異明顯。圖7給出了貝塞爾曲線(xiàn)法推力變化曲線(xiàn)圖，設(shè)計(jì)結(jié)果的累積速度增量為0.2044VU，而文獻(xiàn)[20]方法中的設(shè)計(jì)結(jié)果為0.3231VU，可以看出貝塞爾曲線(xiàn)法的設(shè)計(jì)結(jié)果明顯更優(yōu)，設(shè)計(jì)結(jié)果相比于參考文獻(xiàn)降低了36.7%。需要注意的是文獻(xiàn)[30]方法在文獻(xiàn)[20]方法的基礎(chǔ)上通過(guò)軌道數(shù)值優(yōu)化后的結(jié)果為0.2995VU，該結(jié)果雖略有提升，但仍遠(yuǎn)高于本文貝塞爾曲線(xiàn)法的設(shè)計(jì)結(jié)果。

表2 初始軌道和目標(biāo)軌道參數(shù)Table 2 Parameters of initial orbit and final orbit

表3 貝塞爾曲線(xiàn)法控制點(diǎn)參數(shù)Table 3 Parameters of the control points for Bezier method

當(dāng)轉(zhuǎn)移周期數(shù)N從0到5改變時(shí)，多案例仿真下，貝塞爾曲線(xiàn)法總是可以獲得一個(gè)優(yōu)于文獻(xiàn)設(shè)計(jì)結(jié)果的轉(zhuǎn)移軌跡，累積速度增量降低了15%～40%，具體結(jié)果對(duì)比見(jiàn)表4。極坐標(biāo)系和直角坐標(biāo)系下軌道圖如圖8和圖9所示?？梢钥闯觯D(zhuǎn)移過(guò)程中軌道總是保持在初始軌道和目標(biāo)軌道之間，轉(zhuǎn)移的起止點(diǎn)處，航天器的位置、速度都能與初始軌道和目標(biāo)軌道完美匹配。圖9中任意兩個(gè)終止點(diǎn)間距都為2π，對(duì)應(yīng)一個(gè)軌道周期。

圖4 貝塞爾曲線(xiàn)法設(shè)計(jì)結(jié)果圖Fig.4 The details of the Bezier curves

圖7 推力變化曲線(xiàn)對(duì)比圖Fig.7 The comparison of the thrust profiles

推力變化情況如圖10所示。在N=0時(shí)，貝塞爾曲線(xiàn)法在累積速度增量ΔV和最大推力峰值Tmax方面均優(yōu)于文獻(xiàn)。而當(dāng)N>0時(shí)，貝塞爾曲線(xiàn)法總是得到一個(gè)更小的累積速度增量ΔV，但同時(shí)推力峰值Tmax也更大，然而從圖10可以看出，貝塞爾曲線(xiàn)法僅僅在推力峰值附近，推力更高，其他時(shí)刻推力均更小，因此最終的累積速度增量?jī)?yōu)于文獻(xiàn)[20]。另外，需要注意，當(dāng)N>3時(shí)，IP方法[14]無(wú)解，而本文方法和文獻(xiàn)[20]方法均可獲取有效的解。

圖8 多案例仿真中極坐標(biāo)系下軌道圖Fig.8 Radius profiles when N=0～5

圖9 多案例仿真中直角坐標(biāo)下軌道圖Fig.9 The transfer orbits when N=0～5

圖10 多案例仿真中推力曲線(xiàn)圖Fig.10 The thrust profiles when N=0～5

N貝塞爾曲線(xiàn)法文獻(xiàn)[20]方法ΔVΔtTmaxΔVΔtTmax00.655312.89380.42590.884413.66870.980110.193768.76250.04770.288555.26030.029320.199194.75160.03900.3076110.28790.016730.2044170.47600.04030.3231162.88760.018040.2287191.71620.03660.3590215.35930.015850.2417263.46450.03150.4031267.82110.0132

5.2 多段貝塞爾曲線(xiàn)法轉(zhuǎn)移軌道設(shè)計(jì)仿真

多段貝塞爾曲線(xiàn)法仿真中，參數(shù)設(shè)置均與表2相同，唯一區(qū)別是添加了最大推力峰值約束(仿真中取值為0.018DU/TU2)。多段貝塞爾曲線(xiàn)法設(shè)計(jì)結(jié)果參數(shù)如表5所示，圖11給出了具體的四條貝塞爾曲線(xiàn)形式，仿真結(jié)果中，累積速度增量為0.1535VU，機(jī)動(dòng)持續(xù)時(shí)間為169.8TU，最大推力峰值為0.013DU/TU2。與表4結(jié)果對(duì)比可以看出多段貝塞爾曲線(xiàn)法無(wú)論累積速度增量和最大推力峰值都遠(yuǎn)小于文獻(xiàn)[20]方法，其中累積速度增量?jī)H為文獻(xiàn)的47.5%，這是因?yàn)槎喽畏ㄔ诮档妥畲笸屏Ψ逯档耐瑫r(shí)，由于分段造成的控制點(diǎn)增多，實(shí)際相當(dāng)于階數(shù)也得以提升，因此設(shè)計(jì)過(guò)程也更加靈活，從而可以獲得更優(yōu)的累積速度增量。

圖12、圖13給出多段貝塞爾曲線(xiàn)法的軌道圖。當(dāng)本文在仿真中將進(jìn)化代數(shù)增加至1000時(shí)，累積速度增量可達(dá)到0.135VU(文獻(xiàn)[20]中為0.3231VU)，機(jī)動(dòng)時(shí)間和最大推力峰值分別為178.75TU，0.0136DU/TU2?？梢钥吹竭M(jìn)一步優(yōu)化的累積速度增量?jī)H為文獻(xiàn)[20]的41%，同時(shí)最大推力需求和機(jī)動(dòng)時(shí)間也都更小，設(shè)計(jì)結(jié)果全面優(yōu)于文獻(xiàn)[20]方法。

表5 多段貝塞爾曲線(xiàn)法控制點(diǎn)參數(shù)Table 5 Parameters of the control points for multi-segment Bezier method

圖11 多段貝塞爾曲線(xiàn)法設(shè)計(jì)結(jié)果圖Fig.11 The details of the multi-segment Bezier curves

5.3 貝塞爾曲線(xiàn)法交會(huì)軌道設(shè)計(jì)仿真

當(dāng)任務(wù)約束要求機(jī)動(dòng)時(shí)間Δtfixed=180TU時(shí)，貝塞爾曲線(xiàn)法交會(huì)軌道設(shè)計(jì)與貝塞爾曲線(xiàn)法轉(zhuǎn)移軌道設(shè)計(jì)的仿真結(jié)果對(duì)比如表6所示，可以看出交會(huì)任務(wù)結(jié)果滿(mǎn)足機(jī)動(dòng)時(shí)間約束。圖14為最優(yōu)交會(huì)軌道與最優(yōu)轉(zhuǎn)移軌道對(duì)比圖，可以看出經(jīng)過(guò)梯度下降修正之后，軌道發(fā)生了改變以滿(mǎn)足交會(huì)任務(wù)時(shí)間需求。

圖12 多段貝塞爾曲線(xiàn)法設(shè)計(jì)結(jié)果極坐標(biāo)系下軌道圖Fig.12 Radius profiles of multi-segment Bezier method

圖13 多段貝塞爾曲線(xiàn)法設(shè)計(jì)結(jié)果直角坐標(biāo)系下軌道圖Fig.13 The transfer orbits of multi-segment Bezier method

圖14 轉(zhuǎn)移軌道和交會(huì)軌道直角坐標(biāo)系下軌道對(duì)比圖Fig.14 The comparison of transfer and rendezvous orbits

參數(shù)轉(zhuǎn)移任務(wù)交會(huì)任務(wù)是否滿(mǎn)足時(shí)間約束否是累積速度增量ΔV0.2044VU0.2059VU任務(wù)時(shí)間Δt170.4760TU180.0000TU最大推力需求Tmax0.0403VU/TU20.0401VU/TU2

6 結(jié) 論

本文提出了基于貝塞爾曲線(xiàn)的機(jī)動(dòng)軌道設(shè)計(jì)方法。首先，將貝塞爾曲線(xiàn)的基本特點(diǎn)與軌道機(jī)動(dòng)問(wèn)題相結(jié)合，提出用貝塞爾曲線(xiàn)與軌道形狀方程的復(fù)合函數(shù)描述和設(shè)計(jì)機(jī)動(dòng)軌道的方法，并給出相應(yīng)的約束條件，控制點(diǎn)設(shè)計(jì)和累計(jì)速度增量計(jì)算，所設(shè)計(jì)的控制點(diǎn)和計(jì)算得到的累計(jì)速度增量即為后續(xù)優(yōu)化過(guò)程的優(yōu)化變量與指標(biāo)函數(shù)。仿真結(jié)果表明貝塞爾曲線(xiàn)法的設(shè)計(jì)結(jié)果遠(yuǎn)優(yōu)于其他方法，但仍存在推力峰值過(guò)大的問(wèn)題。進(jìn)而，論文通過(guò)自由分段點(diǎn)的引入提出了多段式貝塞爾曲線(xiàn)法，在解決推力峰值過(guò)大問(wèn)題的基礎(chǔ)上，進(jìn)一步優(yōu)化了最優(yōu)指標(biāo)，最終燃耗僅為對(duì)比方法的41%。在此基礎(chǔ)上，通過(guò)梯度下降算法對(duì)自由控制點(diǎn)進(jìn)行修正來(lái)調(diào)整機(jī)動(dòng)時(shí)間，解決了固定機(jī)動(dòng)時(shí)間約束下的軌道交會(huì)問(wèn)題。仿真結(jié)果證明本文提出的貝塞爾曲線(xiàn)法能更優(yōu)地解決平面連續(xù)推力軌道轉(zhuǎn)移和交會(huì)問(wèn)題。