#### Abstract

A numerical approach for obtaining periodic orbits of satellite relative motion is proposed, based on using the time domain collocation (TDC) method to search for the periodic solutions of an exact nonlinear relative model. The initial conditions for periodic relative orbits of the Clohessy-Wiltshire (C-W) equations or Tschauner-Hempel (T-H) equations can be refined with this approach to generate nearly bounded orbits. With these orbits, a method based on the least-squares principle is then proposed to generate projected closed orbit (PCO), which is a reference for the relative motion control. Numerical simulations reveal that the presented TDC searching scheme is effective and simple, and the projected closed orbit is very fuel saving.

#### 1. Introduction

Relative motion of satellites in earth orbit has attracted a great deal of attentions in the past four decades. Various models accounting for the relative motions of two orbital satellites have been proposed and applied in a variety of space missions. The first developed and perhaps also the most famous mathematical model for the relative motions is the Clohessy-Wiltshire (C-W) equations [1], which have been shown to be useful in solving the rendezvous problems. The C-W equations also provide an approach to determine the initial conditions for bounded or periodic relative motions for satellites in near-circular orbits. However, some of the assumptions used in the C-W equations introduce modeling errors that cannot be ignored. That is, the initial condition based on the C-W equations is valid only in some special cases with circular reference orbits, spherical earth, and linearized gravitational acceleration. Since the periodic relative orbits are important in relative orbit maintenance as well as in the tasks of spacecraft formation flying, more accurate models for relative motions and new approaches for obtaining periodic orbits are necessary.

In literature, many studies have been conducted to generalize the C-W equations. Tschauner and Hempel [2] incorporated the eccentricity for the reference orbit, based on which analytical expressions for relative motions as well as closed form solutions to the periodic relative motion were derived. Multifarious state transition matrices for relative motion near an orbit of arbitrary eccentricity were also developed, including the work of Melton [3] and Broucke [4], which used time as an independent variable, while the others use true anomaly instead. In the work of Inalhan et al. [5], the initial conditions starting from which the T-H equations generate periodic solutions were obtained using the homogenous solutions of the linearized relative equations of motion, which originated from the results by Lawden [6], Carter and Humi [7], and Carter [8]. It can be used to initialize to a closed-form aperture for a large fleet of vehicles with an eccentric reference orbit. However, considering nonlinear differential gravity, perturbation, and large eccentricity, the initial conditions derived by Inalhan et al. do not work anymore.

To take into account the effects of nonlinear terms of gravity potential, eccentricity of reference orbit, and earth oblateness, many efforts have been done to establish more exact models of relative motion for the spacecraft in Earth orbit. The T-H equations perturbed by nonlinear differential gravity with no assumption on eccentricity were first presented by Euler and Shulman [9]. After that, many extended versions dealing with various nonlinearities and perturbations (e.g., the atmospheric drag and third body effects) are proposed. Unfortunately, there are no analytical solutions to these models due to the existence of nonlinearity which nullifies the additivity and homogeneity pertaining to the linear systems.

Xu and Wang [10] presented a much more generalized dynamic model for satellite relative motions, which took perturbation, nonlinear terms of gravity potential, and eccentricity into account in their study. Since the dynamic equations are derived without approximation, it is an exact nonlinear relative model. In this study, we use it to search for some periodic relative orbits using the time domain collocation (TDC) method. The TDC method proposed by Dai et al. [11] is essentially one of the weighted residual methods and has been successfully employed to solve a variety of the nonlinear dynamical problems [12–14]. In the TDC method, the desired periodic solution is preassumed as a truncated Fourier series first, which has been adopted in the previous works by Wang et al. [15] and Kasdin and Gurfil [16]. Then, the approximate solution of the Fourier series is substituted into the equations of nonlinear system, resulting in a residual error function. The residual error function is then enforced to zero at a spectrum of selected collocation points, which leads to a system of nonlinear algebraic equations (NAEs) with Fourier coefficients (plus the natural frequency if an autonomous system is analyzed) being the variables. The NAEs can then be solved by the NAE solvers, such as the classical Newton-Raphson method (NR). In previous studies, the TDC method has been applied to obtain the periodic solutions of the Van der Pol oscillator and the Duffing oscillator. Here, the TDC method is extended to approximate the initial conditions for periodic relative orbits. The initial conditions for generating periodic relative motions of nearly circular orbits without perturbation can be obtained directly by Clohessy-Wiltshire equations. However, in the presence of perturbations, the C-W initial conditions have to be refined to yield closed orbits. Using the TDC method, some adjustments of the initial conditions can be made to accommodate the perturbations.

The structure of this paper is organized as follows. In Section 2, the TDC searching scheme for periodic relative orbits is proposed. In Section 3, the method of selecting initial values for TDC method is introduced. In Section 4, the results of TDC searching scheme are evaluated by simulating the fuel consumption. In Section 5, some numerical results are presented to demonstrate the efficiency of this approach. Finally, we draw some conclusions in Section 6.

#### 2. TDC Searching Scheme for Periodic Relative Orbits

To implement the TDC method, a function in the form of a truncated Fourier expansion is considered as the approximation of the periodic solution, which can be expressed as Here is the number of harmonics included in the approximation, is the supposed frequency of the periodic motion, and are the harmonic coefficient variables. Obviously, the results of TDC method will be inherently periodic.

Collocating at time points equally spaced in a period of , then from (1) we get It can be seen that the coefficients are related to the collocation points via the transformation For clear presentation, we define the transformation matrix Thus, if we can get the values of the time points , the harmonic coefficients can be determined by , where is the inverse () or pseudoinverse () matrix of .

Using expression (1), the first order time derivative of can be written as Collocating at , we have Referring to the derivation of (3), the relationship between the harmonic coefficients and the time derivatives can be written in a matrix formIt can also be expressed in a more subtle way where Note that using (3), can be expressed by . Substituting (3) into (8), we get Thus, the transformation from to is established.

In the following parts, we will use the preceding equations and transformations to derive the TDC algebraic system for a nonlinear relative motion model, thereby approximating the periodic solutions.

The exact nonlinear relative model between two elliptic satellite orbits derived by Xu and Wang is herein employed. Two Cartesian coordinate frames were used as shown in Figure 1. The Earth centered inertial (ECI) coordinate frame is spanned by the unit vectors (). The local vertical local horizontal (LVLH) coordinate is attached on the reference satellite .

With the gravity of the Earth incorporated, the relative motion of the satellite in the LVLH frame is expressed by Remark , , and ; then we get where , , , , , , , , and are periodic time-varying parameters related to the chief orbit. They are described in the Earth centered inertial (ECI) coordinate frame by a set of differential equations which are independent of relative motion. For the reason that the equations of chief orbit are utilized directly in the searching procedure without treatment, they are listed in the appendix. , , and are control forces of the satellite in the LVLH coordinate. , are functions of the relative positions , , and and can be obtained from with Note that although (12) is equivalent to (11) mathematically, it suits TDC method better. By increasing the number of equations in the system, the terms of second order time derivatives are removed. The advantage is that the approximate functions of relative velocities , , and can be selected separately from those of relative positions , , and with (12); that is, the information of both velocities and positions participates in the iteration procedure for solving this system and, thus, the error introduced by approximation will be reduced.

For simplicity, we take on behalf of (12), with , which is the state vector of the relative motion. If there exists any periodic solution for the nonlinear relative model, the components of it can certainly be approximated by the truncated Fourier series expansions similar to .

Considering that the components of relative motion on three directions of LVLH frame can have different frequencies, we assume where , , , , , and take the same form of , except that the frequencies of relative motion around axes , , and are , , and separately. To determine the unknown coefficients and frequencies in (15), we have to derive the TDC algebraic system of (12) first.

Collocate points in a time period of ; then we get Recalling (10), the left hand part of (16) can be converted towith the order of components in (16) being rearranged. For instance, , and so on. Besides, , , and are transformation matrixes corresponding to , , and separately. For brevity, the block diagonal matrix on the right hand of the equation is denoted by in the following parts.

By substituting the preceding equation into (16) and making some arrangements, we get the TDC algebraic system of (12) where . is the rearranged version of the vector .

To solve this algebraic system, the Newton method can be very appropriate. Considering the complexity of (18), we derive the Jacobian matrix in a numerical way. Denote , and assume , with being a small deviation from , which is the th component of . Then Jacobian matrix can be computed by , with being the th column of . Once is obtained, the relationship revealed in (3) can be used to give the approximate periodic solution of the relative model.

With , which is the collocation point of at time , given by , we can use it to initialize the relative orbits. The resulting orbits will be bounded if the approximation is close enough to the accurate periodic solution. In practical terms, the initialization procedure often takes place at a certain relative position. This can be realized by introducing additional constraints on the position of starting point or, in most cases, the first collocation point. The corresponding equations will be shown in Section 3. Here it is also worth noting that different selections of collocation points result in different approximations. Thus, a reasonable amount and spacement of collocation points are significant to reducing the approximation error.

#### 3. Initial Values for TDC Method

To start an iterative process, one has to give some values for initialization. In this paper, the Newton method is applied to solve the resulting TDC algebraic system. In view of the inherent drawback of this method that initial values have to be close enough to the precise solution, some already presented relative dynamic models are used to provide proper relative state vectors for initialization. Under different conditions, the C-W equations and the T-H equations are employed separately.

Before we move on to the following parts, it is necessary to give a brief explanation on the significance of a chief orbit in relative motion. In the literature, a variety of chief orbits with different shapes and sizes were investigated. It was found that a circular or near circular chief orbit brings much convenience to the study of relative motion, while the elliptic one complicates the analysis. When the perturbations are incorporated, the relative motion theories could be more complex, because there is no circular orbit or elliptic orbit under the perturbations and, thus, the osculating orbital elements can hardly describe the satellite orbits conveniently. In the work of Schaub and Alfriend [17], the invariant relative orbits were investigated using mean orbit elements. With the constraints derived in the paper, it reached a conclusion that the size of relative orbit changes as the inclined angle of chief orbit increases and could be extremely large as it approaches a polar orbit. It was also revealed in the work of Inalhan et al. [5] that eccentricity of chief orbit affects the shape of relative motion remarkably, not only on the in-plane motion, but also on the out-of-plane motion. Besides, the height of chief orbit affects the relative motion as well, in terms of the differential air drag. In this study, we focused on the effects of different eccentricities and inclined angles of chief orbit.

##### 3.1. The Clohessy-Wiltshire Equations

As a customary analysis approach of rendezvous mechanics, the C-W equations of relative motion are linearized and can be readily solved to give periodic relative motion. As is known, the constraint on initial conditions gives periodic relative motion, where , with being the constant radius of the chief orbit.

Applied in the exact relative model proposed by Xu and Wang, the preceding initial condition gives unbounded relative orbits that drift away. Suppose that the period of one orbit is time ; the initial values for TDC algebraic system can be chosen by collocating points equally spaced in a period on the unbounded orbits as shown in Figure 2.

Typically, the number of collocation points is chosen to be the same with that of harmonic coefficients in approximate function. Note that the frequencies , and are unknown in the above analysis. In order to solve (18), three more constraints are necessary. Herein, we impose a condition on the starting point; that is, , , and , are predetermined. The additional condition can also be described as follows: Remember that (18) is about the positions and velocities of collocation points. The above equations are further transformed as

With the collocation points settled and the frequencies assumed, a closed orbit can be obtained through the transformations in (1) and (3). Unsurprisingly, this orbit does not satisfy (18) at first. However, with the Newton method, or some other iteration methods, the position and velocity information of the collocation points as well as the values of frequencies can be corrected to minimize the residual error of (18). When the residual error caused by the collocation points is small enough, the corresponding orbit can be used to approximate the periodic relative motion.

Note that the C-W relative model is derived on the assumptions of small relative distances, a circular reference orbit, and a central gravitation field. As the eccentricity of the chief orbit grows, the drift rate of the relative orbits increases, and thus the C-W equations can no longer provide appropriate initial values when eccentricity is too large.

##### 3.2. Tschauner-Hempel Equations

The C-W equations of motion permit a simple analytic description of the relative model, yet it is limited to the circular reference orbit case. When elliptical orbit is considered, we can use the T-H equations instead.

Many researchers have paid efforts on deriving the solution of T-H equations. Various state transition matrices, both those using true anomaly and those using time as the independent variable, are developed. Among them, the work presented by Inalhan et al. [5] gives the periodicity condition for T-H equations.

When the chief is at perigee at the initial time, the bounded relative orbit condition expressed in the time domain is where , with being the semimajor radius of the chief orbit. As we can see, this expression is equivalent to the periodicity condition for C-W equation as .

In the paper of Inalhan et al. [5], (21) was further generalized to other values of ; thus, the initialization procedure can be carried out flexibly. In this paper, however, only the special case when the initialization occurs at apogee is considered for simplicity.

#### 4. Evaluation of the TDC Search Scheme

In practical missions, the maintenance of the satellite formation asks that the deputy satellite flies along a reference orbit. Given that the use of the slowly drifting trajectory obtained with TDC method as reference is impractical in long-term missions, we project it to a periodic closed orbit (PCO). Then a closed loop control strategy based on the projected orbit is introduced to simulate the fuel consumption of the chief/deputy formation flying.

##### 4.1. Projected Closed Orbit

Here, we introduce a straightforward approach to project the nearly bounded relative orbits found with the TDC searching scheme to a PCO.

Suppose the projected closed orbit can be expressed in a system like (15), and , , , , , and are forced to have the same frequency . Then collocate points, normally spaced equally in a time domain, on the slowly drifting trajectory. If the frequency of the projected orbit is predetermined, it would be quite simple to derive the expression of the projected orbit. All we should do is to establish a formula similar to (3) where is a block diagonal matrix , is the vector composed by the harmonic coefficients , and contains the state vectors of the collocation points .

Then multiply the equation with the pseudoinverse of ; we straightly get the least square solution , due to the property of the pseudoinverse operation. Thus, we can derive the projected orbit as , , with .

However, if is unknown, the procedure to determine the projected closed orbit will be laborious. Firstly, define a residual function . Using the least square method, we seek solutions and that minimize the square error . Thus, the formula has to be satisfied.

Substitute the residual function into (23); we get The detailed information for matrices , as well as the Jacobian matrix of (24) is provided in the appendix.

To solve this nonlinear algebraic equation, the globally optimal iterative algorithm (GOIA) is applied. The GOIA method is based on the concept of best decent vector which takes the form of . It is proposed to solve a system of nonlinear algebraic equations (NAEs) without inverting the Jacobian matrix at each step. For more details of GOIA, refer to Liu and Atluri [18].

##### 4.2. Closed Loop Control

The fuel consumption of the deputy satellite when it flies along the projected orbit is evaluated with a control technique called discrete linear quadratic regulator (DLQR). This technique is employed here due to its simplicity and ability to handle this sort of problems. Besides, the actuator is assumed to work impulsively, to approximate the way actual satellite thruster works.

Based on the assumption that the formation keeping is achieved by adjusting the motion of the deputy satellite only, the linearized relative motion dynamics can be obtained as where , is the state vector of the projected orbit, and is the deviation of the actual trajectory from the projected orbit. is given as follows: The detailed information of is provided in the appendix.

Then define the control vector with a sampling interval and a firing duration as follows: where is a sampling instant, , and is the control signal provided by the discrete controller, .

Thus, we get the discrete time model expressed as With , .

Using DLQR control approach, the quadratic performance index below is minimized: And the optimal control law can be achieved simply by the well-known Ricatti solution [19]

#### 5. Numerical Results and Analysis

In the following parts, the numerical simulations demonstrate the effectiveness of the TDC approach on searching for periodic relative orbits, as well as the projected closed orbits and the corresponding fuel consumption.

It is assumed that the initialization procedure occurs at the apogee of the chief orbit, where , km, and km/s, and the relative position of the deputy satellite is set to be km, km, and km. For simplicity, the angular momentum of the chief orbit at the apogee is approximated with that of the Keplerian orbit. In the following simulations, the number of harmonics included in the approximate functions is , while the number of collocation points is . The resulting relative orbits are shown in the LVLH frame.

Figure 3(a) shows the relative orbits obtained using the TDC searching scheme when the eccentricity and the inclined angle of the chief orbit take the values , . Figure 3(b) shows the drifting orbits generated by the periodicity condition of C-W equations. The initial relative orbit is always shown as a bold solid line with a different color to the remaining orbits and the plots obtained with TDC approach always show the data of 10 orbits. In this case, the collocation points are selected in the time period .

**(a)**

**(b)**

Noting that both simulations above are conducted under the same assumptions, the improvement of the result is drastic. In Figure 3(a) the relative orbit is nearly bounded and does not drift away after 10 orbits, while in Figure 3(b), which is plotted for 20 orbits, the substantial drift is apparent.

The following simulations further investigate the performance of TDC approach on searching for periodic relative orbit. The values of eccentricity and inclined angle are changed, so that the TDC searching scheme can be tested in situations where the periodicity condition provided by C-W equations or T-H equations becomes completely inapplicable.

In the case , , the TDC searching scheme gives a fairly satisfying result. is selected to be 7151 s.

Figure 4 makes a comparison between the orbits generated by the periodicity condition of T-H equations and TDC searching scheme. It is obvious that the drift of the orbits is relived greatly with the correction of the TDC method.

**(a)**

**(b)**

As the inclined angle increases up to , the resulting orbit is still nearly bounded and is demonstrated in Figure 5. s.

**(a)**

**(b)**

To illustrate the improvement of the result better, Figure 6 is plotted.

**(a) The initial orbits**

**(b) Enlarged partial view**

The orbit obtained by TDC approach is denoted by D1, which is shown with the solid red line, while the blue line marked with D2 is the orbit result from the T-H equations. Although D1 and D2 start from the same position (10 km, 10 km, and 10 km), the trajectories formed are quite different. In Figure 6(b), it is clearly indicated that the drift rate of orbit D1 is much lower than D2. Besides, the shape of D1 is bounded tightly, making it possible to establish a fuel saving control strategy.

Here it is worth noting that D1 and D2 may appear to be retrograde because of the perspective employed in Figure 6. Actually, the direction of D1 is consistent with D2 in the rough, although there indeed exists an acute angle between these two orbit planes. For further explanation, the intersection angle is caused by the slight differences between the initial velocities provided by the TDC searching scheme and the T-H equations, and the perspective is employed mainly for the purpose of demonstrating the drifting condition of D1 and D2.

In Figure 7, the eccentricity of the reference orbit takes the value , with . The “period” is 7191 s.

**(a)**

**(b)**

In this case, the periodicity condition arising from T-H equations is entirely impractical. As is shown in Figure 7(b), the drift rate on the along-track direction is far from satisfying due to the perturbation and the inclined angle. Nevertheless, the TDC searching scheme still gives a slowly drifting orbit. Although the drift rate gains as the number of orbit increases, the initial orbit and the following several orbits are nearly bounded and can still provide a satisfying approximation to the periodic relative motion.

The initial conditions obtained using T-H equations (IC1) and TDC method (IC2) are presented in Table 1. All of the initialization procedures are carried out at the same relative position (10 km, 10 km, and 10 km). In IC1, the velocities and take the value of zero, for that there is no constraint on and in the periodicity condition developed by Inalhan et al. However, it is clear from the results that these two velocities are no longer negligible when the gravitational perturbation effect is included.

Table 2 shows the approximate frequencies of periodic solutions. As is shown, , , and are slightly different; thus, the approximate solutions are quasiperiodic. Note that the knowledge of the perturbed frequencies due to aids in the design and tuning of the controllers for mitigating its adverse effects. The capability of approximating the in-plane and cross-track frequencies is an advantage of TDC searching scheme. This table also includes the iteration steps for solving the algebraic equations obtained by TDC method (i.e., (18) and (19)). It is shown that the iteration process is quickly convergent and, hence, costs little computer resource.

The drift of the relative orbits on three directions is shown in Figure 8 for the case , . The figures plot the data of 20 orbits, corresponding to roughly 2 days of simulation time. Overall, the drift rates of the first several orbits are small and acceptable. However, as time accumulates, the orbits will eventually drift away. Thus, to generate a reliable projected closed orbit for control, the projection should be carried out in a time period in which the drift is unobvious.

**(a)**

**(b)**

**(c)**

By looking into Figure 8 with some attention, it can be easily found that the drift on the along-track direction is more severe than the others. Thus, the number of orbits used to generate the PCO is often determined by the drift rate on this direction. Normally less than 10 orbits should be up to the requirement.

Figure 9 plots the projected orbits obtained under different situations, as well as the relative trajectories of the satellites which take the projected orbits as references and fly without control for 1 orbit. The tracking errors of different cases are listed in Table 3. As we can see, the tracking error grows as the eccentricity and the inclined angle increase. The reason for it is that the orbit of large eccentricity and inclined angle tends to drift away more quickly under the effect of the perturbations.

**(a)**

**(b)**

**(c)**

**(d)**

The number of points collocated to generate the PCOs is . As is shown in Figure 9, the PCOs (red line) magnificently preserve the characteristics of the nearly bounded orbits presented above. The information of both the in-plane motion and the out-of-plane motion is embodied. The relative trajectories of the uncontrolled satellites (blue line) indicate the efficiency of the PCOs. The tracking errors in some of the figures are not even visible on the scale shown. The frequencies of the PCOs are listed in Table 3 as well. They did not differ much from the original frequencies listed in Table 2. In the formation keeping process, they can help to design and tune the controller.

To further illustrate the efficiency of the projected closed orbit, we take the case where , , for instance. With the corresponding PCO as reference orbit, it is simulated for 40 orbits, almost 4 days, to plot the relative trajectory under LQR control. The sampling time is selected 10,000 s, 1000 s, 100 s, and 10 s separately for each subfigure below. As is presented in Figure 10, all the relative trajectories are bounded, tightly or loosely. In Figure 10(a), the sampling time is obviously too long to keep the relative trajectory close to the PCO, even though we have not lost control of it. So it is indicated explicitly that the PCO is a reliable reference for periodic relative motion control. The total simulated fuel consumption is denoted by Delta . It is shown in Table 4 along with the maximal tracking error.

**(a)**

**(b)**

**(c)**

**(d)**

Here Delta is calculated by the sum , where the magnitude is the norm of control signal , and is the working times of the thruster during the 40-orbit period of relative motion. In this simulation, the firing duration of the thrust pulses is supposed to be and always occurs at the beginning of the sampling time.

As can be seen in Table 4, the sampling time is important to the efficiency of LQR control. The shorter the sampling interval is, the tighter the control could be. As s, the maximal tracking error was reduced to 25.0111 meters, which is rarely satisfying remembering that the magnitude of the relative motion is m. Correspondingly, the fuel consumption Delta , however, did not change much. When s, the Delta needed to maintain the relative motion in a time period of almost 4 days was only 9.0868 m/s, roughly 0.227 m/s per orbit. It can also be easily found that the Delta needed in the cases of s and s is greater than that when sampling time is 10 s. From the numerical simulations it was revealed that there exists a balance between the sampling time and the magnitude Delta ; thus careful consideration must be given to tuning the sampling interval.

#### 6. Conclusions

In this paper, the TDC searching scheme is applied to search for periodic solutions of an exact nonlinear relative model. By assuming that the periodic solutions take the form of truncated Fourier expansion and doing some transformations, we derived the TDC algebraic system for the dynamic model. Solving this system with the Newton-Raphson method, we can obtain some initial positions and velocities for nearly bounded orbits as long as the collocation points for iterative process are selected properly. As is demonstrated above, the drift of the orbits generated by the initial conditions of C-W equations or T-H equations is dramatically reduced with the correction of the TDC method. Then we project these orbits to a closed orbit with a straightforward approach and use the projected orbit as a reference for relative motion control. The numerical results show that the control strategy is very effective and fuel saving.

Advantages of the TDC searching scheme include ability to handle complex equations with nonlinear terms, simplicity of implementation, repeatability of the solutions, and low cost of computational resource. Besides, the TDC approach can be generalized to search for periodic solutions of other dynamic models due to the inherent ability of TDC method to solve nonlinear problems. Actually, with a more exact dynamic model and better approximation of initial conditions for periodic relative motion, the results of TDC searching scheme can be more accurate. It is believed that this approach has potential applications on fuel saving relative orbit maintenance and spacecraft formation flying.

#### Appendix

The motion of the chief satellite can be described by a set of equations as follows:

To generate the projected orbits, matrices , and the Jacobian matrix are necessary. They are expressed as follows: with The matrix in (22) is expressed as with

As is a function of time, it changes along the reference trajectory. However, in the numerical simulation, it is reasonable to assume that is constant for simplicity. The simulation results demonstrate that the error introduced by the assumption is ignorable.

#### Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

#### Acknowledgments

This study is financially supported by the Doctorate Foundation of Northwestern Polytechnical University (CX201305), the Chinese NSF (11172235), and the Doctor Subject Foundation of the Ministry of Education of China (20106102110018).