Modified Chebyshev-Picard Iteration Methods for Station-Keeping of Translunar Halo Orbits

The halo orbits around the Earth-Moon L2 libration point provide a great candidate orbit for a lunar communication satellite, where the satellite remains above the horizon on the far side of the Moon being visible from the Earth at all times. Such orbits are generally unstable, and station-keeping strategies are required to control the satellite to remain close to the reference orbit. A recently developed Modified Chebyshev-Picard Iteration method is used to compute corrective maneuvers at discrete time intervals for station-keeping of halo orbit satellite, and several key parameters affecting the mission performance are analyzed through numerical simulations. Compared with previously published results, the presented method provides a computationally efficient stationkeeping approach which has a simple control structure that does not require weight turning and, most importantly, does not need state transition matrix or gradient information computation. The performance of the presented approach is shown to be comparable with published methods.


Introduction
For the spatial circular restricted three-body problem CR3BP , where the two large bodies move in planar circular orbits about their center of mass and a third body of negligible mass moves under only the 1/r 2 gravitational influence of the large bodies, there exist five stationary points in the rotating reference frame.These points, called the Lagrangian points or libration points, have unique roles for many scientific missions because they lie in the plane of the primary bodies' motion and have fixed positions with respect to these two bodies, similar to the geostationary orbits although a little more complicated.Among the five points, three of them are collinear with the other two bodies, and the motion near these points is unstable; the other two libration points with the primaries form equilateral triangles in the plane of motion of the two large bodies and the motions near these two points are neutrally stable.
There has been a lot of interest in finding and controlling the periodic orbits near the collinear libration points.Farquhar originally proposed using the orbits near the L 2 point of the Earth-Moon system for communication with the far side of the Moon 1 .By appropriate station-keeping and with a modest cost, the satellite is visible from the Earth all the time.However, such orbits are not generally periodic.Later, Farquhar and Kamel found that if the amplitude of the in-plane motion is large enough so that the nonlinear effects become significant, purely periodic three dimensional orbits, or halo orbits, exist that are permanently visible from the Earth 2 .The interest in the Earth-Moon halo orbits later shifted to the Earth-Sun system following the end of the Apollo missions and the launch of International Sun-Earth Explore-3 ISEE3 which is the first libration-point satellite 3 .
The techniques for station-keeping of halo orbits can be classified as 1 using continuous control, or 2 generating impulses at discrete time intervals.Breakwell et al. studied station-keeping for translunar communication station using a continuous feedback control to minimize a cost function, which is a weighted combination of position deviation from the nominal orbit and the control acceleration based on a linearization about the reference orbit 4 .Xin et al. developed a suboptimal closed-form continuous feedback controller by approximately solving a Hamiltonian-Jacobian-Bellman equation 5 .Kulkarni et al. applied the H ∞ approach for station-keeping of the halo orbit around the L 1 of the Earth-Sun system 6 .Cielaszyk and Wie modeled the nonlinearities of the halo orbit problem as persistent disturbances and applied disturbance accommodation with a linear quadratic regulator for the control 7 .Although using continuous control derived from optimal control theory usually generates orbits closer to their reference orbits than using impulses, discrete maneuvers have the advantages of less complexity and risk and higher precision orbit estimation can be obtained during coasting arcs.Dunham and Roberts described the stationkeeping strategies for three Earth-Sun libration point missions, all of which used discrete controls 8 .Two station-keeping techniques for the Earth-Moon libration point orbits were studied by G ómez et al. 9 .One is what they called the target point strategy where a cost function, which is a weighted summation of the position and velocity deviations from several future points and the fuel cost, is minimized to find corrective maneuvers.Another approach they studied is based on Floquet methods which combine invariant manifold theory and Floquet modes to only eliminate the unstable component of the error.Howell and Pernicka also studied using the target point strategy for the Earth-Sun system in a later paper 10 .The optimal spacing time between impulsively applied controls has been studied by Renault and Scheeres 11 , and the optimal time to update control law with continuous control has been recently studied by Gustafson and Scheeres 12 .
A newly developed numerical computation approach, Modified Chebyshev-Picard Iteration MCPI method, is used in this paper for station-keeping to maintain a halo orbit in the Earth-Moon system.This orbit provides an excellent parking orbit for a lunar communication satellite.Details about MCPI methods can be found in the papers by Bai and Junkins 13, 14 and Bai's dissertation 15 .Fusing Chebyshev polynomials with the classical Picard iteration method, MCPI methods iteratively refine an orthogonal function approximation of the entire state trajectory.MCPI methods can solve both initial value problems IVPs and two-point boundary value problems BVPs by constraining the coefficients of the Chebyshev polynomials.A unique characteristic of MCPI is that the Chebyshev coefficients are constrained linearly without approximation on each Picard iteration.As a consequence, shooting techniques to impose boundary conditions are not required.Although perhaps the most striking feature about MCPI methods is their naturally parallel structure because computation of the integrand along each path iteration can be rigorously distributed over many parallel cores with negligible cross communication needed, MCPI methods are computationally efficient even prior to parallelization according to the results reported by Bai and Junkins, 13, 14 where MCPI methods have been compared with several classical methods in solving IVPs and BVPs.In several examples, it has been demonstrated that both greater efficiency and accuracy can be obtained, compared to the pseudospectral method, for example.The MCPI method can be applied to either the impulsive or the continuous control case; in this paper, we will provide only the case of impulsive control.
This paper is structured as follows.We first present the numerical approach to generate the reference halo orbit.We use the third order analytical formulations derived by Richardson 16 to provide a starting guess for the halo orbit, and then use a differential-correction approach to numerically find the accurate initial conditions to generate the periodic halo orbit.The reference orbit is integrated using MCPI method, and the Chebyshev coefficients that precisely represent the nominal orbit are saved for future reference.We then use MCPI methods to generate impulses at fixed discrete time intervals for corrective maneuvers.Several important criteria such as fuel cost and deviations from the reference orbit are analyzed through simulations and we compare our results with previous published results.Conclusions and future directions are presented at the end.

Equations of Motion
As is well known, a good nominal reference orbit reduces the fuel cost for unnecessary maneuvers, so the dynamical model used to generate the reference orbit ideally should include the gravitational perturbation forces from the Sun and other planets as well as solar radiation and other forces.Farquhar and Kamel presented fairly detailed equations of motion where the effects of the Sun's gravitation force, solar pressure, and the Moon's orbital eccentricity were considered 2 .Gómez et al. 9 used JPL DE403 to include the gravitational influences of the Sun, Moon, and all nine planets.To illustrate our novel approach, we have chosen to only consider the CR3BP idealization of the Earth-Moon system.Because the station-keeping technique we present later is a numerical approach to solve for corrective maneuvers and is essentially an open loop control, using a more detailed dynamical model will not change the structure of the algorithm but may moderately affect the performance of the presented technique.For real mission operations, detailed dynamical models should be used to generate the reference orbit as well as the corrective maneuvers.
Most of the development of the equations of motion can be found in the book by Schaub and Junkins 17 .One difference here is that we use the L 2 Lagrangian point as the origin.Figure 1 shows how the rotating reference frame E : { e r , e θ , e 3 } is defined.m 1 , m 2 , and m represent the Earth, Moon, and spacecraft, respectively.Note the vectors denoting the location of the Earth, the Moon, and the center of mass are defined such that r 1 r 1 e r , which results in r 1 , r 2 , and r CM all being negative quantities since the vectors extend from the origin along the negative e r axis.The two central bodies rotate in circles about their center of mass at a constant angular velocity, defined by where r 12 is the distance between m 1 and m 2 .Thus, the E frame in use rotates at a constant velocity ω with respect to the inertial frame, as given by the following equation: ω ω e 3 .

2.2
The inertial position vector of the spacecraft is defined as follows: where r x , r y , and r z are the components of the r vector.The inertial derivative of this equation can be taken using the transport theorem 17 and 2.2 to produce the inertial velocity vector: ṙi ṙx e r ṙy e θ ṙz e 3 ω × r x − r CM e r r y e θ r z e 3 ṙx − r y ω e r ṙy r x − r CM ω e θ ṙz e 3 .

2.4
Finally, the above equation can be differentiated once more to produce the inertial acceleration vector: ri rx − 2 ṙy ω − r x − r CM ω 2 e r ry 2 ṙx ω − r y ω 2 e θ rz e 3 .

2.5
The massless force due to gravity from the two massive bodies is defined in E frame components as

2.6
Mathematical Problems in Engineering 5 with ξ i defined as follows: When combined with 2.5 , the equations of motion can be written as

2.8
The equations of motion are nondimensionalized by first using a new time variable, τ, defined as τ ωt, 2.9 and then using the distance between the Earth and L 2 point |r 2 | , which leads to the new distance variables as and lastly, the masses are eliminated in favor of a nondimensionalized mass ratio μ, defined as Now, rearrange 2.8 and divide by the distance r 2 and also by the factor ω 2 , yielding:

2.12
Replacing m 2 by μ m 1 m 2 and m 1 by 1 − μ m 1 m 2 using 2.11 , and introducing ρ i , the equations of motion become

2.13
Combining 2.13 with 2.1 yields the final, nondimensionalized equations of motion:

2.14
Using 2.15 , the nondimensionalized distance x, which is the distance from the center of mass to the L 2 point, is calculated: The other values, r 1 , r 2 , and r CM , are found using the center of mass equation:

2.16
Table 1 provides the values of the three parameters used in this model.

Finding Initial Conditions
A third-order closed-form solution for the equations of motion is given in Richardson's paper 16 .We use A y 45000 km which is defined as the amplitude of the linearized motion along the y direction in Richardson's paper 16 .This is the same magnitude used in the paper by Breakwell et al. 4 , at which position a satellite is visible from the Earth all the time.The approximate initial conditions obtained from the third order formula quickly diverge when propagated using the full equations of motion due to the inherent instability of the system.A differential corrections method 18 was used to find the initial conditions that lead to a bounded periodic orbit.In order to find a halo orbit one important property of halo orbits was exploited: halo orbits display symmetry about the x-z plane.This means that at the point where the orbit crosses that plane y 0 , the x and z components of velocity are both zero.
As such, choosing a point on the x-z plane as the initial conditions reduces the variable set from three position components x, y, z , three velocity components v x , v y , v z , and the orbit period T to two position components x, z , one velocity component v y , and the orbit period T .Additionally, after a single orbit, or after half an orbit, y, v x , and v y will all be zero again.Table 2 gives the initial conditions we found, where the physical time of the period is T p ≈ 14.5 days.Note that the solution accuracy in the numerical correction process to find these initial conditions has been set as 10 −15 .

Using MCPI Method to Save the Reference Orbit
A nice feature about MCPI methods is that the achieved solutions are represented to high precision in the form of orthogonal Chebyshev polynomials, thus if the coefficients of the states are saved, the state values at arbitrary time can be obtained straightforward by numerically computing the orthogonal polynomials and then a simple inner product matrix multiplication .We integrate the reference orbit using MCPI method for one orbit and save the coefficients of the position and velocity for future station-keeping reference.An important parameter to choose here is the order of the polynomials to use.Generally, higher order approximation leads to higher accuracy solutions but requires larger storage space to save all the coefficients.We use closure error after one period of the reference orbit, defined in 2.17 , as the accuracy criterion to find the optimal order to use: where S T p is the state values position or velocity after one complete orbit and S T 0 is the state values position or velocity at the initial time.where the closure errors of both position and velocity decrease exponentially as the order increases.Thirteen significant digits can be achieved with N ≥ 60, but this is higher precision than required for Earth-Moon halo orbits.In the simulations we present in the next section, we use polynomials of order 30, which still leads to a position accuracy about 0.1 km and velocity accuracy about 1 mm/s as proved by Figure 2 b .As we will discuss in the next section, this obtained accuracy is one magnitude smaller than the range of typical measurement noise.Of course this order can be tuned according to the accuracy requirement, and we simply illustrate the methodology.The reference orbit in three dimension is as well as the projection in the three orthogonal reference planes is shown in Figure 3.

Simulation Studies
We use MCPI method for station-keeping of the halo orbit generated in the last section.The matrix-vector form of MCPI methods for solving this special type of two-point BVPs, where the positions at the two boundaries are constrained, is shown in Figure 4.The matrices shown in the formulas, C x and C B α , are constant once the order of approximation polynomials are specified , and Θ xif is an N 1 × 1 vector depending on the order and the two-point boundary conditions.Detailed derivations of the approach can be found in the paper by Bai and Junkins 14 .Procedures to use the formula and definitions of the matrices and vectors are briefly discussed in the appendix of the present paper for completeness.
For the current problem, at the time when a maneuver is possible, the current position X t i x i , y i , z i and velocity V t i v x , v y , v z are measured.The reference position at the target time X t f x f , y f , z f is computed from the coefficients of the reference orbit through a matrix multiplication and the time can be arbitrary.Then MCPI method is used to solve this two-point BVP using the formula in Figure 4, and the desired velocity V r t i v x r , v y r , v z r at the maneuver time is obtained.The corrective maneuver is computed from As a numerical approach to solve BVPs, MCPI method can be easily adjusted if the dynamical model used for the real mission is more complicated than the model used in the preliminary design.MCPI method does not require neither gradient information nor the computation of state transition matrices, and a prior knowledge can be wisely used to provide a good initial guess to start the interaction, all of which make the method computational very attractive.
The position and velocity of the satellite when it is inserted to the halo orbit will be away from the ideal initial conditions of the reference orbit.The measurement of the position and velocity from Earth will not be perfect.And the implementation of the impulses will also have "execution errors".We include these three kinds of errors in the simulation.All the errors are assumed to be Gaussian distribution with a mean of zero.The standard deviations of these errors are chosen to be consistent with those used by G ómez et al. 9 and are listed as follows i tracking errors: σ x σ y σ z 1 km; σ v x σ v y σ v z 1 cm/s, ii insertion errors: σ x σ y σ z 1 km; σ v x σ v y σ v z 1 cm/s, iii maneuver errors: σ Δv x 0.05; σ Δv y 0.05; σ Δv z 0.02.We look at station-keeping of the halo orbit for 26 revolutions, which corresponds to 377 days, a little more than one year.100 Monte Carlo simulations are run and the results in Table 3 are the averaged numbers from the 100 trials.Here we have used 1/ 7T p as the spacing between two impulsive maneuvers which corresponds to a physical time about 2 days.Later, we will study how this parameter affects the results.Significantly, it might be surprising to notice that MCPI methods only took about three iterations to converge to such accurate solutions.This is because we have used the coefficients of the reference orbit to find the corresponding state values and then use these values to start the MCPI iteration.Thus MCPI method starts from an ideal trajectory then takes only three Picard iterations to converge to the neighboring solutions that restore the halo orbit to within the tolerance adopted.

Mathematical Problems in Engineering
Considering most of the thrusters used to generate impulses will have some limit on the minimum Δv it can provide; we added one parameter Δv min 2 cm/s as a threshold for the control: if the solved impulses are less than the threshold, no correction maneuver will be made.Table 4 shows the results for this setting.As expected, the deviations from the reference become larger and the maximum Δv become larger.The good thing is that both the overall Δv and maneuver numbers are reduced, which in part is because that the small maneuvers that might be in the same magnitude of the measurement noise are avoided.
We also increased this threshold to Δv min 3 cm/s.Unfortunately and not surprisingly, the deviations from the reference become larger and more Δv is required.We notice when we set Δv min 2 cm/s, we have a similar parameter setting as those studied by G ómez et al. 9 , for the case where they set the tracking and minimum maneuver intervals as two days.Compared with their results, our total impulses 1476.5 cm/s are larger than theirs 1111.8 cm/s , and we have more impulsive maneuvers 168 than theirs 80 ; however, our maximum tracking deviations 17.38 km are much smaller than their results 60.9 km and our maximum impulse 29.692 cm/s is also much less than theirs 73.554 cm/s .The reason we have more fuel required than theirs can result from several reasons.First, our reference halo orbit is different from their orbit.Second, they have used two targeting points and have minimized a cost function which is a weighted summation of both position deviations and fuel cost.As they commented in the paper, the weights have been tuned, which very likely emphasized a large penalty on the fuel cost.Instead, we have not optimized anything in our approach.This situation is similar as the one reported by Renault and Scheeres 11 , where the LQR control is more fuel efficient than targeting the stable manifold of the equilibrium point.Additionally, as shown in Figure 4, using the MCPI method only involves evaluating the differential equations on the discrete nodes and doing vector-matrix multiplication and summation.Of special interest for the current case, only three iterations are required for the method to converge.However, the approach by G ómez et al. 9 requires weight tuning beforehand and also a complicated state transition matrix derivation and then computation during the mission when a more realistic dynamical model is used.Thus, we believe our method is computationally very attractive, due to the absence of required tuning, and extremely simple and easy to implement.Figure 5 a shows the history of the magnitude of the impulses from one Monte Carlo simulation.Due to measurement noise and simulated control errors, the computed solution of the two-point BVP will of course contain errors that propagate into an error at the frequent maneuver time.These errors along 26 revolutions 377 days are shown in Figure 5 a .Note the numbers shown in Tables from 3 to 5 are the averaged numbers from 100 Monte Carlo simulations.For the 168 maneuvers along the same orbit of Figures 5 a and 5 b , Figure 5 c shows the number of Picard iterations required to achieve convergence.We conclude, 92.2% of the cases converge in three iterations and 100% converged in four or fewer iterations, a remarkable algorithm!Another important parameter is the time interval to solve for and then execute the maneuver.Usually there exists a minimum maneuver time because of the necessity for orbit determination process and the mission operation capability.On the other hand, if the maneuver intervals become too large, the deviations from the reference can become too large to be easily corrected.We have studied   the performance while changing the interval times from a half day to a little more than three and a half days.The results below show how the position and velocity deviations, number of maneuvers, total ΔV , and maximum ΔV change with respect to the interval time, see Figures 6 a to 6 d .As expected, Figures 6 a and 6 b show that the deviations increase and the number of maneuvers decrease as time interval between corrections becomes larger.It is interesting to see that there exists a broad window of one and half to three days for time intervals between corrections, that gives near-minimum fuel cost, as evident in Figure 6 c .This is similar as a situation reported by Renault and Scheeres 11 , where the optimal spacing for both LQR approach and targeting equisetum method exists for the problems they studied.From Figure 6 d , the fact that the maximum impulse is minimum for 1 to 1.5 days implies that separation is significant.

Approaches to Reduce the Fuel Cost
In this section, we study how to reduce the fuel cost using the presented method while maintaining its advantages such as no requirement for any gradient information, computational efficiency, and the final results in an orthogonal polynomial form.The above-method results in exact satisfaction of terminal boundary conditions on each subinterval if there are no control errors and measurement errors.As a result, the position errors are much smaller than those in the literature 9 , whereas the control required is about 24% higher.To gain some qualitative insight to the tradeoff between allowing increased terminal errors and the associated cost, the current method can be used in a heuristic fashion.We reasoned that applying less than 100% of the control impulse computed would result in increased terminal errors and smaller control, and in fact, we found that using 80% of the computed impulse in the current problem results in a control cost only 4% higher than those in reference 9 , while our terminal errors in fact still remain significantly smaller than those in the references.We emphasize that the relationships between the amount of the control to apply and the performance such as the total fuel cost, the maximum position deviation, and the mean position deviation are not linear.In Figure 7, we show how the total fuel cost, the maximum deviation, and the mean deviation change with respect to the introduced coefficient.The simulation uses the same conditions we studied for the results reported in Table 4. Notice all these three performance variables have been averaged from the 100 Monte Carlo simulations.Especially, we found, by applying 80% of the fuel required for solving the exact two-point boundary value problem, the total fuel cost is 1160 cm/s, the mean position deviation from the reference orbit is 4.84 km, and the maximum position deviation is 22.2 km.Thus with a fuel cost only 4% larger than the reported results 1111.8 cm/s 9 , our position deviations from the reference orbit are much smaller than the reported results, for which the mean and maximum position deviations are 61 km and 74 km, respectively.
Another approach to reduce the fuel cost is to adopt a hybrid propulsion system which uses low thrust during the flight and impulses at some specified points.The hybrid propulsion concept has been studied by Bai et al. for an Earth to Apophis mission design 19 .Since the presented methods can also solve optimal control problems for continuous system 14 , utilizing this hybrid propulsion system is straightforward.We could also reduce the fuel cost by varying the time interval used for the station-keeping instead of using the current constant time interval.However, this will increase the computation requirements, since all the matrices used in the current approach are constant and are precomputed, which have to be recomputed if the time interval changes.We will also explore to efficiently utilize the dynamic information of some specified orbits to reduce the fuel cost in the future study.

Conclusions
Modified Chebyshev-Picard Iteration method is used for station-keeping of a halo orbit around the L 2 libration point of the Earth-Moon system.Compared with other stationkeeping techniques which also generate impulses at discrete times, the presented method is much simpler, since it does not need to compute state transition matrix or gradient information, which can be complicated for the full dynamical model and might be impossible to implement for the on-board system.Previous studies 13, 14 have shown that the MCPI algorithm works well for highly perturbed orbits; it is anticipated that little change will be required in the present algorithm upon including solar and other non-CR3BP effects.There is also no tuning process e.g., for the weight matrices which is required for other approaches.
Although the fuel cost from this approach, for the examples presented herein, might be a little higher than other optimal control results, the simple structure of the method as well as its salient computation efficiency makes this technique very attractive.The method can be used for station-keeping other orbits with the only change being the use of the appropriate  equations of motion.We believe this method can take the same role as the classical Battin's form of the Lambert algorithm for solving orbit transfer problems, with the differences that this method is a numerical approach and can be used for rather general nonlinear systems, including perturbations.We are investigating extending the current approach which solves two-point boundary value problems to solve optimal control problems which include several target points as well as hybrid propulsion mode with both impulsive and continuous thrust.
and the boundary conditions are x t t 0 x 0 and x t t f x f .The first step of MCPI methods is to transform the generic independent variable t to a new variable τ, which is defined on the valid range, the closed interval −1, 1 , of Chebyshev polynomials Introducing this time transformation of A.2 into A.1 , it is rewritten as The position update equation is where X x τ 0 , x τ 1 , x τ 2 , . . ., x τ N T is the vector representing the position trajectory evaluated at all the N 1 Chebyshev-Gauss-Lobatto CGL nodes, which are computed through τ j cos jπ N , j 0, 1, 2, . . ., N. A.5 The components of C x ≡ TW and C B α S B TV are defined as below, and f f τ 0 , f τ 1 , f τ 2 , . . ., f τ N T .Θ xif is defined by the boundary condition as Θ xif x 0 x f , x f − x 0 /2, 0, 0, . . ., 0 T ∈ R N 1 .Notice although the second order formula in A.4 only solves for the position, the velocity solution can be computed conveniently by taking the derivative of the position, which has been obtained in a form of Chebyshev polynomials: and the two diagonal matrices W and V are defined as

A.7
The definition of S B as shown below is more complicated than others where S i, j represents the element in the ith row and jth column.

Figure 4 :
Figure 4: Vector-matrix form of MCPI method for solving second order bvps.
Number of iterations from one Monte Carlo simulation Impulse time (1/7T p ) c Iteration numbers

Figure 5 :
Figure 5: Results from One Monte Carlo Simulation.

Figure 6 :
Figure 6: Effects of the spacing of maneuvers averaged from 100 Monte Carlo simulation .

8
Averaged mean deviation versus coefficient of control c Averaged mean deviation

Figure 7 :
Figure 7: Performance versus coefficient of control.

Table 3 :
Baseline case results.