全電推進衛(wèi)星軌道優(yōu)化的推力同倫解法

中國空間技術(shù)研究院，北京 100094

全電推進(all-electric propulsion)衛(wèi)星是指不僅在入軌后的位置保持采用電推進技術(shù)，而且在火箭與衛(wèi)星分離后到上升至最終工作軌道的過程也采用電推進變軌的衛(wèi)星[1]。由于電推進的比沖相比傳統(tǒng)化學(xué)推進要高近一個量級，因此在相同通信容量的情況下，不配備化學(xué)推進系統(tǒng)的全電推進通信衛(wèi)星的發(fā)射質(zhì)量可以大大減小，或者在相同發(fā)射質(zhì)量的情況下大幅提高有效載荷質(zhì)量。美國波音衛(wèi)星系統(tǒng)公司于2015年發(fā)射的Boeing 702SP平臺是世界上首顆靜止軌道全電推進衛(wèi)星[2]，歐洲空中客車公司制造的全電推進衛(wèi)星EutelSat 172B也于2017年成功發(fā)射[3]。目前世界范圍有數(shù)十顆全電推進衛(wèi)星訂單，中國也開始了全電推進平臺的關(guān)鍵技術(shù)研究及平臺開發(fā)工作[2]。

由于電推進系統(tǒng)的推力非常小，通常只有幾十到數(shù)百毫牛[4]，從同步轉(zhuǎn)移軌道發(fā)射至進入最終靜止軌道的過程通常需要幾個月到半年的時間[5]，需繞行地球數(shù)百圈，此過程中軌道要素和控制量持續(xù)周期性波動變化，這也造成其最優(yōu)軌道求解難度比其它小推力問題更大。

連續(xù)小推力變軌是全電推進平臺最關(guān)鍵的技術(shù)之一，對于連續(xù)推力軌跡優(yōu)化問題的求解，主要可以分為間接法、直接法、混合法，以及狀態(tài)反饋方法[6-8]。間接法可以獲得很高的精度，且尋優(yōu)變量少，但是協(xié)態(tài)變量無實際物理意義，對初值猜測準(zhǔn)確性要求極高，收斂半徑小[6]。直接法對初值要求低，但是尋優(yōu)變量規(guī)模很大，且求解精度比間接法低，容易收斂到局部最優(yōu)解[6]。混合法結(jié)合了間接法和直接法的優(yōu)點，而狀態(tài)反饋法采用反饋控制理論，并不是建立在最優(yōu)控制理論之上，因此只是滿足工程要求的簡化控制方法。由于間接法對初值猜測的敏感性，使其在全電推進入軌的多圈變軌求解中難以應(yīng)用，而直接法因多圈轉(zhuǎn)移時間長而存在大量的待優(yōu)化參數(shù)，因此現(xiàn)有的方法主要是通過軌道平均法或狀態(tài)反饋控制率的方式求解。文獻[7]采用了春分點軌道根數(shù)平均軌道方程的混合法，完成了地球衛(wèi)星多圈轉(zhuǎn)移優(yōu)化。文獻[8]提出了基于Lyapunov反饋函數(shù)的Q-Law方法，給出了簡單實用的反饋控制率，文獻[9-10]的方法也可以歸為這一類。其中平均法的結(jié)果是對多圈取平均得到的軌道要素平根數(shù)，無法獲取在一圈內(nèi)的軌道細節(jié)信息，而近似反饋控制率不是建立在最優(yōu)理論之上的，因此其結(jié)果距離最優(yōu)解往往還有較大改進空間，這些方法都無法像間接法那樣給出精確的最優(yōu)解。

在間接法求解燃料最優(yōu)控制問題方面，文獻[11-12]采用同倫思想，取得了很好的效果。同倫方法從容易求解的問題出發(fā)，通過迭代的方式逼近難求解的問題。需要指出的是，目前在同倫方法的應(yīng)用上，現(xiàn)有文獻主要是在指標(biāo)函數(shù)上進行漸變，從連續(xù)的能量最優(yōu)問題逐漸逼近到最優(yōu)的非連續(xù)Bang-Bang控制問題，而在時間最優(yōu)問題方面還少有應(yīng)用。本文借鑒燃料最優(yōu)問題中的同倫思想，采用推力逐漸縮小的方式完成間接法求解，可獲得全電推進衛(wèi)星最優(yōu)時間變軌問題的理論精確解。

1 數(shù)學(xué)模型

1.1 基于間接法的小推力優(yōu)化

對于偏心率近似為零的圓軌道，以及近似于零傾角的衛(wèi)星而言，近地點幅角、升交點赤經(jīng)將失去意義，此時不宜用傳統(tǒng)的軌道六要素來描述衛(wèi)星的運動，對于靜止軌道衛(wèi)星，可采用春分點軌道根數(shù)[13]即 (p,f,g,h,k,L)，其表達式如下：

(1)

式中：a,e,i,Ω,ω,v為經(jīng)典軌道六要素的半長軸、偏心率、傾角、升交點赤經(jīng)、近地點幅角和真近點角。

令x=[p,f,g,h,k,L]T，衛(wèi)星變軌時的推力為T，質(zhì)量為m，推力器的比沖為Isp，g0為地球海平面加速度，推力方向在軌道的徑向、切向、法向這3個方向的余弦為u=[ur,uτ,un]T，則描述衛(wèi)星運動的微分方程可寫為

(2)

(3)

式(2)中的M表達式為[13]:

由于改進春分點軌道根數(shù)不如傳統(tǒng)軌道根數(shù)和直角坐標(biāo)系那樣直觀，計算結(jié)果需要將改進的春分點軌道根數(shù)轉(zhuǎn)換為傳統(tǒng)軌道六要素，在三維軌跡作圖時需要將改進春分點坐標(biāo)系到地心赤道慣性系的轉(zhuǎn)換矩陣，這兩個計算方法可以參見文獻[14]。

依據(jù)最優(yōu)控制理論[15]，對于時間最優(yōu)指標(biāo)的最優(yōu)控制問題，可以構(gòu)造系統(tǒng)的哈密爾頓函數(shù)H如下：

(4)

式中：λ=[λp,λf,λg,λh,λk,λL]為軌道要素的協(xié)態(tài)變量；λm為質(zhì)量m的協(xié)態(tài)變量。根據(jù)?H/?α=0，可以得到變軌全過程中推力最優(yōu)控制方向為:

(5)

根據(jù)龐特里亞金極值原理，可以得出協(xié)態(tài)變量[λ,λm]的微分方程如下[15]:

(6)

(7)

式(6)中的?D/?x為6×1矩陣，?M/?x為6×3×6矩陣，M可依次取Mij，i=1,2,…,6，j=1,2,3，x依次取6個軌道要素。由于元素較多，這里不給出其具體表達式，詳見文獻[15]。

對于靜止軌道全電推進衛(wèi)星，其最終入軌的軌道根數(shù)應(yīng)滿足靜止軌道的約束條件，靜止軌道半長軸aGEO=42 164 km，則終端約束為:

(8)

式中：tf為轉(zhuǎn)移時間。最優(yōu)控制問題指標(biāo)主要有時間最優(yōu)和燃料最優(yōu)兩種形式。對于全電推進衛(wèi)星，由于其燃料利用率高，但轉(zhuǎn)移時間非常長，因此宜采用時間最優(yōu)作為優(yōu)化目標(biāo)，即使得J=-tf最大。

1.2 推力同倫求解過程

通過電推進完成靜止軌道的入軌過程通常需要數(shù)月至半年，由于圈數(shù)的增加導(dǎo)致軌道根數(shù)呈現(xiàn)振蕩變化，增加了間接法協(xié)態(tài)變量初值的敏感性，降低了猜測值收斂區(qū)間[17]。而對微分方程進行長時間、高精度的積分也極大的增加了求解的計算量和計算時間，增加了協(xié)態(tài)變量初值猜測的難度。為解決這一問題，本文借鑒了小推力軌道燃料最優(yōu)Bang-Bang問題求解中的一種有效策略，即同倫方法[11-12]，此方法主要用在燃料最優(yōu)和能量最優(yōu)兩種指標(biāo)函數(shù)之間的漸變，本文利用該思想將其用在時間最優(yōu)問題的推力大小漸變過程中。

同倫的實質(zhì)是指從一個簡單問題出發(fā)，逐漸改變系統(tǒng)的參數(shù)，使簡單的問題過渡到待求解的一個困難問題。對于最短時間變軌問題，大推力軌道轉(zhuǎn)移的變軌時間更短，協(xié)態(tài)變量收斂域?qū)挘曳e分時間短，求解困難程度將顯著降低。對于兩個不同推力的最優(yōu)變軌問題，只要二者對應(yīng)的推力接近，則可認為二者的協(xié)態(tài)變量初值也相近，因此稍大推力所求解得到的協(xié)態(tài)變量精確解可作為初值代入到稍小推力的問題中，這樣可以避免直接求解微小推力時的初值猜測的盲目性。

得到β0后，可以以此為初值，逐漸減少推力，完成小推力優(yōu)化的求解。本文采用推力大小指數(shù)縮減的方式完成推力同倫過程，定義大于1的縮減系數(shù)s，第j+1次求解時對應(yīng)推力與第j次的關(guān)系為Tj+1=Tj/s。由此構(gòu)成一個大小遞減的系列T0>T1>T2>…>TN，其中TN是要求解的小推力問題對應(yīng)的推力值，在指數(shù)遞減過程中如果推力小于TN，則將推力賦值為TN。將β0最為初始猜測值，代入T1對應(yīng)的優(yōu)化問題，通過SQP算法迭代后得出T1對應(yīng)的精確解β1，重復(fù)以上步驟，即可獲得系列β0，β1，βj，…,βN。其中βN即為最終要求解的小推力問題的協(xié)態(tài)變量初值。以上過程只要縮減系數(shù)s適當(dāng)，就可以保證每次求解是收斂的。

圖1是推力同倫方法求解小推力轉(zhuǎn)移問題的詳細流程。

圖1 推力同倫求解步驟Fig.1 Flow chart of the thrust homotopic

2 仿真算例

考慮一個從地球同步轉(zhuǎn)移軌道(GTO)到地球靜止軌道(GEO)的全電推進轉(zhuǎn)移過程，通過一箭雙星方式發(fā)射兩顆質(zhì)量各為1 800 kg全電推進通信衛(wèi)星。每顆衛(wèi)星依靠兩臺電推進聯(lián)合推進，每臺的推力230.58 mN，比沖為1 800 s。設(shè)星箭分離時衛(wèi)星軌道根數(shù)見表1，目標(biāo)軌道為GEO軌道，目標(biāo)軌道的ω以及v均不作約束。

表1 衛(wèi)星初始和終點軌道要素Table 1 Initial and target orbit elements

使用遺傳算法求解初始值，以及SQP求解非線性規(guī)劃問題需要給出[λ(t0),λm(t0),tf]的上下限，由于協(xié)態(tài)變量是無量綱單位，主要起作用的是相對大小，因此可以根據(jù)經(jīng)驗選取協(xié)態(tài)變量λ(t0),λm(t0)的上限為100，下限為-100。軌道轉(zhuǎn)移時間tf可根據(jù)推力大小和衛(wèi)星質(zhì)量，結(jié)合脈沖轉(zhuǎn)移速度增量ΔV來進行估算，并適當(dāng)擴大得出初始的上下限。

SQP算法采用Matlab中fmincon函數(shù)，大推力初始猜測值求解過程中帶約束遺傳算法采用Matlab中GA工具箱完成。推力初始值為50 N，縮減系數(shù)s=1.25。經(jīng)過21次推力縮減過程，求解出了總推力為461.17 mN的最優(yōu)軌跡，最優(yōu)飛行時間為109.81 d，消耗工質(zhì)為247.87 kg。尋優(yōu)得到的協(xié)態(tài)變量初值為λ(t0)=[-57.86,8.56,30.12,10.12,40.97,0.07]，λm(t0)=-1.56，最終的a=42 160.82 km，e=0.001 2，i=0.000 33°。

圖2是衛(wèi)星半長軸隨時間變化關(guān)系，圖3是偏心率和軌道傾角隨時間變化關(guān)系。在增大半長軸的同時，偏心率和軌道傾角隨之減少，最終均過渡到目標(biāo)軌道。

圖4是衛(wèi)星變軌過程三維軌跡圖，衛(wèi)星一共運行約155圈，坐標(biāo)軸單位是地球半徑Re=6 378 km。

圖5是推力逐次縮減過程中的λ(t0),λm(t0)最優(yōu)解。從圖5中可以看出，最優(yōu)解隨著推力下降逐漸穩(wěn)定，尤其是推力縮減5次之后(推力約16N)，λ(t0),λm(t0)波動幅度較小。由此可見，使用前一次求解得到的λ(t0),λm(t0)精確解作為推力縮減一次之后的尋優(yōu)初始值是合適的。

圖2 半長軸變化過程Fig.2 Semimajor-axis vs time

圖3 偏心率與傾角變化過程Fig.3 Eccentricity and inclination vs time

圖4 衛(wèi)星變軌軌跡Fig.4 Transfer orbit of the satellite

圖5 協(xié)態(tài)變量初值與推力縮減的關(guān)系Fig.5 Initial costate vs thrust reduction times

圖6是最優(yōu)轉(zhuǎn)移時間tf與推力累計縮減倍數(shù)之間的關(guān)系，轉(zhuǎn)移時間基本與推力縮減倍數(shù)呈線性關(guān)系。這與常識相符，在推力同倫逐次求解過程中，也可以利用該規(guī)律來預(yù)計推力縮減一次后的最優(yōu)轉(zhuǎn)移時間，將s倍tf作為下一次SQP算法中轉(zhuǎn)移時間估計值。

圖6 最優(yōu)轉(zhuǎn)移時間與推力縮減量的關(guān)系Fig.6 Minmum transfer time vs thrust reduction

3 結(jié)束語

本文采用推力同倫的思想，從比較容易的大推力問題出發(fā)，通過逐步減少推力的方式，采用間接法求解了全電推進衛(wèi)星最短時間入軌優(yōu)化問題。采用推力同倫的策略大大降低了小推力問題間接法求解過程的初值敏感性問題，可以使問題以較高的概率收斂到最優(yōu)解。以往多圈轉(zhuǎn)移所采用的軌道平均技術(shù)或反饋控制率僅可以獲得近似解或平均解，而該方法可提供另一種途徑，使得采用間接法求解上百圈小推力轉(zhuǎn)移問題的理論最優(yōu)解成為可能。

本文的動力學(xué)方程是基于二體問題的，后續(xù)的進一步研究可以考慮將該方法應(yīng)用到考慮攝動，地影能源約束等非理性情況下的全電推進衛(wèi)星入軌優(yōu)化。