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.