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

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 themethod 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.


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 pointmass bodies under the gravitational influence of a central body. Later, an important contribution was made by Lawden and Tschauner et al. [2][3][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][6][7][8][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][11][12][13][14][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][20][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][25][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 threebody 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 Section 2, 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.

Dynamical Model
In the model of relative motion, the leader satellite 1 moves on a circular orbit around the Earth, and the follower satellite 2 moves under the gravitational field generated by the Earth in an elliptic orbit. The masses of the Earth and satellites 1 and 2 are denoted as 0 , 1 , and 2 , respectively. Compared to the Earth mass, the masses of the satellites 1 and 2 are so small such that their gravitational influences upon each other can be ignored. Consequently, the satellites 1 and 2 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 axis is directed from the Earth towards the leader satellite, the -axis is parallel to the angular momentum vector of the leader satellite, and the -axis completes the right-handed coordinate system. Let us denote this rotating system as − (see Figure 1 for the planar case). 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 by 2 are taken as the units of mass, length, and time, respectively. Under this system of dimensionless units, both the gravitation constant and the spacecraft's mean motion become unitary. In the Earthcentered rotating coordinate system, the position of the leader satellite is fixed at (1, 0, 0).
In the Earth-centered synodic coordinate system − , the equations of motion, describing the motions of the follower satellite, can be written as [28] − 2̇= Ω , where = ( 1 ;̇1) = ( , , ,̇,,) is the state vector of the follower and Ω is the effective potential of the relative motion model, given by with 1 as the distance between the follower and the Earth, given by 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 For the relative motion model, there exist special solutions with zero velocity and acceleration, determined by 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 1 = 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 Mathematical Problems in Engineering 3 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 of − . Denote the rotating coordinate system centered at the equilibrium point of interest by − (refer to Figure 1 Based on the coordinate transformation, the equations of motion in the coordinate system − can be written as [22]̈− 2̇= Ω , where the effective potential of the relative motion model becomes In order to expand the nonlinear term 1/ 1 appearing in (8), the Legendre polynomial defined by is applied, where is the Legendre polynomial of degree , 2 = 2 + 2 + 2 , and 2 = 2 + 2 + 2 . Consequently, the term 1/ 1 can be expressed as To simplify the mathematical expressions, the following system of notation is used: = .
(12) Under this notation system, the equations of relative motion in the coordinate system − can be arranged as where , , and are polynomials about coordinate variables ( , , ) of degree − 1, calculated by the following recursive relations: starting with and is the polynomial of degree about the variables ( , , ) and satisfies by the following recursive relation: which starts with 0 = 1, 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].

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][30][31][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.

Periodic Configurations around Equilibrium Point.
Ignoring the nonlinear terms of the equations of motion, one can get the following linearized equations of motion: Its periodic solution is given by where and are the in-plane and out-of-plane amplitudes; 1 and 2 are the phase angles of motions defined by 1 = + 0 and 2 = + 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 : 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 where , ∈ N and and have the same parity as and , respectively. , , and are the cosine coefficients; , , and are the sine coefficients of coordinates. The phase angles 1 and 2 are defined by 1 = + 10 and 2 = + 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 and and cosine coefficients of 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: where with , ∈ N are the coefficients of frequency and are different from zero only if and are both even numbers.
In the process of solution construction, we denote the order of analytical solutions as = + . For series expansions truncated at order , all the coefficients including coordinate coefficients corresponding to + ≤ and frequency coefficients corresponding to + ≤ − 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 Mathematical Problems in Engineering 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.
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. solution and the analytical solution at time = 2 , denoted by = | − |, can be computed. In this work, we take the base 10 logarithm of the state deviation log 10 ( ) 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.
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 .

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.

Polynomial Expansions of Planar Periodic Motions.
For the planar case, (13) can be written as where and are the known coefficients of the equations of motion. In this work, the state variables in the -direction are taken as a pair of modal coordinates: Polynomial expansion approach implies that the remaining motion components including anḋcan be expressed as where and 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 ( , V) up to order like this: where the coefficients and are functions of the undetermined coefficients.
Taking time derivatives of (25) and (26), we can obtain the following expressions: where the coefficients and 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 and V, we can get the algebraic equations about the undetermined coefficients ( , ) up to order , where , ∈ N and 0 ≤ , ≤ . 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

Polynomial Expansions of Vertical Periodic Configuration.
For the three-dimensional case, the equations of motion represented by (13) can be written as where , , and ℎ are the known coefficients of the equations of motion.
For the three-dimensional case, we take the state components in the -direction as the modal coordinates: Based on polynomial expansion approach, the state components in the -and -directions are expressed as polynomial series, up to order , of the form anḋ( where , , , and are the coefficients to be determined. Substituting (32)- (36) into (31) leads to the threedimensional equations of motion of the form̈= where the coefficients , , and are functions of the undetermined coefficients. Taking time derivatives of (33)-(36), we can derive the following expressions: and̈( Here , , , and 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 and V, we can get the algebraic equations about the undetermined coefficients. By solving them from the first order to the desired order , 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 threedimensional 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.

Results.
As discussed in Sections 4.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 ( 0 , 0 ) is predetermined, and for the vertical case, the pair of state components ( 0 ,̇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.
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 ∈ [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.

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 components ( 0 , 0 ) for the planar case and ( 0 ,̇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 ( 0 coordinate for the planar case and 0 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.

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.

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] − 2̇= Ω , where Ω is the effective potential of system, determined by with 1 and 2 as the distances of the third body from the primaries 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.

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 ( = 3, 4, 5), and Figure 10 gives the vertical periodic orbits around equilibrium points ( = 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.
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 points ( = 1, 2), long-periodic orbits around triangular libration points ( = 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.

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.

Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
The authors declare that they have no conflicts of interest.