韓明仁, 王玉峰,*

(1. 北京控制工程研究所, 北京 100094; 2. 空間智能控制技術(shù)重點實驗室, 北京 100094)

0 引 言

全電推進(jìn)衛(wèi)星目前主要應(yīng)用于地球同步軌道(geostationary orbit, GEO)相關(guān)任務(wù)中,使用電推力器實現(xiàn)全部的軌道轉(zhuǎn)移、位置保持、角動量卸載等功能。得益于電推力器的高比沖和高燃料利用率,全電推進(jìn)衛(wèi)星能夠搭載更多的有效載荷。然而,全電推進(jìn)衛(wèi)星由于推力小,依靠電推力器點火進(jìn)入GEO的時間要長達(dá)幾個月或半年多的時間。因此,優(yōu)化變軌策略,縮短變軌時間,是全電推進(jìn)衛(wèi)星研究中的關(guān)鍵問題之一。

全電推進(jìn)衛(wèi)星軌道優(yōu)化問題屬于最優(yōu)控制問題的范疇,本文主要研究全電推進(jìn)衛(wèi)星在轉(zhuǎn)移軌道段以時間最優(yōu)為目標(biāo)的變軌策略優(yōu)化方法,初始軌道為星箭分離后的同步轉(zhuǎn)移軌道(geostationary transfer orbit,GTO),目標(biāo)軌道為GEO。該問題要求在滿足系統(tǒng)約束和動力學(xué)約束的前提下,在整個轉(zhuǎn)移過程中確定一個推力方向變化策略。高度非線性和非凸動力學(xué)、空間環(huán)境攝動以及局部極小值的存在使優(yōu)化過程十分復(fù)雜。連續(xù)小推力變軌問題目前已存在的研究成果主要包括解析方法和數(shù)值方法兩大類。解析方法通過理論推導(dǎo)獲得最優(yōu)軌道的解析解,但對于大多數(shù)航天器的軌道優(yōu)化問題,可行性得不到保證。大多數(shù)研究者都致力于數(shù)值方法的研究,如直接法、間接法、反饋控制法等。直接法的基本思想是將控制變量離散化,使最優(yōu)控制問題轉(zhuǎn)化成一個非線性規(guī)劃問題,然后使用現(xiàn)有的優(yōu)化方法對其進(jìn)行求解,Ricciardi等采用結(jié)合有限元素全局多目標(biāo)優(yōu)化方法、啟發(fā)式的進(jìn)化算法以及基于梯度的約束非線性規(guī)劃方法開發(fā)了用于空間系統(tǒng)任務(wù)設(shè)計、優(yōu)化的多目標(biāo)直接混合最優(yōu)控制(multi objective direct hybrid optimal control，MODHOC)工具箱。Pritchett等研究了基于配點法的小推力轉(zhuǎn)移軌道的直接優(yōu)化方法,并設(shè)計了軌道鏈疊加法來進(jìn)行初值猜測,將初始軌道和目標(biāo)軌道中間某些弧段的狀態(tài)作為連接點,將初始猜測值引導(dǎo)到特定的解。文獻(xiàn)[9]介紹了空中客車防務(wù)及航天公司基于直接多重打靶法開發(fā)的小推力轉(zhuǎn)移軌道優(yōu)化軟件OptElec,在考慮地球陰影區(qū)和部分?jǐn)z動的情況下能夠設(shè)計從近地軌道(low earth orbit,LEO)或GTO轉(zhuǎn)移至GEO的時間或燃料最優(yōu)的變軌策略。間接法依賴于龐特里亞金極小值原理,將原最優(yōu)控制問題轉(zhuǎn)化成多點邊值問題,再利用微分方程邊值問題的求解方法進(jìn)一步得到最優(yōu)解。Mazzini等在考慮攝動和地影的情況下采用平均軌道根數(shù)法簡化動力學(xué)方程,通過打靶法求解了經(jīng)間接法轉(zhuǎn)化后的兩點邊值問題,得到了GEO衛(wèi)星依靠電推力器提升軌道的最優(yōu)策略。Bastante等研究了基于序列梯度修復(fù)算法的GTO-GEO小推力軌道轉(zhuǎn)移間接優(yōu)化方法,并討論了如何處理推力方向變化的連續(xù)性、最大角速度等約束條件。段傳輝等則提出了推力同倫的方法求解兩點邊值問題,簡化了間接法中對于協(xié)態(tài)變量初值猜測的復(fù)雜度,并以GTO-GEO小推力軌道提升問題為例對算法進(jìn)行了驗證,通過5次推力縮減,最終得到了時間最優(yōu)的轉(zhuǎn)移策略。文獻(xiàn)[12-14]中則采用了反饋控制法,設(shè)計了基于李雅普諾夫函數(shù)的近端商制導(dǎo)律(proximity quotient guidance law,Q-law),作為小推力軌道轉(zhuǎn)移的控制策略。其中,文獻(xiàn)[12]和文獻(xiàn)[13]中對GTO-GEO的小推力轉(zhuǎn)移過程進(jìn)行了仿真,表明了該方法能夠?qū)崿F(xiàn)時間次優(yōu)的軌道轉(zhuǎn)移閉環(huán)控制。

上述研究中所述的方法雖然均能夠得到優(yōu)化的小推力變軌策略,但是在復(fù)雜度、靈活性、最優(yōu)性等方面卻各有優(yōu)劣,直接法的變量規(guī)模巨大,計算復(fù)雜;間接法中協(xié)態(tài)變量初值猜測較為困難;反饋控制法則無法保證策略的最優(yōu)性。近年來有越來越多的研究者都采用強(qiáng)化學(xué)習(xí)方法來解決控制與決策問題,Bertsekas更是全面研究和闡述了強(qiáng)化學(xué)習(xí)與最優(yōu)控制的理論聯(lián)系,同時也有研究者采用強(qiáng)化學(xué)習(xí)方法解決行星際軌道優(yōu)化問題?？梢缘弥捎脧?qiáng)化學(xué)習(xí)方法解決復(fù)雜的最優(yōu)控制問題在簡易性、靈活性和最優(yōu)性方面均具有比較出色的表現(xiàn),能夠處理較為復(fù)雜的約束。因此,本文采用了強(qiáng)化學(xué)習(xí)的方法,求解小推力變軌最優(yōu)策略,在考慮攝動、地影等多種復(fù)雜約束的情況下,訓(xùn)練智能體的神經(jīng)網(wǎng)絡(luò),得到軌道轉(zhuǎn)移決策模型。

本文的主要貢獻(xiàn)如下:①在演員-評論家(Actor-Critic)強(qiáng)化學(xué)習(xí)框架下建立連續(xù)小推力變軌決策模型,結(jié)合廣義優(yōu)勢估計和近端策略優(yōu)化方法訓(xùn)練深度神經(jīng)網(wǎng)絡(luò),以其為非線性逼近器,探索小推力軌道轉(zhuǎn)移的最優(yōu)推力方向變化策略;②以軌道攝動方程、空間環(huán)境攝動、地球陰影區(qū)約束的數(shù)學(xué)模型為基礎(chǔ),建立可與智能體交互的環(huán)境模型;③針對由于狀態(tài)空間過大而導(dǎo)致智能體探索困難的問題,提出了動作網(wǎng)絡(luò)輸出映射和分層獎勵的方法,縮小了探索空間,顯著地提高了智能體探索到目標(biāo)軌道的效率,有效縮短了訓(xùn)練時間。最后,通過幾組不同初始參數(shù)的結(jié)果對比,驗證了求解得到的變軌策略的最優(yōu)性。

1 問題描述

1.1 預(yù)備知識

連續(xù)小推力軌道轉(zhuǎn)移問題中,作用于衛(wèi)星的推力加速度量級一般為10～10km/s,稍大于作用于衛(wèi)星的空間攝動力的大小。因此,若采用經(jīng)典軌道六根數(shù)來描述衛(wèi)星的軌道,那么衛(wèi)星的軌道動力學(xué)模型應(yīng)采用高斯型拉格朗日行星攝動力方程表示:

(1)

式中:

(2)

為半長軸;為偏心率;為軌道傾角;為升交點赤經(jīng);為近地點幅角;為真近點角;為軌道半通徑;為衛(wèi)星質(zhì)心到地心的距離;為地球引力常數(shù);、、分別為作用于衛(wèi)星的合力在衛(wèi)星軌道坐標(biāo)系3個坐標(biāo)軸方向產(chǎn)生的加速度,分別對應(yīng)徑向、軌道面內(nèi)垂直于徑向、軌道面法向3個方向。

本文采用強(qiáng)化學(xué)習(xí)的方法求解小推力變軌問題。強(qiáng)化學(xué)習(xí)方法通過智能體與環(huán)境的交互過程獲得回報,如圖1所示,從而更新智能體的策略,在實現(xiàn)特定目標(biāo)的同時使回報最大化,進(jìn)而得到最優(yōu)策略。

圖1 智能體與環(huán)境交互圖Fig.1 Interaction between agent and environment

在智能體與環(huán)境的交互中,智能體需要多次重復(fù)從初始狀態(tài)到目標(biāo)狀態(tài)的過程以獲得訓(xùn)練數(shù)據(jù),重復(fù)一次稱為一回合,而一回合中又包含很多步,用和分別表示智能體在第步的狀態(tài)和動作,此過程將產(chǎn)生一組狀態(tài)-動作序列(,,…,,),其中(,)表示智能體探索到特定目標(biāo)前的最后一對狀態(tài)-動作,狀態(tài)與動作既可以是標(biāo)量,也可以是矢量,需要根據(jù)實際問題定義。該狀態(tài)-動作序列即代表了智能體在一回合中的軌跡,用表示。可以定義第步的累積回報函數(shù)為

(3)

式中:()表示軌跡的回報;(,)表示在狀態(tài)下采取動作的單步回報。在基于策略的強(qiáng)化學(xué)習(xí)方法中,目標(biāo)函數(shù)可以表示為

(4)

式中：變量泛指智能體從狀態(tài)到動作的映射中涉及到的參數(shù),參數(shù)不同即代表不同的策略；E(·)表示期望；表示在參數(shù)下的策略；(|)表示在參數(shù)下的策略在條件下執(zhí)行的概率。強(qiáng)化學(xué)習(xí)的最終目標(biāo)則是找到最優(yōu)的參數(shù),使得目標(biāo)函數(shù)達(dá)到最大值,用公式描述為

(5)

實際上,強(qiáng)化學(xué)習(xí)問題變成了一個優(yōu)化問題,采用文獻(xiàn)[27]提出的策略梯度優(yōu)化方法,即

(6)

(7)

式中：(,)為評價函數(shù)。在Actor-Critic框架下,(,)為Actor,(,)則為Critic,二者均采用神經(jīng)網(wǎng)絡(luò)來表示,分別稱為動作網(wǎng)絡(luò)和評價網(wǎng)絡(luò)。Actor的實際輸出值為歸一化的值,一般取值范圍為[-1,1],在與環(huán)境交互前需進(jìn)行變換,映射到相應(yīng)的動作空間中。Critic的輸出值為當(dāng)前狀態(tài)及動作的價值。

1.2 問題建模

本文研究全電推進(jìn)衛(wèi)星轉(zhuǎn)移至GEO的時間最短變軌策略,衛(wèi)星的動力學(xué)模型采用式(1)中的微分方程描述,式中的加速度項可以拆分為電推力器提供的加速度和攝動加速度兩項,即

(8)

式中：為在衛(wèi)星軌道坐標(biāo)系中衛(wèi)星的合加速度矢量,=(,,)；為電推力器產(chǎn)生的合力大小；為當(dāng)前時刻的質(zhì)量；為推力在衛(wèi)星軌道坐標(biāo)系下的單位矢量;則為衛(wèi)星所受的攝動加速度,包括地球非球形引力攝動、日月引力攝動、太陽光壓攝動、大氣阻力攝動等。

在變軌過程中,衛(wèi)星質(zhì)量隨著工質(zhì)的消耗而減小,變化率為

(9)

式中：為電推力器的比沖；為重力加速度。

假設(shè)表示推力單位矢量在軌道面內(nèi)的投影與徑向的夾角,表示推力單位矢量與軌道面的夾角,則可以寫為

=(coscos,sincos,sin)

(10)

衛(wèi)星的推力指向即可用和兩個參數(shù)表示,若不考慮工程上的推力指向范圍約束,可取∈[-90°,270°],∈[-90°,90°]。因此,全電推進(jìn)衛(wèi)星的變軌策略優(yōu)化也就是對、兩參數(shù)的變化策略進(jìn)行優(yōu)化。因此,問題可轉(zhuǎn)化為如下優(yōu)化指標(biāo):

(11)

令=(,,,,,),則系統(tǒng)的狀態(tài)方程可以寫為

(12)

式中：為6×3矩陣,為6×1矩陣,具體表達(dá)式見文獻(xiàn)[28],=(,)∈為狀態(tài)空間,=(,)∈為動作空間,優(yōu)化目標(biāo)為時間最短,任務(wù)目標(biāo)軌道為GEO。

2 軌道優(yōu)化方法

2.1 環(huán)境模型

強(qiáng)化學(xué)習(xí)算法要求建立可供與智能體交互產(chǎn)生數(shù)據(jù)的環(huán)境模型,即圖1中的環(huán)境模塊,其輸入是智能體的動作,輸出為采取動作后達(dá)到的狀態(tài)以及當(dāng)前動作的獎勵。由于目標(biāo)軌道是GEO,其傾角為零且為圓軌道。當(dāng)達(dá)到目標(biāo)軌道時,軌道的升交點赤經(jīng)和近地點幅角產(chǎn)生了奇異。為了避免這種奇異,本文采用改進(jìn)的春分點軌道根數(shù)(,,,,,)描述強(qiáng)化學(xué)習(xí)過程中的狀態(tài),與經(jīng)典軌道六根數(shù)轉(zhuǎn)換關(guān)系如下:

(13)

相應(yīng)的狀態(tài)轉(zhuǎn)移方程為

(14)

式中:為輔助量,定義如下:

=1+cos+sin

(15)

本文僅在智能體訓(xùn)練中采用式(14)替代式(12)作為與智能體交互的狀態(tài)轉(zhuǎn)移模型,在分析問題時仍采用式(1)。

為更加貼近真實情況,使得到的策略更具工程意義,需要考慮攝動和地影的影響。考慮的攝動包括地球非球形引力攝動、日月引力攝動、太陽光壓攝動、大氣阻力攝動等,其數(shù)學(xué)模型對應(yīng)文獻(xiàn)中一致,這里不過多贅述。將得到的各種攝動加速度經(jīng)過坐標(biāo)變換并求和,可得到衛(wèi)星軌道坐標(biāo)系下的攝動加速度,代入式(8)中,即可得到衛(wèi)星在空間中的合加速度。

此外,地球陰影區(qū)對全電推進(jìn)衛(wèi)星變軌的影響也需要考慮在內(nèi)。由于電推力器主要靠太陽能轉(zhuǎn)化成的電能驅(qū)動,因此當(dāng)進(jìn)入地影區(qū)時,太陽帆板被遮擋,衛(wèi)星須關(guān)閉電推力器滑行。本文采用圓錐形地影模型,通過太陽矢量和衛(wèi)星位置矢量計算出地球-衛(wèi)星-太陽夾角、衛(wèi)星上觀察太陽的視半徑以及衛(wèi)星上觀察地球的視半徑。通過三者關(guān)系即可得知衛(wèi)星的地影情況,用可視因子表示太陽被地球遮擋的程度,有以下關(guān)系:

(16)

取值范圍為[0,1]。其中,與為中間變量:

(17)

如圖2所示,當(dāng)=1時,衛(wèi)星位于陽照區(qū),地球?qū)μ枱o遮擋;當(dāng)=0時,衛(wèi)星位于本影區(qū),地球完全遮擋太陽;當(dāng)0<<1時,衛(wèi)星位于半影區(qū),地球?qū)μ柕恼趽醭潭入S著的值減小而增加。

圖2 圓錐形地影模型示意圖Fig.2 Schematic diagram of conical earth shadow model

衛(wèi)星位于半影區(qū)內(nèi)太陽可視因子較大的區(qū)域時,電推力器也能夠點火。本文假設(shè)>01時,電推力器可點火,當(dāng)≤01時，電推力器關(guān)機(jī),衛(wèi)星滑行。

2.2 改進(jìn)的近端策略優(yōu)化方法

本文采用結(jié)合廣義優(yōu)勢估計(generalized advantage estimator, GAE)的近端策略優(yōu)化(proximal policy optimization, PPO)(GAE&PPO)方法對所提出的問題進(jìn)行求解。PPO屬于強(qiáng)化學(xué)習(xí)方法的一種,基于Actor-Critic框架而提出,其核心優(yōu)化公式如下:

(18)

其中,

(19)

式中:()為概率比,表示新老策略的差異;clip(·)為剪切函數(shù),用來將()限制在[1-,1+]之間,如果超過邊界值則取邊界值;為剪切率,限制()的剪切范圍,代表了更新策略的探索性,越大,算法越容易探索到新的狀態(tài),但影響穩(wěn)定性。優(yōu)勢函數(shù)定義為

(,)=(,)-()

(20)

(,)表示當(dāng)前步狀態(tài)-動作值函數(shù)(,)與狀態(tài)值函數(shù)()之差,由評價網(wǎng)絡(luò)產(chǎn)生,用來評價在步時采取動作的優(yōu)劣,為值函數(shù)的相關(guān)參數(shù)。

值函數(shù)的參數(shù)屬于評價網(wǎng)絡(luò)中的參數(shù),其更新方法同樣采用策略梯度法：

(21)

(22)

式中:()為值函數(shù)的目標(biāo)函數(shù);為第步的累積獎勵值。

在評價網(wǎng)絡(luò)估計優(yōu)勢函數(shù)時,常常用逼近算法近似計算值函數(shù),因此會引入偏差,文獻(xiàn)[37]中提出了廣義優(yōu)勢函數(shù)估計:

(23)

GAE&PPO算法流程如下。

算法1 GAE&PPO算法1 初始化動作網(wǎng)絡(luò)πφ(s,u)和評價網(wǎng)絡(luò)Ψπ(s,u)的參數(shù)φ0和?0For k=0,1,2,…,M2 策略πφk(s,u)與環(huán)境交互,收集軌跡集合Dk={τi}3 計算每一步累積獎勵{Rt}4 通過式(23)計算每一步優(yōu)勢值{A^t}5 采用策略梯度方法,最小化式(18),更新動作網(wǎng)絡(luò)參數(shù)φk6 采用策略梯度方法,最小化式(22),更新評價網(wǎng)絡(luò)參數(shù)?kEnd for

由于采用了概率分布和采樣的方式確定下一步動作,因此GAE&PPO算法的穩(wěn)定性和收斂性更好,算法具體流程如圖3所示。

圖3 GAE&PPO算法流程圖Fig.3 GAE&PPO algorithm flow diagram

3 訓(xùn)練加速方法

本文中智能體要在七維狀態(tài)空間中探索到目標(biāo)狀態(tài),尋找這一過程中時間最優(yōu)推力指向序列,而目標(biāo)狀態(tài)僅為GEO附近能夠進(jìn)行軌道捕獲的小范圍區(qū)域。智能體在狀態(tài)空間中很難探索到目標(biāo)狀態(tài),導(dǎo)致有效訓(xùn)練數(shù)據(jù)不足,訓(xùn)練過程十分緩慢甚至難以完成訓(xùn)練,這是強(qiáng)化學(xué)習(xí)方法解決小推力變軌問題的難點之一。本文提出了以下方法解決這一難點。

3.1 動作輸出映射

軌道轉(zhuǎn)移的目標(biāo)狀態(tài)限制了半長軸、偏心率和軌道傾角3個參數(shù),衛(wèi)星從初始軌道到目標(biāo)軌道需要完成提高半長軸、減小偏心率、壓低軌道傾角3方面內(nèi)容。本文通過分析軌道動力學(xué)方程,提煉出能夠?qū)崿F(xiàn)這3方面內(nèi)容的推力指向角、的粗略變化規(guī)律。根據(jù)提煉出來的先驗知識,對動作網(wǎng)絡(luò)的輸出進(jìn)一步設(shè)計,使半長軸、偏心率和軌道傾角3個參數(shù)始終向目標(biāo)值靠近,從而縮小探索空間。

由前文可知,動作網(wǎng)絡(luò)的實際輸出值需要做線性變換才能與環(huán)境交互,因此可以在這一過程中引入先驗知識。

暫時忽略式(8)中攝動加速度的影響,將式(10)代入式(8),得到推力加速度矢量三維分量關(guān)于的表達(dá)式

(24)

將式(24)代入式(1)中的半長軸、偏心率和軌道傾角微分方程,整理得到:

(25)

若使、、盡可能向目標(biāo)值靠近,須滿足以下條件:

(26)

由于,∈[-90°,90°],cos≥0,因此,將式(25)代入式(26)可以得到

(27)

由式(27)可知,半長軸、偏心率的增減與∈[0°,180°]有關(guān),軌道傾角的增減與∈[0°,180°]有關(guān),因此可以通過式(27)得到:

(1)當(dāng)∈[0°,180°]時,若使半長軸增大,則∈[0°,90°];若使偏心率減小,則∈[-90°,-]∪[180°-,360°-]。

(2)當(dāng)∈(180°,360°]時,若使半長軸增大,則∈[90°,180°];若使偏心率減小,則∈[180°-,360°-]∪[540°-,270°]。

(3)當(dāng)+∈[0°,90°]∪(270°,360°]時,若使軌道傾角減小,則∈[-90°,0°]。

(4)當(dāng)+∈(90°,270°]時,若使軌道傾角減小,則∈[0°,90°]。

對軌道攝動情況進(jìn)行分析,可以得知在變軌前期軌道高度較低的弧段,地球非球形引力攝動和大氣阻力產(chǎn)生的攝動加速度相對較大,因此變軌前期主要提升半長軸。

綜合以上分析可以得知,能夠?qū)崿F(xiàn)軌道轉(zhuǎn)移任務(wù)的策略具有以下幾個特點:① 變軌前期主要增大半長軸,推力在軌道面內(nèi)分量主要分布在速度方向附近;② 變軌后期進(jìn)一步提升半長軸,同時減小偏心率,推力在軌道面內(nèi)分量隨真近點角程周期變化;③ 整個變軌過程中推力矢量在軌道面法向的分量在升交點側(cè)的一半弧段指向負(fù)法向,在降交點側(cè)的一半弧段指向正法向。

假設(shè)動作網(wǎng)絡(luò)的實際輸出為(′,′),根據(jù)上述特點對其做以下變換:

(28)

(29)

式(28)和式(29)在最優(yōu)策略的必要條件下對、進(jìn)行粗略約束,不影響策略最優(yōu)性的同時使衛(wèi)星更容易探索到目標(biāo)軌道,從而獲得更多的有效訓(xùn)練數(shù)據(jù)。其中,變軌階段臨界點的最優(yōu)值在強(qiáng)化學(xué)習(xí)訓(xùn)練循環(huán)中采用隨機(jī)搜索的方法求解。

3.2 分層獎勵設(shè)計

在強(qiáng)化學(xué)習(xí)中,獎勵通過環(huán)境傳遞給智能體,智能體的訓(xùn)練目標(biāo)是使總收益最大化,因此獎勵函數(shù)的設(shè)計尤為重要。

當(dāng)前的優(yōu)化目標(biāo)是時間最優(yōu),一般的時間最優(yōu)獎勵函數(shù)形式為

(30)

式中：為一個較小值,代表耗時懲罰;為一個較大值,代表完成任務(wù)目標(biāo)的獎勵。

上述獎勵函數(shù)雖然能夠保證優(yōu)化目標(biāo)的實現(xiàn),但是訓(xùn)練過程中可能達(dá)不到任務(wù)目標(biāo),導(dǎo)致很多無效的訓(xùn)練。因此,根據(jù)本文實際問題,設(shè)計獎勵函數(shù)如下：

(31)

式中：、、和、、分別為當(dāng)前時間步和目標(biāo)的半長軸、偏心率和軌道傾角。其他情況中包括耗時懲罰和偏離任務(wù)目標(biāo)懲罰。

前文中獎勵函數(shù)雖然引入了偏離任務(wù)目標(biāo)的懲罰項,但是智能體中神經(jīng)網(wǎng)絡(luò)由于采用了隨機(jī)初始化參數(shù),訓(xùn)練效率受神經(jīng)網(wǎng)絡(luò)初始化影響較大,在訓(xùn)練初期僅通過的獎勵值首次探索到目標(biāo)軌道的用時不穩(wěn)定。因此,在的基礎(chǔ)上,設(shè)計了更加有利于探索到目標(biāo)的獎勵函數(shù):

(32)

式中：為取余符號。中每隔2 000個時間步計算一次當(dāng)前狀態(tài)與目標(biāo)的偏差,其余情況獎勵為0,取消對時間的尋優(yōu),只探索目標(biāo)狀態(tài),提高任務(wù)成功率。

結(jié)合和的優(yōu)勢,本文設(shè)計了分層獎勵方法,將訓(xùn)練過程分為兩個層次,下層訓(xùn)練采用作為獎勵,上層訓(xùn)練采用作為獎勵,訓(xùn)練初期從下層訓(xùn)練開始,優(yōu)先探索能夠達(dá)到目標(biāo)狀態(tài)的策略,當(dāng)某一回合智能體成功達(dá)到目標(biāo)狀態(tài)時,則在此基礎(chǔ)上開始上層訓(xùn)練,進(jìn)行時間最短策略尋優(yōu)。

采用上述方法優(yōu)化后的智能體與環(huán)境交互模型如圖4所示。

圖4 改進(jìn)的智能體與環(huán)境交互模型Fig.4 Improved agent-environment interaction model

4 結(jié)果分析及對比

4.1 訓(xùn)練結(jié)果

采用前文所提到的強(qiáng)化學(xué)習(xí)及訓(xùn)練加速方法對全電推進(jìn)衛(wèi)星軌道轉(zhuǎn)移問題進(jìn)行求解,驗證方法的可行性和最優(yōu)性。

假設(shè)在協(xié)調(diào)世界時(universal time coordinated, UTC)下的歷元時刻2021年1月1日12:00:00,衛(wèi)星位于初始軌道,其初始軌道和目標(biāo)軌道參數(shù)如表1所示,衛(wèi)星初始質(zhì)量為1 600 kg,攜帶兩臺裝有矢量調(diào)節(jié)機(jī)構(gòu)的離子電推力器,可產(chǎn)生推力合力大小為400 m·N,比沖為2 000 s,太陽光壓反射率為1,光壓等效面積為20 m,在低軌道時的大氣阻力系數(shù)為2.2,迎風(fēng)面積為20 m。

表1 衛(wèi)星軌道根數(shù)初始值和目標(biāo)值

當(dāng)衛(wèi)星達(dá)到目標(biāo)軌道附近且半長軸偏差小于0.5 km、偏心率小于0.1、軌道傾角小于0.1°時,即認(rèn)為當(dāng)前回合成功達(dá)到任務(wù)目標(biāo)。當(dāng)衛(wèi)星半長軸大于42 170 km、偏心率大于1、軌道傾角大于90°時,則認(rèn)為探索超出邊界,任務(wù)失敗。

算法中動作網(wǎng)絡(luò)和評價網(wǎng)絡(luò)均采用三隱層平行結(jié)構(gòu),每層神經(jīng)元節(jié)點數(shù)為128個,學(xué)習(xí)率均設(shè)置為10,激活函數(shù)采用ReLu和Hardswish兩種函數(shù),具體結(jié)構(gòu)如圖5所示。

圖5 智能體包含的神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)示意圖Fig.5 Schematic diagram of neural network contained in agent

訓(xùn)練過程中還涉及到前文中提到的一些其他參數(shù),具體參數(shù)值如表2所示。由于算法的穩(wěn)定性和收斂性較好,因此超參數(shù)設(shè)置在合理的范圍內(nèi)即可,選取較為容易。

表2 其余參數(shù)選取

采用Python語言實現(xiàn)網(wǎng)絡(luò)結(jié)構(gòu)并訓(xùn)練,用每回合累積獎勵的變化來表示訓(xùn)練效果,得到的訓(xùn)練曲線如圖6所示。

圖6 學(xué)習(xí)曲線Fig.6 Learning curve

由圖6可知,在訓(xùn)練初期,智能體未能探索到目標(biāo)狀態(tài),獎勵函數(shù)為,此時獎勵為狀態(tài)與目標(biāo)狀態(tài)偏差所產(chǎn)生的懲罰項,因此累積獎勵為負(fù)值。隨后,智能體通過最小化偏差探索到任務(wù)目標(biāo),獲得任務(wù)目標(biāo)獎勵,因此累計獎勵產(chǎn)生了較大的階躍,同時開始對轉(zhuǎn)移時間優(yōu)化。在此階段中,由于策略的探索性,智能體會探索到目標(biāo)狀態(tài)之外,因此會產(chǎn)生累計獎勵突然下降的情況,由于中同樣存在偏差項,智能體能夠重新找到目標(biāo)狀態(tài),從而繼續(xù)優(yōu)化任務(wù)時間,最終累積獎勵穩(wěn)定到一個范圍內(nèi),從中選取最大累積獎勵對應(yīng)的參數(shù),得到動作網(wǎng)絡(luò)模型,即可得到時間最優(yōu)的轉(zhuǎn)移策略。最大累積獎勵為176.84,對應(yīng)的軌道轉(zhuǎn)移時間為142.11 d。

將得到的動作網(wǎng)絡(luò)與環(huán)境交互,即可得到時間最優(yōu)的推力方向變化序列以及相應(yīng)的軌道變化序列,其中半長軸、偏心率、軌道傾角變化曲線如圖7所示。

圖7 衛(wèi)星軌道根數(shù)變化曲線Fig.7 Orbit elements curve of satellite

衛(wèi)星變軌軌跡如圖8所示,其中藍(lán)色部分為陽照區(qū)電推力器點火弧段,紅色部分為地影區(qū)電推力器關(guān)機(jī)滑行階段。

圖8 衛(wèi)星變軌軌跡Fig.8 Trajectory during satellite orbit transfer

4.2 效果對比

在不同的變軌問題以及約束條件下,對比本文方法與其他方法的優(yōu)化效果,采用不同文獻(xiàn)中的初始軌道參數(shù)作為初值,如表3所示,分別求解時間最優(yōu)的軌道轉(zhuǎn)移策略。

表3 參考文獻(xiàn)中的初始軌道參數(shù)

表4中列舉了采用文獻(xiàn)中各自算法求解得出的最優(yōu)轉(zhuǎn)移時間和采用強(qiáng)化學(xué)習(xí)方法求解得到的轉(zhuǎn)移時間。其中,情況1是在開普勒模型下的結(jié)果,忽略了地影和攝動的影響;情況2和情況3只考慮了地球的J2攝動項以及地影區(qū)約束;情況4則是在全引力模型以及地影區(qū)影響下的計算結(jié)果。

表4 不同情況下的軌道轉(zhuǎn)移時間對比

通過表4中的結(jié)果對比,可見在不考慮地影約束和攝動影響的情況下,采用本文所述方法得到的最優(yōu)軌道轉(zhuǎn)移時間與采用間接法得出的最優(yōu)結(jié)果基本一致。在僅考慮J2項非球形引力攝動和地影區(qū)約束的情況下,本文方法得到的最優(yōu)結(jié)果與文獻(xiàn)中的最優(yōu)結(jié)果相近,但略差于最優(yōu)結(jié)果,差距控制在3%以內(nèi)。在考慮全部攝動以及地影區(qū)影響的情況下,本文方法與反饋控制法計算得到的最優(yōu)結(jié)果相比有較大的提升。因此,由上述內(nèi)容可以得出,在不同初始軌道參數(shù)、考慮不同環(huán)境影響的情況下,所設(shè)計的軌道優(yōu)化算法仍然能夠計算出接近最優(yōu)的軌道轉(zhuǎn)移策略。

在求解效率方面,求解GTO-GEO小推力轉(zhuǎn)移軌道問題需要計算每一軌的推力方向變化策略,共包含了幾百甚至上萬個變量,求解計算量較大。直接法求解時涉及到的變量和約束多達(dá)幾十萬個,求解時間長達(dá)幾小時甚至一整天。間接法雖然求解計算量小,但是由于協(xié)態(tài)變量需要猜測初值,初值猜測不準(zhǔn)確將導(dǎo)致求解時間延長甚至無解。本文所提到的方法不需要同間接法一樣進(jìn)行準(zhǔn)確的初值猜測。同時,相比于直接法,本文方法涉及到的變量和約束少,求解時間更短,如圖6所示,經(jīng)過3 000個批次的訓(xùn)練即可達(dá)到最優(yōu)解,采用GPU(NVIDIA GeForce RTX 2080Ti)訓(xùn)練,求解時間為兩小時左右。

5 結(jié) 論

本文研究了將強(qiáng)化學(xué)習(xí)算法引入小推力軌道優(yōu)化問題中的具體實現(xiàn)方法。首先采用高斯型軌道動力學(xué)方程、各種空間攝動方程以及地影模型對強(qiáng)化學(xué)習(xí)所需的環(huán)境進(jìn)行了建模,得到了可以參與訓(xùn)練的環(huán)境模型。隨后,將GAE和PPO方法結(jié)合,得到了改進(jìn)的近端策略優(yōu)化算法,使算法更加穩(wěn)定,同時減小了值函數(shù)的偏差。針對訓(xùn)練中遇到的由于獎勵稀疏而難以收斂的問題,提出了動作輸出映射和分層獎勵等方法,顯著加快了收斂速度。最后,通過數(shù)值仿真,證明了算法的有效性。通過與其他文獻(xiàn)中采用不同方法計算得到的結(jié)果對比,表明了所提算法具有軌道轉(zhuǎn)移時間的最優(yōu)性。同時,所提方法既不需要對初值進(jìn)行猜測,也不需要優(yōu)化大量的變量,相比于其他方法更加簡便、靈活。本文對采用強(qiáng)化學(xué)習(xí)方法解決全電推進(jìn)衛(wèi)星軌道優(yōu)化問題做了初步研究,只考慮了變軌過程中空間環(huán)境的約束,并沒有考慮工程上的約束,如推力指向調(diào)節(jié)范圍、衛(wèi)星姿態(tài)角速度限制、燃料消耗量、同步軌道定點等約束。因此,研究帶有以上工程約束的小推力變軌策略是后續(xù)的研究方向之一。