李玖陽(yáng), 胡 敏, 王許煜, 徐家輝, 李菲菲

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

0 引 言

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

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

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

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

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

1 最優(yōu)控制率

1.1 攝動(dòng)方程

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

(1)

(2)

(3)

(4)

(5)

(6)

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

(7)

變換至an,t,h所在坐標(biāo)系下為

(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ì)量會(huì)隨著小推力器的工作而逐漸減少,小推力器燃料消耗的規(guī)律為

(11)

1.3 最優(yōu)控制率

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

(12)

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

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

(13)

(14)

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

2 ALPSO算法

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

2.1 PSO算法

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

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

(15)

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

(16)

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

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

(17)

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

2.2 增廣拉格朗日方法

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

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

(18)

(19)

(20)

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

2.3 算法流程

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

圖2 外部迭代程序框圖

算法的具體步驟如下。

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

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

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

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

圖3 子程序框圖

3 算例仿真

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

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

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

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

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

表2 算法參數(shù)

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

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

圖4 半長(zhǎng)軸在自然條件下隨時(shí)間變化

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

圖5 整個(gè)迭代過(guò)程半長(zhǎng)軸偏差

圖6 整個(gè)迭代過(guò)程近地點(diǎn)偏差

圖7 第一種處置軌道程序迭代過(guò)程

圖8 第二種處置軌道程序迭代過(guò)程

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

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

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

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

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

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

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

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

4 結(jié) 論

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