連續(xù)推力機動軌道優(yōu)化設(shè)計的貝塞爾曲線法

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

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

0 引 言

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

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

1 貝塞爾曲線

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

(1)

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

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

(2)

(3)

(4)

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

(5)

(6)

(7)

2 貝塞爾曲線法

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

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

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

r=rB1rH1+rB2rH2

(8)

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

(9)

(10)

2.2 約束條件

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

(11)

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

(12)

2.3 控制點設(shè)計

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

首先給出6階貝塞爾曲線B1方程為

(13)

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

(14)

(15)

(16)

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

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

(17)

P17=[rB1|s=1,θB1|s=1]=[0,θf]

(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)于控制點P12，P13，P15，P16的四個中間變量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的控制點參數(shù)可表示為

P12=[r12,θ12]=

(27)

P13=[r13,θ13]=

(28)

P15=[r15,θ15]=

(29)

P16=[r16,θ16]=

(30)

P14作為一個自由移動控制點，其參數(shù)為

P14=[r14,θ14]

(31)

同理，第二條6階貝塞爾曲線B2的7個控制點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]

(38)

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

2.4 相位角歸一化

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

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

(39)

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

(40)

2.5 速度增量計算

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

在小推力機動軌道形狀近似的方法中，方程通常定義在極坐標系下。航天器在機動軌道起止點的速度和位置必須分別與初始軌道和目標軌道一致。

軌道形狀方程為

(41)

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

(42)

(43)

(44)

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

(45)

(46)

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

(47)

(48)

(49)

(50)

(51)

(52)

極坐標系下軌道動力學(xué)方程為

(53)

(54)

(55)

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

(56)

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

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

(57)

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

(58)

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

(59)

(60)

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

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

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

3 多段貝塞爾曲線法

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

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

4 貝塞爾曲線法交會軌道設(shè)計

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

表1 多段貝塞爾曲線法控制點屬性表Table 1 List of the control points for multi-segment method

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

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

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

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

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

(61)

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

(62)

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

(63)

(2)接著，通過式(64)～(65)，利用較小的數(shù)值偏差σ生成第二代個體，通常σ取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ù)計算可得Δt(r)(2)，進而利用如式(66)～(67)進行迭代生成之后的每一代個體，直到滿足如式(68)所示的終止條件后停止，并將最后一代的ΔV(r)作為優(yōu)化算法的指標函數(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)最后，通過SA-DE-RM[29]算法的參數(shù)優(yōu)化即可獲得滿足時間約束的最優(yōu)交會軌道。

5 數(shù)值仿真

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

優(yōu)化過程中，每一代個體數(shù)設(shè)為10，迭代代數(shù)為100。貝塞爾曲線法的目標函數(shù)為minQ=ΔV，這是為了尋找燃耗最優(yōu)的轉(zhuǎn)移軌道。參數(shù)N代表轉(zhuǎn)移所需軌道周期數(shù)，T代表每一時刻所需的推力，是一個沒有任何約束的參數(shù)，完全由設(shè)計結(jié)果確定。貝塞爾曲線的前三個控制點的矢徑參數(shù)r1，r2，r3均為1，而后三個控制點的矢徑參數(shù)r4，r5，r6均為0，這是因為這6個控制點矢徑參數(shù)均受式(12)所示的約束條件限制。設(shè)計結(jié)果中兩條貝塞爾曲線的具體參數(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可以看出，控制點P24的矢徑參數(shù)r24為負數(shù)，這是因為它是一個完全自由的優(yōu)化變量，不受任何約束限制。圖5和圖6分別給出了貝塞爾曲線法設(shè)計結(jié)果與文獻[20]方法在極坐標系和直角坐標系下的軌道對比圖，盡管兩種方法都獲得了一個可行的轉(zhuǎn)移軌道，但設(shè)計結(jié)果差異明顯。圖7給出了貝塞爾曲線法推力變化曲線圖，設(shè)計結(jié)果的累積速度增量為0.2044VU，而文獻[20]方法中的設(shè)計結(jié)果為0.3231VU，可以看出貝塞爾曲線法的設(shè)計結(jié)果明顯更優(yōu)，設(shè)計結(jié)果相比于參考文獻降低了36.7%。需要注意的是文獻[30]方法在文獻[20]方法的基礎(chǔ)上通過軌道數(shù)值優(yōu)化后的結(jié)果為0.2995VU，該結(jié)果雖略有提升，但仍遠高于本文貝塞爾曲線法的設(shè)計結(jié)果。

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

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

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

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

圖7 推力變化曲線對比圖Fig.7 The comparison of the thrust profiles

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

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

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

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

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 多段貝塞爾曲線法轉(zhuǎn)移軌道設(shè)計仿真

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

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

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

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

5.3 貝塞爾曲線法交會軌道設(shè)計仿真

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

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

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

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

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

6 結(jié) 論

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