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

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

0 引 言

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

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

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

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

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

1 最優(yōu)控制率

1.1 攝動方程

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

(1)

(2)

(3)

(4)

(5)

(6)

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

(7)

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

(8)

1.2 推力方向和燃料消耗率

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

圖1 推力方向示意圖

可將推力矢量分解為

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

(9)

推力大小為

(10)

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

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

(11)

1.3 最優(yōu)控制率

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

(12)

式中,λa、λe和λi為協(xié)狀態(tài)變量。

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

(13)

(14)

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

2 ALPSO算法

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

2.1 PSO算法

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

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

(15)

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

(16)

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

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

(17)

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

2.2 增廣拉格朗日方法

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

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

(18)

(19)

(20)

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

2.3 算法流程

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

圖2 外部迭代程序框圖

算法的具體步驟如下。

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

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

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

步驟 5如果不滿足迭代終止條件,重復(fù)步驟2～步驟4,直到滿足迭代次數(shù)和約束閾值要求,輸出優(yōu)化結(jié)果。

圖3 子程序框圖

3 算例仿真

3.1 仿真條件和參數(shù)

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

表1 工作軌道與處置軌道的軌道根數(shù)

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

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

表2 算法參數(shù)

3.2 仿真結(jié)果及分析

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

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

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

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

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

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

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

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

通過仿真分析可以得出以下結(jié)論。

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

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

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

表3 收斂后最優(yōu)解

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

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

4 結(jié) 論

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