J2攝動下多脈沖軌道的主矢量方法求解*

孟雅哲

1. 中國科學院空間應用工程與技術中心太空應用重點實驗室，北京100094 2.中國科學院大學，北京100049

J2攝動二體模型下、任意轉移時間和始末軌道狀態的最小速度脈沖軌道的求解，屬于地球引力場中交會軌道設計的理論研究范疇。地球引力場中需應用交會軌道的航天活動將會越來越多，如大型航天器的在軌組裝、空間站的長期運營、太空碎片處理及衛星捕獲等。這些活動的軌道受J2攝動影響較大，且多為脈沖機動軌道。脈沖軌道優化的一個基礎工具是主矢量方法(Primer Vector Method,PVM)。

PVM被廣泛應用于各種軌道模型。PVM包括主矢量必要條件[1]和不滿足該必要條件的脈沖軌道的改進[2-3]。PVM基于無攝二體模型被提出，被應用于不考慮[4-6]和考慮J2攝動[7-9]的相對運動模型和限制性三體模型[10-11]。不同模型下的PVM理論各有其發展，線性模型中主矢量必要條件被證明為充要條件[5]且被進一步簡化[6]，三體模型中PVM基本公式被給出和證明[10-11]。目前尚無J2攝動二體模型的文章。

J2攝動二體模型下PVM有其發展潛力。該模型下的多脈沖軌道可用基于梯度的數值迭代方法求得。迭代所用初值可以是無攝二體模型多脈沖軌道解[12-13]，也可以是不基于梯度的全局搜索方法直接給出的帶J2攝動的軌道解[14]。但已有文獻中未對所得J2攝動多脈沖解的局部最優性加以驗證。該局部最優性可以用PVM進行驗證和保證。J2攝動二體模型下PVM 的實現也會隨著J2Lambert求解方法[14-15]、并行計算等工具的發展變得容易。本文基于J2攝動二體模型，對PVM的基本公式進行了驗證和整理，并給出了PVM的一種實現方法。算例驗證了該方法的有效性，并給出了一種多解的情況。

1 J2攝動模型及主矢量必要條件

考慮J2攝動時，航天器運動方程為：

(1)

式中：r=[x,y,z]T和v=[vx,vy,vz]T分別為航天器的位置和速度；A和α分別為推進器推力的大小和方向；g(r)=-μr/r3為質點間的萬有引力項；aJ2為式(2)所示J2攝動加速度。上述所有矢量均在地心直角慣性坐標系下表達：

aJ2=1.5μJ2Re(5z2r/r7-Ar/r5)

(2)

式中：μ=398601km3/s2為地球引力常數，Re=6378km為地球半徑，J2=1.0823×10-3為帶諧系數,A=diag(1,1,3)為3×3矩陣。

設航天器轉移軌道的初始和末端時刻分別為t0和tf，對應的航天器位置和速度滿足約束：

r(t0)=r0,v(t0)=v0,r(tf)=rf,v(tf)=vf

(3)

最小速度脈沖性能指標為：

(4)

下面用Pontryagin極小值原理，推導式(1)和(3)約束下、使得式(4)所示J取極小值的A和α所需滿足的必要條件。

首先構造Hamilton函數：

(5)

其中：λr和λv分別為r和v對應的協態變量，滿足微分方程：

(6)

根據Pontryagin極小值原理，使得性能指標J達到極小值的A和α，需使H達到極小值，故有：

(7)

根據動量定理，對于脈沖推進軌道來說，式(4)等價于式(8)

(8)

2)p≤1，其中“=”在且僅在脈沖作用點ti,i=1,2,…,N處成立；

3)推進器施加的脈沖方向和主矢量方向相同，即Δvi/Δvi=p(ti)/p(ti)；

上述主矢量必要條件形式與文獻[1]中給出的相同，且可用文獻[1]中的方法進行證明。

2 J2攝動下PVM公式

下面將給出J2攝動模型下需驗證的PVM基本公式，及PVM中脈沖個數不變和加一時的性能指標J的變化公式及其分析。文中未列出的公式推導過程與文獻[2-3,16]平行。

2.1 基本公式

對于無推力軌道段，式(1)中的A=0，可有：

(9)

注意到G(r)和AJ2均為對稱矩陣，容易驗證：

故在無推力軌道段中，有伴隨方程：

(10)

(11)

ti(ti≠t0且ti≠tf)時刻脈沖作用前后H的變化量為：

(12)

式(10)～(12)將被用于下文主矢量方法公式的推導[2]。

2.2 非最優軌道調整公式

(13)

(14)

(15)

為使新添加脈沖滿足式(16)所示主矢量必要條件，添加脈沖前后tm時刻航天器的位置應從rm變為rm+δrm，其中δrm滿足式(17)：

(16)

(17)

將式(16)代入式(14)可得：

dJ=cm[1-pm]

(18)

可見，當且僅當pm>1時，dJ<0，即J變小；當tm取pm達最大值的時刻，J下降最大。

(19)

其中：a1和a2表達式參看文獻[2-3]。

3 J2攝動下PVM的實現

PVM求解多脈沖軌道的流程如圖 1所示。本節將依圖中流程敘述J2攝動下PVM各步驟的具體實現過程。

圖1 PMV整體流程圖

3.1 脈沖軌道的求解

本文應用表 1所示NLP問題(NLP1)求解某一N值下的J2攝動多脈沖軌道。

表1 求解多脈沖軌道的NLP問題模型(NLP1)

表1中：r(tf)和v(tf)是以r0和v0為初值，用式(1)和各優化變量值求得的tf時刻的位置和速度；rf和vf是tf時刻的位置和速度期望值；不等式約束中：若有始端脈沖，則t0=t1，否則t0

NLP1在PVM中將被多次構造和求解。首次求解NLP1所用初值可以是帶J2的Lambert解[15]、特征點變軌給出的解[19]、含有較多脈沖個數的直接法的解[8]及J2相對模型下的脈沖解[9]等。迭代過程中NLP問題的初值則用3.4節中的方法給出。

3.2 主矢量曲線的求解

3.3 軌道局部最優性判斷

上面求得的p在[t0,tf]上連續，且滿足Δvi/Δvi=p(ti)/p(ti)和p(ti)=1。所以，若該主矢量曲線滿足p(t)≤1，即滿足主矢量必要條件。實際應用中，應首先檢驗有無小脈沖，若無且p(t)≤1，認為已得到了局部最優(滿足主矢量必要條件)的脈沖軌道；否則轉3.4進行處理。

3.4 脈沖作用情況調整

依次檢測如下5種情況是否出現，若某種情況出現，應用本節內容給出脈沖軌道初值，該初值將用于構造和求解表 1所示NLP問題。

情況1：若某個脈沖大小低于設定閾值(如1m/s)，刪除該脈沖點；

情況2：若2個脈沖作用時刻的間隔低于設定值(如1s)，令其中一個速度脈沖為2個速度脈沖的矢量和，刪除另一個速度脈沖；

情況5：若存在使得p(t)>1的時刻t，查找p最大的時刻tm，并在該時刻添加脈沖。ti-1

考慮到NLP求解器具有尋找使J更小的解的能力，在tmtN(其中tN≠tf)時，將tm時刻脈沖初值設為0。

4 算例

本節將給出2個算例，算例一較為詳細地說明了主矢量方法的計算過程，算例二是一個主矢量方法多解情況的例子。

4.1 算例一

設t0=0s，對應經典軌道參數為：半長軸7158.651km、偏心率0.00881、軌道傾角97.689°、升交點赤經301.395°、近地點幅角63.215°和真近點角213.185°，tf=1.241×105s，對應軌道根數為：7168.813km，0.00879，99.234°，285.619°，133.129°，107.167°。應用特征點變軌設定三脈沖軌道初值如式(20)所示，對應的總速度脈沖為2.2673km/s。

0s:[-2.023 -0.809 -0.367]km/s
1.186×105s:[0.005 0.002 0.007]km/s
1.241×105s:[0.017 0.044 0.013]km/s

(20)

以式(20)為初值，構造表 1所示NLP問題，并求得第1個軌道解，而后用PVM對脈沖作用情況進行了4次調整，得到滿足主矢量必要條件的解。詳細調整過程見表 2。滿足必要條件的解如式(21)所示。

表2 主矢量方法迭代過程中各解的情況

圖2 各脈沖解對應的主矢量(Primer Vector，PV)曲線

5.939×103s:[-0.204 -0.057 -0.042]km/s1.155×105s:[0.092 0.112 0.004]km/s1.185×105s:[-1.473 -0.705 -0.318]km/s1.216×105s:[0.162 -0.004 0.020]km/s

(21)

圖2給出了表 2中5個軌道解對應的主矢量曲線。根據該圖，PVM調整前后PV曲線整體形狀變化不大。實際上，本算例PVM調整前后軌道整體形狀和性能指標相對變化也較小。

4.2 算例二

設t0=0s，對應經典軌道參數為：7173.322km，1.773×10-4，99.916°，135.809°，148.415°，182.528°；tf=3.843×105s，對應軌道根數為：7138.306km，0.016，98.760°，310.085°，146.776°，210.694°。應用特征點變軌設定三脈沖軌道初值進行主矢量方法迭代求解。第1個可行解如式(22)所示

0s:[4.536 -5.705 -12.896]km/s
3.804×105s:[1.220 1.297 0.074]km/s
3.843×105s:-[0.043 0.017 0.057]km/s

(22)

0s:[-0.564 2.085 -0.107]km/s
1.985×103s:[-1.028 2.527 -1.745]km/s
3.787×105s:[-1.660 0.798 -2.106]km/s
3.836×105s:[-0.1656 0.805 -0.749]km/s

(23)

9.009×102s:[1.062 -0.381 2.257]km/s
2.893×104s:[0.348 -0.345 2.267]km/s
5.996×104s:[0.102 -0.447 1.084]km/s
3.746×105s:[0.137 -0.484 1.094]km/s
3.809×105s:[1.628 -5.043 10.496]
×10-2km/s

(24)

圖3 第1個解和2個最優解的主矢量曲線

圖4 第1個解和2個最優解對應的軌道

首個可行解和2個優化解對應的主矢量曲線參看圖 3，軌道參看圖 4。可以看出，主矢量方法調整前后主矢量曲線形狀變化較大，3條軌道形狀也相互不同。實際上，當轉移時間較長時，實現同一轉移要求的軌道形式不唯一，每一種形式下都可能有局部最優解存在。

5 總結與展望

PVM實現過程中，NLP和兩點邊值問題初值均選擇已有軌道解附近、使式(8)中J下降的值。若初值應用其他方式給定，如以無攝Lambert解作為初值求解J2Lambert問題，或可得到其他滿足主矢量必要條件的多脈沖解。

轉移時間較長時，滿足主矢量必要條件的軌道常多于1條(如4.2所示)，此時可將主矢量方法與延拓、進化方法和啟發式方法結合起來，以獲得更加接近全局最優的解。

[1] Lawden D F. Optimal Trajectories for Space Navigation[M]. London: Butterworths, 1963: 54-69.

[2] Handelsman M, Lion P M. Primer Vector on Fixed-time Impulsive Trajectories[J]. AIAA Journal, 1967, 6(1): 127-132.

[3] Jezewski D J. Primer Vector Theory and Applications[R]. NASA TR R-454. Washington: NASA, 1975.

[4] Luo Y X, Li H, Tang G J. Interactive Optimization Approach for Optimal Impulsive Rendezvous Using Primer Vector and Evolutionary Algorithms[J]. Acta Astronautica, 2010, 67(3): 396-405.

[5] Prussing J E. Optimal Impulsive Linear Systems: Sufficient Conditions and Maximum Number of Impulses[J]. Journal of the Astronautical Sciences, 1995, 43(2): 195-206.

[6] Carter T E. Necessary and Sufficient Conditions for Optimal Impulsive Rendezvous with Linear Equations of Motion[J]. Dynamics and Control, 2000, 10(3): 219-227.

[7] Arzelier D, Louembet C, Rondepierre A, et al. A New Mixed Iterative Algorithm to Solve the Fuel-Optimal Linear Impulsive Rendezvous Problem[J]. Journal of Optimization Theory and Applications, 2013, 159(1): 210-230.

[8] Roscoe C W, Westphal J J, Griesbach J D, et al. Formation Establishment and Reconfiguration Using Differential Elements in J 2-Perturbed Orbits[J]. Journal of Guidance, Control, and Dynamics, 2015, 38(9): 1725-1740.

[9] Riggi L, D′amico S. Optimal Impulsive Closed-form Control for Spacecraft Formation Flying and Rendezvous[C]. Boston: American Control Conference (ACC), 2016.

[10] Hiday L A. Optimal Transfers Between Libration-point Orbits in the Elliptic Restricted Three-body Problem[D]. West Lafayette: Purdue University, 1992.

[11] Davis K E, Anderson R L, Scheeres D J, et al. Locally Optimal Transfers Between Libration Point Orbits Using Invariant Manifolds[J]. Advances in the Astronautical Sciences, 2009, 135(3): 1623-1641.

[12] Jezewski D J.Optimal Rendezvous Trajectories Subject to Arbitrary Perturbations and Constraints[A]. In: AIAA/AAS ed. Astrodynamics Conference, Guidance, Navigation, and Control and Co-located Conference[C]. Hilton Head Island: American Institute of Aeronautics and Astronautics, 1992:235-245

[13] 羅亞中. 空間最優交會路徑規劃策略研究[D]. 北京： 國防科學技術大學, 2007.

[14] Yang Z, Luo Y Z, Zhang J, et al. Homotopic Perturbed Lambert Algorithm for Long-Duration Rendezvous Optimization[J]. Journal of Guidance, Control, and Dynamics, 2015, 38(11): 2215-2223.

[15] Woollands R M. Regularization and Computational Methods for Precise Solution of Perturbed Orbit Transfer Problems[D]. Texas: Texas A&M University, 2017.

[16] Cheng C W. Optimal Impulsive Trajectories for Orbital Rendezvous Between Elliptic Orbits[D].Ohio: University of Cincinnati, 1992.

[17] Hughes S P, Mailhe L M, Guzman J J. A Comparison of Trajectory Optimization Methods for the Impulsive Minimum Fuel Rendezvous Problem[J]. Advances in the Astronautical Sciences, 2003, 113.

[18] Conway B A. Spacecraft Trajectory Optimization[M]. New York: Cambridge University Press, 2010.

[19] Junkins J L, Schaub H. Analytical Mechanics of Space Systems, Second Edition[J]. AIAA Education, 2003: 512-548.

[20] Glandorf D R. Lagrange Multipliers and the State Transition Matrix for Coasting Arcs[J]. AIAA Journal, 1969, 7(2): 363-365.