基于混合Lie算子辛算法的不变流形计算1)

郑丹丹 罗建军 ,2) 张仁勇 刘 磊

∗(西北工业大学航天学院，西安710072)

†(航天飞行动力学技术重点实验室，西安710072)

∗∗(中国科学院空间应用工程与技术中心，北京100094)

††(北京飞行控制中心，北京100094)

基于混合Lie算子辛算法的不变流形计算1)

郑丹丹∗,†罗建军∗,†,2)张仁勇∗∗刘 磊†,††

∗(西北工业大学航天学院，西安710072)

†(航天飞行动力学技术重点实验室，西安710072)

∗∗(中国科学院空间应用工程与技术中心，北京100094)

††(北京飞行控制中心，北京100094)

平动点附近周期轨道的不变流形因其在低能轨道转移中起着重要作用而受到广泛关注.在设计低能轨道过程中不变流形要实时进行能量匹配，但利用传统数值积分方法进行积分时能量会耗散.显式辛算法具有比隐式辛算法计算效率高的优势，但其要求Hamilton系统必须分成两个可积的部分，而旋转坐标系下的圆型限制性三体问题是不可分的，因而显式辛算法难以用于求解旋转坐标系下的圆型限制性三体问题.本文通过引入混合Lie算子，成功实现了带三阶导数项的力梯度辛算法对圆型限制性三体问题的求解，并将基于混合Lie算子的带三阶导数项的辛算法与Runge-Kutta78算法和Runge-Kutta45算法进行仿真对比，仿真结果表明基于混合Lie算子的含有三阶导数项的辛算法位置精度高、能量误差小且计算效率高.利用基于混合Lie算子的带三阶导数项的辛算法计算不变流形，可以实现低能轨道转移过程中轨道拼接点的能量精准匹配.

圆型限制性三体问题，不变流形，辛算法，混合Lie算子

引言

平动点[1]是限制性三体问题中的动平衡点，在随两个天体一起旋转的参考系中处于静止状态.它们是几何意义上的点，其本身的应用价值十分有限.真正引人关注的是其周围大量存在的周期轨道[2]，以及与之相联的管状不变流形[3](包括稳定流形和不稳定流形).周期轨道是完成某些特殊任务的理想场所，而不变流形则提供了往返于主天体和周期轨道之间的低能耗转移途径[49].太阳系中不同三体系统平动点周期轨道的不变流形还存在相交的情形，从而拓展了系统间的深空转移轨道设计范围，因此寻找航天器从地球停泊轨道到平动点附近周期轨道的转移轨道显得尤为重要.

20世纪90年代，G´omez等[10]引入非线性动力系统理论，利用稳定流形将航天器送到了平动点附近，是转移轨道设计[11]的重大突破.这一发现不但节省了能量也为星际超级公路理论的诞生奠定了基础.但是不论是日--地系统还是地--月系统，不变流形距离地球都比较远，需要在推力或引力辅助[12]作用下将航天器从地球停泊轨道转移到稳定流形上.由于在优化拼接点时需要多次计算不变流形，如何选取不变流形的拼接点是进行低能轨道转移设计的关键.已经有学者对不变流形的计算进行研究，Zhang等[13]利用三次回旋插值法计算圆型限制性三体问题不稳定平动点周期轨道的不变流形，Lei等[1417]利用Lindstedt-Poincar´e方法求解了限制性三体和四体平动点周期轨道不变流形的高阶近似解析解，Dellnitz等[18]提出了方向集数值计算方法，Howell等[19]利用胞映射法[20]描述和存储不变流形的数据.这些计算方法提高了不变流形的计算速度，但是却没有关注不变流形的能量变化.由于圆型限制性三体问题是一个非线性自治哈密顿系统，没有严格的解析解，而且由于强非线性其对初值误差十分敏感，较小的计算误差将导致微分方程的解较大地偏离实际情况.因此需要寻找一种保能量的不变流形计算方法.

普通数值解法具有人为耗散性，长时间计算会导致系统总能量的耗散随时间的平方增长.在研究三体问题的过程中，对辛结构的破坏无疑会对动力学特性造成很大影响，使得系统的长期演化不能真实反映客观事实.辛积分方法[2122]由于具有保辛结构和能量守恒两个重要特性，近年来得到了快速发展[2326].显式辛算法具有比隐式辛算法计算效率高的优势，但其要求Hamilton系统必须可以分成两个可积的部分.旋转坐标系下的圆型限制性三体问题由于受Coriolis力的影响，不能像惯性系下的系统一样可分，因此，从形式上看显式差分格式对其不适用.廖新浩等[27]将圆型限制性三体问题的Hamilton函数分成动能与势能的和以及由坐标旋转产生的非惯性力两部分，然后利用算法合成构造差分结构，计算复杂.Sun等[28]分析了带有三阶导数项的力梯度辛算法在保能量上的优势，但是没有将该算法应用在旋转坐标系下的圆型限制性三体问题中，也没有将辛算法用在不变流形的求解上.

鉴于显式辛算法不适合求解旋转坐标系下圆型限制性三体问题，本文通过重新定义类动能的Lie算子，严格证明基于混合Lie算子的含有三阶导数项的显式辛算法可以求解旋转坐标系下的圆型限制性三体问题.首先，给出文中所用的动力学模型和不变流形的计算方法；然后，建立混合Lie算子并研究显式辛算法在求解Hamilton系统类动能不具有标准二次型的圆型限制性三体问题时的可行性；最后，将本文改进算法与 Runge-Kutta78(RK78)和 Runge-Kutta45(RK45)[29]算法进行仿真对比.

1 空间圆型限制性三体问题

1.1 动力学模型

圆型限制性三体问题描述的是质量为m3的航天器P3在两个绕着共同质心做匀速圆周运动的主天体P1和P2的引力场下运动，其质量分别为m1和m2，且m3≪m1≪m2，以P1和P2的质心为原点，P1指向P2的方向为x轴，y轴在两个主引力体旋转平面上，z轴垂直于xy平面，满足右手法则，如图1所示.

图1 圆型限制性三体问题Fig.1 Circular restricted three-body problem

假设P3在oxy平面运动，为了简化模型，对单位进行无量纲化处理.令引力常量G、主天体P1和P2之间的距离、旋转角速度ω、质量和均为1，因此 P1的质量为1−µ=m1/(m1+m2)，位置坐标为(−µ,0,0)，P2的质量为 µ =m2/(m1+m2)，位置坐标为(1−µ,0,0).文中所有量如无特殊说明都是无量纲的.在此旋转坐标系下航天器的运动方程和Largrange函数[30]分别为

和

Ω(x,y,z)是系统的势函数，通常表示为

其中r1,r2分别表示航天器与P1,P2之间的距离

通过Legendre变换

其中qi,pi(i=1,2,3)分别为航天器的坐标q=(x,y,z)和动量p=(px,py,pz)，可以得到圆型限制性三体问题的Hamilton函数

该系统具有雅可比能量积分

由Hamilton函数(2)，可推导出圆型限制性三体问题的正则方程为

1.2 不变流形

不变流形是与平动点周期轨道光滑连接的一簇空间轨道，航天器在不变流形上可以无动力飞行.因此，不变流形可以作为低能轨道转移的通道，在深空探测轨道转移设计中有着重大应用价值.

其中A(t)为系统(1)的雅可比矩阵，∆¯x为相对于不动点状态的偏移量.由¯x0点出发，积分一个周期后的状态所对应的状态转移矩阵Φ(0,T)称为单值矩阵[31].周期轨道上不同的点对应不同的单值矩阵，每个单值矩阵有一个不稳定特征根λs(λs＞1)及一个稳定特征根 λu=1/λs(λu＞ 1)，它们分别对应着特征向量和，它们包含了稳定流形和不稳定流形的方向信息.在点处沿特征向量方向增加一个微小状态扰动，则可以得到不变流形的积分初值.对于不变流形的计算有

一般都采用数值积分来计算不变流形.不变流形计算的关键是获得积分状态初值，而此初值和Halo轨道上的点有密切关系.由于圆型限制性三体问题的强非线性使得微分方程组(1)的解对初始值十分敏感，可能较小的计算误差都将导致不变流形的能量误差巨大.辛算法因具有保能量的特点而受到广泛关注，而显式辛算法具有比隐式辛算法效率高的优势，因此本文考虑利用含有三阶导数项的力梯度辛算法来计算不变流形.力梯度辛算法是将Hamilton系统分解成动能T(p)(仅关于动量p=(px,py,pz)的正定二次函数)和势能V(q)(关于坐标q=(x,y,z)的函数)两个可积部分，然后分别求解.圆型限制性三体问题的Hamilton函数有多种分解形式，最简单的是分解成类动能T(p)和势能V(q)两部分

其中，ypx−xpy是由于Coriolis力的影响而产生的偏移部分.从式(3)可以看出，类动能T(p)不具有标准的二次型形式，因此从形式上看力梯度辛算法不适合求解圆型限制性三体问题.后文试图通过改进算法来解决这一难题.

2 基于混合Lie算子的力梯度辛算法

考虑可以分解为动能T(p)和势能V(q)的Hamilton系统

其中，动能T(p)仅仅是关于动量p的正定二次函数,即T(p)=p2/2，势能V(q)是关于坐标q的函数.分别定义T(p)和V(q)的Lie算子[33-34]

由于圆型限制性三体问题Hamilton函数分解中类动能T(p)不是动量p的标准二次型，因此具有形式(4)的Lie算子不再适用，考虑将类动能T(p)的Lie算子变为位置与动量的混合型算子

它作用于方程组(3)第一个方程的位置坐标和动量得到微分方程组

该方程组的分析解为

其中,(x0,y0,z0,px0,py0,pz0)为初始状态，t是积分时间，于是新定义的混合算子¯A可积，说明重新定义的混合Lie算子是合理的.容易验证算子B也可积.

利用混合Lie算子¯A和B构造复杂算子

算子D是含有一阶、二阶和三阶导数项的算子，可以与混合Lie算子¯A和Lie算子B一起构造高阶算法

其中,h是积分步长，各积分系数分别为

其差分格式如下

这里，力 f= ∂V/∂q.

由

具体证明过程见附录A,得到

式中，E和F的具体表达形式见附录B.这表明算子D是两主天体的引力梯度，而不是旋转坐标系所产生的非惯性力与主天体引力的混合梯度.

因此，通过引入混合Lie算子，分解成形式(3)的圆型限制性三体问题可以用改进的显式辛算法进行求解.

3 仿真算例

仿真使用Matlab R2008b软件，计算机为Windows 7系统，配置Intel(R)Core(TM)2 Duo CPU E7500处理器，主频2.93 GHz，内存2.00GB，32位操作系统.利用第2节所得到的改进显式辛算法MTGS研究地--月圆型限制性三体问题并与RK78和RK45法进行仿真对比，初始值参数如表1所示.

表1 初始值参数Table 1 The initial value of the parameter

图2 三种算法积分两个周期得到的轨道Fig.2 The orbits of three algorithms integral two orbit periods

积分步长取0.001，积分两个周期得到的轨道如图2所示，可以看出本文改进的算法和RK78算法得到的轨道仍然能够闭合，而通过RK45算法得到的轨道开始偏离周期轨道.当积分时长为3个周期时，轨道如图3所示，改进的显式辛算法、RK78算法与RK45算法都不能保证轨道闭合，RK45算法发散速度最快，改进的算法发散速度最慢.虽然Halo轨道产生偏移的快慢程度与其自身稳定性也有一定关系，但是，在同样时间内，以同样步长计算，不同算法发生偏移的快慢反映了算法的精确度.

图3 三种算法积分3个周期得到的轨道Fig.3 The orbits of three algorithms integral three orbit periods

积分100000步时绝对能量误差如图4所示，可以看出当步长取固定值0.001，积分100000步时改进的显式辛算法的能量误差最大为10−9量级.虽然RK78算法的初始能量误差比较小，但随着时间推移，能量误差积累变得越来越大，出现骤增的情况，而RK45算法的能量误差始终很大.与之相比，显式辛算法的能量误差始终保持在一定范围内，能量误差突然增大是因为航天器与主天体P2的距离突然变得非常小.

图4 三个算法积分100000步的绝对能量误差Fig.4 The absolute energy errors of 100000 steps with three algorithms

利用改进的显式辛算法、RK78和RK45算法分别计算由表1所示初始参数得到的Halo轨道对应的一条不变流形，综合考虑积分时间和积分精度，积分步长取0.0001，积分100000步时三种算法得到的稳定流形如图5所示，不稳定流形与其关于y轴对称，3种算法得到的不变流形从形状来看差别无几，但是不重合.从图6所示能量相对误差分析，可知利用改进的算法得到的不变流形，可以保证在长时间积分过程中能量误差在很小范围内，RK78算法初始能量误差比较小，但呈增大趋势，而RK45的能量误差一直都很大.这3种算法所耗时长和相对能量误差的最大值如表2所示，改进的显式辛算法耗时最短，计算效率约是RK78的12.3084倍，且相对能量误差的最大值最小，因此本文改进的显式辛算法在低能轨道转移过程中可以实现轨道拼接点精准能量匹配.

图5 三种算法得到的一条稳定流形轨道Fig.5 Stable manifold of three algorithms

图6 三种算法积分100000步的相对能量误差Fig.6 The relative energy error of 100000 steps with three algorithms

表2 三种算法的积分效率Table 2 Eefficiency of three algorithms

4 结论

在基于流形拼接法设计航天器低能转移轨道的过程中，需要实时进行不变流形的能量计算.本文针对传统数值积分能量耗散问题，研究了显式辛算法在旋转坐标系下的圆型限制性三体问题不变流形的保能量问题.显式辛算法长时间积分时没有长期的变化，能够保持Hamilton系统的辛结构，本文通过改变圆型限制性三体问题分解形式中类动能(动能部分不是严格的正定二次型)Lie算子的形式，证明含有一阶导数项、二阶导数项和三阶导数项的改进显式辛算法可以求解旋转坐标系下的圆型限制性三体问题，这也说明改进的Lie算子形式合理.最后，应用改进的显式辛算法、RK78和RK45算法分别计算了地--月圆型限制性三体问题L1平动点附近的周期轨道和不变流形，并对其能量误差和计算效率进行了分析.仿真结果显示，本文改进的显式辛算法在精度和能量误差两方面都具有较好的优势.

1 Sezebehely V.Theory of Orbits:The Restricted Problem of Three Bodies.New York:Academic Press，1967

2刘林，侯锡云.深空探测器轨道力学.北京：电子工业出版社，2012(Liu Lin，Hou Xiyun.Orbital Mechanics of the Deep Space Probe.Beijing:Publishing House of Electronics Industry，2012(in Chinese))

3雷汉伦.平动点、不变流形及低能轨道.[博士论文].南京：南京大学，2015(Lei Hanlun.Equilibrium point，invariant manifold and low energy trajectory.[PhD Thesis].Nanjing:Nanjing University，2015(in Chinese))

4 Farquhar RW.The role of Sun-Earth collinear libration points in future space exploration//AAS Annual Meeting,1999

5 Dunham DW，Farquhar RW.Libration point missions，1978-2002//Libration Point Orbits and Applications.Singapore:World Scienti fi c,2003

6 Condon GL,Pearson DP.The role of humans in libration point missionswithspeci fi capplicationtoanEarth-Moonlibrationpointgateway station//AAS 01-307，AAS/AIAA Astrodynamics Specialist Coference,2001

7 Lo MW.The inter-planetary superhighway and the origins program//IEEE Aerospace Conference Proceedings，2002

8 Baoyin HX，McInnes CR.Solar sail Halo orbits at the Sun–Earth arti fi cial L1 point.Celestial Mechanics and Dynamical Astronomy,2006，94:155-17

9 Xu M.Application of Hamiltonian structure-preserving control to formation fl ying on a J2-perturbed mean circular orbit.Celestial Mechanics and Dynamical Astronomy,2012,113(4):403-433

10 G´omez G，Jorba A，Masdenmont J，et al.Study of the transfer from the Earth to a halo orbit around equilibrium point L1.Celestial Mechanics，1993，56:541-562

11袁建平,孙冲,方群.基于虚拟中心引力场方法的航天器转移轨道设计.力学学报,2015,47(1):180-184(Yuan Jianping,Sun Chong,Fang Qun.Spacecraft’s transfer orbit design based on the virtual central gravity?eld method.Chinese Journal of Theoretical and Applied Mechanics,2015,47(1):180-184(in Chinese))

12贾建华,吕敬,王琪.带脉冲的三维引力辅助变轨研究.力学学报,2016,48(2):437-446(Jia Jianhua,L¨u Jing,Wang Qi.Powered gravity assist in three dimensions.Chinese Journal of Theoretical and Applied Mechanics,2016,48(2):437-446(in Chinese))

13 Zhang RY，Luo JJ.Attainable sets approach for low-energy，lowthrustInterplanetarytransfers//InternationalAstronauticalCongress.Beijing，China,2013

14 Lei HL，Xu B，Hou XY，et al.High-order solutions around triangular libration points in the elliptic circular restricted three-body problem and applications to low energy transfers.Communications in nonlinear science and Numerical Simulation，2014，19:3374-3398

15 Lei HL，Xu B.High-order analytical solutions around triangular libration points in the circular restricted three-body problem.Monthly Notices of the Royal Astronomical Society，2013，434(2):1376-1386

16 Lei HL，Xu B，Hou XY，et al.High-order solutions of invariant manifolds associated with libration point in the elliptic circular restricted three-body system.Celestial Mechanics and Dynamical Astronomy，2013，117(4):349-384

17 Lei HL，Xu B.Analytiacl study on the motions around equilibrium points of restricted Four-body problem.Compute Nonlinear Science Number Simulat,2015，29:441-458

18 Dellnitz M，Junge O，Post M.On Target for venus:Set oriented computation of energy efficient low thrust trajectories.Celestial Mechanics and Dynamical Astronomy,2006，95:357-370

19 Howell K，Beckman M，Patterson C，et al.Representations of invariant manifolds for applications in three-body system//14th AAS/AIAA Space Flight Mechanics Conference，Maui，Hawaii，2004

20 Hsu CS.A theory of cell-to-cell mapping dynamical systems.Applied Mechanics，1980,147:931-939

21 Ruth RD．A Canonical integration technique.IEEE Transactions on Nuclear Science，1983，30：2669-2671

22李渊,邓子辰,叶学华等.基于辛理论的载流碳纳米管能带分析.力学学报,2016,48(1):135-139(Li Yuan,Deng Zichen,Ye Xuehua,et al.Analysing the wave scattering in single-walled carbon nanotube conveying fl uid based on the symplectic theory.Chinese Journal of Theoretical and Applied Mechanics,2016,48(1):135-139(in Chinese))

23 Pauli Pihajoki.Explicit methods in extended phase space for inseparable Hamiltonian problems.Celestial Mechanics and Dynamical Astronomy,2015，121:211-231

24 Liu L，Wu X,Huang GQ,et al.Higher order explicit symmetric integrators for inseparable forms of coordinates and momenta MNRAS459.MonthlyNoticesoftheRoyalAstronomicalSociety,2016 25 Hu WP，Deng ZC，Han SM，et al.Generalized multi-symplectic Integrators for a class of Hamiltonian nonlinear wave PDEs.Journal of Computational Physics，2013，235:394-406

26 Mei LJ，Wu X，Liu FY.On preference of Yoshida construction over Forest–Ruth fourth-order symplectic algorithm.The Europe Physical Journal C,2013,73:2413

27廖新浩，刘林.辛算法在限制性三体问题数值研究中的应用.计算物理，1995，12(1):102-108(Liao Xinhao，Liu Lin.Applications of symplectic algorithms to the numerical researches of restricted three-body problem.Physic Computation,1995，12(1):102-108(in Chinese))

28 Sun W，Wu X，Hang GQ.Symplectic integrators with potential derivatives to third order. Research in Astronomy and Astrophysics，2011，11:353-368

29 Hairer E，Wanner G.Solving Ordinary Di ff erential Equation.[s.l.]:Springer Verlag，1991

30李言俊，张科，吕梅柏等.利用拉格朗日点的深空探测技术.西安：西北工业大学出版社，2014(Li Yanjun，Zhang Ke，L¨u Meibo.et al.Deep Space Exploration Using Lagrange Equilibrium Point.Xi’an:Northwestern Polytechnical Chivesity Press,2014(in Chinese))

31 Bernelli-Zazzera F，Topputo F，Massari M.Assessment of Mission Design Including Utilization of Libration Points and Weak Stability Boundaries.[s.l.]:Politecnico di Milano，2004

32 G´omez G，Jorba A，Masdemont J，et al.Study re fi nement of semianalytical halo orbit theory.Final Report，ESOC Contract No:8625/89/D/MD(SC)，1991

33 Omelyan IP，Mryglod IM，Folk R.Optimized Forest-Ruth-and Suzuki-like algorithms for integration of motion in many-body systems.Computer Physics Communications,2002，146:188-202

34 Omelyan IP，Mryglod IM，Folk R.Symplectic analytically integrable decomposition algorithms:Classi fi cation，derivation,and application to molecular dynamics,quantum and celestial mechanics simulations.Computer Physics Communications，2003，151:272-314

COMPUTATION OF INVARIANT MANIFOLD BASED ON SYMPLECTIC ALGORITHM OF MIXED LIE OPERATOR1)

Zheng Dandan∗,†Luo Jianjun∗,†,2)Zhang Renyong∗∗Liu Lei†,††∗(School of Astronautics，Northwestern Polytechnical University，Xi’an 710072，China)
†(National Key Laboratory of Aerospace of Flight Dynamics，Xi’an 710072，China)
∗∗(Technology and Engineering Center for Space Utilization，Chinese Academy of Sciences，Beijing 100094，China)
††(Beijing Aerospace Control Center，Beijing 100094，China)

Invariant manifolds of periodic orbit near the libration points attract a lot of attentions due to their importance in the low-energy orbits transfer problem.In the process of low-energy orbit design,the energy of the invariant manifolds must be matched,but the energy is dissipated when integrating with traditional numerical integration method.The explicit symplectic algorithm with energy conservation is more efficient than the implicit symplectic algorithm,but it requires the Hamiltonian system to be divided into two integral parts,while the circular restricted three-body problem in the rotating coordinate system being inseparable.It is difficult to solve the circular restricted three-body problem in the rotating coordinate system by explicit symplectic algorithm.In this paper,the mixed Lie derivative operator of kinetic energy is used to solve the circular restricted three-body problem in the rotating coordinate system,and the e ff ectiveness of this explicit symplectic algorithm with the third derivation in dealing with this problem has been showed.Compared with the Runge-Kutta45 method and Runge-Kutta78 method,the symplectic algorithm with the third-order derivative term not only has high precision but also the smallest energy error and the highest efficiency.Finally,the invariant manifolds are calculated by the symplectic algorithm with the third derivative term,the patched point can match accurately with this method.

circular restricted three-body problem，invariant manifold，symplectic algorithm，mixed Lie operator

V412.4

A

10.6052/0459-1879-17-079

2017–01–11收稿，2017–07–18 录用,2017–07–27 网络版发表.

1)国家自然科学基金重大项目(61690210，61690211)和航天飞行动力学重点实验室开放基金项目(2015afdl014,2015afd1017)资助.

2)罗建军，教授，主要研究方向：航天飞行动力学与控制.E-mail:jjluo@nwpu.edu.cn

郑丹丹,罗建军,张仁勇,刘磊.基于混合Lie算子辛算法的不变流形计算.力学学报,2017,49(5):1126-1134

Zheng Dandan,Luo Jianjun,Zhang Renyong,Liu Lei.Computation of invariant manifold based on symplectic algorithm of mixed Lie operator.Chinese Journal of Theoretical and Applied Mechanics,2017,49(5):1126-1134

附录A

附录B