基于間接法的最優(yōu)小推力逃逸軌道設(shè)計(jì)

黃鎬 韓潮

（北京航空航天大學(xué)宇航學(xué)院，北京 100191）

1 引言

在小推力星際探測(cè)任務(wù)中，逃逸捕獲階段因其多圈螺旋軌道的特性而在整個(gè)任務(wù)過(guò)程中最為復(fù)雜。目前，基于小推力的軌道優(yōu)化設(shè)計(jì)方法可歸結(jié)為直接法和間接法。直接法采用離散配點(diǎn)法將最優(yōu)控制問(wèn)題轉(zhuǎn)化成高維非線性參數(shù)規(guī)劃問(wèn)題，魯棒性較好，對(duì)于單圈或圈數(shù)較少的軌道轉(zhuǎn)移問(wèn)題比較適用；但對(duì)于多圈（大于50圈）螺旋逃逸軌道設(shè)計(jì)，采用直接法將形成龐大的非線性參數(shù)規(guī)劃問(wèn)題，導(dǎo)致求解難度非常大［1］。間接法利用龐特里亞金極值原理將最優(yōu)控制問(wèn)題轉(zhuǎn)化成兩點(diǎn)邊值問(wèn)題，未知協(xié)狀態(tài)量維數(shù)較少，而且求解的精度高，但兩點(diǎn)邊值問(wèn)題對(duì)協(xié)狀態(tài)量初值非常敏感，對(duì)于多圈長(zhǎng)時(shí)間逃逸問(wèn)題，該問(wèn)題尤甚，要得到收斂解需要提供非常好的初值猜測(cè)，然而協(xié)狀態(tài)變量沒(méi)有實(shí)際物理意義，初值很難猜測(cè)，甚至其量級(jí)都無(wú)法界定，由此給間接法的應(yīng)用帶來(lái)了不小困難。目前，國(guó)內(nèi)外已有一些利用間接法進(jìn)行逃逸螺旋軌跡優(yōu)化的研究［2－5］，也有一些研究在求解整個(gè)星際探測(cè)任務(wù)中考慮了逃逸階段［5－6］。這些研究中，普遍采用了一種協(xié)狀態(tài)控制轉(zhuǎn)換技術(shù)（ACT）［3］來(lái)降低兩點(diǎn)邊值問(wèn)題初值的敏感性，且大部分都只研究了控制無(wú)約束的平面內(nèi)逃逸問(wèn)題，或者采用逃逸時(shí)間［2］或逃逸能量［5］作為優(yōu)化性能指標(biāo)。

本文采用間接法［7］求解小推力燃料最優(yōu)逃逸問(wèn)題，通過(guò)引入一個(gè)待定參數(shù)，求解了復(fù)雜的長(zhǎng)時(shí)間控制有約束的三維燃料最優(yōu)逃逸問(wèn)題。

2 地心球慣性坐標(biāo)系下軌道動(dòng)力學(xué)方程

在地心慣性系中的位置、速度矢量用r和v表示；ex、ey、ez分別表示地心直角慣性系中x，y，z方向上的單位矢量；er、eθ、eφ分別為地心球慣坐標(biāo)系各正交軸單位矢量，如圖1所示。

地心球慣性坐標(biāo)系下航天器軌道動(dòng)力學(xué)系統(tǒng)狀態(tài)量為x＝ ［r，θ，φ，vr，vθ，vφ］T，動(dòng)力學(xué)方程為

圖1 直角坐標(biāo)系與球坐標(biāo)系之間關(guān)系Fig.1 Relation between cartesian and spherical system

式中 r和a 分別為位置和推力加速度幅值；μ為地球引力常數(shù)；u為推力方向矢量，可由角α和角β 來(lái)描述。

對(duì)于變比沖小推力發(fā)動(dòng)機(jī)模型，航天器加速度模型可表示為

式中 F和a分別表示航天器推力和推力加速度矢量；ε為發(fā)動(dòng)機(jī)效率；P為固定的發(fā)動(dòng)機(jī)功率；gn為標(biāo)稱重力加速度；m0和Ispmax分別表示航天器初始質(zhì)量和發(fā)動(dòng)機(jī)最大比沖；m 和Isp則為對(duì)應(yīng)的歸一化處理后的系數(shù)，m ∈［（m0－mfuel）／m0，1］，Isp∈［Ispmin／Ispmax，1］；mfuel為航天器燃料總質(zhì)量；Ispmin為變比沖發(fā)動(dòng)機(jī)能提供的最小比沖。質(zhì)量導(dǎo)數(shù)方程為

3 燃料最優(yōu)控制問(wèn)題

對(duì)于變比沖發(fā)動(dòng)機(jī)燃料最優(yōu)逃逸問(wèn)題，可描述為在初始軌道給定的情況下，在固定時(shí)間通過(guò)控制發(fā)動(dòng)機(jī)比沖Isp的大小以及推力矢量方向，使得航天器軌道能量在達(dá)到逃逸能量的同時(shí)，所消耗的燃料最少。對(duì)應(yīng)的性能指標(biāo)和初始條件的數(shù)學(xué)描述為

式中 t0和tf分別為逃逸始末時(shí)刻；λ0為引入性能指標(biāo)函數(shù)中大于零的待定參數(shù)，若λ0＝1即為原問(wèn)題；若λ0≠1，只要λ0＞0，就不會(huì)改變系統(tǒng)燃料最優(yōu)的本質(zhì)。下標(biāo)0、f分別表示初始時(shí)刻和終點(diǎn)時(shí)刻對(duì)應(yīng)的參量。

由龐特里亞金極小值原理［8］，推導(dǎo)系統(tǒng)Hamilton函數(shù)及最優(yōu)控制為

式中 λ＝［λr，λθ，λφ，λvr，λvθ，λvφ］T為與狀態(tài)量x相對(duì)應(yīng)的協(xié)狀態(tài)量向量；λv＝［λvr，λvθ，λvφ］T為與速度v對(duì)應(yīng)的協(xié)狀態(tài)量向量；λm為與m 對(duì)應(yīng)的協(xié)狀態(tài)量。

由?H／?u＝0，推導(dǎo)系統(tǒng)最優(yōu)控制方向矢量為

將公式（3）帶入Hamilton函數(shù)得：

進(jìn)一步推導(dǎo)最優(yōu)Isp控制為

系統(tǒng)協(xié)狀態(tài)方程為

燃料最優(yōu)逃逸問(wèn)題的終端約束為終點(diǎn)時(shí)刻tf時(shí)的軌道能量εf達(dá)到拋物線逃逸軌道能量值εtar。

拋物線逃逸軌道的能量εtar≡0。由文獻(xiàn) ［8］得到系統(tǒng)增廣約束為

式中 υ為引入的增廣系數(shù)。對(duì)應(yīng)的系統(tǒng)橫截條件為

對(duì)于燃料最優(yōu)逃逸問(wèn)題，由公式（5）、（7）可知˙λθ＝0，且λθf(wàn)＝0，所以λθ≡0。公式（1）、（2）、（5）、（6）、（7）、（8）共同構(gòu)成兩點(diǎn)邊值問(wèn)題。若λ0＝1，即為原問(wèn)題，則由系統(tǒng)Hamilton函數(shù)推導(dǎo)出的兩點(diǎn)邊值問(wèn)題待定參數(shù)向量為［λr0，λφ0，λvr0，λvθ0，λvφ0，λm0，υ］，由于這些變量的取值范圍沒(méi)有任何先驗(yàn)的初始判斷，所以對(duì)于其初值的猜測(cè)難度非常大；若在性能指標(biāo)函數(shù)中引入的正數(shù)λ0為待定參量，則系統(tǒng)待定參數(shù)向量為［λ0，λr0，λφ0，λvr0，λvθ0，λvφ0，λm0，υ］，雖然比原問(wèn)題維數(shù)增加了一維，但卻在不改變?cè)瓋?yōu)化問(wèn)題的前提下使得系統(tǒng)Hamilton函數(shù)中各加減項(xiàng)同質(zhì)［7］，若將系統(tǒng)Hamilton函數(shù)除以大于零的比例系數(shù)K，將在不改變Hamilton函數(shù)的極值屬性前提下得到新系統(tǒng)Hamilton函數(shù)，定義K為待確定參數(shù)向量的模值，即

則由新系統(tǒng)Hamilton函數(shù)推導(dǎo)的兩點(diǎn)邊值問(wèn)題的解向量便為原系統(tǒng)待確定參數(shù)向量除以比例系數(shù)K，即為原系統(tǒng)待確定參數(shù)向量的單位化：

由此可知新系統(tǒng)Hamilton函數(shù)推導(dǎo)的兩點(diǎn)邊值問(wèn)題的待定參數(shù)向量將被約束在一個(gè)多維的單位球里。相對(duì)于原系統(tǒng)Hamilton函數(shù)推導(dǎo)的兩點(diǎn)邊值問(wèn)題，待定參數(shù)向量取值范圍無(wú)任何先驗(yàn)判斷，而由引入未知變量λ0和比例系數(shù)K 構(gòu)成的新系統(tǒng)Hamilton函數(shù)推導(dǎo)的兩點(diǎn)邊值問(wèn)題，每個(gè)待確定參數(shù)的初值猜測(cè)范圍將被約束在區(qū)間 ［－1，1］中，由此將大大提高兩點(diǎn)邊值問(wèn)題的初值猜測(cè)效率，進(jìn)而提高問(wèn)題的收斂性。在求解出由新系統(tǒng)Hamilton函數(shù)推導(dǎo)的兩點(diǎn)邊值問(wèn)題待定參數(shù)向量之后，通過(guò)除以λ＊0，即，將得到原問(wèn)題系統(tǒng)Hamilton推導(dǎo)的兩點(diǎn)邊值問(wèn)題的解向量。

對(duì)于由新系統(tǒng)Hamilton函數(shù)推導(dǎo)的兩點(diǎn)邊值問(wèn)題共有8個(gè)待確定變量，對(duì)應(yīng)的約束條件也為8個(gè)：

在通過(guò)打靶法求解上述待定參數(shù)后，對(duì)公式（1）、（2）、（5）、（6）積分可求得各時(shí)刻軌道狀態(tài)量和協(xié)狀態(tài)，由公式（3）和公式（4）可求得發(fā)動(dòng)機(jī)推力矢量和最優(yōu)Isp控制。

4 仿真算例

設(shè)地球逃逸的初始停泊軌道為半徑為7 500km 的赤道平面內(nèi)圓軌道，由于設(shè)計(jì)參數(shù)與θ無(wú)關(guān)，可以假設(shè)初始軌道的真近點(diǎn)角為0°，航天器總質(zhì)量為100 000kg，發(fā)動(dòng)機(jī)固定功率為5MW，效率為1。

4.1 平面內(nèi)控制有約束長(zhǎng)時(shí)間燃料最優(yōu)地球逃逸

直接采用間接法求解平面內(nèi)控制有約束長(zhǎng)時(shí)間燃料最優(yōu)地球逃逸非常困難，為此，將問(wèn)題分解為兩步：首先由短時(shí)間控制無(wú)約束最優(yōu)逃逸逐步求解長(zhǎng)時(shí)間控制無(wú)約束最優(yōu)逃逸；其次由長(zhǎng)時(shí)間控制無(wú)約束最優(yōu)逃逸逐步求解長(zhǎng)時(shí)間控制有約束最優(yōu)逃逸。

以30天Isp∈［3 600m／s，6 600m／s］區(qū)間的燃料最優(yōu)逃逸問(wèn)題為例，先以0.05天為間隔，分別求解出0.1～5天的燃料最優(yōu)逃逸的收斂解；采用數(shù)值曲線擬合技術(shù)對(duì)0.1～5天的收斂解進(jìn)行擬合，用于外推猜測(cè)5～10 天的未知參量初值，進(jìn)而求解5～10 天燃料最優(yōu)逃逸的收斂解；有了10天的收斂解，擬合形成10～30天的初值猜測(cè)；當(dāng)采用該擬合曲線給出的初值不能收斂的時(shí)候，則利用該發(fā)散時(shí)刻以前的所有已收斂結(jié)果，重新擬合出新的曲線，重新計(jì)算該點(diǎn)初值，然后重新打靶計(jì)算。在仿真試驗(yàn)過(guò)程中，通過(guò)求解最初1～5天燃料最優(yōu)逃逸問(wèn)題發(fā)現(xiàn)，協(xié)狀態(tài)量λvr0變化較大，規(guī)律不明顯，曲線擬合難度較大。為此，借鑒文獻(xiàn)［3］的方法，令在協(xié)狀態(tài)量λvr0隨逃逸時(shí)間tf變化規(guī)律性出現(xiàn)之前，采用擬合b0和λvθ0，然后計(jì)算λvr0的初值猜測(cè)。

根據(jù)0.1～30天的仿真試驗(yàn)數(shù)據(jù)，結(jié)合曲線擬合技術(shù)，得出在采用如下曲線擬合公式時(shí)，能最好地?cái)M合出各待定參數(shù)隨逃逸時(shí)間tf的變化關(guān)系

式中 p1，p2，q1，q2為相應(yīng)的擬合系數(shù)。各待定參數(shù)的擬合系數(shù)如表1所示，各待定參數(shù)對(duì)應(yīng)的擬合曲線與其收斂解之間的關(guān)系如圖2、圖3所示。圖2、3中“◇”為各待定參數(shù)的收斂解，實(shí)線為各待定參數(shù)的擬合曲線。可以看出，采用曲線擬合公式（9）和表1中的參數(shù)，對(duì)各待定參數(shù)擬合效果非常好。

圖2 λr0、λm0和v最優(yōu)解與逃逸時(shí)間關(guān)系Fig.2 Optimal solution and escape time of λr0，λm0and v

圖3 b0、λvr0和λvθ0最優(yōu)解與逃逸時(shí)間關(guān)系Fig.3 Optimal solution and escape time of b0，λvr0andλvθ0

表1 平面內(nèi)Isp無(wú)約束地球逃逸初值猜測(cè)擬合曲線系數(shù)（a＝7 500km，e＝0）Tab.1 Estimation coefficients for earth escape without Ispconstraints（a＝7 500km，e＝0）

采用擬合公式（9）和表1中的系數(shù)，擬合得到tf為30天燃料最優(yōu)逃逸問(wèn)題各待定參數(shù)的初值猜測(cè)，迭代打靶求解平面內(nèi)30天Isp無(wú)約束燃料最優(yōu)逃逸問(wèn)題。在得到平面內(nèi)30天Isp無(wú)約束燃料最優(yōu)逃逸問(wèn)題收斂解的基礎(chǔ)上，以此為初值，利用同倫思想逐步將Isp控制約束從［3 000m／s，10 000 m／s］區(qū)間調(diào)整為實(shí)際的［3 600m／s，5 400m／s］區(qū)間，最終求得Isp屬于［3 600m／s，5 400m／s］區(qū)間的有約束平面內(nèi)逃逸問(wèn)題的收斂解。結(jié)果表面，航天器在30天中繞地球飛行118圈后達(dá)到地球逃逸速度，整個(gè)過(guò)程共消耗燃料13 634kg，占總質(zhì)量的13.6%。Isp無(wú)約束、Isp∈［3 600m／s，6 600 m／s］和Isp∈［3 600m／s，5 400m／s］的30天燃料最優(yōu)平面內(nèi)逃逸軌道如圖4所示，相應(yīng)的Isp隨轉(zhuǎn)移時(shí)間變化的曲線如圖5所示。由圖4可以看出，Isp控制從無(wú)約束到最終約束在［3 600m／s，5 400 m／s］區(qū)間的過(guò)程中，逃逸圈數(shù)沒(méi)有變，依然為118圈，只是逃逸軌跡上有些微小的變化；圖5則清晰地反映了Isp控制由無(wú)約束到約束逐步變嚴(yán)，到最終約束在［3 600m／s，5 400m／s］區(qū)間的過(guò)程。

圖4 平面內(nèi)30天燃料最優(yōu)地球逃逸軌跡Fig.4 Fuel optimal earth escape trajectory for 30days in planar

圖5 Isp隨時(shí)間變化曲線Fig.5 Relations between Ispand transfer time

4.2 三維控制有約束長(zhǎng)時(shí)間燃料最優(yōu)地球逃逸

對(duì)于三維逃逸問(wèn)題，在逃逸點(diǎn)加上軌道傾角約束，假定航天器相對(duì)于地球以特定的軌道傾角逃逸。

軌道傾角在球坐標(biāo)系中可表示為

式中 h為航天器的角動(dòng)量；i為軌道傾角。在逃逸點(diǎn)加上傾角約束之后，燃料最優(yōu)逃逸問(wèn)題的增廣約束方程為

式中 γ為由加入傾角約束而引入的待確定增廣系數(shù)；itar為三維逃逸問(wèn)題的目標(biāo)傾角要求。

利用Isp無(wú)約束平面內(nèi)燃料最優(yōu)逃逸問(wèn)題的解作為初值猜測(cè)，可逐步求得Isp無(wú)約束三維燃料最優(yōu)逃逸問(wèn)題的解；有了Isp無(wú)約束三維燃料最優(yōu)逃逸問(wèn)題的解，逐步加入Isp控制約束，進(jìn)而最終求解Isp有約束三維燃料最優(yōu)逃逸問(wèn)題。仿真同樣以半長(zhǎng)軸為7 500km 的地球赤道圓軌道為初始停泊軌道，假定要求30天后達(dá)到逃逸能量，Isp約束范圍為 ［2 700m／s，6 500m／s］，且逃逸點(diǎn)軌道傾角為30°。以Isp無(wú)約束平面內(nèi)30天燃料最優(yōu)逃逸的最優(yōu)解為初值，將終端軌道傾角約束由3°開(kāi)始，逐步加大到最終設(shè)計(jì)的逃逸傾角，即30°，上一步得到的最優(yōu)解作為下一步計(jì)算的初值猜測(cè)，最終求解Isp無(wú)約束、末端30°傾角的燃料最優(yōu)逃逸問(wèn)題。在得到Isp無(wú)約束30°燃料最優(yōu)逃逸收斂解之后，以此為初值，逐步加嚴(yán)Isp控制約束，同樣利用上一步得到的最優(yōu)解作為下一步計(jì)算的初值猜測(cè)，最終求得Isp∈［2 700m／s，6 500m／s］的30°燃料最優(yōu)逃逸軌道。圖6給出了逃逸傾角為30°、Isp約束范圍為 ［2 700m／s，6 500m／s］的30天燃料最優(yōu)逃逸軌跡，航天器逃逸過(guò)程共繞地球飛行118圈，消耗燃料14489.6kg，約占總質(zhì)量的14.5%；圖7給出了30°傾角逃逸，Isp無(wú)約束和有約束情況下Isp隨轉(zhuǎn)移時(shí)間變化的關(guān)系曲線，可以看出，Isp有約束和無(wú)約束隨轉(zhuǎn)移時(shí)間變化規(guī)律基本相同，從另一個(gè)側(cè)面說(shuō)明了以無(wú)約束問(wèn)題的收斂解作為有約束問(wèn)題的初值猜測(cè)的合理性。

6 30天Isp有約束地球30°逃逸燃料最優(yōu)逃逸軌跡Fig.6 30days earth escape trajectory under 30°

圖7 地球30°燃料最優(yōu)逃逸Isp隨時(shí)間變化曲線Relations between Ispand transfer time under 30°

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

本文通過(guò)在燃料最優(yōu)性能指標(biāo)中引入待定參數(shù)，使由此推導(dǎo)的兩點(diǎn)邊值問(wèn)題的各待定參數(shù)初值約束在 ［－1，1］區(qū)間內(nèi)，初值猜測(cè)范圍縮小，大大提高了間接法的收斂性。數(shù)值仿真結(jié)果表明，改進(jìn)的間接法結(jié)合同倫思想和曲線擬合技術(shù)，能提供較好的初值猜測(cè)，從而減少迭代次數(shù)，提高收斂速度，能有效地求解復(fù)雜的長(zhǎng)時(shí)間控制有約束三維逃逸問(wèn)題，是一種高效的小推力逃逸軌道設(shè)計(jì)方法。此外，本文以小于0.5天的時(shí)間間隔，計(jì)算并繪制了逃逸時(shí)間從0.1～30天的控制無(wú)約束平面內(nèi)燃料最優(yōu)逃逸問(wèn)題的最優(yōu)解，其呈現(xiàn)出的規(guī)律趨勢(shì)不僅可以用曲線擬合技術(shù)外推產(chǎn)生較好的長(zhǎng)時(shí)間最優(yōu)逃逸初值猜測(cè)，也將為進(jìn)一步探索其中的規(guī)律提供一定的參考。

［1］BETTS J T.Very low－thrust trajectory optimization using a direct SQP method［J］.Journal of Computational and Applied Mathematics，2000，120：27－40.

［2］嚴(yán)輝，吳宏鑫，吳新珍.小推力軌道優(yōu)化研究 ［J］.中國(guó)空間科學(xué)技術(shù)，1998，18（2）：8－13 YAN HUI，WU HONGXIN，WU XINZHEN.Study on optimization of low－thrust trajectories ［J］.Chinese Space Science and Technology，1998，18（2）：8－13.

［3］RANIERI C L，OCAMPO C A.Indirect optimization of spiral trajectories［J］.Journal of Guidance，Control，and Dynamics，2006，29（6）：1360－1366.

［4］LEE D H，BANG H.Efficient initial costates estimation for optimal spiral orbit transfer trajectories design［J］.Journal of Guidance，Control，and Dynamics，2009，32（6）：1943－1947.

［5］NAH R S，VADALI S R.Fuel－optimal，low－thrust，three－dimensional earth－mars trajectories［J］.Journal of Guidance，Control，and Dynamics，2001，24（6）：1100－1107.

［6］RANIERI C L，OCAMPO C A.Indirect optimization of three dimensional finite－burning interplanetary transfers including spiral dynamics［J］.Journal of Guidance，Control，and Dynamics，2009，32（2）：444－454.

［7］JIANG F H，BAOYIN H X，LI J F.Practical techniques for low－thrust trajectory optimization with homotopic approach ［J］.Journal of Guidance，Control，and Dynamics，2012，35（1）：245－258.

［8］HULL D G.Optimal control theory for applications［M］.New York：Springer，2003：199－247.