基于ALPSO算法的低軌衛星小推力離軌最優控制方法

李玖陽, 胡 敏, 王許煜, 徐家輝, 李菲菲

(航天工程大學航天指揮學院, 北京 101416)

0 引 言

隨著人類航天發射任務的不斷增加,低軌及中高軌軌道空間逐漸擁擠,使在軌運行的衛星碰撞概率逐漸增大。根據機構空間碎片協調委員會(Inter-Agency Space Debris Coordination Committee,IADC)編訂的《IADC空間碎片減緩指南》[1],在自然軌道超過25年的情況下,低軌到壽衛星需要在25年內進入處置軌道進而在大氣層內燒毀[2]。同時,小衛星技術的發展和發射成本的降低使得低軌衛星星座中衛星數量劇增。OneWeb公司計劃發射由2 620顆衛星組成的星座[3]、Samgsung公司計劃發射由4 600顆衛星組成的星座[4]、Boeing計劃發射2 956顆衛星組成的星座,Starlink計劃發射4.2萬顆衛星組成星座,這些星座規模巨大,其衛星到壽后將會影響低軌的空間環境[5]。因此,采取合理方式對低軌衛星星座中到壽衛星進行離軌控制變得尤為重要。

針對低地球軌道(low earth orbit,LEO)衛星,離軌方法一般分為電動力電纜離軌、大氣阻力離軌和衛星自身推力器離軌[6-8],衛星自身推力器離軌相較于其他離軌方式具有離軌時間短、控制靈活等特點。電推力器作為一種典型的小推力器已經在深空探測、軌道維持和軌道機動中得到廣泛的應用[9],Starlink星座中所有衛星均帶有電推力器。Fromm等人[10]分別仿真不同小推力器作用下衛星的離軌時間,并得出了適于仿真對象衛星離軌的小推力器推力范圍。Huang等人[11]基于兩種不同的小推力離軌策略對單星和OneWeb星座離軌過程進行了仿真分析,并得出最優的單星及星座離軌策略。Trofimov等人[12]研究了被動穩定衛星的小推力離軌問題,對比了降低半長軸和降低近地點兩種離軌方式。以上離軌控制研究主要集中在控制率推導方面,對算法實現少有提及,粒子群優化(particle swarm optimization, PSO)算法在各類連續空間優化、神經網絡訓練等領域中取得良好效果[13],已經被多位學者用于求解轉移軌道最優控制問題[14-18]。文獻[19]基于多目標PSO算法解決燃料與時間同時最優的地-木和地-土轉移軌道。沈如松等人[20]采用多鄰域PSO算法對同步軌道入軌問題進行了仿真,但約束處理較為簡單。文獻[21]采用PSO算法與直接法結合的方法,但未考慮攝動影響。Wang等人[22]將PSO算法與序列二次規劃法相結合,提高了算法的全局和局部搜索能力,但所求解的精度不高。

本文主要研究了基于增廣拉格朗日PSO(augmented Lagrangian PSO, ALPSO)算法的低軌衛星小推力離軌最優控制問題。

小推力離軌控制需要將衛星從壽命末期軌道轉移至處置軌道,對這一過程的優化是一種小推力轉移軌道優化問題,主要通過數值解法解決,通常分為兩種[23]:① 采用龐特里亞金極大值原理,推導出狀態方程,結合打靶法進行求解,該方法的收斂對初值非常敏感,初值求解困難,可分為直接打靶法、直接配點法、偽譜法和微分包含法[24];② 將控制率參數化,通過懲罰函數等方式將約束融入目標函數,采用非線性規劃的方式求解,可結合遺傳算法、模擬退火和PSO算法求全局最優解,該方法避免了上一種方法初值求解困難的問題。深空探測軌道轉移、高軌衛星變軌等小推力轉移軌道優化設計問題的特點是約束復雜,使用間接法求解難度大,受攝動影響較小。低軌衛星離軌問題必須要考慮地球非球形攝動和大氣阻力攝動影響,使軌道演化過程的復雜性大大增加,采用遺傳算法、模擬退火等可能會出現罰因子過大的問題,無法得出正確結果,而ALPSO算法具有設置參數少,計算代價低,目標函數不易因為罰因子過大而陷入病態的特點,有效地避免了這一問題。

本文采用第二種方法得出了最優控制率,結合運動方程列出哈密爾頓函數,得出含協狀態參數的最優控制率;分別闡述PSO算法和增廣拉格朗日方法,得出具體的算法流程;利用兩種處置軌道的仿真算例驗證算法,并與遺傳算法優化結果對比,得出適于離軌的處置軌道。

1 最優控制率

1.1 攝動方程

小推力器推力大小為毫牛量級,相較于中心天體引力為小量,因此可將其作為攝動力處理,其對軌道根數的影響采用高斯型攝動方程[25],可表示為

(1)

(2)

(3)

(4)

(5)

(6)

低軌衛星在離軌過程中所受攝動力主要為地球非球形J2攝動和大氣阻力攝動,兩者對軌道根數的影響均可通過將攝動加速度代入式(1)～式(6)求得,由兩者造成的總攝動加速度為

(7)

變換至an,t,h所在坐標系下為

(8)

1.2 推力方向和燃料消耗率

設推力矢量在軌道平面內投影與速度方向的夾角為α,推力矢量與軌道平面的夾角為β,如圖1所示。

圖1 推力方向示意圖

可將推力矢量分解為

[an,at,ah]T=|Fthrust|[sinαcosβ, cosαcosβ, sinβ]T

(9)

推力大小為

(10)

式中,η為小推力器效率;P為小推力器功率;m為航天器質量;g為重力加速度;Isp為小推力器比沖。

航天器的總質量會隨著小推力器的工作而逐漸減少,小推力器燃料消耗的規律為

(11)

1.3 最優控制率

軌道轉移的目標軌道由半長軸、偏心率和軌道傾角所確定,因此由變分法可以得到哈密爾頓函數[26]為

(12)

式中,λa、λe和λi為協狀態變量。

由哈密爾頓方程達到最優控制條件:?H/?α=0,?H/?β=0,可得到噴射角度的最優控制率[26]分別為

(13)

(14)

將式(13)和式(14)代入式(9),即為最優控制量u*,最優控制率中含有3個協狀態變量,為了確定最優控制率,需要借助ALPSO算法在一定范圍內尋找滿足約束條件的最優解。

2 ALPSO算法

本文使用文獻[27]提出的ALPSO算法進行低軌衛星小推力離軌問題中最優控制的計算。文獻[27]對該算法進行了系統的闡述,在處理帶約束優化問題中具有較快的收斂速度和較高的精度。

2.1 PSO算法

Kennedy等人[28]于1995年在模擬鳥群尋找棲息地這一行為的基礎上提出了PSO算法,該算法可以在方程梯度信息未知的情況下解決不可微分方程的連續非凸性問題[29]。因此,PSO算法被廣泛應用在非線性函數優化領域。

設向量λi=[λa,λe,λi],該向量代入式(13)和式(14)代表某一種最優控制率,在PSO算法中λi代表一個粒子,向量λi的集合即為粒子群。粒子群中任意粒子都具有位置和速度兩個特性,位置和速度更新式[30]為

(15)

上述方法為標準PSO算法迭代方法,該方法在初期搜索速度較快,后期搜索速度較慢,當某粒子當前位置與個體歷史最佳位置和種群歷史最佳位置相等時,如果速度很小,會造成算法提前收斂。針對這些問題,粒子的第k+1次位置和速度[31]可以變為

(16)

式中,m1為[0,1]之間的隨機數;ρk為隨機方向參數。

為了避免過早收斂的問題,慣性權重ω采用線性微分遞減的方法,即

(17)

慣性權重ω在初期變化緩慢,使粒子具有較好的全局探索能力,而在后期變化速度較快,使粒子能快速收斂于全局最優值。隨機方向參數ρk使粒子在種群歷史最優解附近進行隨機搜索,如式(17)所示,如果連續多次種群歷史最佳位置相同并超出閾值,搜索空間會減小,反之會擴大搜索空間,該方法進一步提高了粒子的搜索能力,避免在粒子速度較低時陷入局部最優解。

2.2 增廣拉格朗日方法

上述基本PSO算法只能解決無約束優化問題,而離軌最優控制問題為帶約束的優化問題,解決帶約束優化問題通常采用罰函數法。但隨著迭代次數的增加,罰函數中的罰因子會逐漸趨于無窮,使目標函數病態逐漸加重,影響算法的收斂性。理論研究結果表明,增廣拉格朗日函數法有效地克服了罰函數法的缺點[32]。

增廣拉格朗日函數法是在拉格朗日函數上增加一個與約束相關的二次罰函數項。因此,結合KKT條件,在目標函數中加入二次罰函數項,保證約束可行性[27],可得

(18)

(19)

(20)

為了加快算法的運行速度并減少目標函數的約束,一些約束條件通過相應的處理可以在軌道積分程序中實施。由于粒子的初始位置值是隨機的,某些粒子在軌道積分時可能會造成軌道高度增加和偏心率異常等異常狀況,針對這一問題,可以在軌道積分程序中設置中斷條件,并將離軌時間返回為無窮大值。這些粒子的適應度值也就變為無窮大,從而在適應度篩選中淘汰。

2.3 算法流程

通過上述分析,ALPSO算法可將帶約束優化問題有效轉化為多個無約束優化的子問題,進而采用PSO算法求解子問題,算法的外部迭代程序和子程序框圖分別如圖2和圖3所示。

圖2 外部迭代程序框圖

算法的具體步驟如下。

步驟 1參數值初始化。首先將粒子分為3個維度,分別代表協狀態變量λa、λe和λi,而后在可行域范圍內,隨機初始化種群中所有粒子的位置和速度,并分別將拉格朗日乘子和懲罰因子初始化為0和r0。

步驟 3更新粒子參數。依據式(16)更新粒子的位置和速度,再次進行循環,直到達到最大迭代次數kmax,輸出本次子程序優化結果。

步驟 4更新增廣拉格朗日函數參數。利用子程序輸出的優化結果帶入式(20)更新拉格朗日乘子和懲罰因子,生成新的目標函數。

步驟 5如果不滿足迭代終止條件,重復步驟2～步驟4,直到滿足迭代次數和約束閾值要求,輸出優化結果。

圖3 子程序框圖

3 算例仿真

3.1 仿真條件和參數

本文的仿真算例考慮將衛星軌道高度從821 km的工作軌道降低至150 km的處置軌道或增大偏心率使衛星近地點軌道高度降至150 km,工作軌道和兩種處置軌道的半長軸、偏心率和軌道傾角如表1所示。

表1 工作軌道與處置軌道的軌道根數

衛星質量1 000 kg,大氣阻力系數為2,受曬面質比為0.02,小推力發動機推力為3 mN,效率為24.5%,比沖為4 660 s,功率為115 W。衛星在離軌過程中所受攝動力為地球非球形J2攝動和大氣阻力攝動,大氣密度模型采用美國標準大氣模型SA76。

PSO算法參數和初始懲罰因子設置如表2所示。針對高維問題,s、f和隨機方向參數初值ρ0推薦選取表中數值[31]。迭代終止條件為半長軸偏差絕對值或近地點偏差絕對值小于約束閾值εf。

表2 算法參數

3.2 仿真結果及分析

軌道積分程序仿真了衛星在僅受大氣阻力和地球非球形J2攝動下的軌道半長軸演化結果。半長軸為7 200 km的衛星軌道高度在800 km以上,大氣密度極低,在自然條件下軌道半長軸的變化如圖4所示,離軌進入處置軌道的時間為140年,遠遠超出了IADC規定的25年內離軌的標準。

圖4 半長軸在自然條件下隨時間變化

在其他條件相同的情況下,針對第一種和第二種處置軌道采用遺傳算法的迭代過程分別如圖5和圖6所示,在達到100次最大迭代次數時,偏差分別為66 m和9 608 m,并未達到約束閾值。針對第一種處置軌道,半長軸偏差絕對值隨迭代次數變化如圖7所示,從圖7中可以看出,程序總共迭代60次,在第57次時滿足收斂條件,在第20次至第30次迭代之間程序并未滿足收斂條件,造成后續偏差值小幅提高。第二種處置軌道近地點偏差隨程序迭代次數變化如圖8所示,程序總共迭代25次,在第12次時滿足收斂條件,程序收斂速度較快。

圖5 整個迭代過程半長軸偏差

圖6 整個迭代過程近地點偏差

圖7 第一種處置軌道程序迭代過程

圖8 第二種處置軌道程序迭代過程

兩種處置軌道的最優解、離軌時間和燃料消耗如表3所示,從表3數據得知,第二種處置軌道的離軌時間比第一種長,且燃料消耗比較大。以兩種處置軌道為目標,衛星在最優控制率作用下軌道參數隨時間變化分別如圖9和圖10所示。從圖9和圖10可以看出,兩種情況下軌道傾角基本不變,以第一種處置軌道為目標軌道時,偏心率有小幅度的上升,但對半長軸和軌道傾角的變化幾乎無影響,以第二種處置軌道為目標軌道時,偏心率隨時間變化呈近似線性,半長軸隨時間整體變化幅度較小,存在較小幅度波動,但不影響近地點高度隨時間的變化。

通過仿真分析可以得出以下結論。

(1) ALPSO算法可以有效地在給定協狀態量范圍的情況下找出最優的協狀態量,進而得出最優控制率。

(2) 相較于遺傳算法,ALPSO算法可以在較少的迭代次數下使最終的偏差更小。

(3) 在僅受大氣阻力和地球非球形J2攝動下,衛星不能在25年內實現離軌,而在星上小推力器的作用下,衛星以兩種處置軌道為目標均可滿足25年內離軌要求。其中,以第一種處置軌道為目標的情況下離軌時間較短,適合作為軌道高度821 km衛星的離軌處置軌道。

表3 收斂后最優解

圖9 最優控制率作用下衛星半長軸、偏心率、軌道傾角變化(第一種處置軌道)

圖10 最優控制率作用下衛星半長軸、偏心率、軌道傾角和近地點高度變化(第二種處置軌道)

4 結 論

本文將ALPSO算法應用于低軌衛星星座離軌最優控制中,相較于傳統打靶法,ALPSO算法對初值不敏感,在給定搜索范圍后就能進行搜索,大大減少了初值確定的困難。針對PSO算法的位置、速度迭代公式和慣性權重迭代公式進行了改進,減少了算法陷入局部最優解的可能性,提高了算法的整體搜索能力和全局收斂性。采用增廣拉格朗日方法引入了二次懲罰項,進而通過迭代校正拉格朗日乘子,避免了罰函數法因參數值過大引起的病態問題,并將一些約束通過軌道積分程序轉化到適應度函數中,減少了目標函數的約束并提高了程序的計算速度。針對兩種處置軌道,仿真對比了遺傳算法和ALPSO算法的優化結果,相較于遺傳算法,ALPSO算法迭代收斂更快,收斂時偏差較小。仿真對比了以兩種處置軌道為目標的離軌時間,以第一種處置軌道為目標離軌時間較短,適合作為軌道高度821 km衛星離軌的處置軌道。