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

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

#### 2. Dynamical Model

In the model of relative motion, the leader satellite moves on a circular orbit around the Earth, and the follower satellite moves under the gravitational field generated by the Earth in an elliptic orbit. The masses of the Earth and satellites and are denoted as , , and , respectively. Compared to the Earth mass, the masses of the satellites and are so small such that their gravitational influences upon each other can be ignored. Consequently, the satellites and 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 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 Earth-centered rotating coordinate system, the position of the leader satellite is fixed at .

In the Earth-centered synodic coordinate system , the equations of motion, describing the motions of the follower satellite, can be written as [28]where is the state vector of the follower and is the effective potential of the relative motion model, given bywith as the distance between the follower and the Earth, given byThe 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 byFor the relative motion model, there exist special solutions with zero velocity and acceleration, determined bySimilar 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 , 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 to describe its coordinate like this: . It is observed that the position of the leader corresponds to an equilibrium point specified by . 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 , 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 of . Denote the rotating coordinate system centered at the equilibrium point of interest by (refer to Figure 1 for the planar case), and denote the state of the follower satellite in this coordinate system by . The coordinate transformation between and can be realized byBased on the coordinate transformation, the equations of motion in the coordinate system can be written as [22]where the effective potential of the relative motion model becomeswithIn order to expand the nonlinear term appearing in (8), the Legendre polynomial defined byis applied, where is the Legendre polynomial of degree , , and . Consequently, the term can be expressed asTo simplify the mathematical expressions, the following system of notation is used:Under this notation system, the equations of relative motion in the coordinate system can be arranged aswhere , , and are polynomials about coordinate variables of degree , calculated by the following recursive relations:starting withand is the polynomial of degree about the variables and satisfies by the following recursive relation:which starts withBy 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 , 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:Its periodic solution is given bywhere and are the in-plane and out-of-plane amplitudes; and are the phase angles of motions defined by and with and 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 and :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 formwhere and and have the same parity as and , respectively. , , and are the cosine coefficients; , , and are the sine coefficients of coordinates. The phase angles and are defined by and , with as the motion frequency. For the equilibrium points specified by 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 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 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 and , (21) can describe the planar periodic configuration; if and , it can describe the vertical periodic configuration; and if and , 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 , , , and 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.2, 0.3, and 0.4. The vertical cases are presented in Figure 3, where the motion amplitudes are taken as , 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 and that of vertical amplitude as and evenly distribute the pair of in the domain . Totally, there are grids of to be evaluated. The corresponding amplitude pairs are denoted by with and . For a certain grid , the corresponding state of periodic configuration at time 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 to time ( corresponds to the period), and the integrated state at time is denoted by . At the same time, the series expansions constructed can directly produce the state of the follower at time and is denoted by . Hence, the state deviation between the numerical solution and the analytical solution at time , denoted by , can be computed. In this work, we take the base 10 logarithm of the state deviation 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 , , , and . 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 .

#### 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 aswhere 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 and can be expressed as polynomial series of modal coordinates up to order of the formwhere 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 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:andwhere 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 , we can get the algebraic equations about the undetermined coefficients up to order , where and . 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 aswhere , , 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 formandwhere , , , and are the coefficients to be determined. Substituting (32)–(36) into (31) leads to the three-dimensional equations of motion of the formwhere the coefficients , , and are functions of the undetermined coefficients. Taking time derivatives of (33)–(36), we can derive the following expressions:andHere , , , 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 , 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 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 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 is predetermined, and for the vertical case, the pair of state components is predetermined.

For the equilibrium point specified by 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 , , , and . Starting from these initial states provided by polynomial expansions, the equations of motion are integrated over time interval , 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.

##### 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 components for the planar case and 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 , the accuracy performances are shown in Figure 8 as a function of the motion amplitude ( coordinate for the planar case and 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.

#### 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 . In this part, with the aid of continuation concept, we endeavour to identify periodic orbits in the CRTBP by continuing the mass parameter from to its real value.

##### 5.1. Parameter-Dependent Dynamical Model

Let us introduce a parameter 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]where is the effective potential of system, determined bywith and as the distances of the third body from the primariesObserving (42)–(44), one can get that the case of 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 corresponds to the real CRTBP, and the case of corresponds to the intermediate dynamical model. By continuing from to , 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 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 , , , and 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 , and Figure 10 gives the vertical periodic orbits around equilibrium points . 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 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 , long-periodic orbits around triangular libration points , 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.

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

#### Acknowledgments

This work was carried out with financial support of the National Natural Science Foundation (no. 11603011, 41774038), the Natural Science Foundation of Jiangsu province (no. BK20160612), the National Basic Research Program 973 of China (2015CB857100), the Satellite Communication and Navigation Collaborative Innovation Center (no. SatCN-201409), and National Defense Scientific Research Fund (no. 2016110C019).