基于脈沖初值的小推力轉(zhuǎn)移軌道優(yōu)化研究

王旭生，施偉璜，王 偉，彭玉明

(1.上海衛(wèi)星工程研究所，上海 201109； 2.上海市深空探測(cè)技術(shù)重點(diǎn)實(shí)驗(yàn)室，上海 201109)

0 引言

深空探測(cè)是當(dāng)今世界航天活動(dòng)的熱點(diǎn)研究領(lǐng)域，是一個(gè)國(guó)家綜合國(guó)力和創(chuàng)新能力的體現(xiàn)。我國(guó)已將深空探測(cè)納入國(guó)家重大科技專項(xiàng)，正在實(shí)施首次自主火星探測(cè)任務(wù)，并規(guī)劃了多目標(biāo)小行星探測(cè)、火星取樣返回、木星系及行星穿越等任務(wù)。其中，小行星探測(cè)任務(wù)的總速度增量要求大于8 km/s，由于傳統(tǒng)低比沖的化學(xué)推進(jìn)系統(tǒng)無(wú)法使用，因此小推力、高比沖、高可靠的電推進(jìn)系統(tǒng)成為工程實(shí)施的優(yōu)選動(dòng)力形式。

在小推力作用下，深空探測(cè)器的軌道優(yōu)化在本質(zhì)上是一個(gè)最優(yōu)控制問(wèn)題。目前求解方法主要有直接法、間接法和混合法[1]，它們都不可避免地需要進(jìn)行初值猜測(cè)。初值猜測(cè)不僅在理論上作為算法的啟動(dòng)值必不可少，而且在工程應(yīng)用上也是粗略優(yōu)化的過(guò)程，可提供任務(wù)可行性、性能指標(biāo)、各階段大致時(shí)間節(jié)點(diǎn)等關(guān)鍵信息，為進(jìn)一步優(yōu)化軌道提供參考。脈沖軌道作為小推力轉(zhuǎn)移軌道的初值，具有形式簡(jiǎn)單、計(jì)算方便的優(yōu)點(diǎn)，被廣泛應(yīng)用于小推力轉(zhuǎn)移軌道優(yōu)化的求解，如：美國(guó)噴氣推進(jìn)實(shí)驗(yàn)室(JPL)基于圓錐曲線拼接法開(kāi)發(fā)了MALTO軟件包的初值，采用一系列脈沖替代小推力軌道[2]，該方法已應(yīng)用于木星冰衛(wèi)星軌道器任務(wù)的軌道設(shè)計(jì)中[3]；陳揚(yáng)等[4]基于脈沖估算結(jié)果設(shè)計(jì)了電推進(jìn)軌道，采用間接法求解燃料最優(yōu)控制問(wèn)題，得到了小推力的最優(yōu)軌跡；李九天[5]提出了基于小推力可實(shí)現(xiàn)性的脈沖軌道約束，研究了基于該約束的燃料最優(yōu)脈沖轉(zhuǎn)移軌道設(shè)計(jì)方法；蔣小勇等[6]提出了一種基于多沖量能耗估算的任務(wù)窗口搜索方法，可用于搜索小推力任務(wù)的發(fā)射窗口。以上研究均基于單一形式的脈沖初值，在工程應(yīng)用中未考慮脈沖軌道的多樣性，缺少對(duì)不同形式的脈沖初值與小推力優(yōu)化結(jié)果之間影響關(guān)系的分析。

本文基于小行星探測(cè)任務(wù)，以燃料最優(yōu)為性能指標(biāo)，提出了一種采取粒子群優(yōu)化(PSO)和序列二次規(guī)劃(SQP)的組合優(yōu)化算法，可適應(yīng)于不同形式的脈沖初值輸入，并保證收斂性；設(shè)計(jì)了地球1∶1共振近地小行星2016HO3交會(huì)任務(wù)的小推力轉(zhuǎn)移軌道，給出了軌道設(shè)計(jì)參數(shù)；對(duì)比分析了可用于小推力轉(zhuǎn)移軌道優(yōu)化初值結(jié)果的3種典型脈沖軌道，討論了不同脈沖初值對(duì)小推力轉(zhuǎn)移軌道優(yōu)化結(jié)果的影響。

1 小推力轉(zhuǎn)移軌道優(yōu)化問(wèn)題

在日心黃道慣性坐標(biāo)系下，探測(cè)器典型的二體動(dòng)力學(xué)方程為

(1)

式中：r，v為探測(cè)器的位置矢量和速度矢量；U為發(fā)動(dòng)機(jī)的推力矢量，U=[FxFyFz]；m為探測(cè)器的瞬時(shí)質(zhì)量；Isp為發(fā)動(dòng)機(jī)比沖；f為除發(fā)動(dòng)機(jī)推力加速度之外的攝動(dòng)加速度；g0為海平面重力加速度。r，v，m為狀態(tài)變量，X=[rvm]。將式(1)簡(jiǎn)寫(xiě)為

(2)

將小推力軌道優(yōu)化問(wèn)題轉(zhuǎn)換為尋找最優(yōu)控制函數(shù)U(t)，并使性能指標(biāo)

J=J(tf,Xf)

(3)

達(dá)到最小，同時(shí)滿足動(dòng)力學(xué)方程約束條件、初末狀態(tài)約束條件、路徑約束條件，即

(4)

X(t0)=X0,X(tf)=Xf

(5)

ψ(U(t))≤0

(6)

2 優(yōu)化數(shù)學(xué)模型

小推力軌道優(yōu)化問(wèn)題的最優(yōu)控制描述形式會(huì)引入無(wú)明顯物理意義的協(xié)態(tài)變量，導(dǎo)致收斂域窄，造成初值猜測(cè)困難[7]。本文基于直接法的離散思想，建立非線性規(guī)劃(NLP)形式的小推力轉(zhuǎn)移軌道優(yōu)化數(shù)學(xué)模型。以離散節(jié)點(diǎn)的中點(diǎn)作為強(qiáng)制匹配點(diǎn)，將微分方程約束轉(zhuǎn)化為殘差形式的等式約束。

將整個(gè)飛行時(shí)間分為N個(gè)時(shí)間段，各節(jié)點(diǎn)處的時(shí)刻為

t0

(7)

(8)

則S∈[0,1]。對(duì)狀態(tài)變量進(jìn)行三階Hermite插值，控制變量采用線性插值。以x代表狀態(tài)變量的任一分量，則

x=C0+C1S+C2S2+C3S3

(9)

(10)

(11)

(12)

對(duì)于狀態(tài)變量中任一分量均可執(zhí)行上述操作。至此，連續(xù)的動(dòng)力學(xué)方程約束轉(zhuǎn)化為離散的7N個(gè)等式約束。

在離散形式下，初始狀態(tài)與終端狀態(tài)分別與出發(fā)天體和到達(dá)天體的位置速度相同，即

(13)

給定出發(fā)時(shí)間和到達(dá)時(shí)間，通過(guò)查找星歷獲取出發(fā)、到達(dá)天體的位置速度。末端質(zhì)量自由，初始質(zhì)量為發(fā)射質(zhì)量，假設(shè)發(fā)射質(zhì)量為m0，則

m(t0)=m0

(14)

離散形式的路徑約束為各節(jié)點(diǎn)處的推力值小于等于最大推力值Fmax，即

‖U(ti)‖≤Fmax

(15)

以燃料最優(yōu)為指標(biāo)，將離散形式表示為

(16)

至此，小推力轉(zhuǎn)移軌道問(wèn)題優(yōu)化的數(shù)學(xué)模型完成建立。尋優(yōu)變量Z={t0，tf，X0，…，XN，U0,…,UN}，優(yōu)化指標(biāo)為式(16)，約束條件為式(12)～(15)。其中，約束條件和性能指標(biāo)可根據(jù)具體情況進(jìn)行相應(yīng)修改。

3 基于PSO和SQP的組合優(yōu)化算法

3.1 PSO算法搜索脈沖初值

傳統(tǒng)PSO算法具有全局優(yōu)化和并行計(jì)算的能力，相比于經(jīng)典的基于梯度的優(yōu)化算法具有更好的全局收斂性[9]，但PSO算法出現(xiàn)早熟現(xiàn)象的概率較大，以致收斂于局部最優(yōu)解。為此，本文通過(guò)線性減小慣性權(quán)重系數(shù)和局部學(xué)習(xí)因子，增大全局學(xué)習(xí)因子，使算法在初期具備較強(qiáng)探索能力，在后期又有較好的收斂性，使其在最優(yōu)解附近精細(xì)搜索。本文算法主要包括以下步驟。

1) 初始化粒子群位置p和速度v，粒子位置信息包含脈沖作用時(shí)刻和脈沖作用矢量。

2) 對(duì)每個(gè)粒子用構(gòu)造的位置函數(shù)Fitness進(jìn)行評(píng)價(jià)，此時(shí)Fitness為總速度增量。

3) 更新每個(gè)粒子的歷史最優(yōu)位置pb和群體的全局最優(yōu)位置gb。

4) 分別按照式(17)，(18)更新每個(gè)粒子的速度和位置，其表達(dá)式為

vi+1=w·vi+c1·rand()·(pb-p)+

c2·rand()·(gb-pb)

(17)

pi+1=pi+vi+1

(18)

式中：w為粒子的慣性權(quán)重系數(shù)；c1為局部學(xué)習(xí)因子；c2為全局學(xué)習(xí)因子；rand()表示0～1的隨機(jī)數(shù)。同時(shí)設(shè)置粒子學(xué)習(xí)速度在[-vmax，vmax]內(nèi)，若更新后的速度超出邊界，則取相應(yīng)的邊界值。w，c1，c2的更新依據(jù)為

(19)

(20)

(21)

式中：c為運(yùn)行代數(shù)；cmax為最大運(yùn)行代數(shù)。

5) 判斷終止循環(huán)條件。若條件不滿足，則循環(huán)執(zhí)行步驟2～4；若滿足，則輸出結(jié)果，生成初值Z0={t0，tf，X0，…，XN，U0，…，UN}0。其中，{t0，tf}和{X0,X1,…,XN}中的位置和速度分量可直接從結(jié)果中生成，質(zhì)量可基于其消耗特性按等差遞減生成；{U0,U1,…,UN}0可在[0，F(xiàn)max]內(nèi)隨機(jī)猜測(cè)。為避免隨機(jī)性對(duì)結(jié)果帶來(lái)的影響，Ui0=Fmax[0.5, 0.5, 0.5]。算法相關(guān)參數(shù)設(shè)置見(jiàn)表1。

表1 粒子群算法參數(shù)

3.2 SQP算法求解NLP問(wèn)題

SQP算法是求解NLP問(wèn)題的有效算法[10]，但SQP算法是基于可行方向搜索的一種約束優(yōu)化方法，其只具有局部?jī)?yōu)化能力，因此將SQP與PSO算法結(jié)合可充分發(fā)揮兩者優(yōu)勢(shì)。算法的主要步驟如下。

1) 給定PSO算法搜索后生成的初值Z0，選定正定矩陣H0。

2) 求解的二次規(guī)劃問(wèn)題為

mindT

s.t. con(Zk)+dTcon(Ζk)≤0

(22)

式中：obj為目標(biāo)函數(shù)；con為約束；dk為搜索方向，當(dāng)‖dk‖<10-6時(shí)，迭代結(jié)束。

3) 更新Zk+1=Zk+αkdk，其中αk由線性搜索確定。

4) 修正Hk，使Hk+1保持正定。

5) 判斷是否滿足迭代終止條件：若迭代次數(shù)超過(guò)2 000次或‖dk‖<10-6，則迭代結(jié)束；否則重復(fù)步驟2～4。

為便于計(jì)算，初值生成過(guò)程和優(yōu)化求解過(guò)程的尋優(yōu)變量均進(jìn)行歸一化處理。組合優(yōu)化算法的流程如圖1所示。

圖1 組合優(yōu)化算法流程Fig.1 Combinatorial optimization algorithm flow

4 近地小行星2016HO3交會(huì)任務(wù)

2016HO3是2016年發(fā)現(xiàn)的1顆阿波羅型地球1∶1共振近地小行星，其日心軌道周期與地球公轉(zhuǎn)周期呈現(xiàn)1∶1的關(guān)系。雖然2016HO3為繞日小行星，但是其相對(duì)地球的軌道較穩(wěn)定，因此被稱為地球的第2個(gè)月亮。該小行星具有獨(dú)特的探測(cè)價(jià)值，是我國(guó)小行星探測(cè)任務(wù)的首要備選目標(biāo)，其軌道參數(shù)見(jiàn)表2。表中：AU為天文單位；MJD為約化儒略日。

分析不同初值對(duì)優(yōu)化結(jié)果的影響。假設(shè)探測(cè)器的發(fā)射質(zhì)量為1 500 kg，發(fā)射日期在2022年1月1日至2024年1月1日之間，飛行時(shí)間小于700 d。發(fā)射C3(離開(kāi)地球引力影響球時(shí)與地球相對(duì)速度的平方)不大于20 km2/s2，2016HO3交會(huì)任務(wù)的發(fā)射C3等高線如圖2所示，發(fā)射時(shí)間從2022年1月1日開(kāi)始。由圖可見(jiàn)：若僅從發(fā)射C3的角度考慮，則每年有4個(gè)發(fā)射窗口，平均每季度出現(xiàn)1次波谷，分為長(zhǎng)轉(zhuǎn)移和短轉(zhuǎn)移，兩者交替出現(xiàn)。其中，長(zhǎng)轉(zhuǎn)移時(shí)間約為380～550 d，短轉(zhuǎn)移時(shí)間約為50～200 d。

表2 2016HO3小行星軌道參數(shù)(MJD=58 000.0)

圖2 發(fā)射C3等高線Fig.2 Contour line of launch C3

假設(shè)小推力大小為120 mN，發(fā)動(dòng)機(jī)比沖為3 500 s，推力大小和方向可調(diào)。選取單圈雙脈沖、三脈沖、多圈雙脈沖3種典型的脈沖初值軌道。將轉(zhuǎn)移軌道等時(shí)間間隔離散為150個(gè)軌道段，經(jīng)PSO算法搜索后的初值軌道參數(shù)及經(jīng)SQP算法優(yōu)化后的小推力軌道參數(shù)見(jiàn)表3，由表可得：

1) 3種初值軌道優(yōu)化后得到了2022年11月29日和2022年5月28日這2個(gè)發(fā)射窗口，其中三脈沖初值和多圈雙脈沖的初值總速度增量差距最大，但優(yōu)化結(jié)果最為接近。

2) 優(yōu)化軌道和初值軌道的發(fā)射時(shí)間與飛行時(shí)間具有強(qiáng)耦合關(guān)系，提供的初值飛行時(shí)間越長(zhǎng)，優(yōu)化得到的結(jié)果飛行時(shí)間也越長(zhǎng)。SQP算法只是局部?jī)?yōu)化，因此優(yōu)化軌道與初值軌道的飛行時(shí)間非常接近。

表3 2016HO3小行星交會(huì)任務(wù)軌道優(yōu)化結(jié)果

3) 3種形式的脈沖轉(zhuǎn)移軌道總速度增量差距最大達(dá)58%，而3種小推力轉(zhuǎn)移方案燃料消耗差距不超過(guò)6%。由此可見(jiàn)，本文算法對(duì)初值的敏感度小，具有一定的全局優(yōu)化能力。

各初值對(duì)應(yīng)的小推力轉(zhuǎn)移軌道的控制曲線如圖3所示。圖中，推力方向以黃經(jīng)、黃緯表示。

圖3 不同初值的最優(yōu)控制Fig.3 Optimal control curves of different initial values

由圖3可見(jiàn)：

1) 小推力發(fā)動(dòng)機(jī)基本符合滿推或關(guān)機(jī)的形式，即為Bang-Bang控制的小推力軌道最優(yōu)控制形式[11]。

2) 最優(yōu)小推力轉(zhuǎn)移軌道總體上存在2個(gè)主推進(jìn)段，不同的脈沖軌道初值對(duì)應(yīng)的最優(yōu)控制曲線的差異主要體現(xiàn)在推力器開(kāi)關(guān)機(jī)時(shí)刻的不同上。

3) 對(duì)于不同初值對(duì)應(yīng)的優(yōu)化結(jié)果，其推力方向的變化在相位上有明顯差別，可看出2個(gè)發(fā)射窗口對(duì)應(yīng)2種推力變化形式。

4) 從飛行時(shí)間來(lái)看，多圈雙脈沖初值優(yōu)化得到的小推力轉(zhuǎn)移方案在第2個(gè)推進(jìn)段關(guān)機(jī)后有約10 d的自由飛行段。由此可見(jiàn)，三脈沖初值具有更優(yōu)的性能。

圖4 燃料最優(yōu)小推力轉(zhuǎn)移軌道根數(shù)變化Fig.4 Time histories of fuel-optimal low-thrust transfer orbital elements

單圈雙脈沖初值和三脈沖初值對(duì)應(yīng)的小推力轉(zhuǎn)移軌道根數(shù)變化過(guò)程如圖4所示。由圖可見(jiàn)：對(duì)應(yīng)于日期為2022年11月29日的發(fā)射窗口，發(fā)動(dòng)機(jī)的工作序列為“先長(zhǎng)后短”；對(duì)應(yīng)于日期為2022年5月28日的發(fā)射窗口，發(fā)動(dòng)機(jī)的工作序列為“先短后長(zhǎng)”。11月的發(fā)射窗口與5月的發(fā)射窗口相比，當(dāng)探測(cè)器從地球出發(fā)時(shí)，前者為小偏心率大軌道傾角，后者為大偏心率小軌道傾角。從半長(zhǎng)軸變化趨勢(shì)來(lái)看，前者為先減小后增大，后者為先增大后減小。

5 結(jié)束語(yǔ)

為解決小推力軌道優(yōu)化未考慮脈沖初值多樣性的問(wèn)題，利用不同形式的脈沖軌道初值對(duì)小推力轉(zhuǎn)移軌道優(yōu)化結(jié)果的影響進(jìn)行了研究，提出了能適應(yīng)不同脈沖初值的組合優(yōu)化算法，以我國(guó)深空探測(cè)規(guī)劃背景任務(wù)中的小行星探測(cè)任務(wù)為例進(jìn)行了仿真分析。針對(duì)近地小行星2016HO3交會(huì)任務(wù)，以3種典型的脈沖軌道優(yōu)化得到了3種具體的小推力轉(zhuǎn)移軌道方案和2個(gè)發(fā)射窗口。其中，三脈沖軌道初值與多圈雙脈沖軌道初值對(duì)應(yīng)的優(yōu)化結(jié)果主要體現(xiàn)在開(kāi)關(guān)機(jī)時(shí)刻的不同，單圈雙脈沖初值與三脈沖初值、多圈雙脈沖初值的優(yōu)化結(jié)果主要體現(xiàn)在推力方向變化的差異，3種初值軌道對(duì)應(yīng)最優(yōu)小推力軌道燃料消耗量差距不超過(guò)6%。研究結(jié)果表明：本文算法對(duì)初值的敏感度低，能收斂于Bang-Bang最優(yōu)控制形式；不同初值對(duì)小推力軌道的整體性能指標(biāo)影響較小，但開(kāi)關(guān)機(jī)時(shí)刻和推力方向的變化會(huì)產(chǎn)生較大差異，從而得到不同的最優(yōu)控制曲線。該研究成果對(duì)小推力軌道優(yōu)化初值的選取、我國(guó)未來(lái)小行星探測(cè)任務(wù)的工程實(shí)施具有一定的參考價(jià)值。由于本文未對(duì)2個(gè)窗口的控制策略進(jìn)行優(yōu)劣性分析，且優(yōu)化模型僅以燃料消耗作為評(píng)價(jià)指標(biāo)，未考慮軌道設(shè)計(jì)中飛行時(shí)間、控制誤差、多體模型等其他重要參考因素，因此后續(xù)將綜合考慮多種評(píng)價(jià)指標(biāo)，進(jìn)一步優(yōu)化設(shè)計(jì)。