張 冉,李小娟,韓 潮,張亞航,于俊慧
(1.北京空間飛行器總體設(shè)計部·北京·100094; 2.北京航空航天大學(xué) 宇航學(xué)院·北京·100191)
電推進(jìn)系統(tǒng)以其比沖大、推力小的技術(shù)特點在深空探測、近地軌道阻力補(bǔ)償、軌道位置保持等領(lǐng)域得到應(yīng)用[1]。近年來,電推進(jìn)器在GEO軌道轉(zhuǎn)移中嶄露頭角。歐洲先進(jìn)通信衛(wèi)星ARTEMIS(2001)和美國先進(jìn)極高頻軍事通信衛(wèi)星AEHF-1(2010)原計劃以化學(xué)推進(jìn)和電推進(jìn)結(jié)合的混合方式進(jìn)行GTO-GEO軌道的轉(zhuǎn)移,但因不同的故障導(dǎo)致化學(xué)推進(jìn)失效,電推進(jìn)系統(tǒng)成為GEO軌道轉(zhuǎn)移的主推力器,超時工作挽救了任務(wù)。波音公司基于702SP平臺發(fā)展了全電推進(jìn)通信衛(wèi)星,迅速投入商業(yè)市場;截至目前,已有五顆702SP平臺全電推進(jìn)衛(wèi)星發(fā)射入軌,單星質(zhì)量在1950 kg~2300 kg之間,每顆衛(wèi)星配備了4臺XIPS-25離子推力器,總推力大小為165 mN,GEO入軌時間為6個月;其中兩次發(fā)射采用獵鷹九號“一箭雙星”的發(fā)射方式,極大降低了發(fā)射成本??罩锌蛙嚪绖?wù)與航天基于E3000EOR平臺也發(fā)展了全電推通信衛(wèi)星,單星質(zhì)量3550 kg~5300 kg之間,每顆衛(wèi)星配備一臺SPT140D穩(wěn)態(tài)等離子推力器,推力大小約290 mN,GTO-GEO軌道轉(zhuǎn)移時間為4~6個月;目前已有三顆E3000EOR平臺的通信衛(wèi)星發(fā)射入軌。全電推進(jìn)通信衛(wèi)星采用比沖更高的電推進(jìn)系統(tǒng)進(jìn)行GTO-GEO的軌道轉(zhuǎn)移,可以大幅減少燃料質(zhì)量,有效降低衛(wèi)星總重或提高有效載荷的質(zhì)量比,延長任務(wù)的壽命,降低發(fā)射成本[2-4]。因此,在GTO-GEO軌道轉(zhuǎn)移中采用電推進(jìn)技術(shù)已經(jīng)成為中小型通信衛(wèi)星發(fā)展的趨勢,研究電推進(jìn)GEO軌道轉(zhuǎn)移優(yōu)化和制導(dǎo)方法具有重要的工程意義。
全電推進(jìn)GEO衛(wèi)星的軌道轉(zhuǎn)移設(shè)計主要分為軌跡優(yōu)化和制導(dǎo)控制兩個方面。求解航天器軌跡優(yōu)化最優(yōu)控制問題的數(shù)值方法分為直接法、間接法和結(jié)合二者的混合法[5-6]。直接法通過離散變量和參數(shù)將最優(yōu)控制問題轉(zhuǎn)換為參數(shù)優(yōu)化問題,然后通過非線性規(guī)劃算法獲得最優(yōu)解,直接法的計算量比較大。間接法通過變分法和龐特里亞金極大值原理構(gòu)造兩點邊值問題(Two-Point-Boundary-Value-Problem,TPBVP),然后通過打靶法求解靶函數(shù);間接法的困難在于靶函數(shù)的收斂域很小,靶函數(shù)對于初值特別敏感。一些學(xué)者提出同倫法[7-9]和非線性濾波參數(shù)估計算法[10]等解決初值敏感的數(shù)值問題?;旌戏▌t結(jié)合了直接法和間接法各自的優(yōu)點,通常假定推力結(jié)構(gòu),由間接法的必要條件保證控制的最優(yōu)性。在對控制結(jié)構(gòu)適當(dāng)假設(shè)和簡化的前提下,一些學(xué)者獲得了軌跡優(yōu)化的解析求解方法。Kechichian[11]提出了保持偏心率為零的前提下使半長軸和軌道傾角最快改變的分段常值控制律;田百義等[12]在小推力GEO入軌問題中假設(shè)平面內(nèi)外控制解耦,平面內(nèi)保持遠(yuǎn)地點高度不變,平面外在升降交點采取控制的解析策略。高揚(yáng)[13]提出三種分別控制軌道半長軸、偏心率和軌道傾角的切向、慣性和偏航控制律。解析法的控制律只能達(dá)到次優(yōu)的控制結(jié)果,與最優(yōu)轉(zhuǎn)移相比有一定差別。在小推力軌道轉(zhuǎn)移制導(dǎo)方法研究方面,許多文獻(xiàn)提出了跟蹤最優(yōu)軌跡的控制思路[14-15],通過設(shè)計一組控制參數(shù),修正實際軌道與參考軌道之間的偏差;另外,基于李雅普諾夫函數(shù)理論的制導(dǎo)律也得到利用,主要問題是李雅普諾夫函數(shù)的定義方式和增益的計算[16-17]。
本文主要研究全電推進(jìn)GEO衛(wèi)星的入軌轉(zhuǎn)移制導(dǎo)策略問題,提出了一種跟蹤最優(yōu)參考軌跡的分段常值制導(dǎo)策略,該策略形式簡單,易于星上執(zhí)行,并且控制結(jié)果與最優(yōu)控制律相差不大,可為我國全電推進(jìn)通信衛(wèi)星平臺的研究提供參考。
本文采用經(jīng)典軌道要素的高斯攝動方程作為軌道轉(zhuǎn)移優(yōu)化的控制模型,形式如下:
(1)
式中,μ為地球引力常數(shù),a為軌道半長軸,e為軌道偏心率,i為軌道傾角,Ω為升交點赤經(jīng),ω為近地點幅角,θ為真近點角;p為軌道半通徑,h為軌道角動量大小,r為衛(wèi)星地心距,v為衛(wèi)星速度大?。豢刂萍铀俣萢在衛(wèi)星速度方向、垂直速度方向和動量矩方向的三個分量分別為ft,fn,fh。a是除二體加速度外的其他外力加速度矢量,包含發(fā)動機(jī)推力和攝動加速度矢量兩部分;本文攝動加速度僅考慮了地球非球形引力場J2項的影響。
設(shè)衛(wèi)星質(zhì)量為m,全電推進(jìn)發(fā)動機(jī)具有上限推力大小Tmax和常值比沖Isp,g0為赤道重力加速度,式(1)改寫為矩陣的形式:
(2)
對于全電推進(jìn)GEO衛(wèi)星的燃料最優(yōu)軌道轉(zhuǎn)移問題,性能指標(biāo)為
(3)
轉(zhuǎn)移時間tf是相應(yīng)時間最短問題轉(zhuǎn)移時間tfmin的倍數(shù),用系數(shù)Ctf表示,Ctf>1。
利用變分法和極大值原理,小推力燃料最優(yōu)GEO軌道轉(zhuǎn)移問題可以表述為下面的兩點邊值問題:
(4)
其中,λ是與軌道要素對應(yīng)的協(xié)態(tài)變量,λm是與質(zhì)量對應(yīng)的協(xié)態(tài)變量,λθ,λΩ,λω分別是與真近點角、升交點赤經(jīng)、近地點幅角對應(yīng)的協(xié)態(tài)變量,Ctf是飛行時間系數(shù);下標(biāo)“0”表示初始時刻的狀態(tài),下標(biāo)“f”表示變軌終止時刻的狀態(tài)。
上述TPBVP問題需要利用打靶法求解協(xié)態(tài)變量的初值λ(t0)、λm(t0)。對于電推進(jìn)GEO軌道的轉(zhuǎn)移問題,軌道轉(zhuǎn)移過程包含幾百圈軌道,轉(zhuǎn)移時間長達(dá)半年左右,靶函數(shù)對于初值非常敏感。為了解決數(shù)值求解困難,本文采用了文獻(xiàn)[10]提出的UKF非線性濾波參數(shù)估計算法。
本文仿真中采取的初始軌道是GTO轉(zhuǎn)移軌道,軌道要素如表1所示。飛行時間系數(shù)Ctf設(shè)為1.5,使用UKF濾波參數(shù)估計算法求解不同推力值Tmax情況下的燃料最優(yōu)GEO軌道轉(zhuǎn)移問題。為避免結(jié)果陷入局部極值,本文采取逐步縮小推力幅值Tmax的延拓思路,每一步計算時以上一步計算結(jié)果為初值,并根據(jù)控制結(jié)構(gòu)篩查全局最優(yōu)的極值。
表1 軌道轉(zhuǎn)移初始與目標(biāo)軌道要素
圖1(a)顯示了燃料最優(yōu)GEO軌道轉(zhuǎn)移問題不同推力幅值Tmax條件下軌道半長軸、偏心率和軌道傾角隨時間的變化歷程,為便于比較,橫軸的時間做了去單位化處理(t=t/tf)。Tmax較大時,軌道要素的變化呈現(xiàn)出“階梯”形,這是由Bang-bang控制結(jié)構(gòu)造成的。隨著Tmax的減小,軌道要素的變化變得平滑,局部放大后仍存在“階梯”特性。不同Tmax條件下,軌道要素呈現(xiàn)出相似的變化趨勢,可以考慮用一條平滑曲線去擬合軌道要素的變化規(guī)律。
第1節(jié)獲得了燃料最優(yōu)GEO轉(zhuǎn)移的最優(yōu)軌跡和最優(yōu)控制曲線,這些軌道和控制信息可以作為實際工程中飛行制導(dǎo)的參考,然而直接存儲上述軌道與控制信息有以下幾個缺點:①軌道圈數(shù)多,轉(zhuǎn)移時間長,存儲數(shù)據(jù)量大;②最優(yōu)控制律的表達(dá)式與協(xié)態(tài)變量相關(guān),而協(xié)態(tài)變量沒有其他物理意義;③狀態(tài)變量和控制變量與時間關(guān)聯(lián),攝動力和控制偏差會造成相應(yīng)狀態(tài)在指定時間達(dá)到預(yù)定數(shù)值,使存儲信息失效;④軌道要素隨時間的小幅度振蕩影響制導(dǎo)律的設(shè)計。
為解決以上問題,本文利用擬合方式獲得簡單解析形式的曲線,更適合在星上計算機(jī)存儲和利用。根據(jù)最優(yōu)軌跡的形狀,本文選擇簡單的多次多項式對最優(yōu)軌跡進(jìn)行擬合,數(shù)值試驗表明,采用四次及以上的多項式對最優(yōu)軌跡的擬合效果比較好。
圖1(b)為使用四次多項式對推力幅值Tmax為2.5 N、1.0 N、0.5 N和0.2 N條件下的燃料最優(yōu)GEO軌道轉(zhuǎn)移軌跡進(jìn)行擬合時獲得的參考軌道,四條擬合曲線非常接近。Ctf固定時,燃料最優(yōu)軌道轉(zhuǎn)移軌跡的軌道要素變化趨勢與推力大小關(guān)系不大,在設(shè)計小推力燃料最優(yōu)轉(zhuǎn)移問題的參考軌道時,僅需計算較大推力條件下的軌道轉(zhuǎn)移工況,通過曲線擬合方式可以快速得到小推力條件下的參考軌道。
(a)最優(yōu)軌跡
(b)擬合軌跡
第一節(jié)的最優(yōu)軌跡對應(yīng)的最優(yōu)控制是一組開環(huán)的控制,這些數(shù)據(jù)中控制力的大小和方向等數(shù)據(jù)量非常大,不適合在星上計算機(jī)存儲;控制力的方向與當(dāng)前協(xié)態(tài)變量和狀態(tài)變量以及時間相關(guān),協(xié)態(tài)變量本身沒有具體的物理含義,這給制導(dǎo)律的設(shè)計帶來了很多不便。在衛(wèi)星的實際變軌過程中,不確定的擾動因素、未建模的攝動力、導(dǎo)航與控制偏差會導(dǎo)致衛(wèi)星偏離參考軌道,必須設(shè)計含閉環(huán)反饋的制導(dǎo)控制律。在獲得解析化的參考軌道之后,多圈長時的軌道轉(zhuǎn)移問題可以分解為多個單圈轉(zhuǎn)移設(shè)計問題,根據(jù)參考軌道的形式,將整個軌道改變量以一種近乎優(yōu)化的方式分配到每一個軌道周期內(nèi),本節(jié)基于此基礎(chǔ)設(shè)計單圈內(nèi)的軌道制導(dǎo)律。
電推進(jìn)控制加速度[ft,fn,fh]T與引力加速度相比是小量,使用偏近點角E代替真近點角θ作為第六個軌道要素,dE/dt可以近似為:
(5)
式中,n為軌道角速度。將式(5)代入式(1),可得到以偏近點角為獨立變量的微分方程:
(6)
在推力加速度[ft,fn,fh]T作用下,衛(wèi)星軌道偏近點角由E1到E2的過程中軌道要素的變化量具有解析表達(dá)式,軌道要素X的改變量具有以下形式:
(7)
根據(jù)式(7),根據(jù)不定積分的公式,可以獲得TX,NX,HX的表達(dá)式,具體形式可以參考文獻(xiàn)[18]。
衛(wèi)星在偏近點角由E1到E2的過程中滿推力工作時,質(zhì)量的改變量為:
(8)
圖2顯示了燃料最優(yōu)GEO軌道轉(zhuǎn)移的最優(yōu)控制的曲線和設(shè)計的常值推力控制結(jié)構(gòu),其中ut,un,uh是推力加速度單位矢量在速度方向、垂直速度方向和動量矩方向上的三個分量。細(xì)線表示的是最優(yōu)控制曲線,黑粗實線是根據(jù)最優(yōu)控制曲線形狀設(shè)計的分段常值控制結(jié)構(gòu)。
(a)軌道控制前中期
(b)軌道控制末期圖2 最優(yōu)和分段常值推力控制結(jié)構(gòu)Fig.2 Optimal and piecewise constant thrust control structure
在軌道轉(zhuǎn)移的前中期,平面外有兩個推力弧段:[E1,E2]和[E2,E3]段,平面外推力方向分別用β1和β2描述;平面內(nèi)的推力方向角用α1和α2描述。常值推力結(jié)構(gòu)見表2所示,在軌道轉(zhuǎn)移的前中期,每個軌道周期內(nèi),在三段常值推力的作用下,軌道要素的變化量為:
(9)
軌道轉(zhuǎn)移的末期,平面外有四個推力弧段,弧段用五個設(shè)計變量描述:E1,E2,E3,E4,E5;平面內(nèi)根據(jù)E3和π的大小關(guān)系,分為五個推力弧段;為減少變量數(shù)目,分別用四個方向角α1,α2,β1,β2描述推力矢量的方向,軌道轉(zhuǎn)移末期的常值推力結(jié)構(gòu)如表3所示。在軌道轉(zhuǎn)移末期,在五段常值推力作用下,軌道要素的變化量為:
ΔX=
(10)
表3 軌道末期的常值推力弧段
在考慮J2攝動時,單圈內(nèi)的控制結(jié)構(gòu)如圖3所示。在軌道轉(zhuǎn)移的前期和中期,控制序列包含三個連續(xù)的推力弧段和一個漂移弧段,設(shè)計變量為
σ=[E1,E2,E3,α1,α2,β1,β2]
(11)
性能指標(biāo)為:
J=(E3-esinE3)-(E1-esinE1)
(12)
約束條件為:
0≤α1,α2,β1≤π/2,-π/2≤β2≤π/2
0 (13) (a)軌道控制前中期 (b)軌道控制末期圖3 單圈內(nèi)的控制序列示意圖Fig.3 Schematic diagram of control sequence in a single loop 在軌道轉(zhuǎn)移的末期,單圈內(nèi)控制序列包括五個推力弧段和兩個漂移弧段,設(shè)計變量為: σ=[E1,E2,E3,E4,E5,α1,α2,β1,β2] (14) 性能指標(biāo)為: J=(E4-esinE4)-(E2-esinE2)+(E1-esinE1)-(E5-esinE5) (15) 約束條件為: (16) 其中ΔX為根據(jù)式(10)所計算的軌道要素改變量。 以上完成了控制單圈內(nèi)數(shù)學(xué)規(guī)劃問題的數(shù)學(xué)建模,這是一個具有簡單邊界約束的優(yōu)化問題,該類問題有較成熟的優(yōu)化算法;控制律充分利用了最優(yōu)控制的變化規(guī)律,給出了較好的初值;每一圈的優(yōu)化結(jié)果可以作為下一圈的設(shè)計變量初值,提高了優(yōu)化求解速率。 選取表1所示的初末軌道要素,最大推力值為Tmax=0.2N,在最優(yōu)軌跡計算時考慮J2項攝動,最優(yōu)軌跡的飛行時間為265.53天,燃料消耗123.84 kg;采用四次多項式擬合生成參考軌道。圖4顯示了本文提出的分段常值推力跟蹤制導(dǎo)策略的控制結(jié)果,半長軸、偏心率和軌道傾角的跟蹤效果非常好。在考慮J2項攝動的情況下,變軌結(jié)束后軌道半長軸控制誤差為0.83 km,偏心率誤差為0.00145,軌道傾角控制誤差為0.0025°。由近地點幅角的變化規(guī)律可以看出,跟蹤控制并不會抑制近地點幅角ω的漂移。分段常值跟蹤控制的飛行時間為266.96天,燃料消耗132.72 kg。本文提出的分段常值跟蹤制導(dǎo)方法與最優(yōu)軌跡相比,飛行時間延長0.54%,燃料消耗增加了7.17%。燃料消耗的主要原因是簡化的分段常值的推力結(jié)構(gòu)控制效率低于連續(xù)變化的最優(yōu)控制結(jié)構(gòu)。如果能根據(jù)最優(yōu)控制結(jié)構(gòu)的形狀,找出一組簡化的數(shù)學(xué)解析式來模擬推力值與偏近點角的關(guān)系,燃料消耗能夠進(jìn)一步優(yōu)化。 圖4 GTO-GEO轉(zhuǎn)移分段常值制導(dǎo)軌跡、擬合軌跡和最優(yōu)軌跡Fig.4 GTO-GEO transfer piecewise constant guidance trajectory, fitting trajectory and optimal trajectory 本文針對全電推進(jìn)GEO衛(wèi)星的入軌轉(zhuǎn)移制導(dǎo)策略進(jìn)行了研究分析,提出了一種跟蹤參考軌道的制導(dǎo)策略,并結(jié)合算例分析了制導(dǎo)策略的效果。根據(jù)最優(yōu)軌跡的變化規(guī)律,提出采用多次多項式擬合最優(yōu)軌跡生成參考軌道的方法,參考軌道形式簡單,能夠表征最優(yōu)轉(zhuǎn)移軌道的變化特點,易于在星上計算機(jī)存儲和使用。跟蹤參考軌道的閉環(huán)反饋控制策略將多圈軌道轉(zhuǎn)移的控制量分配到每一個軌道單圈內(nèi),保證了制導(dǎo)的燃料近優(yōu)性。在軌道單圈內(nèi),本文提出一種分段常值推力控制結(jié)構(gòu),這種控制結(jié)構(gòu)形式簡單,在此假設(shè)下軌道要素改變量有解析解,簡化了制導(dǎo)策略設(shè)計;并且常值推力結(jié)構(gòu)易于工程實施。仿真結(jié)果顯示,在考慮J2攝動條件下,分段常值的制導(dǎo)策略比最優(yōu)控制解多消耗7.17%的燃料,對參考軌道的跟蹤效果好,控制精度高。上述設(shè)計可為我國全電推進(jìn)GEO通信衛(wèi)星的入軌控制提供參考。為進(jìn)一步提高跟蹤制導(dǎo)策略的效率,后續(xù)可研究一類模擬最優(yōu)控制結(jié)構(gòu)的解析連續(xù)變化的曲線,并進(jìn)一步考慮更多的工程約束條件。4 結(jié) 論