黃鎬 韓潮

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

1 引言

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

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

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

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

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

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

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

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

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

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

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

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

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

式中 λ＝［λr，λθ，λφ，λvr，λvθ，λvφ］T為與狀態(tài)量x相對應(yīng)的協(xié)狀態(tài)量向量；λv＝［λvr，λvθ，λvφ］T為與速度v對應(yīng)的協(xié)狀態(tài)量向量；λm為與m 對應(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)逃逸問題的終端約束為終點時刻tf時的軌道能量εf達(dá)到拋物線逃逸軌道能量值εtar。

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

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

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

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

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

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

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

4 仿真算例

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

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

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

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

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

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

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

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

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

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

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

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

對于三維逃逸問題，在逃逸點加上軌道傾角約束，假定航天器相對于地球以特定的軌道傾角逃逸。

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

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

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

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

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

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

5 結(jié)束語

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