# Semianalytical Solutions of Relative Motions with Applications to Periodic Orbits about a Nominal Circular Orbit

Authors: Qiwei Guo; Hanlun Lei; Bo Xu
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

In the dynamical model of relative motion with circular reference orbit, the equilibrium points are distributed on the circle where the leader spacecraft is located. In this work, analytical solutions of periodic configurations around an arbitrary equilibrium point are constructed by taking Lindstedt-Poincaré (L-P) and polynomial expansion methods. Based on L-P approach, periodic motions are expanded as formal series of in-plane and out-of-plane amplitudes. According to the method of polynomial expansions, a pair of modal coordinates is chosen, and the remaining state variables are expressed as polynomial series about the modal coordinates. In order to check the validity of series solutions constructed, the practical convergence is evaluated. Considering the fact that relative motion model is a special case of restricted three-body problem, the periodic configurations constructed in the model of relative motion are taken as starting solutions to numerically identify the periodic orbits in restricted three-body problem by means of continuation technique with the mass of system as continuation parameter.


---

## Body

## 1. Introduction

Studying the relative motion among satellites is of significance for configuration design of formation flying, constellation, and distributed spacecraft systems, where the mode of working in cooperation with several probes is required in many space missions.In the relative motion problem, the commonly adopted model is the Hill-Clohessy-Wiltshire (HCW) equations [1], which are used to describe the relative motion of two point-mass bodies under the gravitational influence of a central body. Later, an important contribution was made by Lawden and Tschauner et al. [2–4], who generalized the HCW equations with an elliptic reference orbit. The results in the case of elliptic reference orbit were studied and extended by many researchers [5–9]. Some researchers studied the relative motions, starting from the solutions of HCW equations and latter moving to consider some perturbations, such as oblateness [10–15], drag [16, 17], and third-body effects [18].In recent decades, Lindstedt-Poincaré (L-P) method has been applied to the problem of relative motion. High-order solutions have been constructed by several authors [19–21]. In particular, Gómez and Marcote [22] took the L-P method to obtain series solution of periodic configurations up to an arbitrary order about in-plane and out-of-plane amplitudes, and Ren et al. [23] presented the third-order expression for the solutions of elliptic relative problem.In addition, Shaw et al. [24–26] proposed an invariant manifold approach based on polynomial expansions to perform nonlinear modal analysis. Recently, Qian et al. [27] provided the third-order polynomial relations among the three motion directions for vertical periodic orbits around triangular equilibrium points in the circular restricted three-body problem based on the polynomial expansion method.To the authors’ knowledge, the periodic motions around an arbitrary equilibrium point in the relative motion model have not been studied. In this work, we aim to investigate analytical expressions for them by taking both the L-P and polynomial expansion approaches. As an application, the periodic objects obtained in the relative motion model are continued to periodic orbits in the circular restricted three-body problem (CRTBP) by means of numerical continuation technique with the mass parameter continuing from zero to the real value.The remaining part of this paper is structured as follows. In Section2, the dynamical model of the relative motion is briefly introduced and the equations of motion in the rotating coordinate system centered at the equilibrium point of interest are provided. Sections 3 and 4 discuss how to construct high-order analytical solutions of periodic configurations in the relative motion model by means of the L-P and polynomial expansion techniques, respectively. In Section 5, a continuation methodology is given for identifying periodic orbits in the CRTBP. Finally, the conclusions are drawn in Section 6.


## 2. Dynamical Model

In the model of relative motion, the leader satelliteS1 moves on a circular orbit around the Earth, and the follower satellite S2 moves under the gravitational field generated by the Earth in an elliptic orbit. The masses of the Earth and satellites S1 and S2 are denoted as m0,  m1, and m2, respectively. Compared to the Earth mass, the masses of the satellites S1 and S2 are so small such that their gravitational influences upon each other can be ignored. Consequently, the satellites S1 and S2 move in Keplerian orbits around the Earth.In this work, we focus on the relative motions of the follower with respect to the leader satellite. In order to formulate the equations of motion, we introduce a rotating coordinate system, where the origin is centered at the Earth, the starting axisX is directed from the Earth towards the leader satellite, the Z-axis is parallel to the angular momentum vector of the leader satellite, and the Y-axis completes the right-handed coordinate system. Let us denote this rotating system as O-XYZ (see Figure 1 for the planar case).Figure 1
The rotating coordinate system centered at the Earth and the coordinate system centered at the equilibrium point specified byθ. The position vectors of the follower satellite in these two coordinate systems are denoted by R1 and r.To ensure the computation accuracy, the following normalized units are used. The mass of the Earth, the distance between the Earth and the leader satellite, and the leader satellite's period divided by2π are taken as the units of mass, length, and time, respectively. Under this system of dimensionless units, both the gravitation constant G and the spacecraft’s mean motion n become unitary. In the Earth-centered rotating coordinate system, the position of the leader satellite is fixed at 1,0,0.In the Earth-centered synodic coordinate systemO-XYZ, the equations of motion, describing the motions of the follower satellite, can be written as [28](1)X¨-2Y˙=ΩX,Y¨+2X˙=ΩY,Z¨=ΩZ,where X=R1;R˙1=X,Y,Z,X˙,Y˙,Z˙T is the state vector of the follower and Ω is the effective potential of the relative motion model, given by(2)Ω=12X2+Y2+1R1with R1 as the distance between the follower and the Earth, given by(3)R1=X2+Y2+Z2.The dynamical model represented by (1) is an autonomous system and possesses an integral of motion which is similar to the Jacobi constant existing in the CRTBP, given by(4)C=2Ω-X˙2+Y˙2+Z˙2=2Ω-V2.For the relative motion model, there exist special solutions with zero velocity and acceleration, determined by(5)ΩX=X1-1R13=0,ΩY=Y1-1R13=0,ΩZ=-ZR13=0.Similar to the case of the CRTBP, we call these special solutions ‘equilibrium points’. Observing (5), we can conclude that the equilibrium points are located on the motion plane of the leader satellite and satisfy the relation R1=1.0, indicating that all the points located on the planar unitary circle are equilibrium solutions. Around these equilibrium points, there are periodic orbits including planar and vertical periodic motions, which are important and useful in the configuration design of formation flying and satellite constellation.In order to specify a certain equilibrium point on the unitary circle, we introduce an angle parameterθ∈[0,2π] to describe its coordinate like this: (cos⁡θ,sin⁡θ,0). It is observed that the position of the leader corresponds to an equilibrium point specified by θ=0. For this special case, the periodic configuration around the equilibrium point (i.e., the position of the leader satellite) can be applied to satellite formation design. For the other cases corresponding to θ≠0, the periodic configurations around the equilibrium point can be applied to satellite constellation design.To study the motions around an equilibrium point, it is necessary to move the origin of the rotating coordinate system from the Earth to the equilibrium point of interest, and the coordinate axes have the same directions as the ones ofO-XYZ. Denote the rotating coordinate system centered at the equilibrium point of interest by o-xyz (refer to Figure 1 for the planar case), and denote the state of the follower satellite in this coordinate system by x=r;r˙=x,y,z,x˙,y˙,z˙T. The coordinate transformation between O-XYZ and o-xyz can be realized by(6)x=X-cos⁡θ,y=Y-sin⁡θ,z,x˙,y˙,z˙=Z,X˙,Y˙,Z˙.Based on the coordinate transformation, the equations of motion in the coordinate system o-xyz can be written as [22](7)x¨-2y˙=Ωx,y¨+2x˙=Ωy,z¨=Ωz,where the effective potential of the relative motion model becomes(8)Ω=12x+cos⁡θ2+y+sin⁡θ2+1r1with(9)r1=x+cos⁡θ2+y+sin⁡θ2+z2.In order to expand the nonlinear term 1/r1 appearing in (8), the Legendre polynomial defined by(10)1x-A2+y-B2+z-C2=1D∑n≥0ρDnPnAx+By+CzDρis applied, where Pn is the Legendre polynomial of degree n, D2=A2+B2+C2, and ρ2=x2+y2+z2. Consequently, the term 1/r1 can be expressed as(11)1r1=1x+cos⁡θ2+y+sin⁡θ2+z2=∑n≥0ρnPn-xcos⁡θ+ysin⁡θρ.To simplify the mathematical expressions, the following system of notation is used:(12)Tn=ρnPn-xcos⁡θ+ysin⁡θρ,Rnx=∂Tn∂x,Rny=∂Tn∂y,Rnz=∂Tn∂z.Under this notation system, the equations of relative motion in the coordinate system o-xyz can be arranged as(13)x¨-2y˙-3xcos2⁡θ-3ycos⁡θsin⁡θ=∑n≥3Rnx,y¨+2x˙-3xcos⁡θsin⁡θ-3ysin2⁡θ=∑n≥3Rny,z¨+z=∑n≥3Rnz,where Rnx, Rny, and Rnz are polynomials about coordinate variables x,y,z of degree n-1, calculated by the following recursive relations:(14)Rnx=-2n-1ncos⁡θTn-1-2n-1nxcos⁡θ+ysin⁡θRn-1x-2n-1nxTn-2-n-1nρ2Rn-2x,Rny=-2n-1nsin⁡θTn-1-2n-1nxcos⁡θ+ysin⁡θRn-1y-2n-1nyTn-2-n-1nρ2Rn-2y,Rnz=-2n-1nxcos⁡θ+ysin⁡θRn-1z-2n-1nzTn-2-n-1nρ2Rn-2z,starting with(15)R1x=-cos⁡θ,R1y=-sin⁡θ,R1z=0,R2x=x3cos2⁡θ-1+32ysin⁡2θ,R2y=32xsin⁡2θ+y3sin2⁡θ-1,R2z=-z,and Tn is the polynomial of degree n about the variables x,y,z and satisfies by the following recursive relation:(16)Tn=1-2nnxcos⁡θ+ysin⁡θTn-1-n-1nρ2Tn-2,which starts with(17)T0=1,T1=-xcos⁡θ+ysin⁡θ.By now, we have formulated the equations of motion in the rotating coordinate system centered at the equilibrium point specified by the parameter θ. For the equilibrium point specified by θ=0, Gómez and Marcote [22] have constructed high-order analytical solutions of periodic configurations by means of L-P technique. As an extension, in this work we aim to construct the periodic motions including planar and vertical periodic objects around an arbitrary equilibrium point in the relative motion model by means of L-P as well as polynomial expansion techniques. As we consider a general case of relative motion in this work, our results can cover the ones discussed in [22].


## 3. High-Order Solutions Based on L-P Method

By taking advantage of L-P method, many researchers have constructed analytical solutions around equilibrium points in the CRTBP [29–32]. In this part, we take L-P method to construct the planar and vertical periodic configurations around an arbitrary equilibrium point in the relative model of motion.


### 3.1. Periodic Configurations around Equilibrium Point

Ignoring the nonlinear terms of the equations of motion, one can get the following linearized equations of motion:(18)x¨-2y˙-3xcos2⁡θ-32ysin⁡2θ=0,y¨+2x˙-32xsin⁡2θ-3ysin2⁡θ=0,z¨+z=0.Its periodic solution is given by(19)xt=αcos⁡θ1,yt=κ-1αcos⁡θ1+κ-2αsin⁡θ1,zt=βcos⁡θ2,where α and β are the in-plane and out-of-plane amplitudes; θ1 and θ2 are the phase angles of motions defined by θ1=t+θ0 and θ2=t+θ20 with θ10 and θ20 as the initial phase angles. Substituting the linear expression given by (19) in the linearized equations of motion represented by (18) leads to the expressions for computing the coefficients κ-1 and κ-2:(20)κ-1=-3sin⁡2θ6sin2⁡θ+2,κ-2=-23sin2⁡θ+1.Considering the perturbation of the nonlinear part, the motions around an equilibrium point in the relative motion model are expressed as formal series about the in-plane parameter α and out-of-plane parameter β of the form(21)xt=∑i+j=1N∑m≤i,n≤jxijmncos⁡mθ1+nθ2+x-ijmnsin⁡mθ1+nθ2αiβj,yt=∑i+j=1N∑m≤i,n≤jyijmncos⁡mθ1+nθ2+y-ijmnsin⁡mθ1+nθ2αiβj,zt=∑i+j=1N∑m≤i,n≤jzijmncos⁡mθ1+nθ2+z-ijmnsin⁡mθ1+nθ2αiβj,where i,j∈N and m and n have the same parity as i and j, respectively. xijmn, yijmn, and zijmn are the cosine coefficients; x-ijmn, y-ijmn, and z-ijmn are the sine coefficients of coordinates. The phase angles θ1 and θ2 are defined by θ1=ωt+θ10 and θ2=ωt+θ20, with ω as the motion frequency. For the equilibrium points specified by θ=0 and θ=π in the relative motion model, the symmetry properties of equations of motion in the equilibrium point-centered synodic system imply that the sine coefficients of x and z and cosine coefficients of y are zero. For these two special cases, the series expansions are in line with the ones given in [22].Due to the perturbation of nonlinear terms, the frequencies are also expanded as formal series ofα and β as follows: (22)ω=∑i+j=0Nωijαiβj,where ωij with i,j∈N are the coefficients of frequency and are different from zero only if i and j are both even numbers.In the process of solution construction, we denote the order of analytical solutions asN=i+j. For series expansions truncated at order N, all the coefficients including coordinate coefficients corresponding to i+j≤N and frequency coefficients corresponding to i+j≤N-1 need to be determined. Starting with the linear solutions given by (19), the coefficients of high-order series expansions can be constructed iteratively. The details on the construction process are ignored.The series expansions represented by (21) can be used to describe the periodic configurations around the equilibrium point specified by θ. If α≠0 and β=0, (21) can describe the planar periodic configuration; if α=0 and β≠0, it can describe the vertical periodic configuration; and if α≠0 and β≠0, it corresponds to the general periodic configuration.


### 3.2. Results

According to the method mentioned above, we accomplished the computer program on constructing series expansions of motions around the equilibrium points up to an arbitrary order. The program runs fast on our personal computer with the hardware environment Intel(R) core(TM) i5-4590 CPU, 3.30GHz, and 8.0GB RAM. It is about 0.03 sec for the series expansions truncated at order 5, 1.9 sec for the ones truncated at order 10, and 26 sec for the ones truncated at order 15.As an example, the series expansions truncated at order 15 are used to produce the periodic configurations. In practice, four representatives of equilibrium points specified byθ=0, θ=π/3, θ=π, and θ=5π/3 are taken into consideration. The periodic configurations around these considered equilibrium points are reported in Figure 2, where the motion amplitudes are assumed as α=0.1, 0.2, 0.3, and 0.4. The vertical cases are presented in Figure 3, where the motion amplitudes are taken as β=0.1, 0.2, 0.3, and 0.4. From Figures 2 and 3, we can reach the conclusion that the series expansions constructed by means of L-P method up to the 15th order are accurate enough to describe the planar periodic and vertical periodic configurations. Therefore, the series expansions represented by (21) provide high-accuracy parametrization expressions for the motions around equilibrium points in the relative motion model, which can be directly applied to the configuration design of satellite formation flying and constellation, especially in the step of preliminary design. It should be noted that, in all figures presented in this work, the symbol ‘adim’ represents the dimensionless unit of length, unless otherwise stated.Figure 2
Planar periodic configurations with amplitudesα=0.1,0.2,0.3,0.4 around the equilibrium points specified by θ=0, θ=π/3, θ=π, and θ=5π/3 in the Earth-centered rotating coordinate system. These periodic orbits are determined by the analytical solutions truncated at the 15th order.Figure 3
Vertical periodic configurations with amplitudesβ=0.1,0.2,0.3,0.4 around the equilibrium points corresponding to θ=0, θ=π/3, θ=π, and θ=5π/3 in the Earth-centered rotating coordinate system. These periodic orbits are determined by the 15th-order analytical solutions.In order to check the validity of the analytical solutions constructed, the corresponding convergence is evaluated by means of numerical approach. In the practical computation, we assume the maximum value of planar amplitude asαmax=0.5 and that of vertical amplitude as βmax=0.5 and evenly distribute the pair of (α,β) in the domain [0,αmax]×[0,βmax]. Totally, there are 100×100 grids of (α,β) to be evaluated. The corresponding amplitude pairs are denoted by (αi,βi) with i=1,2,…,100 and j=1,2,…,100. For a certain grid i,j, the corresponding state of periodic configuration at time t=0 can be directly obtained by means of the series expansions represented by (21). Starting with the initial state, the equations of motion can be integrated from time t=0 to time t=2π (T=2π corresponds to the period), and the integrated state at time t=2π is denoted by XN. At the same time, the series expansions constructed can directly produce the state of the follower at time t=2π and is denoted by XA. Hence, the state deviation between the numerical solution and the analytical solution at time t=2π, denoted by δij=XA-XN, can be computed. In this work, we take the base 10 logarithm of the state deviation log10(δij) as the performance index of accuracy. Figure 4 shows the practical convergence of the 15th-order series expansions of periodic configuration around the equilibrium points specified by θ=0, θ=π/3, θ=π, and θ=5π/3. It is observed that (a) the convergence domain is different for different equilibrium points and (b) with the decreasing of amplitude, the series expansions constructed perform better.Figure 4
Practical convergence of the 15th-order series expansions of periodic configuration around the equilibrium points specified byθ=0 (up-left panel), θ=π/3 (up-right panel), θ=π (bottom-left panel), and θ=5π/3 (bottom-right panel). The index shown in the color bar represents the performance of accuracy.It should be noted that the numerical integrator we take is a classic eighth-order Runge-Kutta solver with a seventh order for automatic step-size control [33], and the absolute tolerance is taken as 1×10-14.


## 3.1. Periodic Configurations around Equilibrium Point

Ignoring the nonlinear terms of the equations of motion, one can get the following linearized equations of motion:(18)x¨-2y˙-3xcos2⁡θ-32ysin⁡2θ=0,y¨+2x˙-32xsin⁡2θ-3ysin2⁡θ=0,z¨+z=0.Its periodic solution is given by(19)xt=αcos⁡θ1,yt=κ-1αcos⁡θ1+κ-2αsin⁡θ1,zt=βcos⁡θ2,where α and β are the in-plane and out-of-plane amplitudes; θ1 and θ2 are the phase angles of motions defined by θ1=t+θ0 and θ2=t+θ20 with θ10 and θ20 as the initial phase angles. Substituting the linear expression given by (19) in the linearized equations of motion represented by (18) leads to the expressions for computing the coefficients κ-1 and κ-2:(20)κ-1=-3sin⁡2θ6sin2⁡θ+2,κ-2=-23sin2⁡θ+1.Considering the perturbation of the nonlinear part, the motions around an equilibrium point in the relative motion model are expressed as formal series about the in-plane parameter α and out-of-plane parameter β of the form(21)xt=∑i+j=1N∑m≤i,n≤jxijmncos⁡mθ1+nθ2+x-ijmnsin⁡mθ1+nθ2αiβj,yt=∑i+j=1N∑m≤i,n≤jyijmncos⁡mθ1+nθ2+y-ijmnsin⁡mθ1+nθ2αiβj,zt=∑i+j=1N∑m≤i,n≤jzijmncos⁡mθ1+nθ2+z-ijmnsin⁡mθ1+nθ2αiβj,where i,j∈N and m and n have the same parity as i and j, respectively. xijmn, yijmn, and zijmn are the cosine coefficients; x-ijmn, y-ijmn, and z-ijmn are the sine coefficients of coordinates. The phase angles θ1 and θ2 are defined by θ1=ωt+θ10 and θ2=ωt+θ20, with ω as the motion frequency. For the equilibrium points specified by θ=0 and θ=π in the relative motion model, the symmetry properties of equations of motion in the equilibrium point-centered synodic system imply that the sine coefficients of x and z and cosine coefficients of y are zero. For these two special cases, the series expansions are in line with the ones given in [22].Due to the perturbation of nonlinear terms, the frequencies are also expanded as formal series ofα and β as follows: (22)ω=∑i+j=0Nωijαiβj,where ωij with i,j∈N are the coefficients of frequency and are different from zero only if i and j are both even numbers.In the process of solution construction, we denote the order of analytical solutions asN=i+j. For series expansions truncated at order N, all the coefficients including coordinate coefficients corresponding to i+j≤N and frequency coefficients corresponding to i+j≤N-1 need to be determined. Starting with the linear solutions given by (19), the coefficients of high-order series expansions can be constructed iteratively. The details on the construction process are ignored.The series expansions represented by (21) can be used to describe the periodic configurations around the equilibrium point specified by θ. If α≠0 and β=0, (21) can describe the planar periodic configuration; if α=0 and β≠0, it can describe the vertical periodic configuration; and if α≠0 and β≠0, it corresponds to the general periodic configuration.


## 3.2. Results

According to the method mentioned above, we accomplished the computer program on constructing series expansions of motions around the equilibrium points up to an arbitrary order. The program runs fast on our personal computer with the hardware environment Intel(R) core(TM) i5-4590 CPU, 3.30GHz, and 8.0GB RAM. It is about 0.03 sec for the series expansions truncated at order 5, 1.9 sec for the ones truncated at order 10, and 26 sec for the ones truncated at order 15.As an example, the series expansions truncated at order 15 are used to produce the periodic configurations. In practice, four representatives of equilibrium points specified byθ=0, θ=π/3, θ=π, and θ=5π/3 are taken into consideration. The periodic configurations around these considered equilibrium points are reported in Figure 2, where the motion amplitudes are assumed as α=0.1, 0.2, 0.3, and 0.4. The vertical cases are presented in Figure 3, where the motion amplitudes are taken as β=0.1, 0.2, 0.3, and 0.4. From Figures 2 and 3, we can reach the conclusion that the series expansions constructed by means of L-P method up to the 15th order are accurate enough to describe the planar periodic and vertical periodic configurations. Therefore, the series expansions represented by (21) provide high-accuracy parametrization expressions for the motions around equilibrium points in the relative motion model, which can be directly applied to the configuration design of satellite formation flying and constellation, especially in the step of preliminary design. It should be noted that, in all figures presented in this work, the symbol ‘adim’ represents the dimensionless unit of length, unless otherwise stated.Figure 2
Planar periodic configurations with amplitudesα=0.1,0.2,0.3,0.4 around the equilibrium points specified by θ=0, θ=π/3, θ=π, and θ=5π/3 in the Earth-centered rotating coordinate system. These periodic orbits are determined by the analytical solutions truncated at the 15th order.Figure 3
Vertical periodic configurations with amplitudesβ=0.1,0.2,0.3,0.4 around the equilibrium points corresponding to θ=0, θ=π/3, θ=π, and θ=5π/3 in the Earth-centered rotating coordinate system. These periodic orbits are determined by the 15th-order analytical solutions.In order to check the validity of the analytical solutions constructed, the corresponding convergence is evaluated by means of numerical approach. In the practical computation, we assume the maximum value of planar amplitude asαmax=0.5 and that of vertical amplitude as βmax=0.5 and evenly distribute the pair of (α,β) in the domain [0,αmax]×[0,βmax]. Totally, there are 100×100 grids of (α,β) to be evaluated. The corresponding amplitude pairs are denoted by (αi,βi) with i=1,2,…,100 and j=1,2,…,100. For a certain grid i,j, the corresponding state of periodic configuration at time t=0 can be directly obtained by means of the series expansions represented by (21). Starting with the initial state, the equations of motion can be integrated from time t=0 to time t=2π (T=2π corresponds to the period), and the integrated state at time t=2π is denoted by XN. At the same time, the series expansions constructed can directly produce the state of the follower at time t=2π and is denoted by XA. Hence, the state deviation between the numerical solution and the analytical solution at time t=2π, denoted by δij=XA-XN, can be computed. In this work, we take the base 10 logarithm of the state deviation log10(δij) as the performance index of accuracy. Figure 4 shows the practical convergence of the 15th-order series expansions of periodic configuration around the equilibrium points specified by θ=0, θ=π/3, θ=π, and θ=5π/3. It is observed that (a) the convergence domain is different for different equilibrium points and (b) with the decreasing of amplitude, the series expansions constructed perform better.Figure 4
Practical convergence of the 15th-order series expansions of periodic configuration around the equilibrium points specified byθ=0 (up-left panel), θ=π/3 (up-right panel), θ=π (bottom-left panel), and θ=5π/3 (bottom-right panel). The index shown in the color bar represents the performance of accuracy.It should be noted that the numerical integrator we take is a classic eighth-order Runge-Kutta solver with a seventh order for automatic step-size control [33], and the absolute tolerance is taken as 1×10-14.


## 4. High-Order Solutions Based on Polynomial Expansion Method

To construct the polynomial solutions of periodic motions, we take the position and velocity coordinates in one direction as a pair of modal coordinates [25, 27] and express the components in the other directions as polynomial expansion of the modal coordinates.


### 4.1. Polynomial Expansions of Planar Periodic Motions

For the planar case, (13) can be written as(23)x¨-2y˙=∑k=1n∑i+j=kcijxiyj,y¨+2x˙=∑k=1n∑i+j=kdijxiyj,where cij and dij are the known coefficients of the equations of motion. In this work, the state variables in the x-direction are taken as a pair of modal coordinates:(24)x=u,x˙=v.Polynomial expansion approach implies that the remaining motion components including y and y˙ can be expressed as polynomial series of modal coordinates up to order n of the form(25)yu,v=a1,0u+a0,1v+a2,0u2+a1,1uv+a0,2v2+…=∑k=1n∑i+j=kai,juivi,(26)y˙u,v=b1,0u+b0,1v+b2,0u2+b1,1uv+b0,2v2+…=∑k=1n∑i+j=kbi,juivi,where aij and bij are the coefficients to be determined.By substituting (24)–(26) into (23), the planar equations of motion can be expressed as polynomial series with respect to the modal coordinates u,v up to order n like this:(27)x¨=2y˙+∑k=1n∑i+j=kcijxiyj=∑k=1n∑i+j=klijuivj,(28)y¨=-2x˙+∑k=1n∑i+j=kdijxiyj=∑k=1n∑i+j=kpijuivj,where the coefficients lij and pij are functions of the undetermined coefficients.Taking time derivatives of (25) and (26), we can obtain the following expressions:(29)y˙u,v=∂y∂uu˙+∂y∂vv˙=∂y∂uv+∂y∂vx¨=∑k=1n∑i+j=kqijuivj,and(30)y¨=∂y˙∂uu˙+∂y˙∂vv˙=∂y∂uv+∂y∂vx¨=∑k=1n∑i+j=krijuivj,where the coefficients qij and rij are also functions of the coefficients to be determined.Comparing (26) with (29), (28) with (30), and equating the coefficients of the terms with the same powers of u and v, we can get the algebraic equations about the undetermined coefficients aij,bij up to order n, where i,j∈N and 0≤i,j≤n. By solving these equations from the first order to the desired order, all the unknown coefficients can be determined.Following the above-mentioned procedure, we can construct the polynomial expressions represented by (25) and (26), which can be used to generate the states of periodic configurations. For the planar motions, if any two state components are known, the remaining state components can be directly identified by solving the polynomial expansion constraints. The simulation results are to be reported in Section 4.3.


### 4.2. Polynomial Expansions of Vertical Periodic Configuration

For the three-dimensional case, the equations of motion represented by (13) can be written as(31)x¨-2y˙=∑m=1N∑i+j+k=meijkxiyjzk,y¨+2x˙=∑m=1N∑i+j+k=mfijkxiyjzk,z¨=∑m=1N∑i+j+k=mhijkxiyjzk,where eijk, fijk, and hijk are the known coefficients of the equations of motion.For the three-dimensional case, we take the state components in thez-direction as the modal coordinates:(32)z=u,z˙=v.Based on polynomial expansion approach, the state components in the x- and y-directions are expressed as polynomial series, up to order n, of the form(33)xu,v=a1,0u+a0,1v+a2,0u2+a1,1uv+a0,2v2+…=∑k=1n∑i+j=kai,juivj,(34)yu,v=b1,0u+b0,1v+b2,0u2+b1,1uv+b0,2v2+…=∑k=1n∑i+j=kbi,juivj,(35)x˙u,v=c1,0u+c0,1v+c2,0u2+c1,1uv+c0,2v2+…=∑k=1n∑i+j=kci,juivj,and(36)y˙u,v=d1,0u+d0,1v+d2,0u2+d1,1uv+d0,2v2+…=∑k=1n∑i+j=kdi,juivj,where aij, bij, cij, and dij are the coefficients to be determined. Substituting (32)–(36) into (31) leads to the three-dimensional equations of motion of the form(37)x¨=2y˙+∑m=1n∑i+j+k=meijkxiyjzk=∑k=1n∑i+j=klijuivj,y¨=-2x˙+∑m=1n∑i+j+k=mfijkxiyjzk=∑k=1n∑i+j=koijuivj,z¨=∑m=1n∑i+j+k=mhijkxiyjzk=∑k=1n∑i+j=kpijuivj,where the coefficients lij, oij, and pij are functions of the undetermined coefficients. Taking time derivatives of (33)–(36), we can derive the following expressions:(38)x˙u,v=∂x∂uu˙+∂x∂vv˙=∂x∂uv+∂x∂vz¨=∑k=1n∑i+j=kqijuivj,(39)y˙u,v=∂y∂uu˙+∂y∂vv˙=∂y∂uv+∂y∂vz¨=∑k=1n∑i+j=krijuivj,(40)x¨u,v=∂x˙∂uu˙+∂x˙∂vv˙=∂x˙∂uv+∂x˙∂vz¨=∑k=1n∑i+j=ksijuivj,and(41)y¨u,v=∂y˙∂uu˙+∂y˙∂vv˙=∂y˙∂uv+∂y˙∂vz¨=∑k=1n∑i+j=kwijuivj.Here qij, rij, sij, and wij are functions of the coefficients of polynomial expansions.Similar to the planar case, comparing (35)–(37) with (38)–(41) and equating the coefficients of the terms with the same powers of u and v, we can get the algebraic equations about the undetermined coefficients. By solving them from the first order to the desired order n, the polynomial expansions of vertical periodic configurations in the relative motion model can be constructed.In fact, the polynomial expansions represented by (33)–(36) provide four constraint functions, indicating that there are only two independent state variables for the three-dimensional periodic objects. That is to say, given that any two state components are known, all the remaining four motion components can be generated by the polynomial relations. Therefore, the polynomial expansions constructed can be used to describe the three-dimensional periodic motion around the equilibrium points in the model of relative motion.


### 4.3. Results

As discussed in Sections4.1 and 4.2, the polynomial expansions represented by (25) and (26) provide two polynomial constraints for the planar periodic motions, and the ones represented by (33)–(36) provide four polynomial constraints for vertical periodic configurations around equilibrium points in the relative motion model. With considerations of the polynomial constraints, for both the planar and vertical periodic motion modes, there are only two independent state variables. Given that two of the state components are provided, the remaining state variables can be directly generated by solving the nonlinear system determined by the polynomial expansions. Usually, for the planar case, the pair of state components x0,y0 is predetermined, and for the vertical case, the pair of state components (z0,z˙0) is predetermined.For the equilibrium point specified byθ=0 corresponding to the position of the leader satellite, the integrated orbits for the planar and vertical cases are compared in Figure 5, where the initial states are provided by the polynomial expansions truncated at different orders. For the planar case, the first- and fourth-order polynomial expansions are considered, and for the vertical case, the fourth- and eighth-order polynomial expansions are taken into account. Figure 5 tells us that the order of the polynomial expansion is higher; the corresponding integrated trajectory performs better in terms of periodicity.Figure 5
Comparison of the integrated orbits for the planar (left pane) and vertical (right panel) configurations around the equilibrium point specified by θ=0, where the initial states are provided by the polynomial expansions truncated at different orders.Next, the polynomial expansions of planar and vertical periodic configurations are constructed up to order eight, and they are used to generate the initial states of periodic configurations around the equilibrium points specified byθ=0, θ=π/3, θ=π, and θ=5π/3. Starting from these initial states provided by polynomial expansions, the equations of motion are integrated over time interval t∈[0,2π], and the resulting integrated orbits are reported in Figure 6 for the planar case and in Figure 7 for the vertical case. Figures 6 and 7 indicate that the polynomial expansions constructed are accurate enough for describing the periodic configurations within a certain amplitude range in the relative motion model.Figure 6
The integrated planar periodic configurations around the equilibrium points corresponding toθ=0, θ=π/3, θ=π, and θ=5π/3, where the initial states are generated by the 8th-order polynomial expansions.Figure 7
The integrated vertical periodic configurations around the equilibrium points corresponding toθ=0, θ=π/3, θ=π, and θ=5π/3, where the initial states are provided by the 8th-order polynomial expansions.


### 4.4. Accuracy Analysis

In order to check the validity, it is necessary to perform accuracy analysis for the polynomial expansions constructed. As discussed above, for a given pair of state componentsx0,y0 for the planar case and (z0,z˙0) for the vertical case, the remaining state variables of periodic motions can be produced by the polynomial expansions. By taking the state generated by polynomial expansions as initial guesses, the accurate state of periodic configuration can be identified by means of numerical techniques, such as shooting method [34]. Then the accuracy performance is defined as the Euclidian norm of the state deviation between the initial state generated by polynomial expansions and the numerically corrected state.For the periodic motions around the equilibrium point specified byθ=0, the accuracy performances are shown in Figure 8 as a function of the motion amplitude (x0 coordinate for the planar case and z0 for the vertical case) for the polynomial expansions truncated at the sixth and eighth orders. The following conclusions can be reached: (a) the performance increases with the increasing of the motion amplitude, and (b) the higher-order polynomial expansions perform better than the lower-order ones in terms of accuracy.Figure 8
The state deviation between the initial state provided by polynomial expansions and the numerically corrected states as a function of the motion amplitude for the planar periodic configurations (left panel) and vertical periodic configurations (right panel) around the equilibrium point corresponding to θ=0 in the relative motion model.


## 4.1. Polynomial Expansions of Planar Periodic Motions

For the planar case, (13) can be written as(23)x¨-2y˙=∑k=1n∑i+j=kcijxiyj,y¨+2x˙=∑k=1n∑i+j=kdijxiyj,where cij and dij are the known coefficients of the equations of motion. In this work, the state variables in the x-direction are taken as a pair of modal coordinates:(24)x=u,x˙=v.Polynomial expansion approach implies that the remaining motion components including y and y˙ can be expressed as polynomial series of modal coordinates up to order n of the form(25)yu,v=a1,0u+a0,1v+a2,0u2+a1,1uv+a0,2v2+…=∑k=1n∑i+j=kai,juivi,(26)y˙u,v=b1,0u+b0,1v+b2,0u2+b1,1uv+b0,2v2+…=∑k=1n∑i+j=kbi,juivi,where aij and bij are the coefficients to be determined.By substituting (24)–(26) into (23), the planar equations of motion can be expressed as polynomial series with respect to the modal coordinates u,v up to order n like this:(27)x¨=2y˙+∑k=1n∑i+j=kcijxiyj=∑k=1n∑i+j=klijuivj,(28)y¨=-2x˙+∑k=1n∑i+j=kdijxiyj=∑k=1n∑i+j=kpijuivj,where the coefficients lij and pij are functions of the undetermined coefficients.Taking time derivatives of (25) and (26), we can obtain the following expressions:(29)y˙u,v=∂y∂uu˙+∂y∂vv˙=∂y∂uv+∂y∂vx¨=∑k=1n∑i+j=kqijuivj,and(30)y¨=∂y˙∂uu˙+∂y˙∂vv˙=∂y∂uv+∂y∂vx¨=∑k=1n∑i+j=krijuivj,where the coefficients qij and rij are also functions of the coefficients to be determined.Comparing (26) with (29), (28) with (30), and equating the coefficients of the terms with the same powers of u and v, we can get the algebraic equations about the undetermined coefficients aij,bij up to order n, where i,j∈N and 0≤i,j≤n. By solving these equations from the first order to the desired order, all the unknown coefficients can be determined.Following the above-mentioned procedure, we can construct the polynomial expressions represented by (25) and (26), which can be used to generate the states of periodic configurations. For the planar motions, if any two state components are known, the remaining state components can be directly identified by solving the polynomial expansion constraints. The simulation results are to be reported in Section 4.3.


## 4.2. Polynomial Expansions of Vertical Periodic Configuration

For the three-dimensional case, the equations of motion represented by (13) can be written as(31)x¨-2y˙=∑m=1N∑i+j+k=meijkxiyjzk,y¨+2x˙=∑m=1N∑i+j+k=mfijkxiyjzk,z¨=∑m=1N∑i+j+k=mhijkxiyjzk,where eijk, fijk, and hijk are the known coefficients of the equations of motion.For the three-dimensional case, we take the state components in thez-direction as the modal coordinates:(32)z=u,z˙=v.Based on polynomial expansion approach, the state components in the x- and y-directions are expressed as polynomial series, up to order n, of the form(33)xu,v=a1,0u+a0,1v+a2,0u2+a1,1uv+a0,2v2+…=∑k=1n∑i+j=kai,juivj,(34)yu,v=b1,0u+b0,1v+b2,0u2+b1,1uv+b0,2v2+…=∑k=1n∑i+j=kbi,juivj,(35)x˙u,v=c1,0u+c0,1v+c2,0u2+c1,1uv+c0,2v2+…=∑k=1n∑i+j=kci,juivj,and(36)y˙u,v=d1,0u+d0,1v+d2,0u2+d1,1uv+d0,2v2+…=∑k=1n∑i+j=kdi,juivj,where aij, bij, cij, and dij are the coefficients to be determined. Substituting (32)–(36) into (31) leads to the three-dimensional equations of motion of the form(37)x¨=2y˙+∑m=1n∑i+j+k=meijkxiyjzk=∑k=1n∑i+j=klijuivj,y¨=-2x˙+∑m=1n∑i+j+k=mfijkxiyjzk=∑k=1n∑i+j=koijuivj,z¨=∑m=1n∑i+j+k=mhijkxiyjzk=∑k=1n∑i+j=kpijuivj,where the coefficients lij, oij, and pij are functions of the undetermined coefficients. Taking time derivatives of (33)–(36), we can derive the following expressions:(38)x˙u,v=∂x∂uu˙+∂x∂vv˙=∂x∂uv+∂x∂vz¨=∑k=1n∑i+j=kqijuivj,(39)y˙u,v=∂y∂uu˙+∂y∂vv˙=∂y∂uv+∂y∂vz¨=∑k=1n∑i+j=krijuivj,(40)x¨u,v=∂x˙∂uu˙+∂x˙∂vv˙=∂x˙∂uv+∂x˙∂vz¨=∑k=1n∑i+j=ksijuivj,and(41)y¨u,v=∂y˙∂uu˙+∂y˙∂vv˙=∂y˙∂uv+∂y˙∂vz¨=∑k=1n∑i+j=kwijuivj.Here qij, rij, sij, and wij are functions of the coefficients of polynomial expansions.Similar to the planar case, comparing (35)–(37) with (38)–(41) and equating the coefficients of the terms with the same powers of u and v, we can get the algebraic equations about the undetermined coefficients. By solving them from the first order to the desired order n, the polynomial expansions of vertical periodic configurations in the relative motion model can be constructed.In fact, the polynomial expansions represented by (33)–(36) provide four constraint functions, indicating that there are only two independent state variables for the three-dimensional periodic objects. That is to say, given that any two state components are known, all the remaining four motion components can be generated by the polynomial relations. Therefore, the polynomial expansions constructed can be used to describe the three-dimensional periodic motion around the equilibrium points in the model of relative motion.


## 4.3. Results

As discussed in Sections4.1 and 4.2, the polynomial expansions represented by (25) and (26) provide two polynomial constraints for the planar periodic motions, and the ones represented by (33)–(36) provide four polynomial constraints for vertical periodic configurations around equilibrium points in the relative motion model. With considerations of the polynomial constraints, for both the planar and vertical periodic motion modes, there are only two independent state variables. Given that two of the state components are provided, the remaining state variables can be directly generated by solving the nonlinear system determined by the polynomial expansions. Usually, for the planar case, the pair of state components x0,y0 is predetermined, and for the vertical case, the pair of state components (z0,z˙0) is predetermined.For the equilibrium point specified byθ=0 corresponding to the position of the leader satellite, the integrated orbits for the planar and vertical cases are compared in Figure 5, where the initial states are provided by the polynomial expansions truncated at different orders. For the planar case, the first- and fourth-order polynomial expansions are considered, and for the vertical case, the fourth- and eighth-order polynomial expansions are taken into account. Figure 5 tells us that the order of the polynomial expansion is higher; the corresponding integrated trajectory performs better in terms of periodicity.Figure 5
Comparison of the integrated orbits for the planar (left pane) and vertical (right panel) configurations around the equilibrium point specified by θ=0, where the initial states are provided by the polynomial expansions truncated at different orders.Next, the polynomial expansions of planar and vertical periodic configurations are constructed up to order eight, and they are used to generate the initial states of periodic configurations around the equilibrium points specified byθ=0, θ=π/3, θ=π, and θ=5π/3. Starting from these initial states provided by polynomial expansions, the equations of motion are integrated over time interval t∈[0,2π], and the resulting integrated orbits are reported in Figure 6 for the planar case and in Figure 7 for the vertical case. Figures 6 and 7 indicate that the polynomial expansions constructed are accurate enough for describing the periodic configurations within a certain amplitude range in the relative motion model.Figure 6
The integrated planar periodic configurations around the equilibrium points corresponding toθ=0, θ=π/3, θ=π, and θ=5π/3, where the initial states are generated by the 8th-order polynomial expansions.Figure 7
The integrated vertical periodic configurations around the equilibrium points corresponding toθ=0, θ=π/3, θ=π, and θ=5π/3, where the initial states are provided by the 8th-order polynomial expansions.


## 4.4. Accuracy Analysis

In order to check the validity, it is necessary to perform accuracy analysis for the polynomial expansions constructed. As discussed above, for a given pair of state componentsx0,y0 for the planar case and (z0,z˙0) for the vertical case, the remaining state variables of periodic motions can be produced by the polynomial expansions. By taking the state generated by polynomial expansions as initial guesses, the accurate state of periodic configuration can be identified by means of numerical techniques, such as shooting method [34]. Then the accuracy performance is defined as the Euclidian norm of the state deviation between the initial state generated by polynomial expansions and the numerically corrected state.For the periodic motions around the equilibrium point specified byθ=0, the accuracy performances are shown in Figure 8 as a function of the motion amplitude (x0 coordinate for the planar case and z0 for the vertical case) for the polynomial expansions truncated at the sixth and eighth orders. The following conclusions can be reached: (a) the performance increases with the increasing of the motion amplitude, and (b) the higher-order polynomial expansions perform better than the lower-order ones in terms of accuracy.Figure 8
The state deviation between the initial state provided by polynomial expansions and the numerically corrected states as a function of the motion amplitude for the planar periodic configurations (left panel) and vertical periodic configurations (right panel) around the equilibrium point corresponding to θ=0 in the relative motion model.


## 5. Applications to Periodic Orbits in the CRTBP

The problem of relative motion can be regarded as a special case of the CRTBP with system parameter asμ=0. In this part, with the aid of continuation concept, we endeavour to identify periodic orbits in the CRTBP by continuing the mass parameter from μ=0 to its real value.


### 5.1. Parameter-Dependent Dynamical Model

Let us introduce a parameterε∈[0,1] to control the mass parameter of the CRTBP. The controlled equations of motion in the barycentric synodic coordinate system, describing the motions of the third body, can be written as [28, 35](42)X¨-2Y˙=ΩX,Y¨+2X˙=ΩY,Z¨=ΩZ,where Ω is the effective potential of system, determined by(43)Ω=12X2+Y2+1-εμR1+εμR2,0≤ε≤1,with R1 and R2 as the distances of the third body from the primaries(44)R1=X+εμ2+Y2+Z2,R2=X-1+εμ2+Y2+Z2.Observing (42)–(44), one can get that the case of ε=0 represents the model of relative motion about circular reference orbit, which has been studied in Sections 3 and 4 in this paper, the case of ε=1 corresponds to the real CRTBP, and the case of ε∈(0,1) corresponds to the intermediate dynamical model. By continuing from ε=0 to ε=1, the dynamical model can transit from the relative motion model to the real CRTBP. Naturally, the solutions in the relative motion model could be continued to the ones in the CRTBP.


### 5.2. Results

In this part, the Earth-Moon CRTBP is adopted as an example dynamical system with mass parameter asμ=0.012150585609624 obtained by JPL planetary ephemeris DE405 [36]. Based on the continuation process discussed in Section 5.1, the periodic solutions in the CRTBP are identified by continuing from the periodic configurations obtained by high-order analytical solutions. In the process of continuation, differential corrector is used [34].Periodic configurations around the equilibrium points specified byθ=0, θ=π/3, θ=π, and θ=5π/3 in the dynamical model of relative motion are generated by the analytical solutions constructed in Sections 3 or 4. Starting from these initial solutions, the periodic orbits in the Earth-Moon CRTBP can be identified by means of the continuation technique discussed in Section 5.1. Figure 9 shows the planar periodic orbits around the Moon and equilibrium points Lii=3,4,5, and Figure 10 gives the vertical periodic orbits around equilibrium points Lii=3,4,5. It is observed that the orbit shown in the up-left panel of Figure 9 corresponds to the distant retrograde orbits around the Moon, which is continued from the periodic configuration around the equilibrium point specified by θ=0 in the relative motion model.Figure 9
Planar periodic orbits around the Moon and equilibrium pointsLii=3,4,5 in the Earth-Moon CRTBP, continued from the planar periodic configurations around the equilibrium points corresponding to θ=0,π,π/3, and 5π/3 in the model of relative motion.Figure 10
Vertical periodic orbits around the libration pointsLii=3,4,5 in the Earth-Moon CRTBP, continued from the planar periodic configurations around the equilibrium points corresponding to θ=π,π/3, and 5π/3 in the model of relative motion.It should be pointed out that not all the periodic orbits in the CRTBP can be continued from the solutions in the relative motion model, such as periodic orbits around collinear libration pointsLii=1,2, long-periodic orbits around triangular libration points Lii=4,5, and halo orbits. This is because we cannot find the corresponding periodic objects in the model of relative motion for these kinds of periodic motions. However, the continuation method taking the mass parameter as a continuing variable provides an alternative approach to determining the periodic orbits in the CRTBP.


## 5.1. Parameter-Dependent Dynamical Model

Let us introduce a parameterε∈[0,1] to control the mass parameter of the CRTBP. The controlled equations of motion in the barycentric synodic coordinate system, describing the motions of the third body, can be written as [28, 35](42)X¨-2Y˙=ΩX,Y¨+2X˙=ΩY,Z¨=ΩZ,where Ω is the effective potential of system, determined by(43)Ω=12X2+Y2+1-εμR1+εμR2,0≤ε≤1,with R1 and R2 as the distances of the third body from the primaries(44)R1=X+εμ2+Y2+Z2,R2=X-1+εμ2+Y2+Z2.Observing (42)–(44), one can get that the case of ε=0 represents the model of relative motion about circular reference orbit, which has been studied in Sections 3 and 4 in this paper, the case of ε=1 corresponds to the real CRTBP, and the case of ε∈(0,1) corresponds to the intermediate dynamical model. By continuing from ε=0 to ε=1, the dynamical model can transit from the relative motion model to the real CRTBP. Naturally, the solutions in the relative motion model could be continued to the ones in the CRTBP.


## 5.2. Results

In this part, the Earth-Moon CRTBP is adopted as an example dynamical system with mass parameter asμ=0.012150585609624 obtained by JPL planetary ephemeris DE405 [36]. Based on the continuation process discussed in Section 5.1, the periodic solutions in the CRTBP are identified by continuing from the periodic configurations obtained by high-order analytical solutions. In the process of continuation, differential corrector is used [34].Periodic configurations around the equilibrium points specified byθ=0, θ=π/3, θ=π, and θ=5π/3 in the dynamical model of relative motion are generated by the analytical solutions constructed in Sections 3 or 4. Starting from these initial solutions, the periodic orbits in the Earth-Moon CRTBP can be identified by means of the continuation technique discussed in Section 5.1. Figure 9 shows the planar periodic orbits around the Moon and equilibrium points Lii=3,4,5, and Figure 10 gives the vertical periodic orbits around equilibrium points Lii=3,4,5. It is observed that the orbit shown in the up-left panel of Figure 9 corresponds to the distant retrograde orbits around the Moon, which is continued from the periodic configuration around the equilibrium point specified by θ=0 in the relative motion model.Figure 9
Planar periodic orbits around the Moon and equilibrium pointsLii=3,4,5 in the Earth-Moon CRTBP, continued from the planar periodic configurations around the equilibrium points corresponding to θ=0,π,π/3, and 5π/3 in the model of relative motion.Figure 10
Vertical periodic orbits around the libration pointsLii=3,4,5 in the Earth-Moon CRTBP, continued from the planar periodic configurations around the equilibrium points corresponding to θ=π,π/3, and 5π/3 in the model of relative motion.It should be pointed out that not all the periodic orbits in the CRTBP can be continued from the solutions in the relative motion model, such as periodic orbits around collinear libration pointsLii=1,2, long-periodic orbits around triangular libration points Lii=4,5, and halo orbits. This is because we cannot find the corresponding periodic objects in the model of relative motion for these kinds of periodic motions. However, the continuation method taking the mass parameter as a continuing variable provides an alternative approach to determining the periodic orbits in the CRTBP.


## 6. Conclusions

In this work, high-order analytical solutions are constructed for periodic configurations around equilibrium points in the relative motion model with circular reference orbit by means of L-P and polynomial expansion approaches. The series expansions constructed by these two methods could be used to describe the planar and vertical periodic configurations. Particularly, the solution constructed by L-P method provides an explicit parametrization expression, while the solution constructed by polynomial expansion technique determines several polynomial constraints for the state variables of periodic motions. To check the validity, the practical convergence is evaluated by means of numerical technique. The results indicate that the analytical solutions constructed by these two approaches are available for describing periodic objects in the relative motion model. Subsequently, considering the fact that the relative motion model corresponds to a special case of the CRTBP with zero mass parameter, a control parameter is introduced to continue the dynamical system from the relative motion model to the real CRTBP. Following the continuation concept, the periodic orbits in the CRTBP are identified by continuing from the periodic configurations in the relative motion model.



---

*Source: 4183941-2018-07-09.xml*