晁 寧,李言俊
(西北工業(yè)大學航天學院,西安710072)
眾所周知,平動點附近的暈軌道具有指數(shù)不穩(wěn)定性,因此對運行在其上的探測器進行軌道保持是十分必要的。而自從開始了各種平動點空間任務,暈軌道保持控制就一直是研究的熱點。Farquhar[1]和Breakwell[2]最早提出了target控制模式。隨后Giamberardino基于線性模型設計了漸近穩(wěn)定的非線性控制器,實現(xiàn)了Halo軌道的漸近跟蹤和擾動補償[3]。Howell等[4]同樣基于線性化模型的思想,給出了控制器精度與能耗均滿足要求的折中策略。David Cielaszyk和BongWie[5]提出了用LQR線性二次型最優(yōu)控制方法來保持暈軌道的穩(wěn)定。Rahmani[6]利用最優(yōu)控制理論中極值曲線的變分成功求解了兩點邊值問題,實現(xiàn)了Halo軌道維持??紤]到最優(yōu)控制求解的難度,Ming Xin等[7]提出了用于近似求解HJB方程的θ-D方法,這種方法具有很高的實時性。胡少春等[8]將序優(yōu)化理論與微分修正法相結合,優(yōu)化了暈軌道的入軌機動問題。這種方法在優(yōu)化過程中具有收斂速度快、對初值不敏感與計算量小的優(yōu)點。
基于最優(yōu)控制思路的方法在Halo軌道維持中具有較高控制精度,但在每個控制點求解系統(tǒng)微分方程組時需要進行大量計算且具有較多控制點,系統(tǒng)實時性能受到影響。而θ-D法雖具有很高實時性但屬于次最優(yōu)方法。本文在三體模型誤差線性模型的基礎上,利用有限時間調(diào)節(jié)器問題推導了暈軌道周期內(nèi)的連續(xù)小推力控制方案。針對整周期控制方式在超調(diào)后狀態(tài)量收斂速度慢的問題,通過分段連續(xù)推力模式來近似瞬時脈沖推力控制模式,給出了最短分段控制時間的計算方法。方法基于小偏差假設對非線性模型進行了線性近似,以最優(yōu)控制方案推導線性反饋控制律。同時在整個Halo軌道周期僅需1~2個控制點,減小了計算量且使得系統(tǒng)結構簡化。仿真實驗表明,控制方法能夠根據(jù)實際軌道與標稱軌道偏差的大小,調(diào)整控制時間區(qū)間長度,以盡可能低的能耗快速消除探測器入軌狀態(tài)偏差。
天體力學中多數(shù)情況都可描述為一個質(zhì)量忽略不計的小天體P在兩個相互環(huán)繞運動的大天體P1、P2(P1>P2)引力作用下的運動狀態(tài),這是比二體模型更精確的一種合理近似。對于圓限制性三體問題,通常利用以兩大天體質(zhì)心為圓心旋轉的會合坐標系或稱旋轉坐標系來研究。假定P不影響P1、P2的運動,兩大天體共同繞其質(zhì)心做角速度為ω的圓周運動。P1指向P2的方向為會合坐標系的x軸,ω方向為z軸,y軸與x、z軸成右手系。這種假設下的動力學模型不考慮攝動因素,即無攝圓限制性三體模型,小天體的無量綱化運動對應一個二階常微初值問題
分量形式為
其中,Ω(x,y,z)=(x2+y2)/2+U(r1,r2),引力勢函數(shù)U(r1,r2)=(1- μ)/r1+μ/r2,μ=m2/(m1+m2)是航天器到大天體的距離是航天器到小天體的距離。
對于該CRTBP目前僅找到5個特解和一個Jaccobi積分[9]。前者對應5個引力平衡點,其中三個共線平衡點不穩(wěn)定,位于兩大天體連線上,記作Li(i=1,2,3);兩個三角平衡點Lyapunov穩(wěn)定。由于前兩個點的應用價值較大,引發(fā)了學者們的深入研究。五個平衡點處對應的Jacobi積分常數(shù)Ci(μ)有如下關系:
曲面2Ω(x,y,z)=C即為零速度面。
對于三體模型來說,地月系統(tǒng)數(shù)據(jù)量級一般都比較大。為了簡化計算,需要利用無量綱化方法對計量單位進行省略并將數(shù)據(jù)縮小相應倍數(shù)。地月系數(shù)據(jù)計算中無量綱化時對應的長度時間和質(zhì)量單位分別為:
[M]=m1+m2——兩天體質(zhì)量
[L]=3.84401×105km,[M]=6.0477×1024kg,[T]=3.7519×105s
因此,月球的無量綱質(zhì)量為μ=0.012153,地球的無量綱質(zhì)量為1-μ。若按照1年為365.25個平太陽日計算,則1年和1天分別是3.15576×107s和8.64×104s,無量綱化后的時間為84.1110和0.2303。無量綱速度V00和實際速度(km/s)的轉化關系為
式(4)表示了地月系中一個單位的無量綱速度與實際速度的關系。其中q、m分別表示無量綱系統(tǒng)中的距離和時間。同樣可以得到無量綱加速度a00和實際加速度的轉換關系為
3.1.1不穩(wěn)定性
暈軌道是存在于共線平動點附近的一類周期軌道。這種軌道具有指數(shù)不穩(wěn)定性,發(fā)散快,并且對初始值十分敏感。即使不考慮攝動因素,并且入軌精度均達到10-8(10m)的量級,無控的暈軌道最多也只能維持3個周期左右(以THalo=13.5808天為例),之后就開始發(fā)散,如圖1所示。利用Richardson三階近似解獲得的軌道雖然能夠較準確地表現(xiàn)軌道的三維周期運動,但它是理想化的解模型,不能夠體現(xiàn)出暈軌道的弱穩(wěn)定性??紤]到入軌誤差和各種實際因素引發(fā)的模型變異,需要對在軌探測器按照偏離規(guī)律進行軌控。
圖1 地月系L1點某暈軌道Fig.1 Earth-Moon system L1 point some halo orbit
3.1.2周期性
原系統(tǒng)線性化模型中的系數(shù)矩陣A是時變的。但是在理想入軌狀況下,經(jīng)歷了n個暈軌道周期的時間后,A中的元素具有不變性,即
于是,線性化后的時變系統(tǒng)就可以化為單個周期內(nèi)的定常系統(tǒng)。系統(tǒng)矩陣元素依積分起始點而定。因此暈軌道上所有點的狀態(tài)都能夠作為積分起始點,只要積分周期為軌道周期的整數(shù)倍,時變系統(tǒng)就能夠化為定常系統(tǒng),從而簡化計算。
3.1.3可控性
其中
取定系統(tǒng)矩陣A(t)為多元向量函數(shù)f(X)的雅克比矩陣,控制矩陣b為增廣單位陣[03I3]T。從式(9)可以看出,系統(tǒng)僅需要對后3階進行控制,因此分析可控性前首先將系統(tǒng)降為3階。控制矩陣b取為I3,則從上式中容易看到可控性矩陣滿秩,線性時變系統(tǒng)可控。
系統(tǒng)矩陣A(t)是以狀態(tài)向量X(t)為變量的函數(shù)。以暈軌道上任意點為研究對象時,時變系統(tǒng)又能夠轉換為定常系統(tǒng)。線性定常系統(tǒng)完全能控的充要條件為
其中n為系統(tǒng)階數(shù)。線性系統(tǒng)可控性特征為:
(1)系統(tǒng)在平動點處可控;
(2)一個周期內(nèi),暈軌道上所有點對應的可控性矩陣序列K(t)均滿秩,t∈[t0,t0+THalo]。
將線性時變系統(tǒng)式在各軌道機動點處做定常變換后,容易驗證對應的(8)式是成立的。因此誤差線性系統(tǒng)在小偏差范圍內(nèi)是完全可控的,其上所有點都能夠作為系統(tǒng)控制作用u的施加點。
用于深空探測的無攝圓限制性三體模型表現(xiàn)出很強的非線性,在考慮攝動因素后這種非線性特性就更加明顯。而當前比較成熟的控制律多數(shù)以線性系統(tǒng)為研究對象,并且具有良好的控制效果與魯棒性。當實際軌道與標準軌道的狀態(tài)誤差不大時,誤差線性系統(tǒng)模型具有較高的準確度。通常用到的線性化方法有兩種:
(1)將誤差狀態(tài)作為新的狀態(tài)變量,利用多元向量函數(shù)的雅克比矩陣作為線性系統(tǒng)的系數(shù)矩陣;
(2)仍以原位置速度作為狀態(tài)量,將動力學系統(tǒng)轉換關系作為系數(shù)矩陣。
本文以前者作為研究重點進行分析,不對后者詳細敘述。
對上式微分并在平衡位置進行一階泰勒展開,可得
忽略高階項,并取向量函數(shù)的雅克比矩陣為系數(shù)矩陣,即
雅克比矩陣為
計算可知,處于暈軌道狀態(tài)時,交叉項均為0。于是,非線性系統(tǒng)可近似化為線性系統(tǒng)
這種方法在小擾動情況下具有較高準確性。在此假設下,可以通過選擇能量函數(shù)通過使其非正定獲得控制律,或者選擇線性調(diào)節(jié)器的指標函數(shù)利用極小值原理或動態(tài)規(guī)劃法求解到控制律u(t)。同時可以近似認為擾動周期與軌道周期相同,即TR=T。
上面A(t)為系統(tǒng)矩陣,其中變量元素的值會隨探測器在會合坐標系中的位置而變化。對誤差線性系統(tǒng)來說,控制器設計的目標即推導出合適的控制向量u,使得誤差系統(tǒng)狀態(tài)量~X<ε,ε為正常數(shù)。
軌道保持過程中要求以最小的速度增量實現(xiàn)最高的位置精度,這樣的控制要求類似于最優(yōu)控制中對性能指標的表述。因此考慮利用極小值原理及Riccati方程來推導一個THalo內(nèi)的小推力連續(xù)控制方案。
以一個暈軌道周期T為控制區(qū)間,取有限時間線性調(diào)節(jié)器性能指標為
終端固定,其狀態(tài)約束為X(tf)=0。因此,單周期暈軌道控制問題屬于終端固定的有限時間狀態(tài)調(diào)節(jié)器。式(14)中Q>0、R>0分別為加權矩陣,用來控制各分量的比重。通常Q為對角陣,各個元素取狀態(tài)向量對應量綱數(shù)量平方的倒數(shù)。
Hamilton函數(shù)取為
根據(jù)線性調(diào)節(jié)器問題相關結論[11],要求K(t)陣滿足對稱正定和逆Riccati方程,即
線性動力學誤差系統(tǒng)對應Riccati方程階數(shù)較高,通常利用數(shù)值積分進行求解來提高運算速度。積分中需要注意到:如果按照時間反向推演的方向進行積分,數(shù)值積分的初值就是Riccati方程解的終值,即K-1(tf)=0。因此設置K(0)=[0,…,0],K∈Rk(n),其中為系統(tǒng)矩陣階數(shù)。因此可以在以每個暈軌道周期為單位求解Riccati方程來構造反饋控制律。
于是有最優(yōu)控制
設總的控制量為
式中u0(t)為變軌需要的加速度。當飛行器通過一次推力由L1點穩(wěn)定流形進入暈軌道運行時,u0(t)=0,于是控制總量u(t)=u*(t)。
最優(yōu)軌跡為下面一階線性微分方程的解
小推力控制模式具備一定的優(yōu)勢,如推力穩(wěn)定,線性特性好。但是這種方式控制周期長,施控前后時間系統(tǒng)較難統(tǒng)一,在某些需要付出較高能耗來實現(xiàn)軌道快速跟蹤的場合就不合適使用。SCTC模式的主要方法就是縮短控制區(qū)間,以此來近似瞬時脈沖推力模式,使其兼具連續(xù)推力控制與脈沖控制的優(yōu)點。軌道修正時,每一個分段連續(xù)控制區(qū)間被模擬為瞬時推力模式的一個脈沖。通常每條暈軌道都會事先規(guī)劃2~3個機動點,這些機動點處就是每個分段控制器開始工作的位置;而控制器的工作時間就由控制區(qū)間tf-t0確定。
但是SCTC模式存在一個限制:如果入軌誤差比較大,并且要求航天器在較短的時間內(nèi)收斂至標準軌道,就需要提高機動加速度,而這樣又會造成振蕩加劇,這樣也會影響航天器的控制性能。因此,不同的狀態(tài)誤差就對應了不同的連續(xù)推力模式能夠接受的最短控制時間min(Tcon)。
為使得控制向量在規(guī)定的時間內(nèi)到達零點附近,結合小推力模式提供的機動能力,能夠獲得一個基本保持不變的閾值。具體來講,利用距離誤差除以控制時間,能夠獲得一個速度閾值;而利用速度誤差除以控制時間,能夠獲得一個加速度閾值。結合小推力模式下飛行器能夠提供的變軌能力和實驗數(shù)據(jù),在地月三體系統(tǒng)中,飛行器進行軌道修正能夠承受速度和加速度選擇如下經(jīng)驗值:
無量綱速度:V00=7.4006×10-4;
無量綱加速度:a00=0.001365。
當給定了入軌的初始狀態(tài)誤差,通過上面兩個閾值能夠確定兩個控制時間,即
為了保證修正軌道的收斂性,當t1≠t2時,控制區(qū)間Tcon=max(t1,t2)。
由于地月系單位長度為地月均距約為3.8×105千米,因此初始狀態(tài)包含了方差為10-6的隨機高斯噪聲擾動,意味著存在100m級的偏差。
圖2 誤差等級σ2=10-6時軌控前后情況Fig.2 Orbit control situation under the error grade ofσ2=10-6
可以看出,通過線性最優(yōu)控制,狀態(tài)曲線在經(jīng)過若干個減幅振蕩周期后,能夠在有限時間內(nèi)(THalo=12.33天)逐漸收斂至標準暈軌道。同時,入軌誤差等級大,軌道發(fā)散嚴重,受控軌道振幅也就大。另外,小推力控制的優(yōu)點是能量消耗低,所有能耗被分配整個軌道,使得每時刻的能量需求很小。以圖2的情況為例,整個暈軌道周期中三軸消耗的總速度增量共4.7782×10-6m/s。
以10-6入軌誤差等級的某個初值為例,比較了SCTC和TCTC模式下的控制效果,如圖3。從圖3中比較可以看出,分段連續(xù)軌控下實際軌道與標準軌道的貼近度比較高,并且通過1個平太陽日(一個平太陽日的無量綱數(shù)為0.2303)左右的時間實際軌道基本收斂至標準軌道。一個暈軌道周期實施1~2次狀態(tài)反饋控制即可。
圖3 兩種軌控效果對比Fig.3 Contrast of two kinds of orbit control effect
表1 不同eX對應的最佳控制時間及跟蹤方差和(一個平太陽日的無量綱數(shù)為0.2303)Table 1 Best control time and sum of tracking variance corresponding to different eX
圖4 一個T Halo內(nèi)的加速度曲線Fig.4 Acceleration curves in a T Halo
下面以10-7、10-6和10-5三個不同入軌誤差等級為例,對比分段連續(xù)推力模式的控制情況。表1是不同誤差eX對應的t1、t2及跟蹤結果,跟蹤方差和(TVS)表示三軸分量方差在單個THalo內(nèi)的總和。TVS按照式(21)計算
本文在小偏差假設基礎上將原限制性三體模型問題轉換為誤差線性模型進行控制方法研究,并利用最優(yōu)控制方法推導了暈軌道周期內(nèi)的連續(xù)小推力控制方案。結合減少能耗和提高狀態(tài)量收斂速度的綜合考慮,提出分段連續(xù)推力控制模式來改善超調(diào)后的不良品質(zhì),其線性反饋控制律也使得系統(tǒng)結構得到簡化。方法能夠根據(jù)實際軌道與標稱軌道偏差的大小,調(diào)整控制時間區(qū)間長度,以盡可能低的能耗快速消除探測器入軌狀態(tài)偏差。仿真結果驗證了方法的有效性。
[1]Farquhar R W.The control and use of libration point satellite[R].NASA TR R-346.
[2]Breakwell J V,Kamel A A,Ratner M J.Stationkeeping of a translunar communication station[J].Celestial Mechanics,1974,10(3):357-373.
[3]Giamberardino P,Monaco S.On halo orbits spacecraft stabilization[J].Acta Astronautica,1996,38(12):903-925.
[4]Howell K C,Pernicka J.Stationkeeping method for libration point trajectories[J].Journal of Guidance,Control and Dynamics,1993,16(1):151-159.
[5]Cielaszyk D,Wie B.New approach to halo orbit determination and control[J].Journal of Guidance,Control and Dynamics,1996,19(2):266-273.
[6]Rahmani A,JalaliM A,Pourtakdoust S H.Optimal approach to halo orbit control[C].AIAA Guidance,Navigation,and Control Conference and Exhibit,Austin,2003.
[7]Xin M,Dancer M W.Station-keeping of an L2Libration point satellite withθ-D technique[C].The 2004 American Control Conference,Boston,2004:1037-1042.
[8] 胡少春,孫承啟,劉一武.基于序優(yōu)化理論的暈軌道轉移軌道設計[J].宇航學報,2010,31(3):662-668.[Hu Shaochun,Sun Cheng-qi,Liu Yi-wu.Transfer trajectory design for halo orbit based on ordinal optimization theroy[J].Journal of Astronautics,2010,31(3):662-668.]
[9] 俞輝,寶音賀西,李俊峰.雙三體系統(tǒng)不變流形拼接成的低成本探月軌道[J].宇航學報,2007,28(3):637-642.[Yu Hui,BaoYin He-xi,Li Jun-feng.Low energy transfer to the Moon using the patching of invariant manifolds of two there-body systems[J].Journal of Astronautics,2007,28(3):637-642.]
[10] 闕志宏,周鳳岐,羅健,等.線性系統(tǒng)理論[M].西安:西北工業(yè)大學出版社,1994.
[11] 程國采.彈道導彈制導方法與最優(yōu)控制[M].長沙:國防科技大學出版社,1987.