李 鑒，韓 潮

(北京航空航天大學(xué)，北京100191)

0 引言

小推力推進系統(tǒng)具有比沖高、消耗工質(zhì)少，有效載荷比高，體積小，工作時間長，推力精確可控等特點，所以在空間任務(wù)中得到越來越多的應(yīng)用。

小推力軌道轉(zhuǎn)移具有機動時間長，非線性強的特點，其優(yōu)化問題一直是空間軌跡優(yōu)化領(lǐng)域的研究熱點和難點，目前主要的研究內(nèi)容是時間最短軌道轉(zhuǎn)移問題和燃料最優(yōu)軌道轉(zhuǎn)移問題。就數(shù)值解法而言，連續(xù)推力軌跡優(yōu)化的求解方法可分為直接法和間接法。直接法通過不同的離散化方法對控制量進行參數(shù)化，把軌道優(yōu)化問題轉(zhuǎn)化為帶有約束的參數(shù)優(yōu)化問題，然后采用非線性規(guī)劃算法進行求解。這類方法目前的研究熱點是偽譜法［1－2］，該方法與傳統(tǒng)方法相比具有求解精度較高，尋優(yōu)參數(shù)少，收斂性好的優(yōu)點，但它仍然需要猜測初值，而且嚴重依賴非線性規(guī)劃軟件的求解能力。目前所見文獻中，用偽譜法求解小推力最優(yōu)軌道轉(zhuǎn)移問題多局限于同平面圓軌道之間。對于大尺度長時間小推力軌道機動問題，其求解效果尚不明確。

間接法主要通過變分法和極大值原理將原始的連續(xù)推力軌跡優(yōu)化問題轉(zhuǎn)化為兩點邊值問題，然后利用各種數(shù)值方法來求解該問題，常用的方法為打靶法。在應(yīng)用打靶法求解兩點邊值問題時，有3個難點:(1)需要獲取目標(biāo)函數(shù)或者終值誤差函數(shù)對于狀態(tài)量和協(xié)態(tài)變量的梯度矩陣，這些矩陣推導(dǎo)復(fù)雜，并且可能出現(xiàn)病態(tài)。一般采用數(shù)值微分方法避免推導(dǎo)，但精度有限。(2)打靶法對于兩點邊值問題中的協(xié)態(tài)變量初值很敏感，收斂半徑小，而協(xié)態(tài)變量又缺乏物理意義，難以估計初值。于是一部分學(xué)者研究了協(xié)態(tài)變量的初值猜測技術(shù)，如魯棒性算法［3］，進化算法［4］和粒子群算法［5］。隨機搜索算法雖然可以得到全局最優(yōu)解，但是計算效率太低，不適合需要長時間積分的小推力軌跡優(yōu)化問題。Lee［6］提出了一種新的螺旋形軌道協(xié)態(tài)變量初值猜測方法，但僅適用于平面內(nèi)軌道機動。(3)控制切換時刻不易確定，文獻［7］提出了控制切換點的搜索算法，但無法保證一定能獲得全部控制切換點。

文獻［8－10］通過構(gòu)造兩層同倫，以控制連續(xù)、求解相對容易的能量最優(yōu)問題解作為初值，逐步逼近燃料最優(yōu)問題的解。這樣使每一層優(yōu)化問題都獲得了良好的協(xié)態(tài)變量初值，克服了協(xié)態(tài)變量初值猜測的困難。文獻［11］在文獻［9］算法基礎(chǔ)上對協(xié)態(tài)變量進行歸一化處理，并加入控制切換點搜索算法和粒子群算法以獲取全局最優(yōu)解。同倫法具有較高的穩(wěn)定性、精確性和求解效率，是目前小推力軌跡優(yōu)化文獻中最有效的算法。但是它引入了兩個新的優(yōu)化問題，需要另外構(gòu)造其最優(yōu)控制律和相關(guān)的梯度矩陣，增加了求解的繁瑣程度，并且還需要使用同倫算法程序包，這些都增加了該算法對于工具軟件的依賴程度。

針對上述研究現(xiàn)狀，本文提出一種基于UKF參數(shù)估計算法的小推力最優(yōu)軌道轉(zhuǎn)移問題求解方法，將軌道機動最優(yōu)控制問題所對應(yīng)的兩點邊值問題轉(zhuǎn)化為參數(shù)估計問題，然后用UKF濾波器進行參數(shù)估計。該方法流程簡單，不依靠強大的非線性規(guī)劃軟件包。其基本原理是估計理論，不依賴梯度信息，故不需要良好的協(xié)態(tài)變量初值，也不用推導(dǎo)梯度矩陣，同時又具有較好收斂性和精確性，可以很好地解決小推力軌道機動優(yōu)化問題。

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

1. 1 原始優(yōu)化問題描述

小推力軌道轉(zhuǎn)移的時間歷程通常很長，飛行圈數(shù)可達幾百甚至上千圈，如果采用直角坐標(biāo)描述航天器的運動會導(dǎo)致狀態(tài)量的時間歷程振蕩劇烈，不利于快速積分。本文采用改進春分點軌道要素來描述航天器的軌道運動，這樣既有利于提高積分速度和精度，又避免了可能出現(xiàn)的奇點問題。改進春分點軌道要素［p，ex，ey，hx，hy，L］與經(jīng)典 Kepler 軌道要素的關(guān)系如公式(1)所示:

式中:p為軌道半通徑，a為軌道半長軸，e為偏心率，i為軌道傾角，Ω為升交點赤經(jīng)，ω為近地點幅角，θ為真近點角。同時，定義如下輔助變量:

以改進春分點軌道要素表示的航天器軌道動力學(xué)方程為:

式中:

式(6)中，a為除二體引力加速度之外的其他外力加速度矢量，［ar，at，ah］T為 a 在軌道徑向、垂直于徑向和動量矩方向的分量，aeng為發(fā)動機推力加速度矢量，adis為其他攝動加速度矢量。本文的討論范圍為只有二體引力加速度和發(fā)動機推力加速度情況下的小推力軌道機動優(yōu)化問題，其他攝動加速度的影響將在后續(xù)工作中研究。為簡化起見，在后文公式中用a代替aeng作為發(fā)動機推力加速度矢量。

設(shè)航天器質(zhì)量為m，發(fā)動機推力大小具有上限Tmax和常值比沖Isp，g0=9.780 4m/s2為赤道重力加速度，則方程(3)～(6)可改寫為

1. 2 最優(yōu)控制兩點邊值問題

要求解MF問題必須先求解TF問題，因為MF問題的飛行時間必須大于求解TF問題所得的最短飛行時間tfmin。本文基于極大值原理求解最優(yōu)控制問題，由文獻［9－10］可知，TF問題與MF問題所滿足的狀態(tài)變量與協(xié)態(tài)變量微分方程形式一致，兩種問題的最優(yōu)控制律也具有相似的形式。應(yīng)用本文提出的算法求解TF問題與MF問題的過程完全相同，并且TF問題更容易求解，因此將MF問題作為討論重點，對TF問題只給出求解結(jié)果以便與文獻的算例作比較。

引入與改進春分點軌道要素對應(yīng)的協(xié)態(tài)變量λ= ［λp，λex，λey，λhx，λhy，λL］T和與質(zhì)量對應(yīng)的協(xié)態(tài)變量λm。MF問題的Hamilton函數(shù)為:

其中α為矢量u與BTu之間的夾角。與終端狀態(tài)~x(tf)和終端時刻tf相關(guān)的性能指標(biāo)K為:

狀態(tài)約束為:

系統(tǒng)定常，根據(jù)極大值原理，最優(yōu)控制的必要條件為:

狀態(tài)與協(xié)態(tài)方程:

初始條件:

引入狀態(tài)約束g的拉格朗日乘子ν，則末值條件為:

Hamilton函數(shù)對最優(yōu)控制u*有極小值，即:

由式(8)和(14)可知，控制取最優(yōu)時α=π。系統(tǒng)定常，Hamilton函數(shù)不顯含時間，故當(dāng)控制取最優(yōu)時H在［t0，tf］上為常數(shù)。

考察λm的方程

易知當(dāng)α=π時，λm的導(dǎo)數(shù)非正。又由式(13)可得當(dāng)控制為最優(yōu)控制時，tf時刻λm應(yīng)為0，故λm在［t0，tf］上非負，這一結(jié)論對TF問題也成立。令開關(guān)函數(shù)ψ為:

由文獻［9］可知最優(yōu)控制律應(yīng)具有如下形式:

下面討論如何將上述兩點邊值問題改寫成參數(shù)估計問題，然后利用UKF濾波估計算法求解。

2 UKF參數(shù)估計及其求解方法

參數(shù)估計問題，又被稱為系統(tǒng)辨識或者機器學(xué)習(xí)問題，其目的是為了確定一個非線性映射:

該映射的輸入為xk，輸出為yk，w為非線性映射的參數(shù)。一般來說，映射輸入xk和期望輸出dk是不變的，輸出誤差定義為ek=dk－G(xk，w)，求解參數(shù)估計問題就是要估計w的均值，使映射G(xk，w)的輸出誤差最小。用濾波器求解參數(shù)估計問題的相關(guān)理論，文獻［12］作了詳細闡述，本文只作簡要介紹。

將原始參數(shù)估計問題寫成狀態(tài)空間表達式:

上式代表一個狀態(tài)轉(zhuǎn)移矩陣為單位陣的靜態(tài)過程，rk為過程噪聲，而期望輸出dk則與對wk的非線性觀測相對應(yīng)，ek為觀測噪聲。這樣，原參數(shù)估計問題就可以用EKF，UKF等濾波器求解?；赨KF濾波器的參數(shù)估計算法流程［12］如圖1所示。

圖1 UKF參數(shù)估計算法Fig.1 UKF parameter estimation

UKF參數(shù)估計算法公式中，

Rr和Re分別是過程噪聲和觀測噪聲的協(xié)方差陣。N是w的維數(shù)，η是尺度參數(shù)，常量ε決定了UT變換的σ點相對于w當(dāng)前均值的分布范圍，一般設(shè)為小量，取值范圍為［10－4，1］。常量 κ 也是個尺度參數(shù)，一般取為0或者3－N。β是與w的先驗分布相關(guān)的常量，對于高斯分布，β=2是最優(yōu)的。ρRLS是遺忘因子，用于防止因模型誤差較大造成的濾波發(fā)散，其取值范圍為(0，1］。關(guān)于這些常量選取方法的詳細討論可以參看文獻［12］。

參數(shù)估計問題相當(dāng)于求解如下以w為優(yōu)化變量的優(yōu)化問題:

式(23)中

定義觀測誤差為qk=dk－G(xk，wk)，將式(21)重寫為如下觀測誤差形式:

上式表明觀測誤差的期望值為0，求解原始參數(shù)估計問題最終轉(zhuǎn)化為求解式(24)所代表的非線性過程的參數(shù)估計問題。

3 UKF參數(shù)估計算法求解小推力最優(yōu)軌道轉(zhuǎn)移問題

將式(19)代表的燃料最優(yōu)兩點邊值問題改寫為參數(shù)估計問題，選擇待估計參數(shù)w為:

觀測誤差為q:

其中G為如下微分方程組初值問題的解:至此，TPBVPMF問題(19)與式(24)形式上一一對應(yīng)，已轉(zhuǎn)化為參數(shù)估計問題，可通過UKF參數(shù)估計算法求解。

4 數(shù)值算例與結(jié)果分析

首先求解TF問題。以表1給出的初末軌道參數(shù)和發(fā)動機參數(shù)為例，求解不同最大推力值下的時間最短小推力軌道轉(zhuǎn)移問題，并將計算結(jié)果(由PE法表示)與同倫法的計算結(jié)果［10］作比較，如表2所示。

表1 初始計算參數(shù)Table 1 Initial parameters

PE法可以直接用隨機初值得到Tmax=0.01N時TF問題的解，文獻［10］只給出了直至Tmax=0.14N時的計算結(jié)果，所以在表中沒有數(shù)據(jù)。

文獻［8］證明了當(dāng)最大推力值趨近于0時，最小轉(zhuǎn)移時間和最大推力值之間存在如下簡單關(guān)系:

式中C為常數(shù)，僅與機動過程的初末軌道參數(shù)和發(fā)動機參數(shù)相關(guān)。表2中的PE數(shù)據(jù)是完全符合上述關(guān)系式的，也說明本文算法能夠正確求解時間最短小推力軌道轉(zhuǎn)移問題。

表2 最短軌道轉(zhuǎn)移時間Table 2 Minimum transfer times

下面以表1給出的初末軌道參數(shù)、發(fā)動機參數(shù)為例，取飛行時間參數(shù)ctf=1.5，求解不同最大推力值下的燃料最優(yōu)小推力軌道轉(zhuǎn)移問題，并在表3中與同倫法的計算結(jié)果作比較。

文獻［9］所給出的計算條件為SUN-Blade 1 000，本文所得數(shù)據(jù)的計算條件為Intel i3 3.10GHz，由于沒有確切的計算能力參數(shù)，只能根據(jù)CPU的浮點計算能力估計后者的計算能力是前者的10倍。由表3可見，兩種方法的求解效率大致相當(dāng)。

表3 不同T max情況下的燃料最優(yōu)軌道轉(zhuǎn)移問題控制切換次數(shù)和計算時間Table 3 Switches and computation times of fuel-minimum orbit transfer problem with various T max

可以看出在推力較大的情況下兩者的計算結(jié)果是非常吻合的。當(dāng)推力小于1N時，結(jié)果出現(xiàn)了差別，但由表中數(shù)據(jù)算出的比值仍然是吻合得比較好的，基本上保證了一圈內(nèi)有2次控制切換，這與文獻［9］中得出的結(jié)論一致。數(shù)據(jù)的差別可以解釋為兩種方法得到了不同的局部最優(yōu)解。文獻［9］明確指出，求解小推力TPBVPMF問題(19)時極易陷入局部最優(yōu)解，其具體表現(xiàn)就是相同的飛行時間對應(yīng)不同的Lf，即不同的飛行圈數(shù)。事實上，文獻［9］得到的最優(yōu)解最終剩余質(zhì)量大約為1 382kg，而本文算法得到的最終剩余質(zhì)量約為1 378kg，從實際應(yīng)用角度出發(fā)，這個偏差是可以接受的。

為了驗證本文方法所得結(jié)果的精確性，下面以Tmax=10N為例求解TPBVPMF問題(19)，結(jié)果如圖2～圖3所示。對比圖2～圖3和文獻［9］的圖示結(jié)果，兩者幾乎是完全一致的。圖3中Hamilton函數(shù)值在整個機動過程中近似保持不變則驗證了所得結(jié)果的最優(yōu)性。

5 問題與經(jīng)驗總結(jié)

圖2 燃料最優(yōu)軌道機動問題最優(yōu)解的軌道要素隨時間變化歷程(最大推力值為10N)Fig.2 Modified equinoctial elements vs time for minimum-fuel orbit transfer problem(T max=10N)

圖3 燃料最優(yōu)軌道機動問題的最優(yōu)控制量和Hamilton函數(shù)隨時間變化歷程(最大推力值為10N)Fig.3 Optimal control and Hamilton function vs time for minimum-fuel orbit transfer problem(T max=10N)

(1)應(yīng)用UKF參數(shù)估計方法，協(xié)態(tài)變量初值可以隨機選取，但λm的初值一定要為正。初值的選取范圍與構(gòu)造問題時使用的量綱是緊密相關(guān)的。本文所有算例的長度量綱為106m，時間量綱為小時，以便與文獻［9］對比結(jié)果，不一定是最優(yōu)的選擇。針對不同的軌道轉(zhuǎn)移問題，選取合適的量綱，使各個協(xié)態(tài)變量量級盡量接近，可以提高求解效率。

(2)UKF濾波算法公式中的過程噪聲Rr、常量ε和遺忘因子ρRLS對濾波收斂過程的影響很大。這些量的選取和更新方法屬于UKF濾波器算法改進范疇，不是本文的討論重點，可以作為下一步工作。本文通過大量數(shù)值仿真，總結(jié)常量ε和遺忘因子ρRLS的選取和更新方法如下:對于常量ε，濾波初始階段ε應(yīng)盡量取較大值(如0.8)，隨著濾波迭代的進行，ε應(yīng)該隨觀測誤差的減小而減小，這樣有利于算法的快速收斂。本文的算例中，ε最終可減小至0.001。對于遺忘因子ρRLS，當(dāng)遺忘因子ρRLS較小時，濾波過程的振蕩幅度很大，反之則振蕩幅度較小。推力較大時(1N/1 500kg左右)，采用較小的 ρRLS(如 0.2 ～0.3)，濾波可以很快收斂。當(dāng)推力較小時，則應(yīng)采用較大的 ρRLS(如0.5 ～ 0.7)。

6 結(jié)論

本文提出的基于UKF參數(shù)估計算法的小推力軌道轉(zhuǎn)移優(yōu)化設(shè)計方法具有簡潔、精確、高效的優(yōu)點，可作為求解小推力軌道機動優(yōu)化問題的一個有力工具。

［1］ Ross I M，Gong Q，Sekhavat P.Low-thrust，high-accuracy trajectory optimization［J］.Journal of Guidance and Control，2007，30(4):921－933.

［2］ 尚海濱，崔平遠，徐瑞，等.基于高斯偽光譜的星際小推力轉(zhuǎn)移軌道快速優(yōu)化［J］.宇航學(xué)報，2010，31(4):1005－1011.［Shang Hai-bin，Cui Ping-yuan，Xu Rui，et al.Fast Optimization of interplanetary low-thrust transfer trajectory based on Gauss pseudospectral algorithm［J］.Journal of Astronautics，2010，31(4):1005 －1011.］

［3］ 劉滔，何兆偉，趙育善.持續(xù)推力時間最優(yōu)軌道機動問題的改進魯棒算法［J］.宇航學(xué)報，2008，29(4):1216－1221.［Liu Tao，He Zhao-wei，Zhao Yu-shan.Continuous-thrust orbit maneuver optimization using modified robust algorithm［J］.Journal of Astronautics，2008，29(4):1216 －1221.］

［4］ Igarashi J，Spencer D B.Optimal continous thrust orbit transfer using evolutionary algorithms［J］.Journal of Guidance，Control，and Dynamics，2005，28(3):547 －549.

［5］ Pontani M，Conway B A.Particle swarm optimization applied to space trajectories［J］. Journal of Guidance， Control， and Dynamics，2010，33(5):1429 －1441.

［6］ Lee D H，Bang H C.Efficient initial costates estimation for optimal spiral orbit transfer trajectories design［J］.Journal of Guidance，Control，and Dynamics，2009，32(6):1943 －1947.

［7］ Jamison B R，Coverstone V.Analytical study of the primer vector and orbit transfer switching function［J］.Journal of Guidance，Control，and Dynamics，2010，33(1):235 －245.

［8］ Bombrun A，Pomet J B.Asymptotic behavior of time optimal orbital transfer for low thrust 2-body control system［J］.Discrete and Continuous Dynamical Systems，2007(Supplement):122 －129.

［9］ Haberkorn T，Martinon P，Gergaud J.Low-thrust minimum-fuel orbital transfer:a homotopic approach［J］.Journal of Guidance，Control，and Dynamics，2004，27(6):1046 －1060.

［10］ Caillau J，Gergaud J，Noailles J.3D geosynchronous transfer of a satellite:continuation on the thrust［J］.Journal of Optimization theory and Applications，2003，118(3):541 －565.

［11］ Jiang F H，Baoyin H X，Li J F.Practical techniques for lowthrust trajectory optimization with homotopic approach［J］.Journal of Guidance，Control，and Dynamics，2012，35(1):245－258.

［12］ Haykin S.Kalman Filtering and neural networks［M］.USA:John Wiley＆ Sons Inc，2002.