遲進(jìn)梓 余紅英 張子雄

中北大學(xué)，太原 030051

0 引言

目前，化學(xué)燃料推進(jìn)和電推進(jìn)是航天器的主要推力來源[1]。由于化學(xué)燃料推進(jìn)方式比沖小，能量密度低，在執(zhí)行任務(wù)時(shí)，攜帶大量推進(jìn)劑，增加了衛(wèi)星發(fā)射成本[2]。而電推進(jìn)屬于小推力推進(jìn)方式，由于其具有比沖大、控制精度高等特點(diǎn)，提高了探測(cè)器有效載荷的比重，具有優(yōu)秀的研究應(yīng)用前景[3]。電推進(jìn)技術(shù)[4]是利用電力技術(shù)使推進(jìn)劑膨脹，向后高速噴發(fā)，產(chǎn)生航天器所需推力的技術(shù)，適用于速度矢量較大以及保持軌道精度等變軌情況。

使衛(wèi)星等航天器從原軌道轉(zhuǎn)移至期望軌道上，此過程稱之為軌道機(jī)動(dòng)[5]，本文針對(duì)小推力類軌道機(jī)動(dòng)的軌道轉(zhuǎn)移問題進(jìn)行研究。針對(duì)衛(wèi)星等航天器模型的變軌問題，Kluever[6]提出了一種電推進(jìn)軌道轉(zhuǎn)換權(quán)衡的算法，主要對(duì)化學(xué)-電混合推進(jìn)的GEO軌道轉(zhuǎn)移方案進(jìn)行研究，實(shí)現(xiàn)了快速開展大推力軌道變換。Fakoor等[7]將解析方法與基于種群的人工蜂群算法相結(jié)合，通過人工蜂群算法搜索控制變量的初值和目標(biāo)函數(shù)的終止時(shí)間，給出了軌道轉(zhuǎn)移過程中狀態(tài)和控制變量的變化率，所得結(jié)果具有相當(dāng)?shù)氖諗啃院妥銐虻木?。Caillau等[8]對(duì)從GEO軌道到地月系L1點(diǎn)和月球軌道的航天器轉(zhuǎn)移軌道進(jìn)行相關(guān)研究，雖然能實(shí)現(xiàn)航天器的軌道轉(zhuǎn)移，但此方法并未對(duì)小推力轉(zhuǎn)移軌道優(yōu)化問題進(jìn)行研究。付磊等[9]建立了多沖量變軌模型，得到最優(yōu)沖量的時(shí)刻和大小，解決了遠(yuǎn)程導(dǎo)引多沖量變軌問題。黎桪等[10]創(chuàng)造性地使用兩次偽譜法優(yōu)化小推力軌道，較好克服一次偽譜法在求解燃料最優(yōu)時(shí)產(chǎn)生的震蕩問題。Pontani等[11]聯(lián)合使用VTD-NOG和PD-RM控制完成了低推力月球軌道轉(zhuǎn)移。

本文提出一種新的軌道轉(zhuǎn)移和規(guī)劃算法，將小推力衛(wèi)星變軌過程轉(zhuǎn)化為最優(yōu)控制中兩點(diǎn)邊值問題，然后針對(duì)該問題引入混合遺傳算法進(jìn)行解算，完成小推力航天器軌道轉(zhuǎn)移，實(shí)現(xiàn)了對(duì)小推力航天器由低軌向高軌的飛行軌道規(guī)劃和優(yōu)化設(shè)計(jì)。

1 動(dòng)力學(xué)模型

小推力航天器質(zhì)量隨燃料消耗而減小，在飛行時(shí)主要受到發(fā)動(dòng)機(jī)推力與地球重力作用，而春分點(diǎn)軌道根數(shù)與春分點(diǎn)坐標(biāo)系OE-xyz相聯(lián)系，坐標(biāo)系變換見文獻(xiàn)[2]。經(jīng)過計(jì)算，以春分點(diǎn)軌道根數(shù)和衛(wèi)星的燃料質(zhì)量[r′,xh,yh,n1,n2,λ,m]為狀態(tài)量，不考慮攝動(dòng)，改進(jìn)春分點(diǎn)根數(shù)下的動(dòng)力學(xué)模型為：

(1)

2 軌道轉(zhuǎn)移優(yōu)化問題

軌道轉(zhuǎn)移優(yōu)化是在已知航天器姿態(tài)以及初始狀態(tài)情形下，并在航天器滿足以春分點(diǎn)根數(shù)變分方程為基礎(chǔ)的航天器運(yùn)動(dòng)方程條件下，通過優(yōu)化算法得到變軌過程中燃料消耗和時(shí)長(zhǎng)的最優(yōu)解。

2.1 小推力航天器運(yùn)動(dòng)方程

本文用春分點(diǎn)根數(shù)變分方程表示航天器的運(yùn)動(dòng)方程：

(2)

以k=[r,xh,yh,n1,n2,λ]T表示t時(shí)刻軌道狀態(tài)向量，k可表示為位置和速度的函數(shù)：

(3)

則有：

(4)

其中有：

(5)

P與m分別表示推力矢量與衛(wèi)星的燃料質(zhì)量，R為從地心到在軌服務(wù)對(duì)象的位置矢徑。式(4)中的春分點(diǎn)根數(shù)相對(duì)于位置和速度的偏導(dǎo)數(shù)通過春分點(diǎn)根數(shù)間的泊松括號(hào)與位置速度對(duì)于春分點(diǎn)根數(shù)的偏導(dǎo)數(shù)聯(lián)系起來，求得春分點(diǎn)軌跡根數(shù)矢量狀態(tài)變分表達(dá)式：

(6)

由式(6)可知，求出春分點(diǎn)根數(shù)關(guān)于速度的偏導(dǎo)數(shù)，代入式(5)得到最終的軌道根數(shù)變分方程。

2.2 最優(yōu)控制問題描述

小推力變軌機(jī)動(dòng)的軌跡優(yōu)化，本質(zhì)上是求燃料消耗最少且用時(shí)最短的軌跡。由于研究問題的背景不同，約束條件不同，導(dǎo)致最優(yōu)控制問題不盡相同，需利用數(shù)值方法得解析解。本文選取的優(yōu)化量是推力方向矢量。而目標(biāo)函數(shù)分為燃料消耗以及時(shí)間消耗兩類。但如果變軌過程中持續(xù)推力一直為最大值，則兩類目標(biāo)函數(shù)效果是一致的。根據(jù)飛行任務(wù)，目標(biāo)衛(wèi)星最終狀態(tài)要符合前期需求中的限定因素。在最短時(shí)間的約束條件下，其模型性能指標(biāo)可表示為：

(7)

tf為最終時(shí)間，t0為初始時(shí)間。根據(jù)極大值原理以及式(1)中的燃料質(zhì)量方程和式(2)，哈密頓函數(shù)(H函數(shù))為：

(8)

式(8)中，λxT為協(xié)態(tài)變量矩陣，λλ和λm為協(xié)態(tài)變量。由最優(yōu)控制理論可知，若要達(dá)到最優(yōu)控制目標(biāo)，需使H函數(shù)可在某區(qū)間內(nèi)取到極小值。而由式(8)，航天器最佳推力方向應(yīng)與矢量λxTN反向，即：

(9)

P應(yīng)該滿足約束：

0≤P≤Pmax

(10)

P應(yīng)該在這一區(qū)間內(nèi)使H函數(shù)取到極小值，將H對(duì)P求偏導(dǎo)，得：

(11)

為H取極小值，選取合適點(diǎn)火開關(guān)函數(shù)[12]。協(xié)態(tài)方程為：

(12)

δ和η分別為推力的方位角和高低角，由δ和η確定推力方向角。協(xié)態(tài)變量在函數(shù)中還需要受到截?cái)嘁蛩氐挠绊懀?/p>

(13)

式(13)中，拉格朗日乘子用γ進(jìn)行表示。由于tf自由，因此在最佳軌道時(shí)刻，哈密頓方程理應(yīng)滿足以下表達(dá)：

H[k*(tf*)，λ*(tf*)，P*(tf*),tf*]=

(14)

將式(8)、式(13)第二式代入式(14)可得：

(15)

式(7)～(13)一起組成最短時(shí)間約束條件下最佳目標(biāo)參數(shù)函數(shù)方程組。本文求解該方程組思路是將狀態(tài)變量和協(xié)態(tài)變量轉(zhuǎn)化為兩點(diǎn)邊值問題，然后選用混合遺傳算法進(jìn)行求解。

3 算法流程

遺傳算法具有較強(qiáng)的全局遍歷能力，使時(shí)間效率得到有效提高，但其局部查找能力不強(qiáng)，而經(jīng)典優(yōu)化算法恰好具有與之相反的特點(diǎn)。因此本文取兩者進(jìn)行混合使用，使整個(gè)算法有良好的時(shí)間復(fù)雜度。

遺傳算法是一種通過分析系統(tǒng)整體特征，產(chǎn)生具有自我優(yōu)化參數(shù)能力的算法模型。由于當(dāng)前模型的特殊性，懲罰函數(shù)進(jìn)行如式(16)定義：

(16)

p為整數(shù)，通常取1或2。βj是表示約束權(quán)重的因子。dj(Z)=0(j=1,2,…,n)表示等式約束。廣義適應(yīng)度函數(shù)在懲罰函數(shù)約束下為：

F(Z)=J(Z)-φQ(Z)

(17)

式(17)中，φ為懲罰因子，當(dāng)φ>0，φ→+∞時(shí)，原函數(shù)收斂得到最終結(jié)果。但因子系數(shù)較高或較低會(huì)導(dǎo)致結(jié)果與真實(shí)結(jié)果值相去甚遠(yuǎn)，甚至破壞種群更迭條件。因此，本文使用具有自我調(diào)節(jié)能力退火篩選函數(shù)約束：T=T0，種群中懲罰因子為：

(18)

式(18)中，μ∈[0,1]表示溫度制冷速率常數(shù)；G代表種群繁衍層數(shù)，Tf是退火凝結(jié)溫度。式(16)中的βj選擇可以通過隨機(jī)權(quán)重的方法進(jìn)行獲取：

(19)

式(19)中，lj是非負(fù)隨機(jī)數(shù)。

4 仿真實(shí)驗(yàn)及結(jié)果

本章基于電推式小推力衛(wèi)星軌道規(guī)劃的理論基礎(chǔ)，在開源科學(xué)工程計(jì)算軟件SCILAB6.0.1上的SiROS下開發(fā)了仿真模塊。并利用仿真模塊搭建了GTO-GEO軌道轉(zhuǎn)移仿真模型，對(duì)電推式小推力衛(wèi)星變軌機(jī)動(dòng)方案的優(yōu)化算法進(jìn)行驗(yàn)證。

4.1 仿真模型搭建

本文4個(gè)仿真模塊是在SiROS下的宇航函數(shù)庫(kù)中開發(fā)的，分別是模塊1：初始軌道規(guī)劃模塊，模塊2：軌道優(yōu)化模塊，模塊3：軌道動(dòng)力學(xué)模塊，模塊4：變軌制導(dǎo)控制算法模塊。初始軌道規(guī)劃模塊封裝及連接如圖2所示。其是基于標(biāo)稱軌道法設(shè)計(jì)的，輸出為小推力軌道優(yōu)化結(jié)果，表示多圈小推力標(biāo)稱軌道上各分段節(jié)點(diǎn)的軌道根數(shù)初值。

圖1 混合遺傳算法流程圖

圖2 初始軌道規(guī)劃模塊封裝及連線

如圖3所示為軌道優(yōu)化模塊封裝及連接圖，軌道優(yōu)化設(shè)計(jì)模塊是根據(jù)本文優(yōu)化方法開發(fā)的，給出基于退火遺傳算法得到的小推力軌道優(yōu)化結(jié)果，其結(jié)果表示小推力軌道上各分段節(jié)點(diǎn)的軌道根數(shù)初值，參數(shù)值與初始標(biāo)稱軌道規(guī)劃模塊一致。

圖3 軌道優(yōu)化模塊封裝及連線

如圖4所示為軌道動(dòng)力模塊封裝及連接圖，將該模塊數(shù)據(jù)輸入模塊4，隨仿真時(shí)間的推進(jìn)不斷計(jì)算新的軌道根數(shù)、飛行器質(zhì)量、推力，直至仿真任務(wù)完成。該模塊的攝動(dòng)加速度設(shè)置為常量。

圖4 軌道動(dòng)力模塊封裝及連線

而變軌制導(dǎo)控制模塊用來對(duì)規(guī)劃好的軌道進(jìn)行跟蹤。變軌制導(dǎo)控制模塊封裝及連接如圖5所示。

圖5 變軌制導(dǎo)模塊封裝及連線

各模塊需要配置參數(shù)如表1，輸出參數(shù)如表2。

表1 各模塊參數(shù)配置

表2 各模塊輸出參數(shù)

4.2 軌道轉(zhuǎn)移仿真與分析

變軌開始時(shí)間為2020年1月1日0時(shí)0分，初始質(zhì)量為1500 kg，電推力比沖為3800，推力方位角上限為6.28 rad，推力上限為1 N。GTO-GEO軌道轉(zhuǎn)移仿真模型如圖6所示，各參數(shù)如表3所示。

表3 模型各參數(shù)初值

圖6 GTO-GEO軌道轉(zhuǎn)移仿真模型圖

為驗(yàn)證該算法的正確性和有效性，將本文所提方法解與最優(yōu)解進(jìn)行性能比較。圖7為半長(zhǎng)軸隨時(shí)間變化圖。圖8為偏心率矢量隨時(shí)間變化圖，圖9為傾角隨時(shí)間變化圖，圖10為推力隨時(shí)間變化圖。虛擬小推力飛行器模型在2020年01月01日00時(shí)00分出發(fā)，攜帶的初始燃料質(zhì)量為1500 kg，在2020年02月 13日21時(shí)36分到達(dá)，所用的時(shí)間是44.9天，消耗燃料為171.2 kg。

圖7 半長(zhǎng)軸隨時(shí)間變化圖

圖8 偏心率矢量隨時(shí)間變化圖

圖9 傾角隨時(shí)間變化圖

圖10 推力隨時(shí)間變化圖

圖11所示為推力示意圖，其中Ux，Uy和Uz分別為徑向、切向和法向小推力加速度，各推力分量變化如圖12～14所示，推力方位角變化如圖15～16所示。

圖11 推力示意圖

圖12 推力分量1示意圖

圖13 推力分量2示意圖

圖14 推力分量3示意圖

圖15 推力方位角1變化示意圖

圖16 推力方位角2變化示意圖

在航天器由GTO軌道轉(zhuǎn)移到目標(biāo)軌道時(shí)，由公式(1)可知，推力器產(chǎn)生的推力使軌道半長(zhǎng)軸增大，軌道偏心率降低，軌道傾角減小。由圖7～10可知，在小推力航天器進(jìn)行軌道變換期間，偏心率矢量與軌道傾角單調(diào)下降，而半長(zhǎng)軸單調(diào)上升，推力的方向和大小都成周期性變化，與理論分析相符。推力方位角1在[-45°，45°]范圍內(nèi)變化，推力方位角2在[-90°，90°]范圍內(nèi)變化，直至變軌結(jié)束。通過對(duì)GTO-GEO軌道轉(zhuǎn)移的過程進(jìn)行仿真與分析，可以得到仿真結(jié)果如表4所示。

其中，最優(yōu)解結(jié)果是不考慮對(duì)軌道產(chǎn)生顯著影響的攝動(dòng)因素的結(jié)果，混合動(dòng)力是本算法下化學(xué)燃料和電推進(jìn)混合推進(jìn)方式。由表4可知，本文所提算法消耗了占總質(zhì)量的11.4%的燃料，在燃料消耗上，僅比最優(yōu)解多1.2%，而相較于傳統(tǒng)遺傳算法，本文算法節(jié)省了3.7%的燃料，相較于混合動(dòng)力模型，本文算法節(jié)省了10.7%的燃料。本文所提算法完成變軌用時(shí)44.9天，時(shí)間消耗上，僅比最優(yōu)解多13.9%，而相較于傳統(tǒng)遺傳算法，本文算法節(jié)省了30.9%的時(shí)間，相較于混合動(dòng)力模型，本文算法節(jié)省了19.8%的時(shí)間。由此可見，本文所提算法能順利完成GTO軌道至GEO軌道的轉(zhuǎn)移任務(wù)，各模塊運(yùn)行正常，能以接近于最優(yōu)控制解的結(jié)論完成變軌，輸出預(yù)期的結(jié)果值，具有燃料消耗次優(yōu)的特點(diǎn)，降低了發(fā)射成本。

表4 GTO-GEO變軌燃料和時(shí)間消耗

5 結(jié)論

在連續(xù)小推力衛(wèi)星變軌以及制導(dǎo)控制過程中，本文將電推進(jìn)式小推力衛(wèi)星模型的變軌過程轉(zhuǎn)化為最優(yōu)控制中兩點(diǎn)邊值問題，然后針對(duì)該問題使用混合遺傳算法進(jìn)行解算。在SCILAB6.0.1上進(jìn)行實(shí)驗(yàn)，驗(yàn)證了本文提出的小推力衛(wèi)星變軌算法優(yōu)化設(shè)計(jì)方法。實(shí)驗(yàn)結(jié)果表明，變軌實(shí)際消耗燃料較最優(yōu)解僅多1.2%，消耗時(shí)間較最優(yōu)解僅多13.9%，且相較于傳統(tǒng)遺傳算法以及傳統(tǒng)動(dòng)力模型，燃料和時(shí)間消耗都大大減少，完成了小推力航天器由低軌向高軌的飛行軌道規(guī)劃及優(yōu)化。