Periodic Solutions of Nonlinear Relative Motion Satellites

: The relative motion of an outline of the rendezvous problem has been studied by assuming that the chief satellite is in circular symmetric orbits. The legitimacy of perturbation techniques and nonlinear relative motion are investigated. The deputy satellite equations of motion with respect to the ﬁxed references at the center of the chief satellite are nonlinear in the general case. We found the periodic solutions of the linear relative motion satellite and for the nonlinear relative motion satellite using the Lindstedt–Poincaré technique. Comparisons among the analytical solutions of linear and nonlinear motions and the obtained solution by the numerical integration of the explicit Euler method for both motions are investigated. We demonstrate that both analytical and numerical solutions of linear motion are symmetric periodic. However, the solutions of nonlinear motion obtained by the Lindstedt–Poincaré technique are periodic and the numerical solutions obtained by integration by using explicit Euler method are non-periodic. Thus, the Lindstedt–Poincaré technique is recommended for designing the periodic solutions. Furthermore, a comparison between linear and nonlinear analytical solutions of relative motion is investigated graphically.


Importance of Relative Motion Satellites
Relative motion is considered one of the fundamental problems in orbital mechanics, it has different applications that arise in the literature of astrodynamics for rendezvous and formation of flying satellites. The first idea concerning this problem was explained by the author of [1] in the late 19th century, who proposed analyzing the Moon's motion. His objective was to establish a more mathematically sound method for developing tables of lunar motion, this was at the time based on "practical astronomy rather than of mathematics" in his talk.
The first practical space of relative motion was in the field of intercept and space rendezvous throughout the late 1950s and extends up to the current time. In the intercept model, a chase satellite is compelled to move in such a way that its orbit intersects the trajectory of the chief satellite at a certain time. In the rendezvous problem the relative velocities of both vehicles must be pushed to zero at the intersection instant; thereby, a docking or berthing procedure or other activities may be carried out.
in a circular orbit, and acts as a control center. It sends messages with relative location and velocity data to slave satellites. Based on that data, the satellites perform a docking maneuver and rendezvous by using an on-board propulsion system.
The comprehension of relative motion satellites is substantial for two significant applications of orbital mechanics: • The space rendezvous problem. • The formation flight problem.
The first problem investigates the designing of centralization and decentralization algorithms to dock a satellite successfully, for instance, the space shuttle with the International Space Station. These algorithms demand the modeling of relative motion satellites, without needing permanent communication with ground-based stations. Rendezvous orbital dynamics and control is one of the key elements of space flight technology for operating space rendezvous and docking missions. Rendezvous demands a specific congruent of the orbital velocities and position vectors of the two vehicles, permitting them to stay at a constant separation during orbital station-keeping. It may or may not be followed by docking or berthing, procedures which bring the spacecraft into physical connection and generate a link between them. The same rendezvous method could be used for a spacecraft "landing" on a natural celestial body with a weak gravitational field, for example, landing on one of the Martian moons will demand the same conveniency of orbital velocities, followed by a "descent" that shares some similarities with docking.
The second problem is related to the formation of a satellite system, where the dynamics of the system are exploited to place a number of satellites in space. Formation of satellites is a novel concept that distributes a large spacecraft function among several small satellites, which are less expensive and cooperative. Most applications to relative motion are in the formation satellites, which have a considerable significance, because the use of a large number of satellites with low impact satellites running, in a concrete style, could introduce a better outcome than a single high implementation satellite. Furthermore, these formations do not require a big cost, which means a great opportunity for space mission success with more resiliency. Satellite formations have good benefits for Earth observation of space missions, where the group distribution of low resolution devices, operating in conjunction with each other, may provide a higher overall information quality than a single high resolution device.

Mathematical Models of Relative Motion Satellites
The relative motion between two satellites is described by several mathematical models [3][4][5][6]. The "Hill-Clohessy-Wiltshire model (HCW)" is as a classical approach described by a set of linearized, time-invariant ordinary differential equations. However, this model assumes a circular chief orbit and it is valid only for small initial separations between the satellites.
An assumption of an elliptical shape of the chief satellite orbit results in time-varying differential equations. In this work such a model is represented by the Tschauner-Hempel (TH) equations. However, since the applicability of this model is also limited to small initial separations between the satellites, this investigation deploys it only for comparative purposes.
In cases where a chief satellite (target or passive vehicle) is moving in a circular orbit, the orbital relative motion can be described using HCW equations [2]. A principal assumption for the application of HCW equations is a distance of separation between chief and deputy (chase or active vehicle) satellites of less than 1 km.
The authors of [7] were the first to find a closed-form solution for linearized relative motion within the frame of elliptic orbits. The authors of [8] found a singularity in Lawden's solution and then constructed a state-transition matrix that depends explicitly on the true and the eccentric anomaly, overcoming the singularity. We find a closed-form solution for a linearized relative motion within the frame of elliptic orbits in [7] for the first time. The author of [8] found a singularity in Lawden's solution. Afterwards, they proposed a new state-transition matrix without a singularity by explicitly using the true and the eccentric anomaly data. They derived the state transition matrix (STM) of the (TH) equations for orbits with small eccentricity, based mainly on the solutions of the HCW equations using a perturbation theory. In the framework of improving the obtained solution, approximation solutions of the second-order relative motion equations set in spherical coordinates are constructed [9].
Several versions of the relative motion dynamics have been developed to involve the perturbations effect [10][11][12][13]. The inhomogeneous gravitational field generated by the presence of the non-sphericity of objects, such as the Earth's oblateness, is based on an assumption that the gravitational field is a sum of zonal harmonics coefficients; the authors of [10] derived a set of relative dynamics time-variant equations under the effect of the zonal harmonic parameter. An analytical solution of relative motion subject to Lorentzforce perturbations is constructed in [14]. Within the frame of developing the analysis of relative motion using state transition, matrices for the elliptical case are found, by the authors in [15,16], to be under the perturbation effect.
Furthermore, an outstanding piece of work studied the reformulation of relative motion based on non-linear system dynamics; see [17]. Another excellent piece of work proposed to analyze relative motion dynamics, control techniques and specific applications for formation flying, deployment, station keeping, and reconfiguration; see [18][19][20][21][22].
The aim of the proposed paper is to find the Periodic solutions of nonlinear relative motion using the perturbation technique of the Lindstedt-Poincaré method, which enable us to remove the secular terms, which will be growth with increasing time and will go to infinity with increasing time to infinity. However, there are some numerical methods, which can be used to estimate such solutions; for example: The Verrlet method, the step size control algorithm and the explicit Runge-Kutta method; for details, comparisons and effective methods, see [23,24]. However, we emphasize that the numerical methods are not used to verify the analytical solutions in most cases, and vice versa. In addition, there is no full warranty for convergence of numerical solutions in most cases. Thus, we are interested to find the periodic solution of non-linear relative motion, which is more realistic in practical applications by using the perturbation technique of Lindstedt-Poincaré.

Formulation of Motion Model
We consider a set of two satellites model, the first is called target satellite (or chief or passive vehicle) and the second is called chaser (or deputy or active vehicle). In general, this model is implemented to accomplish some space mission such as space rendezvous. The mission between two space crafts or satellites is successful when the two satellites have the same velocity and position vectors at a certain time. In the initiates of rendezvous time the two parts in a rendezvous sequence are:

•
Phasing of maneuvers to accomplish a rendezvous in the timing sequence, which will bring the two satellites into close proximity. • Maneuvers of finalizing rendezvous that involve the relative motion between the two satellites for rendezvous and docking. Now we will assume that the inertial references frame with origin, which coincides with the center of the main body, say Earth's center; the orthogonal unit vector is denoted by {i x , i y , i z }. The vector i x in the direction of the equinoxes and the vectors i y is perpendicular to it, while i z passes through the North Pole. To describe a relative motion satellite in Local-Vertical-Local-Horizontal (LVLH) coordinates, we proposed that the orientation of these frame by the unit vector triad {i r , i θ , i h }, where the vector i r lies in the chief radial direction, i h lies in the direction of the chief angular momentum, and i θ completes the right-handed orthogonal triad, such that i θ = i h ∧ i r ; see Figure 1. According to Figure 1, the deputy satellite relative position vector in LVLH coordinates can be written as: where ξ, ζ, and ϑ are the components of the relative position vector along the radial, long track, and out of plane, respectively. Since where r c = ri r represents the current radius of the chief's satellite orbit. The LVLH frame rotation with angular velocity ω is given by where ω r = ω θ = 0 and ω = ω h i h , while the attached frame rotates with the angular velocity ω h =θ =ḟ . We have to note that θ = Ω + f in which Ω and f are the argument of perigee and the true anomaly, respectively, within the frame of two-body motion. The angular momentum magnitude of the chief satellite is expressed by where r is the trajectory of the chief (target) satellite, which may be circular or elliptical. On the other hand, the expression of the chief's true anomaly rateθ can be written in terms of a semi-latus rectum asθ wherein the elliptical case p = a(1 − e 2 ), a are the semi-major axes and e is the eccentricity, µ is the gravitational constant and r = p/(1 + e cos f ). The angular velocity and the acceleration of the the LVLH reference frame are given by Now, the general vectorial equations that model the relative motion in the LVLH reference frame areR where T is the sum of total external forces acting on the deputy and chief, respectively; (.) and (..) denote the first and second derivative with respect to time. Substituting Equations (1) and (2) into Equation (3), we obtain Now we can write the right-hand side of Equation (4) in the form where T d is the external forces acting on the deputy and T c is the external forces acting on the chief satellite, which can be written in the case of the gravitational field being generated by a spherical object, such as Hence, Substitution Equation (5) into Equation (4); then, the vector form or relative motion equations will be controlled by The vector Equation (6) can be degraded tö Equation (7) represents the general equations of the satellite relative motion in LVLH coordinates.

General Circular Relative Motion Satellites
Now, we denote v c and v d as the potential functions for gravitational forces T c and T d , acting on the chief and deputy satellite, respectively; hence, Note that (ξ/r) = (ξ/R).(R/r); Equation (8) is a generating function for Legendre polynomials with argument −(ξ/R); consequently, where P n is the Legendre polynomial with degree n, where P 0 (ξ/R) = 1, where v h is the higher-order of a Legendre polynomial given by Then the total external forces can be rewritten as Since ∆v d is the gravitational acceleration acting on the deputy and ∆v c is the gravitational acceleration acting on the chief's satellite, while ∆v is the total relative gravitational acceleration, then and Utilizing Equations (9) and (10) with Equation (4), the vector equation of relative motion can be rewritten as Hence, Equation (11) can be written as three differential equations of second order Equation (12) represents the general dynamical system of the relative motion satellite. However, we have to observe that the quantities in the right-hand sides, ∂v h /∂ξ, ∂v h /∂ϑ, ∂v h /∂ζ, represent the components of perturbing gravitational accelerations, which are considered the main contribution for higher order of nonlinearity to a dynamical system of relative motion.

Linear Circular Relative Motion Satellites
In the case where the chief's orbit is moving in circular orbits, then r = a andθ = n are constants. Thereby, the rotation frame with these properties is referred to as Hill's frame and rotates with angular velocity n = µ/a 3 , and the resulting equations are referred to as HCW equations, as in [2]. In addition, if the perturbing acceleration or the higher order on nonlinearity terms (v h = 0) are neglected. Then the linearized relative motion equations will be governed byξ Let ξ = x, ϑ = y, ζ = z; thereby, the above equations can be written as The above dynamical system is famous and is called the Hill-Clohessy-Wiltshire equations, and represents the linearity of relative satellite motion. Furthermore, the dynamical system of relative motion (HCW) is invariant under the symmetry (x, y, z) → (−x, −y, −z).

Solutions of HCW Equations
In order to construct the solutions of the dynamical system which represent the linear relative motion satellite, we will assume that the deputy satellite starts its motion under the following set of the initial conditions z(0) = z 0 0 ,ż(0) =ż 0 0 , Thereby, after integration of the second equation in System (13), we obtaiṅ where X 0 is a constant which can be determined from the initial conditions, so one obtains and Equation (15) will take the following forṁ Substitutingẏ into the first equation in System (13), we geẗ Let ψ = x − 2X 0 /n, then dψ/dt = dX/dt and d 2 ψ/dt 2 = d 2 X/dt 2 and, the Equation (17) becomes Hence, Equation (18) has the solution Then the final form solution is Substituting Equation (19) into Equation (16) and integrating it, we get where Y 0 the constant of the integration.
To complete our steps, we have to integrate the third equation in System (13), which represents the out-plane motion; hence, we get z = Z sin(nt + ϕ 0 ).
Using the obtained solutions of Equations (19)- (21), the solutions of System (13) are summarized in the following formulas where X, Y, Z, θ 0 , and ϕ 0 are arbitrary constants that must match the initial conditions. Here, we observe that the second formula in Solutions (22) contains an unbounded term, which will grow to infinity with increasing time. To get bounded solutions, in particular, periodic solutions, the coefficient of the unbounded term must equal zero; these will give an additional condition to get periodic and bounded solutions; thereby, we will obtain the new conditions on the initial conditions aṡ So we rewrite the periodic solutions of System (13) in the following formulas Applying the first derivative to Equations (24), we obtain the relative velocities in the formẋ If we apply the initial conditions on Equation (24) and the condition in Equation (23), then the constants of integration can be determined by Equations (24) and (25) represent the relative positions and velocities of the chaser satellite for space rendezvous within the framework of linearized motion.
With the help of Equations (24)- (27), the initial amplitude of the chaser satellite for space rendezvous, i.e., (x 0 0 , y 0 0 , z 0 0 ) is taken as (0.13, 0.12, 0.10), with different values of phase angle. The projection of the periodic solution in the xy-plane is shown in Figure 2, whereas the orbit in three dimensions is shown in Figure 3 at different values of the phase angle. It is clear from Figure 2a that the analytical solution of linear relative motion in two dimensions is periodic, and the obtained solution by a numerical integration of the explicit Euler method using step size 10 −3 at a set of initial condition x 0 0 , y 0 0 ,ẋ 0 0 ,ẏ 0 0 as (0, 1.2, 1.0, 0) is also periodic; this solution is shown in Figure 2b. In addition, the two solutions are symmetric due to the y-axis (x = 0).
The linear motion in three dimensions is also periodic due two both solutions. We remark that both analytical and numerical solutions of linear relative motion in three dimensions are also periodic. The analytical solution is given in Figure 3a, while the solution by a numerical integration of the explicit Euler method using step size 10 −3 at a set of initial condition x 0 0 , y 0 0 ,ẋ 0 0 ,ẏ 0 0 as (0, 1.2, 1.0, 0) is given in Figure 3b.

Nonlinear Circular Relative Motion Satellites
In the case of considering the perturbed potential v h = 0, and keeping the nonlinearity up to second order terms, the relative motion satellite equations of System (12) are approximated to the following equations The above dynamical system consists of the three second-order differential equations with second degree, which represent the relative motion of the deputy satellite.
Now we intend to find the periodic solutions of these equations to find the relative position and velocity of the deputy satellite to accomplish space rendezvous. The solutions of these equations will be generalized or extended for the linear solutions; then, we have to find the solutions under the previous conditions of linear relative motion. Unfortunately, there are no exact analytical solutions for these equations, and the exact solutions are not possible. Thereby, we must introduce techniques or perturbation methods, which can be used to evaluate approximate or semi-analytical solutions; in particular, we must construct periodic solutions. So, the next section will be devoted to analysis of a perturbation technique, which enables us to find an appropriate approximated analytical solution to perform a space rendezvous mission.

Periodic Solution of Nonlinear Relative Motion Satellites
This section is devoted to constructing the periodic solutions of the nonlinear relative motion satellite. From physical point of view, it is sufficient to find the solutions of the dynamical system within the frame of second degree or third in maximum. This view can be applied for the model of relative motion satellite, because the relative distance (distance between chief and deputy satellites) is very small compared to the initial position of the target satellite.

Legitimacy of Perturbation Techniques
The classical perturbation theory, which is called the straightforward expansion technique, the Lindstedt-poincaré method, and the multiple scales method are used when the dynamical systems are either weakly nonlinear or weakly non-autonomous; this means that the effects of either nonlinearity or non-autonomy are very small. Alternatively, the dynamical systems with these properties can be considered either linear or quasi-linear. In most cases, the systems of differential equations have linear expressions and small nonlinear or non-autonomous expressions in such a way that they are separated from each other.
In the context of either weak nonlinearity or weak non-autonomy, more approximations may still be required by finding a parameter in the dynamical system , which can be used as a very small perturbation effect, but its value will be replaced by one at the end. This technique is reasonable and accepted when the terms of the dynamical system are quasi-linear or small in themselves. The obtained parameters in this case are called "place keeping parameters".

Legitimacy of Relative Motion Equations
In line with the reasons outlined in the previous subsection, we can neither apply the multiple scales technique nor the Lindstedt-Poincaré method to System (28). Thereby, we will use the place keeping parameters method, after carrying out the relative equations of motion in dimensionless variables.
In order to set the equations in dimensionless variables, we assume that (ξ, ϑ, ζ) = r(x, y, z); then, System (28) can be written as Since R 0 << r, we can write = R 0 r << 1, where R 0 is the measure of relative orbit size. Regarding this idea, we can use the place keeping parameters method, and replace the coordinate (x, y, z) by ( x, y, z). Then System (29) can be written as The above equations represent the nonlinear relative motion satellite in dimensionless variables, which are convenient for applying some perturbation techniques, such as multiple scales technique [25][26][27] and KBM methods or Lindstedt-Poincaré [28][29][30][31]. However, we will examine the provided solution by the latter method with the initial conditions of the linear relative motion satellite; these conditions are describe in Equation (14).

Periodic Solutions by the Lindstedt-Poincar é Technique
Now we will apply the Lindstedt-Poincaré technique on the relative satellite motion to find the periodic solution of the equations of motion of the second degree, say System (30). For simplicity, we suppose that θ = nt; then, System (30) can be written in the form To account for the nonlinear dependence of the frequency, we will introduce the stretch variable, τ = ωθ; we denote dΓ/dτ = Γ and d 2 Γ/dτ 2 = Γ ; then, System (31) can be written in the form Now we look for approximated solutions in the form of a power series as and The initial conditions of position components can be rewritten in the form Substituting Equation (33) into Equation (32) with the help of Equation (34) and keeping the terms up to first order in , we obtain where the frequency is The set of initial conditions in Equations (35) and (36) may be rewritten as and x (0, ) = x 0 (0) + x 1 (0) = x 0 0 , y (0, ) = y 0 (0) + y 1 (0) = y 0 0 , z (0, ) = z 0 (0) + z 1 (0) = z 0 0 , Equating each of the coefficients of that have the same power in System (37), using Conditions (38) and (39), we get the unperturbed equations system in the form and the set of initial conditions is while the equations which represent the first corrections are with the initial conditions To accomplish our procedure and find the approximated periodic solutions, we have to find solutions for both the dynamical systems in Equations (40) and (42)

Zero order solutions
According to the solutions of linear motion, the general periodic solutions of the unperturbed relative motion will be controlled by With a help of Conditions (41), the constants of integration are given by

First order solutions
Substituting Equation (44) into Equation (42), we obtain x 1 − 2y 1 − 3x 1 = X 11 + X 12 cos(τ + θ 0 ) + X 13 sin(τ + θ 0 ) + X 14 sin(τ + ϕ 0 ) + X 15 cos 2(τ + θ 0 ) + X 16 cos 2(τ + ϕ 0 ), where X 1i , i = 1, 2, 3, . . . 6, Y ij and Z 1j , j = 1, 2, 3 are constants, which are calculated by the following relations Equation (45) consists of a dynamical system of second-order non-homogeneous differential equations. After integrating this system with the help of Conditions (43), we observe that the terms with coefficients X 1i = 1, 2, 3, 4, Y 11 , Y 12 , and Z 12 represent secular terms, as well as the extra term with coefficient S = 2X 2 + Z 2 + 3X 2 cos 2θ 0 . Therefore, if these coefficients are nonzero, then the correction solutions will involve secular terms in the scale τ. However, that is exactly what we need to avoid. Hence we will remove these terms by equating their coefficients by zeros. These procedures ensure that the solutions will not include secular terms and at least no secularities appear in the first two terms in the perturbation series. Furthermore, the terms with coefficients X 1i = 1, 2, 3, 4, Y 11 , Y 12 , and Z 12 represent solutions for the associated homogeneity of Equation (45); hence, we set the coefficients to zeros before integration.
In this context, the removal of the secular terms yields Then the first correction will be controlled by Therefore, the periodic solutions of the nonlinear relative motion satellite can be written in the form where τ = nωt.
Since we use the method of place keeping parameters, we have to take = 1; then, Solutions (47) will be rewritten as x(τ, ) = x 0 (τ) + x 1 (τ), y(τ, ) = y 0 (τ) + y 1 (τ), z(τ, ) = z 0 (τ) + z 1 (τ). With the help of Equation (48), the initial amplitude of the chaser satellite for space rendezvous, i.e., (x 0 0 , y 0 0 , z 0 0 ), works similar to in the previous section (0.13, 0.12, 0.10) with different values for the phase angle. The projection of the periodic solution in the xy-plane is shown in Figure 4, whereas the orbit in three dimensions is shown in Figure 5 with different values of phase angle.
A comparison between the nonlinear relative motion solutions of the Lindstedt-Poincaré technique obtained by the numerical integration of the explicit Euler method are shown in Figure 4, in two dimensions, and in Figure 5, in three-dimensional motion. While numerical solutions for linear motion are periodic, the numericality of nonlinear motion is not periodic; this can be observed in Figures 4b and 5b The projections in the xy-plane of the periodic solution of linear and nonlinear motion are shown in Figure 6 at different phases, whereas the differences between orbits in three dimensions of the linear and nonlinear motion at the same different phases angle with ϕ 0 = π/4 are shown in Figure 7.    Figures 6 and 7 show that the analytical solutions of linear and nonlinear motions are periodic; however, there are differences between the style of the two motions. The periodicity of the nonlinear relative motion solutions gives a considerable verification due to the solution obtained with the Lindstedt-Poincaré technique.

Conclusions
In this paper, we studied the relative motion of an outline of the space rendezvous problem, assuming that the chief satellite was in circular symmetric orbits. It was found that the equations that describe the motion of the deputy satellite with respect to the references fixed at the center of the chief satellite is nonlinear in the general case. Then the periodic solutions first for linear relative motion satellite were evaluated. Here a key role is the symmetry of the periodic solution joined with the fact that the system HCW is invariant under negative symmetry.
The Lindstedt-Poincaré technique was used to get the approximated periodic solutions for nonlinear differential equations of the second degree. This method is referred to by many sources as a method of successive approximations. In our study, we use it to find the first corrections of a relative motion satellite. In this context, the obtained solutions are more general compared to the linear solutions. The analysis of a linear relative motion satellite and its periodic solutions are constructed. The approximate periodic solutions of the nonlinear relative motion satellite are constructed using the Lindstedt-Poincaré technique.
Comparisons among the analytical solution of linear motion and the obtained solution by the numerical integration of the explicit Euler method are investigated, and we show that bth solutions are periodic in two and three dimensions. Moreover, comparisons among the obtained solutions of nonlinear relative motion using the Lindstedt-Poincaré technique and the numerical solution by integration of the explicit Euler method were studied. However, the numerical solutions are not periodic in both motions in two and three dimensions. Thus, the Lindstedt-Poincaré technique is recommended for designing the periodic solutions.