武 迪，程 林，王 偉，李俊峰

（1. 清華大學(xué) 航天航空學(xué)院，北京 100084；2. 北京航空航天大學(xué) 宇航學(xué)院，北京 100191）

引 言

以電推進(jìn)為代表的小推力推進(jìn)系統(tǒng)近年來受到了廣泛的關(guān)注，相較于傳統(tǒng)化學(xué)推進(jìn)，電推進(jìn)具有更高的推進(jìn)效率，即更能節(jié)省燃料[1]。自“深空1號”（Deep Space 1）探測器成功驗(yàn)證電推進(jìn)技術(shù)后，“黎明號”（Dawn）、“隼鳥號”（Hayabusa）和正在進(jìn)行的“貝皮科倫布”探測器（BepiColombo）等深空任務(wù)均采用了電推進(jìn)發(fā)動機(jī)。而在近地任務(wù)中，波音公司的Boeing 902SP和中國的“實(shí)踐十七”“實(shí)踐二十”將電推進(jìn)技術(shù)應(yīng)用于地球同步軌道的入軌和軌道保持任務(wù)[2]。因此，電推進(jìn)技術(shù)由于其出色的推進(jìn)效率，在未來航天任務(wù)設(shè)計中將逐漸獲得更普遍的應(yīng)用。

在軌道設(shè)計方面，由于電推進(jìn)所產(chǎn)生的推力較小，需要長時間的開機(jī)工作才能產(chǎn)生足夠的速度增量，以電推進(jìn)驅(qū)動的航天器小推力最優(yōu)軌跡將和脈沖推進(jìn)的最優(yōu)軌跡產(chǎn)生顯著的差別。圓錐曲線拼接[1]、基于形狀的標(biāo)稱軌道方法[3]等相對解析的設(shè)計方法難以實(shí)現(xiàn)對小推力最優(yōu)軌跡的準(zhǔn)確近似。而對于一般的深空探測軌道設(shè)計，在連續(xù)推力作用下的軌道無法采用類似解析的方法得到最優(yōu)解，而是采用較為耗時的數(shù)值優(yōu)化的方法（主要為直接法和間接法）[4]。

對于探測任務(wù)尤其是多目標(biāo)任務(wù)的全局規(guī)劃問題[5-6]，需要快速且準(zhǔn)確地確定小推力軌跡的燃料消耗、轉(zhuǎn)移時間等參數(shù)。因此，實(shí)現(xiàn)對小推力軌跡的快速優(yōu)化具有重要意義。

目前，主要的數(shù)值優(yōu)化方法分為直接法和間接法[7]。

直接法將小推力優(yōu)化問題進(jìn)行離散，轉(zhuǎn)化為非線性規(guī)劃問題求解，通常需要較多的離散點(diǎn)以提高解的精度和最優(yōu)性[8]。相較于直接法，間接法利用極大值原理構(gòu)建兩點(diǎn)邊值問題求解，能夠保證最優(yōu)性的一階必要條件[9]。但是這些數(shù)值方法均需要提供初值進(jìn)行優(yōu)化，直接法需要提供離散點(diǎn)的狀態(tài)和控制初值，比較容易利用設(shè)計經(jīng)驗(yàn)或標(biāo)稱軌道方法得到，而且直接法的收斂域相對較大，其求解對初值敏感性較弱[8]，因此其計算性能主要取決于離散點(diǎn)數(shù)。間接法需要提供協(xié)態(tài)的初值，相較于具有實(shí)際物理意義較容易猜測的狀態(tài)，協(xié)態(tài)并無明確的物理意義，其取值范圍也難以確定，因此較難提供初值[1]。同倫方法[9]將小推力燃料最優(yōu)問題通過同倫參數(shù)聯(lián)系轉(zhuǎn)換為求解相對較容易的能量最優(yōu)問題求解，提高了初值猜測成功的概率。協(xié)態(tài)歸一化方法[10]通過引入歸一化參數(shù)將協(xié)態(tài)的猜測范圍限制在高維單位球面上，提供了初值猜測的范圍，進(jìn)一步提高了初值猜測的效率?；跇?biāo)稱軌道的協(xié)態(tài)初值估計方法[1]能夠提供能量最優(yōu)問題的初值猜測，但是其猜測的成功率依賴于標(biāo)稱軌道的設(shè)計，對于較復(fù)雜的任務(wù)，標(biāo)稱軌道與最優(yōu)軌跡偏離較大，協(xié)態(tài)估計的效果較差。因此針對小推力間接法軌跡優(yōu)化，解析的保證收斂的協(xié)態(tài)初值能夠顯著降低求解失敗的概率，以提高小推力軌跡優(yōu)化的效率。

本研究提出了一種解析的協(xié)態(tài)初始化方法，引入切換系統(tǒng)，將小推力燃料最優(yōu)問題同一個具有解析初值的軌跡優(yōu)化問題相聯(lián)系，通過迭代方法求解得到小推力燃料最優(yōu)軌跡。本文首先將小推力優(yōu)化模型嵌入到更為廣義的切換系統(tǒng)模型中，給出了切換系統(tǒng)的構(gòu)建形式，并設(shè)計了相應(yīng)的切換函數(shù)，從而將小推力優(yōu)化問題轉(zhuǎn)化為更簡單優(yōu)化模型的優(yōu)化問題；其次，給出了具有解析協(xié)態(tài)初值的簡單系統(tǒng)優(yōu)化模型，能夠通過切換系統(tǒng)與小推力優(yōu)化模型相聯(lián)系，進(jìn)一步給出了計算迭代算法，并通過數(shù)值仿真驗(yàn)證了該方法相較于傳統(tǒng)同倫方法，具有更高的求解效率。

1 切換系統(tǒng)及切換函數(shù)

針對常見的小推力軌跡優(yōu)化問題，本節(jié)首先給出了一般指標(biāo)形式下的小推力軌跡優(yōu)化模型，其中動力學(xué)模型也采用了一般的矩陣形式描述。在此優(yōu)化模型的基礎(chǔ)上，進(jìn)一步介紹了切換系統(tǒng)的概念，從而將小推力優(yōu)化模型嵌入進(jìn)更廣義的切換系統(tǒng)優(yōu)化模型中，同時給出了切換函數(shù)設(shè)計。

1.1 小推力軌跡優(yōu)化模型

一般而言，考慮如下小推力軌跡優(yōu)化模型，優(yōu)化指標(biāo)為

其中：指標(biāo)J1通常取為總的燃料消耗[9]、能量消耗[1]、飛行時間[8]或混合指標(biāo)；表示狀態(tài)矢量，通常由航天器的位置、速度（或軌道要素）和質(zhì)量組成，即表示控制量，通常包括控制量幅值及其方向。動力學(xué)方程為

此優(yōu)化模型可以轉(zhuǎn)化為兩點(diǎn)邊值問題并采用間接法求解。相比于直接法，間接法能夠保證最優(yōu)的一階必要條件，因此一旦求解成功即可得到局部最優(yōu)的軌道。但另一方面，在間接法的求解過程中需要人為提供無物理意義的協(xié)態(tài)初值猜測作為求解的初始化，而且間接法的收斂域也較直接法更小，難以對估計不準(zhǔn)的協(xié)態(tài)初值實(shí)現(xiàn)成功求解。因此，優(yōu)化時的初值猜測情況會顯著影響間接法的求解成功率和效率，當(dāng)前主要采用隨機(jī)猜測的方式多次猜測求解，猜測成功率較低且求解效率較慢。

基于協(xié)態(tài)初值估計的方法，本文將小推力優(yōu)化模型嵌入進(jìn)切換系統(tǒng)中，同時設(shè)計切換系統(tǒng)形式使其具有解析的協(xié)態(tài)初值，避免了初值猜測困難。以此解析協(xié)態(tài)初值作為初始化，通過切換函數(shù)迭代實(shí)現(xiàn)了小推力最優(yōu)軌跡的快速求解。

1.2 切換系統(tǒng)設(shè)計

基于小推力軌跡優(yōu)化模型，引入如下一般形式的切換系統(tǒng)[11]優(yōu)化模型。優(yōu)化指標(biāo)為

動力學(xué)方程為

一般而言，傳統(tǒng)的從能量最優(yōu)問題同倫求解燃料最優(yōu)問題的算法，將同倫參數(shù)添加到優(yōu)化指標(biāo)中，而動力學(xué)方程保持不變。此時，首先求解能量最優(yōu)問題的協(xié)態(tài)初值，然后逐漸改變同倫參數(shù)值，以得到的協(xié)態(tài)初值作為猜測值，求解對應(yīng)同倫參數(shù)下的優(yōu)化問題，最終得到燃料最優(yōu)解。

與傳統(tǒng)方法不同，本文的優(yōu)化模型將同倫參數(shù)引入到兩個系統(tǒng)之間的切換函數(shù)中，系統(tǒng)1的優(yōu)化指標(biāo)為動力學(xué)方程的右函數(shù)為。系統(tǒng)2的優(yōu)化指標(biāo)為動力學(xué)方程的右函數(shù)為采用同樣的同倫過程即可得到燃料最優(yōu)解。在本文中，系統(tǒng)1設(shè)計為小推力燃料最優(yōu)軌跡優(yōu)化模型，系統(tǒng)2及其解析協(xié)態(tài)初值的設(shè)計將在下一節(jié)中詳細(xì)給出，切換函數(shù)設(shè)計為

根據(jù)極大值原理，系統(tǒng)的最優(yōu)控制使哈密頓函數(shù)

2 燃料最優(yōu)軌跡優(yōu)化

基于標(biāo)稱軌道線性化的方法，本節(jié)首先給出了系統(tǒng)2的模型，并進(jìn)一步給出了確定小推力燃料最優(yōu)軌跡的求解算法。

2.1 解析協(xié)態(tài)初值

系統(tǒng)2考慮如下線性優(yōu)化模型。優(yōu)化指標(biāo)為

另外，對于標(biāo)稱軌道的選取，在本文中僅考慮的簡單線性的標(biāo)稱軌道，具體形式為

由于系統(tǒng)2設(shè)計即為簡單線性的優(yōu)化模型，此種標(biāo)稱軌道能夠進(jìn)一步簡化解析解求解的難度，快速給出解析協(xié)態(tài)初值。本文的數(shù)值仿真證明此標(biāo)稱軌道能夠處理較為復(fù)雜的燃料最優(yōu)問題。

根據(jù)極大值原理可得最優(yōu)控制律為

因此，得到了一個線性時變的動力學(xué)微分方程，此方程中系數(shù)僅和標(biāo)稱軌道以及當(dāng)前時刻有關(guān)，而和狀態(tài)、協(xié)態(tài)無關(guān)。此方程的解可以解析得到

對于深空探測中的交會問題，通常考慮如下末端條件

由于此系統(tǒng)中質(zhì)量為常數(shù)，質(zhì)量協(xié)態(tài)的初值并不影響此系統(tǒng)的解，因此在計算中將質(zhì)量協(xié)態(tài)的初值設(shè)為0，即將式（16）代入式（15），其它協(xié)態(tài)初值的解析表達(dá)式為

2.2 燃料最優(yōu)軌跡優(yōu)化算法

對于燃料最優(yōu)問題，系統(tǒng)1的指標(biāo)取為燃料消耗，系統(tǒng)2的指標(biāo)中取值和系統(tǒng)1保持一致，即

將式（18）代入切換系統(tǒng)的表達(dá)式，可得切換系統(tǒng)的優(yōu)化指標(biāo)為

3）以上一步收斂的協(xié)態(tài)初值作為打靶量的初值進(jìn)行打靶；

4）若打靶收斂，則轉(zhuǎn)到步驟2繼續(xù)求解過程，若打靶不收斂，則以更小的步長重新轉(zhuǎn)到步驟3計算；

5）當(dāng)同倫參數(shù)為0時，輸出收斂的小推力燃料最優(yōu)軌跡的結(jié)果。

3 數(shù)值仿真

本節(jié)仿真試驗(yàn)了自地球出發(fā)前往不同目標(biāo)[1]的深空探測軌跡，以驗(yàn)證方法的有效性。在仿真中，采用了ODE45積分器，其相對誤差和絕對誤差設(shè)置為打靶算法采用Minpack-1，其參數(shù)設(shè)置參考相關(guān)文獻(xiàn)[10]。為加快積分并降低數(shù)值敏感性，歸一化參數(shù)設(shè)置為長度單位1 AU（天文單位），時間單位1年，質(zhì)量單位歸一化為航天器初始質(zhì)量。此時，太陽的引力常數(shù)為所有算例均以C++語言編寫（Microsoft Visual Studio 2012編譯），并在臺式機(jī)使用單核運(yùn)行（Intel Core i7-7 700 CPU，主頻3.6 GHz，內(nèi)存8.0 GB）。

算例任務(wù)參數(shù)如表1所示，由于本文采用春分點(diǎn)軌道根數(shù)描述，因此各個算例的航天器軌道圈數(shù)同時給定。相應(yīng)時刻的地球和目標(biāo)的位置速度可根據(jù)JPL系統(tǒng)查詢得到。

表1 數(shù)值仿真任務(wù)參數(shù)表Table 1 Parameters for the numerical simulation missions

3.1 計算效率對比

本文方法的計算效率和傳統(tǒng)能量到燃料最優(yōu)同倫方法的計算效率進(jìn)行了對比。這兩種方法均應(yīng)用了前文所述的迭代過程。由于能量最優(yōu)求解需要進(jìn)行隨機(jī)猜測初值，因此其每個算例進(jìn)行了1 000次隨機(jī)猜測并記錄收斂個數(shù)。而本文方法采用線性的標(biāo)稱軌道，且由于具有解析初值僅進(jìn)行單次計算即可。相應(yīng)的計算收斂率、所需的迭代步數(shù)和得到單個解的平均計算時間如表2所示。本文方法對3個算例均能夠?qū)崿F(xiàn)100%收斂，而傳統(tǒng)同倫方法的收斂率分別為86.7%、83.6%和39.3%，尤其是對第3個算例，傳統(tǒng)方法收斂率較低，計算效率差別較大。

表2 該方法與傳統(tǒng)同倫方法計算效率對比Table 2 Table1 Comparisons of proposed method with the traditional homotopic method

3.2 燃料最優(yōu)解

對每個任務(wù)，迭代過程中的控制曲線如圖1所示。對于地球至金星的任務(wù)，飛行時間為700 d，燃料最優(yōu)控制包含6個關(guān)機(jī)段和5個開機(jī)段，總的開機(jī)時長為275.1 d，消耗燃料210.45 kg。對于地球至“坦普爾1號”（Tempel 1）任務(wù)，燃料最優(yōu)控制包含2個關(guān)機(jī)段和3個開機(jī)段，最短的關(guān)機(jī)段為初始的弧段，時長0.63 d。任務(wù)中總開機(jī)時間為197.7 d，消耗燃料348.36 kg。地球至小行星3 671（Dionysus）的任務(wù)總時長為3 534 d，燃料最優(yōu)控制包含7個關(guān)機(jī)段和6個開機(jī)段，總計開機(jī)時間1 361.5 d，總計燃料消耗1 279.93 kg。3個任務(wù)場景的燃料最優(yōu)軌跡如圖2所示，地球至金星任務(wù)為2圈近圓軌跡，軌道傾角改變較小，求解相對較為容易，因此本文方法和傳統(tǒng)同倫方法均具有較高的求解效率。Tempel 1的軌道為橢圓軌道，因此，此任務(wù)的控制曲線較為復(fù)雜，而計算所需的迭代步數(shù)較多。Dionysus的軌道為傾角較大的橢圓軌道，本文方法求解此問題所需時間僅為傳統(tǒng)方法的1/5。

圖1 各任務(wù)迭代過程的推力曲線圖Fig. 1 History of the thrust magnitude of different missions

圖2 各任務(wù)場景的燃料最優(yōu)軌跡圖Fig. 2 Fuel-optimal trajectories of different missions

4 結(jié) 論

該研究提出了基于切換系統(tǒng)的解析初始化方法，將小推力燃料最優(yōu)軌跡優(yōu)化問題同具有解析協(xié)態(tài)初值的線性優(yōu)化問題聯(lián)系起來，能夠解析初始化，并提高間接法小推力燃料最優(yōu)軌跡的求解效率。仿真表明：本方法和從能量最優(yōu)同倫至燃料最優(yōu)軌跡的方法所需的迭代步數(shù)基本相同，而其協(xié)態(tài)初值能夠解析得到，因而具有更高的求解效率。