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

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

0 引言

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

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

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

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

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

(1)

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

(2)

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

J=J(tf,Xf)

(3)

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

(4)

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

(5)

ψ(U(t))≤0

(6)

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

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

將整個飛行時間分為N個時間段，各節(jié)點處的時刻為

t0

(7)

(8)

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

x=C0+C1S+C2S2+C3S3

(9)

(10)

(11)

(12)

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

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

(13)

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

m(t0)=m0

(14)

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

‖U(ti)‖≤Fmax

(15)

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

(16)

至此，小推力轉(zhuǎ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īng)典的基于梯度的優(yōu)化算法具有更好的全局收斂性[9]，但PSO算法出現(xiàn)早熟現(xiàn)象的概率較大，以致收斂于局部最優(yōu)解。為此，本文通過線性減小慣性權(quán)重系數(shù)和局部學(xué)習(xí)因子，增大全局學(xué)習(xí)因子，使算法在初期具備較強探索能力，在后期又有較好的收斂性，使其在最優(yōu)解附近精細(xì)搜索。本文算法主要包括以下步驟。

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

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

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

4) 分別按照式(17)，(18)更新每個粒子的速度和位置，其表達(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的隨機數(shù)。同時設(shè)置粒子學(xué)習(xí)速度在[-vmax，vmax]內(nèi)，若更新后的速度超出邊界，則取相應(yīng)的邊界值。w，c1，c2的更新依據(jù)為

(19)

(20)

(21)

式中：c為運行代數(shù)；cmax為最大運行代數(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)隨機猜測。為避免隨機性對結(jié)果帶來的影響，Ui0=Fmax[0.5, 0.5, 0.5]。算法相關(guān)參數(shù)設(shè)置見表1。

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

3.2 SQP算法求解NLP問題

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

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

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

mindT

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

(22)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

由圖3可見：

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

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

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

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

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

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

5 結(jié)束語

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