An improved numerical method for constructing Halo/Lissajous orbits in a full solar system model

Yingjing QIAN,Xiodong YANG,*,Wuxing JING,Wei ZHANG

aBeijing Key Laboratory of Nonlinear Vibrations and Strength of Mechanical Structures College of Mechanical Engineering,Beijing University of Technology,Beijing 100124,China

bDepartment of Aerospace Engineering,Harbin Institute of Technology,Harbin 150001,China

1.Introduction

Libration points are the five equilibrium solutions in the Circular Restricted Three-Body Problem(CRTBP)founded by Euler and Lagrange.Three of the five equilibrium points,L1,L2,and L3,are called collinear libration points,and the remaining two,L4and L5,are triangular libration points.Studies about probes moving around orbits in the vicinity of the Earth-Moon collinear libration points have shown potential of being the communication spot for exploring the far side of the Moon,1being established as a space hub for the exploration of asteroids and serving,for instance,as a construction and repair facility.2

The linearized motions around collinear libration points in the CRTBP are typical center× center×saddle type.There exist center manifolds,stable and unstable,arriving at or departing from the corresponding Lyapunov orbit exponentially with time.3Hence,the bounded orbits in the vicinity of the collinear Earth-Moon libration points are very sensitive to any variation of initial conditions or external perturbations.

In literature,systems like the Sun-Earth-probe and Sun-Jupiter-probe systems are always considered as circular restricted three-body problems due to the small eccentricity of the secondary primary and neglectable gravitational influences from other planets.However,whether the real Earth-Moon-probe system can be modeled as simple as the CRTBP model depends on the problems of concern.Specifically,Qi and Xu4pointed out that,outside the sphere of influence of the Earth-Moon system(the radius is about 579734.2 km),the real Earth-Moon-probe system can be modeled as the CRTBP model with the Sun and the Earth-Moon barycenter as two primaries,while inside the sphere of influence of the Earth-Moon system,the eccentricity of the Moon,and other planets’direct gravitational influences on the probe and indirect gravitational influences on the Earth-Moon system both cannot be neglected.5Therefore,the libration mission trajectory design in the Earth-Moon system requires a high- fidelity model and an improved understanding of external perturbations.It is worth to mention that exact periodic orbits in the CRTBP do not exist in the vicinity of the Earth-Moon collinear libration points.Instead,perturbed Halo/Lissajous orbits may dominate considering the two center manifolds and strong external perturbations.6–10This research effort focuses on the development of an efficient construction method to best satisfy future mission requirements in these complex dynamical regimes.

The construction method for Halo/Lissajous orbits was firstly investigated by perturbation methods.11Although exact solutions cannot be derived for the CRTBP,perturbation methods are available to construct approximate analytical solutions of the motions around the libration points,which can provide deep understanding of the dynamics and can present explicit parameter contribution from the global perspective.12The mostrepresentative work was successfully demonstrated by Richardson13and Farquhar and Kamel,14who derived the third-order and fourth-order approximate solutions of halo orbits around the collinear libration points in the CRTBP with the Lindstedt-Poincaré(L-P)method.Masdemont15introduced the solutions of Lissajous orbits and quasi-halo orbits by considering perturbations.

Differential correction methods have been introduced into the construction of periodic orbits in the RTBP since the 1980s.Based on the differential correction idea,a multiple shooting method was proposed by Howell and Pernicka.16The multiple shooting method demonstrates a powerful tool,in which an initially-guessed trajectory is discretized into n‘patch points’connected by n–1 trajectory arcs.Newton’s method is incorporated in the differential correction process and converges to a solution that is continuous in position and velocity at all internal patch points and satisfies any additional problem-dependent end constraints.17Later,a Fourier series correction method was proposed by Jorba18and Gomez and Mondelo.19By overcoming the time-consuming problem,Ren and Shan20presented a generic numerical algorithm that can generate long-term libration point orbits in the CRTBP.Lian and Gomez21combined the multiple shooting method and a re fined Fourier analysis to study the dynamics of a massless particle around the collinear libration points of the Earth-Moon system.Zhang and Li22studied numerical techniques to accurately compute and robustly extend the libration point orbits.

Based on the combination of quantitative analytical solutions and numerical methods, final periodic orbits can be constructed systematically.To be detailed,the aforementioned analytical solutions from Richardson13and Farquhar and Kamel14have been implemented as the guessed-initial value.Numerical methods which utilize the desired symmetry of the trajectory or the patch points as the terminal constraints were applied to derive the whole trajectory.This procedure has been successfully served as the design technique for most of the libration mission trajectories in the Sun-Earth system,such as ISEE-3,GENESIS,GAIA and PLANCK,23and binary asteroid systems.24However,when the libration mission in the Earth-Moon system in the full solar system model is considered,the eccentricity of the Moon and the influence from the Solar system are not neglectable.Firstly,the Sun always lies in the same side of the Earth-Moon plane during one Halo/Lissajous revolution,and its gravitational force makes the trajectory certainly asymmetry.The orbit loses its periodicity under the influences of the Sun and the eccentricity of the Moon.Secondly,the distance between the Earth and libration points is no longer a constant but a time-dependent trajectory moving along the line connecting the Earth and the Moon.For large time-spans,the convergence of multiple shooting strongly depends on the accuracy of the initial conditions of the patch points.Therefore,the guessed-initial condition from approximation of the CRTBP may not result in a converged solution in the high- fidelity model.Obtaining more accurate patch points becomes the major concern in constructing Halo/Lissajous orbits in the full solar system model.

The Poincarésection technique offers a good way to identify Halo and Lissajous features,which presents an alternative option to determine the initial conditions of the patch points.Based on the Poincarésections technique,Smith and Szebehely25and Winter26studied a family of simply symmetric retrograde periodic orbits around the Moon.The investigations about the dynamical behavior with Poincarésections of the bounded orbits in the vicinity of the Earth-Moon libration points motivate the new idea for an improved design method.

In this paper,we make contributions to obtain Halo/Lissajous trajectories in the Earth-Moon system around the collinear libration points.An ephemeris model involving all the solar system’s gravitations and the Moon’s influence on the rotating coordinate system is proposed as well as the angular velocity of the rotating coordinate system relative to the inertial coordinate system.Patch points in the multiple shooting method are obtained based on a dynamical analysis with Poincarésections.The methodology of constructing the quasiperiodic orbits are presented in Section 4.Results and some examples are given in Section 5,and conclusions are drawn in Section 6.

2.Equation of motion

The motions in the vicinity of libration points shown in Fig.1 are defined in the synodical frames.In the full solar system gravitational model,a spacecraft is considered ‘massless’,and the geometry of this model is conveniently described in a Geocentric Rotating Coordinate(GRC)system Oxyz,centered at the center of the Earth O.The x axis direction is parallel to the Earth-Moon line,and the positive z axis is also oriented perpendicular to the orbital plane of the Earth-Moon system and parallel to the Moon’s instantaneous orbital angular momentum.The y axis then completes the righthanded triad.In Fig.1,EM-Lidenote the Earth-Moon system Lipoint.

Fig.1 Coordinate systems.

An inertial reference frame,J2000,is represented by three orthogonal unit vectors ex,ey.ezof X,Y,and Z axis.The X axis points to the vernal equinox at noon on January 1st,2000,and the Z axis points to the North Pole.The Y axis then completes the right-handed coordinate system.Generally,standard ephemerides are documented in this coordinate system.

To investigate the motion around the libration points,it is necessary to move the origin of the coordinate system to the libration points,namely,Li(i=1,2,3)Synodic Coordinate(LiSC)systems.Li-xi,Li-yi,and Li-ziare parallel to Ox,Oy,and Oz,respectively.It is worth to mention that the coordinate is a time-dependent coordinate system since the Lipoint is not an inertial point in the Earth-Moon system.

The equation of motion in the vicinity of the Earth-Moon libration points in the J2000 coordinate system is known as

According to the relative differential relationship of the moving vector between the inertial coordinate system and the rotating coordinate system,the following relation is obtained:

where ω denotes the angular velocity of the GRC with respect to the J2000,and,,andare the position,velocity,and acceleration from the center of mass of the Earth to the probe in the GRC,respectively.

Considering Eqs.(1)and(3),the equation of motion for a probe in the vicinity of the Earth-Moon libration points in the GRC is derived as

where μe, μm,and μsiare the gravitational constants for the Earth,Moon,Sun,and other planets including Mercury,Venus,Mars,Jupiter,Saturn,Uranus,and Neptune in the solar system.rpis the position vector from the mass center of the Earth to the probe,and rmand rsidenote the Moon’s position vector,the Sun’s position vector,and other planets’position vectors relative to the center of the Earth,respectively.rm,rp,rsidenote the scalar form of rm,rpand rsi

where

Clearly,Gedenotes the gravitational forces from the Earth.Gmdenotes the gravitational forces from the Moon.Gsidenotes the gravitational influence from the Sun and the other 7 planets.Note that the gravitational influence of the Sun and other 7 planets also acts on the motion of the Earth-Moon system,which results in a variation of the angular velocity ω.

According to the definition of the GRC,it is noticed that ω is actually the orbital angular velocity of the Moon’s motion relative to the Earth.In order to derive the expression for ω,the equation of motion for the Moon is firstly analyzed as follows:

in which rmsi= ‖rm-rsi‖ .

According to the definition of the GRC,the unit vectors i,j,and k of x axis,y axis,and z axis directions are expressed as27

where the Moon’s orbital angular momentum

and hmdenotes the value of hm.

In addition,the unit vectors i,j,and k are not fixed in space,but are continuously changing directions;therefore,their time derivatives are not zero.They obviously have a constant magnitude(unity),and being attached to the xyz frame,they all have the same angular velocity ω.Hence,

and

in which i× (ω × i)is a triple vector product,‘×” denotes the cross product for vectors, ‘·” denotes the dot product for vectors,and （ω ·i）i is dyad.The multiplication symbols in j× (ω × j)and k × (ω × k)are the same as those in i× (ω × i).The angular velocity can be expressed in the GRC as

Considering Eqs.(12)and(13),the following relation is obtained:

Then,substituting Eq.(9)into Eq.(14),the expression for ω is derived as

The angular acceleration can be derived by the first derivative of Eq.(15)as

where

and

Then,substituting the expression of ω into Eq.(4),the equation of motion for the probe in the vicinity of the Earth-Moon libration points in the GRC is derived as

It is noticed that most ephemeris models involving all the solar system’s gravitations in Refs.17,28are presented in the inertial frames,while the model proposed in this paper is a full solar system gravitational model in the GRC with a clear presentation of the angular velocity of the GRC relative to the inertial coordinate system.In this model,the Sun,Moon,and other 7 planets’orbital parameters,such as rsi,rm,and hm,can be directly calculated with the ephemerides(DE421 ephemerides).

If the direct/indirect influence of the Sun and other 7 planets and the Moon’s orbital eccentricity are neglected,i.e.,μsi=0,ωx=0,˙ωx=0,and˙ωz=0,Eq.(18)recovers to the traditional circular restricted three-body problem model.

3.Dynamical behavior analysis with Poincarésections

In this section,Poincarésections are used to identify Halo and Lissajous orbits in the full solar system model,which provides an alternative way to generate patch points for multiple shooting.

The use of a Poincarémap allows a continuous-time system of n dimensions to be reduced to a discrete-time system of n–1 dimensions,which may aid with the visualization of higherdimensional systems.In the CRTBP model,a Poincarémap is usually chosen with a certain Jacobi constant.As examples,typical Poincarésections for Halo and Lissajous orbits in LiSC are shown in Fig.2.Per the Kolmogorov-Arnold-Moser(KAM)theory,the point represents a Halo orbit in the rotating frame(see Fig.2(a)),and the closed curves around the point correspond to Lissajous orbits(see Fig.2(b)).29However,when the Sun’s motion is considered,the relative relations of the Sun,Earth,and Moon could provide various potential energy.30Even starting from the same initial position and velocity,the total energy still depends on the epoch time.There is no Jacobi constant in the full solar system gravitational model.Based on this fact,we choose the same epoch time for mapping Poincarésections.Two types of Poincarésections are defined.The x-y plane in the GRC is defined as the first Poincarésection to analyze the dynamical behavior of the Halo/Lissajous orbits in the GRC.The xi-yiplanes in the LiSC are defined as the second-type section.

Five types of perturbed periodic orbits are studied in this investigation.Since the linearized in-plane frequency of the bounded orbits in the vicinity of the L3point is very close to the linearized out-of-plane frequency,only a Halo orbit is assumed in the investigation.Taking the Earth-Moon system for example,the in-plane non-dimensional linear frequency is 1.01041 and the out-of-plane non-dimensional linear frequency is 1.00533 with dimensionless time and space units.31Hence,the bounded orbits in the vicinity of the L3point always present the Halo-like type.For the L1and L2points,the difference between linearized in-plane and out-of-plane frequencies has potential to present a Lissajous or Halo orbit in a linearized form.Therefore,both perturbed Halo and Lissajous orbits are considered in the study of L1and L2points.

Fig.2 Typical Poincarésections for Halo and Lissajous orbits in LiSC for CRTBP.

Due to the existence of the stable/unstable manifold,it is not possible to locate the exact bounded orbits.A small deviation may cause a path to depart the bounded orbit after only a few revolutions.However,a few revolutions are enough for revealing the feature for the orbits about the Earth-Moon collinear libration points.For this purpose, five groups of initial conditions for Earth-Moon perturbed Halo and Lissajous orbits are adopted from STK/Astrogator and integrated with Eq.(18)for about 3–4 revolutions.The initial conditions are listed in Table 1,the epoch time is 1 Oct 2022,12:00:00.000.The Poincarésections about perturbed Halo and Lissajous orbits in the GRC and the LiSC are shown on Figs.3–7.It is noticed that the solid lines are the projections of the perturbed Halo and Lissajous orbits on the x-y plane,and the red stars are the intersection points of trajectories with the x-y plane,which are the points on the Poincarésections.

The Poincarésections of perturbed halo orbits in the LiSC coordinate system,shown in Figs.3(b),5(b),and 7(b),are still similar to the typical periodic orbit type as presented in Fig.2(a).Although the points on a section are not overlapping,they still locate closely to each another and are symmetrically distributed abouttheyi-axis,which demonstrate approximate periodic features.However,the Poincarésections of the same perturbed Halo orbits in the GRC coordinate system presented in Figs.3(a),5(a),and 7(a)lose their features of periodicity.The major factor that accounts for the huge different phenomenon in the GRC and the LiSC is the movements of the libration points.In the LiSC,the origin of the coordinate system is the libration points which are no longer inertial points but time-dependent points changing along with the radius of the Moon.

On the other hand,the Poincarésections about the perturbed Lissajous orbits in the LiSC coordinate system in Figs.4(b)and 6(b)show points located along the closed curves around the libration points,which are also similar to the Poincarésection presented in Fig.2(b).However,irregular motions dominate the features of Figs.4(a)and 6(a)as the points on the Poincarésections are dispersively distributed.Similarly,the main reason for Lissajous orbits losing their periodicity in the x-y plane is due to the movements of the libration points.The libration points in the full solar system gravitational model are no longer inertial points but time-dependent points moving along the line connecting the Earth and the Moon.

Based on the phenomenon revealed in the Poincarésection analysis,a novel approach to locate more accurate conditions is proposed.Obviously,the new patch points consist of two parts,the inherited periodic part and the perturbed part.The inherited periodic part can be adopted from the orbits in the CRTBP model.When orbits are expressed in the Lisynodic coordinate systems,only the periodic part shows,since the origins of the Lisynodic coordinate systems represent the movements of the libration points.The perturbed part can be achieved by solving the instantaneous state of the Li(i=1,2,3)point in the GRC.The instantaneous state of the Lipoint contains perturbation results from the Moon eccentricity and the gravitational influence from the solar system.The aforementioned analysis converts complicated determinations of the initial conditions problem into obtaining the instantaneous state of Liand the basic initial conditions in the CRTBP model.Therefore,the problem to determine the initial conditions problem can be solved as the instantaneous state of Liis obtained.In the next section,the methodology to deriveinitial conditions by the superposition in detail and its application in orbit computations are described.

Table 1 Initial conditions for Earth-Moon perturbed Halo/Lissajous orbits in the GRC.

Fig.3 Poincarésections of Halo orbit in GRC and L1SC.

Fig.4 Poincarésections of Lissajous orbit in GRC and L1SC.

Fig.5 Poincarésections of Halo orbit in GRC and L2SC.

4.Numerical procedure

The accuracy of solutions in the full solar system gravitational model will be determined by appropriate initial conditions.Inspired by the dynamical behaviors analysis of Poincarésections,more accurate conditions for patch points are obtained and employed.Then,the multiple-shooting method is used to derive perturbed Halo/Lissajous trajectories.The general approach is presented as follows:

Step 1.Determination of the initial condition for each trajectory arc in the CRTBP model

Fig.6 Poincarésections of Lissajous orbit in GRC and L2SC.

Fig.7 Poincarésections of Halo orbit in GRC and L3SC.

The trajectory generated from the CRTBP model in LiSC(i=1,2,3)is discretized into m ‘patch points’connected by m–1 trajectory arcs.Each patch point is represented by a six-dimensional vector,xk(k=1,2,...,m),the integration time associated with each arc is τk(k=1,2,...,m),and the time interval between each point is Tk-1= τk-τk-1(k=2,3,...,m).Then,the variable vector~X is comprised as~X=[x1,x2,...,xm,T1,T2,...,Tm-1,τ1,τ2,...,τm].

We locate those patch points by the phase angle η which is defined in the x-y plane since halo and Lissajous orbits are all periodic in the x-y plane in the CRTBP model.η is a nongeometric angle describing the positions of a spacecraft in the orbits in a manner similar to mean anomaly.A phase angle of zero occurs when the spacecraft intersects the x-z plane and travels in the positive y direction.The phase angle increases in the direction of the orbital motion and is defined as

where p is the x-y period for halo and Lissajous orbits.Fig.8 shows some values of η which are not evenly spaced.

Step 2.Determination of the instantaneous state of the

libration points

Geometrically,the libration points belong to the instantaneous plane of motion of the Moon around the Earth so that the distances to the Earth and to the Moon ful fil the same relationships as in the CRTBP.The effects of the remaining solar system bodies,as well as the non-circular motion of the Moon around the Earth,prevent these points to be relative equilibrium ones.In the CRTBP,the relation parameter λi(i=1,2,3)between the position Li(i=1,2,3)and the Moon are parameters only relevant to the mass parameter.Due to the influence of the Moon’s eccentricity and the gravitation from the solar system,λibecome slowly time-varying parameters.32

However,variations of λiwith respect to the initial epoch and the time-span are very small that can still be considered as a constant for each short trajectory arc.Meanwhile,the relation between the position Liand the Moon is assumed as3

where Rmis the position vector from the center of the Earth to that of the Moon in the GRC.

Since the libration points are stationary solutions for the equation of motion such as Eq.(18),with consideration of the relation Eq.(20),the solutions should subject to RLi=[λirm(t0),0,0]and˙RLi=[λi˙rm(t0),0,0]Tin the GRC.The Newton iteration process is incorporated to determine the value of λifor each short trajectory arc.The initial-guessed value of λi0is chosen from the CRTBP.

Fig.8 Positions of various values of phase angle η,on a Halo/Lissajous orbit in the Earth-Moon system projected onto x-y plane.

For the kth arc,since the initial-guessed solutions of the libration points[λi0rm(t0),0,0,λi0˙rm(t0),0,0]Tare not accurate enough,they cannot maintain stationary after being integrated for a time interval τk(k=1,2,...,m).P is denoted as the solution after integration.The following equation should be minimized to ensure solutions subjecting to[λirm(t0),0,0]T:

where Tk–1is the initial time for the kth arc,and τkdenotes the time interval for the kth arc.Therefore,the updated λi,j+1can be obtained as

in which Δλi,jis a small value for the iterative least squares algorithm.Generally,over six times of iterations could offer satisfying results.

Step 3.Determination of initial condition for each trajectory arc with consideration of the instantaneous state of LiAfter choosing the epoch of time and the patch points from Step 1,according to the corresponding time-dependent relationship between the LiSC and the GRC,each patch point can be transformed into the corresponding point in the GRC.The exact initial conditions can be obtained by the superposition of the instantaneous state of Liand the inherited initial conditions in the CRTBP model with fixed Lipoints.

The schematic diagram of transformation is shown in Fig.9,and the following relationship is obtained:

where Rkdenotes the position vector of the kth patch point in the GRC,RLidenotes the position vector from the center of the Earth to Liin the GRC,and rkdenotes the position vector in the GRC from Lito the probe.

Differentiating Eq.(23)with respect to time yields

wheredenotes the kth patch point’s information xkin Step 1.ωGCis the angular velocity vector of the GRC with respect to the LiSC.Obviously,according to the definitions of these two coordinate systems,we have

Substituting Eqs.(25)and (20)into Eq.(24),and considering

Fig.9 Schematic diagrams of geometrical relationships for LiSC and GRC.

in whichis the transformation matrix from the J2000 to the GRC and can be expressed as Eq.(9),the following relationship is obtained:

Herein,rmand˙rmdenote the position and velocity vectors of the Moon in the J2000 coordinate system,respectively.Since the dot product of two vectors is independent of the coordinate system,the dot product of˙Rmand Rmcan be computed directly based on ephemeris data in the J2000.

Step 4.Connecting patch points

Here we assume that the time along each segment is not constant and the constraints should satisfy that the updated set of states Φk(xk)=xk+1for k=1,2,....,m-1.Besides,for most missions,the final epoch is given as τm.Therefore,the constraint vector F（）is written as

Sequential quadratic programming is used to obtain a solution that is continuous in position and velocity at all internal patch points,as shown in Fig.10.

5.Numerical simulations

Based on the discussions presented in Sections 2–4,it is concluded that we made contributions in two aspects.(A)We proposed a full solar system gravitational model in the GRC with a clear presentation of the angular velocity of the GRC relative to the inertial coordinate system.In this model,the Sun,Moon,and other 7 planets’orbital parameters rsi,rm,and hmcan be directly calculated with ephemerides,such as the DE421 ephemerides.(B)We provided an alternative way to search the patch points in the multiple shooting method based on the dynamical analysis with Poincarésections.Finally,those patch points were used to construct quasi-periodic orbits in the vicinity of the libration points with the proposed model.With the contributions above,we are able to improve the numerical method for constructing quasi-periodic orbits,which is listed as follows:

Fig.10 Schematic diagrams for improved multiple-shooting method to govern perturbed Halo/Lissajous trajectories.

(1)To simplify the computational process:In the literature,the multiple shooting method is always implemented with an ephemerides model in the inertial coordinate system,such as Eq.(1)in this manuscript.Patch points are firstly selected from a periodic orbit in the CRTBP model.Then,the patch points need to be transformed to the inertial coordinate system according to ephemeris data.During the coordinate transformation,the rotational angular velocity is obtained based on the twobody assumption for the Earth and the Moon.This assumption would result in rounding errors.In this study,we establish the dynamic model in the GRC.In the angular velocity and angular acceleration that we present in Eqs.(15)–(16),the Sun,Moon,and other 7 planets’orbital parameters rsi,rm,and hmcan be directly calculated with ephemerides,such as the DE421 ephemerides.We implement the multiple shooting method directly with the solar system gravitational model in the GRC without a two-body assumption of the angular velocity.

(2)To locate the patch points:In the literature,the shooting process is started with roughly guessed patch points from the CRTBP,and then the patch points are with re fined trigonometric approximations during the iteration.We locate patch points differently.Through the dynamical behaviors analysis with Poincarésections,we can conclude that the main reason for halo and Lissajous orbits losing their periodicity in the x-y plane is due to the movements of the libration points.Patch points can be described by the superposition of both the periodic and perturbed parts.Based on the Poincarésection analysis,we propose a new approach to locate patch points for the numerical procedure.Simulations verify the proposed method.Therefore,we start the shooting process with more accurate patch points from an ephemerides model and then update the patch points with sequential quadratic programming.Therefore,we improve the method by simplifying the computational process and locating patch points for the numerical procedure.

Fig.11 Halo and Lissajous orbits in CRTBP(viewed from L1SC).

Fig.12 Halo and Lissajous orbits in CRTBP(viewed from L2SC).

In this section,perturbed Halo/Lissajous orbits for the L1/L2points and a perturbed Halo orbit for the L3point in the Earth-Moon system are generated as examples to illustrate the proposed method.The initial epoch time is 12:00 Oct.1,2030(Greenwich Mean Time).

Firstly,Halo and Lissajous trajectories are generated from the CRTBP model in LiSC(i=1,2)for the L1/L2point case as well as a Halo orbit for the L3point in L3SC,which are shown in Figs.11–13.The parameters for those orbits are listed in Table 2.

Fig.13 Halo orbit in CRTBP(viewed from L3SC).

In order to describe patch points explicitly,the state xiwhere the phase angle η is 45°,135°,225°,or 315°for each revolution are chosen as the initial conditions for patch points.(Based on our simulation experience,four patch points for one revolution are enough for a convergence.)Besides,the initial state x0and the final state xmare also chosen for patch points.For instance,in the L2Halo case,the orbit goes through the pre-set phase angles(45°,135°,225°,or 315°)28 times.The initial state and the final state are also selected as patch points.Therefore,the whole trajectory is discretized into 30 ‘patch points’connected by 29 trajectory arcs.The total number of patch points for all cases are listed in Table 2 as well.Since the x-y frequency in the L3point case is only half of the x-y frequency in the L2point case,during the same time interval,only 14 ‘patch points’are chosen.

So far,the initial conditions for all patch points in the CRTBP model are determined.

Secondly,based on the method presented in Section 4 Step 2,each λicorresponding to the trajectory arcs for all five cases are obtained.The initial-guessed value of λi0is chosen from the CRTBP as λ10=0.849065766935798,λ20=1.16783268238542,and λ30=0.99291206846683 in Eq.(22).Then,according to Eq.(27),new patch points with consideration of the Moon’s real motion are obtained.Hence,the time-dependent positions of Liare determined by Eq.(20)through values of λi0.

Table 2 Parameters for five types of periodic orbits in LiSC.

Fig.14 Patched L1perturbed Halo orbit.

Fig.15 Patched L1perturbed Lissajous orbit.

Lastly,sequential quadratic programming is used to obtain a solution that is continuous in position and velocity at all internal patch points,which are shown in Figs.14-18 for five types of trajectories.

6.Conclusions

Fig.16 Patched L2perturbed Halo orbit.

Fig.17 Patched L2perturbed Lissajous orbit.

Fig.18 Patched L3perturbed Halo orbit.

Perturbed Halo/Lissajous orbits in the Earth-Moon system are investigated in this paper by an improved numerical method.Since perturbations such as the Moon’s eccentricity and the solar system gravity may cause unfavorable excursions of a probe,a high- fidelity dynamical model is developed.Through a dynamical behaviors analysis with Poincarésections,it is concluded that the main factor leading to aperiodicity of Halo and Lissajous orbits in the x-y plane is the movements of the libration points which change along with the radius of the Moon.To obtain more accurate bounded orbits,better initial conditions have been studied in the numerical integrations.Initial conditions are described by the superposition of both the periodic and perturbed parts.Inspired by the Poincarésection analysis,an improved multiple-shooting method that utilizes the initial conditions with consideration of the instantaneous state of libration points for each integral arc is proposed.By applying the proposed model and the improved multipleshooting method, five types of perturbed Halo/Lissajous orbits are obtained.Furthermore,in a real mission design process,the solar radiation must be considered since it is one of the contributing sources of perturbation.However,obtaining the solar radiation requires the geometrical size and attitude information of a probe,which is beyond the purpose of the paper.We believe that the method proposed in this study can be extended to a case with consideration of the solar radiation.

Acknowledgements

The authors gratefully acknowledge the supports of the National Natural Science Foundation of China (Nos.11772009,11402007 and 11672007),and the Funding Project for Academic Human Resources Development in Institutions of Higher Learning under the Jurisdiction of Beijing Municipality.

1.Farquhar RW.The control and use of libration point satellites[dissertation].Stanford:Stanford University;1969.

2.Kakoi M,Howell KC,Folta D.Access to Mars from Earth-Moon libration point orbits:manifold and direct options.Acta Astronaut 2014;102:269–86.

3.Szebehely V.Theory of orbits—The restricted problem of three bodies.Salt Lake:Academic Press;1976.p.5–10.

4.Qi Y,Xu SJ.Study of lunar gravity assist orbits in the restricted four-body problem.Celest Mech Dyn Astr 2016;125(3):333–61.

5.Heritier A,Howell KC.Dynamical evolution of natural formations in libration point orbits in a multi-body regime.Acta Astronautica 2014;102:332–40.

6.Meng YH,Zhang YD,Dai JH.Floquet-based design and control approach to spacecraft formation flying in libration point orbits.Sci China Technol Sci 2011;54(3):758–66.

7.Qian YJ,Liu Y,Zhang W,Yang XD,Yao MH.Stationkeeping strategy for quasi-periodic orbit around Earth-Moon L-2 point.Proc Inst Mech Eng G-J Aerospace Eng 2016;230(4):760–75.

8.Peng HJ,Chen BS,Wu ZG.Multi-objective transfer to librationpoint orbits via the mixed low-thrust and invariant-manifold approach.Nonlinear Dyn 2014;77(1–2):321–38.

9.Zeng H,Zhang JR.Design of impulsive Earth-Moon Halo transfers:lunar proximity and direct options.Astrophys Space Sci 2016;361(10):1–10.

10.Zeng H,Zhang JR,Qi R,Li MT.Study of time-free transfers into libration point orbits with multiple constraints.J Guidance Control Dyn 2017;40(11):2752–70.

11.Qian YJ,Yang XD,Zhai GQ,Zhang W.Analytical and numerical construction of vertical periodic orbits about triangular libration points based on polynomial expansion relations among directions.Astrophys Space Sci 2017;362(8):1–14.

12.Zhang JR,Zhao SG,Yang YZ.Characteristic analysis for elliptical orbit hovering based on relative dynamics.IEEE Trans Aerospace Electron Syst 2013;49(4):2742–50.

13.Richardson DL.Analytic construction of periodic-orbits about the collinear points.Celestial Mech 1980;22(3):241–53.

14.Farquhar RW,Kamel AA.Quasi-periodic orbit about the translunar libration.Celestial Mech 1972;7:458–73.

15.Masdemont JJ.High-order expansions of invariant manifolds of libration point orbits with applications to mission design.Dyn Syst 2005;20(1):59–113.

16.Howell KC,Pernicka HJ.Numerical determination of Lissajous trajectories in the restricted 3-body problem.Celestial Mech 1987;41(1–4):107–24.

17.Pavlak TA,Howell KC.Evolution of the out-of-plane amplitude for quasi-periodic trajectories in the Earth-Moon system.Acta Astronautica 2012;81(2):456–65.

18.Jorba A.Numerical computation of the normal behaviour of invariant curves of n-dimensional maps.Nonlinearity 2001;14(5):943–76.

19.Gomez G,Mondelo JM.The dynamics around the collinear equilibrium points of the RTBP.Physica D 2001;157(4):283–321.

20.Ren Y,Shan JJ.A novel algorithm for generating libration point orbits about the collinear points.Celest Mech Dyn Astr 2014;120(1):57–75.

21.Lian YJ,Gomez G,Masdemont JJ,Tang GJ.A note on the dynamics around the Lagrange collinear points of the Earth-Moon system in a complete solar system model.Celest Mech Dyn Astr 2013;115(2):185–211.

22.Zhang HQ,Li S.A general method for the generation and extension of collinear libration point orbits.Celest Mech Dyn Astr 2016;126(4):339–67.

23.Hechler M,Cobos J.HERSCHEL,PLANCK and GAIA orbit design.Proceedings of the conference on libration point orbits and applications;2002 Jun 10–14;Aiguablava,Spain.2003.p.115–35.

24.Bu SC,Li S,Yang HW.Artificial equilibrium points in binary asteroid systems with continuous low-thrust.Astrophys Space Sci 2017;362(8):1–14.

25.Smith RH,Szebehely V.The onset of chaotic motion in the restricted problem of three bodies.Celest Mech Dyn Astr 1993;56(3):409–25.

26.Winter OC.The stability evolution of a family of simply periodic lunar orbits.Planet Space Sci 2000;48(1):23–8.

27.Jing WX,Qian YJ,Gao CS,Li JQ.Autonomous navigation to quasi-periodic orbits near translunar libration points.Chin J Aeronaut 2013;26(5):1259–68.

28.Haapala AF,Howell KC.Representations of higher-dimensional Poincare maps with applications to spacecraft trajectory design.Acta Astronaut 2014;96:23–41.

29.Dutt P,Sharma RK.Analysis of periodic and quasi-periodic orbits in the Earth-Moon system.J Guid Control Dyn 2010;33(3):1010–7.

30.Qian YJ,Zhang W,Yang XD,Yao MH.Energy analysis and trajectory design for low-energy escaping orbit in Earth-Moon system.Nonlinear Dyn 2016;85(1):463–78.

31.Hou XY,Liu L.Bifurcating families around collinear libration points.Celest Mech Dyn Astr 2013;16(3):241–63.

32.Erdi B.3-Dimensional motion of Trojan asteroids.Celestial Mech 1978;18(2):141–61.