Periodic Model Predictive Control for Tracking Halo Orbits in the Elliptic Restricted Three-Body Problem

A periodic model predictive control (MPC) scheme is proposed for tracking halo orbits. The problem is formulated and solved in the elliptic restricted three-body problem (ER3BP) setting. The reference trajectory to be tracked is designed by using eccentricity continuation techniques. The MPC design exploits the periodicity of the tracking model and guarantees exponential stability of the linearized closed-loop system, through a suitable choice of the terminal set and weight matrices. A sum-of-norms cost function is adopted to promote fuel saving. The proposed control scheme is validated on two simulated missions in the Earth–Moon system, which, respectively, involve station keeping on a halo orbit near the L1 Lagrange point and rendezvous to a halo orbit near the L2 Lagrange point. Results illustrate the advantage of designing the reference trajectory and the periodic control directly in the ER3BP setting versus approximate solutions based on the circular restricted three-body problem (CR3BP).


Periodic Model Predictive Control for Tracking
Halo Orbits in the Elliptic Restricted Three-Body Problem Renato Quartullo , Andrea Garulli , Fellow, IEEE, and Ilya Kolmanovsky , Fellow, IEEE Abstract-A periodic model predictive control (MPC) scheme is proposed for tracking halo orbits. The problem is formulated and solved in the elliptic restricted three-body problem (ER3BP) setting. The reference trajectory to be tracked is designed by using eccentricity continuation techniques. The MPC design exploits the periodicity of the tracking model and guarantees exponential stability of the linearized closed-loop system, through a suitable choice of the terminal set and weight matrices. A sumof-norms cost function is adopted to promote fuel saving. The proposed control scheme is validated on two simulated missions in the Earth-Moon system, which, respectively, involve station keeping on a halo orbit near the L1 Lagrange point and rendezvous to a halo orbit near the L2 Lagrange point. Results illustrate the advantage of designing the reference trajectory and the periodic control directly in the ER3BP setting versus approximate solutions based on the circular restricted three-body problem (CR3BP).
Index Terms-Aerospace control, control design, halo orbits, model predictive control (MPC), Space vehicles.

I. INTRODUCTION
L OW-THRUST station-keeping and orbital rendezvous in cis-lunar space play a key role for long-term solar system exploration missions as well as lunar landing [1]. In particular, parking orbits in the lunar vicinity are receiving increasing attention from several space agencies [2]. Near rectilinear halo orbits (NRHOs) are limit cycles typically found close to Lagrange points in the three-body problem of orbital mechanics. Thanks to their properties, halo orbits near L1 and L2 Lagrange points in the Earth-Moon system are deemed promising candidates for parking orbits in cis-lunar space missions. In particular, they benefit from the existence of Manuscript  Long-term station-keeping and trajectory design in the circular restricted three-body problem (CR3BP) have been extensively treated in the literature (see [6] and references therein). The CR3BP describes the motion of a satellite attracted by the gravitation of two massive bodies orbiting their center of mass on circular orbits and maintaining a constant distance between them during the motion. Choosing a rotating reference frame that keeps the position of the primaries fixed, the dynamics is represented by a set of autonomous ordinary differential equations, i.e., a time-invariant model; such a model has been extensively employed for station keeping [7], [8], [9] low-thrust trajectory design [10], [11], formation flight [12], and other applications.
However, the lunar orbit around the Earth is elliptic with a nonnegligible eccentricity (≃0.055). For this reason, the CR3BP represents only an approximation of the three-body problem for the Earth-Moon system. The elliptic restricted three-body problem (ER3BP) takes into account the eccentricity of the orbit of the primaries, and it is, therefore, a more accurate model than the CR3BP. However, the time-varying distance between the primaries renders the equations of motion nonautonomous; thus, the model is periodically time-variant. Furthermore, the generation of a reference halo orbit to be tracked becomes more involved. A number of methods have been developed for periodic orbit design in the ER3BP, see, e.g., [13], [14], [15], [16], [17], [18]. Control problems in the ER3BP setting have also been addressed [19], [20], although typically the reference trajectory to be tracked is designed in the CR3BP setting [21], [22].
In the last decade, model predictive control (MPC) has emerged as a promising technology for enhancing autonomy of the flight control systems in space applications [23]. The ability of MPC to handle state and input constraints and to optimize suitable performance indexes has made this technique attractive, especially for low-thrust operations and proximity maneuvers, see, e.g., [24], [25], [26], [27], [28]. Most popular MPC schemes are based on the minimization of a cost function which is quadratic in both state and input vectors. However, it has been observed that the use of alternative performance indexes may be convenient to achieve specific control requirements. In particular, sum-of-norms cost functions have been recognized to provide desirable properties in terms of control sparsity and fuel saving [27], [29], [30]. MPC can also be adapted to deal with inherently periodic systems or to track periodic references (see, e.g., [31], [32], [33] and references therein). In [34], an MPC strategy has been derived for periodic systems involving a sum-of-norms objective function. In this article, periodic MPC solutions based on the sum-ofnorms objective function are derived for tracking halo orbits.
In recent years, several MPC schemes have been proposed for problems involving halo orbits. The use of linear MPC for station keeping while tracking a halo orbit was considered in [35]. Nonlinear MPC is adopted in [36] for halo orbits in the Sun-Earth CR3BP. In [37], a quadratic MPC approach is proposed to stabilize a multirevolution halo orbit in the elliptic Sun-Mercury model. The problem is formulated directly in the ER3BP setting, but the periodicity of the model is not directly exploited and stability of the control scheme is not discussed. A chance constrained MPC for spacecraft rendezvous in NRHO has been proposed in [38], to ensure robustness with respect to probabilistic disturbances. In [39], a nonlinear continuous-time control law is coupled with a sampled-data MPC to perform station keeping of quasi-halo orbits near L2. The reference model is that of the CR3BP and the effect of eccentricity is treated as a disturbance.
In this article, a constrained stabilizing control law for halo orbit tracking in the ER3BP is presented. By exploiting the periodicity of the ER3BP model, a periodic MPC controller based on the application of the methodology proposed in [34] is developed to control a spacecraft involved in cis-lunar space missions. The novelty of the contribution with respect to the literature is twofold. Firstly, the control design problem is formulated as a periodic MPC with a sum-of-norms objective function, instead of the usually employed quadratic performance index. The control input is computed via convex optimization. The second key feature of the proposed approach is that the reference halo orbit is generated directly in the ER3BP setting via eccentricity continuation techniques [15], [17]. The resulting control scheme is validated on two simulated space missions, by employing a high-fidelity model based on nonlinear ER3BP spacecraft dynamics, affected by several disturbance sources, such as localization errors, thrust imperfections, and fourth body perturbation. Simulation results demonstrate the potential for successful application of the periodic MPC control scheme to station keeping and rendezvous on NRHOs. In particular, it is shown that formulating and solving the orbit-tracking problem directly in the ER3BP setting leads to a remarkable reduction in control effort, and thus fuel consumption, with respect to tracking a halo orbit designed under the CR3BP assumption. Moreover, a comparison with a classical quadratic MPC controller is presented, showing the advantages of adopting a sum-of-norms cost function in terms The article is organized as follows. In Section II, a description of the ER3BP equations is provided and a model suitable for the considered orbit-tracking problem is described. The control problem and the proposed MPC are presented in Section III. Section IV details the generation of the periodic orbit in ER3BP used as reference trajectory. The validation of the proposed method through numerical simulations is reported in Section V. Section VI contains some concluding remarks.

II. DYNAMIC MODEL
The general restricted three-body problem describes the motion of a body with negligible mass m 3 , under the gravitation attraction of two massive bodies m 1 and m 2 (with m 1 > m 2 ≫ m 3 ), namely the primaries, whose mass ratio is defined as ρ = m 2 /(m 1 + m 2 ). In the ER3BP, the primaries move in elliptic orbits, with eccentricity e and semi-major axis a, around their center of mass, according to Kepler's law. In this article, the particle represents a controlled satellite and the primaries are the Earth and the Moon.
The motion of the satellite is described in a rotating frame centered in the Earth-Moon center of gravity, where the position of the primaries is fixed on the x-axis (also known as syzygy-axis), the z-axis is normal to the Earth-Moon orbital plane, i.e., in direction of their angular momentum, and y-axis completes a right-handed triad. In this frame, the distances are instantaneously normalized, i.e., divided, by the primaries separating distance, which changes with the small primary true anomaly θ as where p = a(1 − e 2 ) is the semiparameter of the primaries' orbit. In this way, the position of the bodies is expressed in nondimensional units, the distance between Earth and Moon is equal to 1 and their coordinates are (−ρ, 0, 0) and (1−ρ, 0, 0), respectively. Since the normalization factor is not constant, the frame is also called pulsating [40].
Let r = [x, y, z] T be the nondimensional satellite position in the described rotating and pulsating frame. Its dynamics are then described by the following system of ordinary differential equations, where the independent variable is the true anomaly of the primaries θ: In (2), ( ) ′ indicates the differentiation with respect to θ is the pseudo-potential, andū x ,ū y , andū z are the nondimensional forcing accelerations. For completeness and in order to elucidate the definition and meaning of scaled inputs, we include the derivation of (2) in Appendix A. Fig. 1 depicts the rotating and pulsating frame with respect to an Earthcentered inertial (ECI) frame.
T be the state vector collecting the nondimensional position and velocity components of the satellite andū = [ū x ,ū y ,ū z ] T be the control input vector. In this notation, system (2) can be written as As shown in Appendix A, the nondimensional acceleration u is related to the actual (dimensional) acceleration exerted by the propulsion system u d by the equation where h = (G(m 1 + m 2 ) p) 1/2 is the magnitude of the angular momentum of the primaries and G is the universal gravitational constant. Note that the control law needs to ultimately govern the physical acceleration of the satellite; hence, we define and B(θ ) = (1/(1 + e cos θ) 3 )B, so that (4) becomes Since the right-hand side of system (7) is periodic in θ with period T = 2π, (7) is a nonlinear periodic system.

III. PERIODIC MPC
In this article, the control objective is to track a reference trajectory, representing a close periodic orbit, within the family of halo orbits [4]. In the following sections, a solution relying on the receding horizon periodic MPC [34] is developed. Before presenting the results, the following definition is given.
Let ξ r be an uncontrolled reference trajectory, obtained as an unforced solution of (7), so that By linearizing system (7) around ξ r , one obtains a linear time-varying system where x = ξ − ξ r represents the deviation from the reference trajectory .
In order to design a digital control scheme, system (9) is ZOH-discretized with sampling interval θ s . Thus, the resulting periodic discrete-time linear system has the form where k ∈ N and the matrix sequences A k ∈ R n×n and B k ∈ R n×m are N -periodic with period N = T /θ s . The use of a linearized model for control design has a twofold motivation. On one hand, typical halo orbit missions involve transfer to a neighborhood of the desired reference orbit, at which point the controller is engaged to stay on the halo orbit itself. In such a scenario, a linearized model provides a sufficiently accurate representation of the dynamics for such a controller. Moreover, the use of (10) as a prediction model allows one to formulate a computationally feasible MPC problem, which can be solved by standard convex optimization tools. In this work, in the controller design phase, we assume that a reliable state estimate is available at each time instant k from a localization system. This assumption is in line with the current state of art, see, e.g., [41], [42], while measurement noise will be considered during closed-loop simulations. Moreover, in order to meet recent mission technology requirements, the satellite is assumed to be equipped with a single electric propulsion system. In particular, maneuvering is achieved by firing a single thruster and steering the thrust vector via attitude control. Therefore, constraints on the maximum deliverable thrust can be modeled as ||u(k)|| ≤ 1 (a more general input constraint ∥u(k)∥ ≤ u max can be recast as ∥u(k)∥ ≤ 1 by scaling (10) by u max ). In this article, the attitude control problem is not addressed, thus the orientation of the thruster is assumed to be accurately realized during the design, while thrust errors will be considered during closed-loop simulations. This is a reasonable assumption in practice, because the attitude control authority has typically a much higher bandwidth than the translational one [43].
Let us consider the following optimization problem: where H is a given time horizon length,x k ( j) denotes the predicted state j steps ahead of k, and the decision variables are the elements of the control sequencê The objective function J (Û k ) is chosen as in which Q is a full-rank matrix, while W k+H and S k+H are full-rank matrices belonging to N -periodic sequences W k and S k = S T k , respectively. Matrix Q can be adjusted to trade-off tracking performance and fuel consumption. The proposed MPC design is based on the solution, at each time instant k ∈ N, of problem (11). Then, as is common in the receding horizon control, the first element of the optimal solution is applied to the system, i.e., Note that problem (11) is a second-order cone program (SOCP), thus its solution is computationally affordable with convex optimization tools [44].
It is worth stressing that the proposed MPC scheme is characterized by the cost (13), which is different from standard quadratic performance indexes, being instead a sum-of-norms of states and inputs. It has been observed that this choice is useful to promote control sparsity and fuel saving (see, e.g., [27], [29], [30]). Sum-of-norms MPC schemes have been studied for both time-invariant and periodic time-varying systems [30], [34]. Hereafter, their main theoretical properties are briefly recalled.
In order to ensure closed-loop stability of system (10) with the control law (11)- (15), it is crucial to suitably design the terminal set and the terminal cost, defined respectively by the N -periodic matrix sequences S k and W k . Due to the structure of matrices A k and B k , system (10) is stabilizable via N -periodic linear feedback. Hence, consider an auxiliary asymptotically stabilizing control law where the feedback gain K k ∈ R m×n is N -periodic and can be computed, for instance, by solving a periodic Riccati equation [45]. The resulting closed-loop system is given by which is clearly N -periodic. The following results establish the desired theoretical properties of the proposed MPC scheme. Proposition 1: Let S k = S T k ∈ R n×n be a N -periodic matrix sequence such that it is the solution of the following set of periodic linear matrix inequalities (LMIs): for k = 0, 1, . . . , N . Then, if problem (11) is feasible at time k 0 , then it is also feasible for all k > k 0 . Proof: See [34].
where λ m (D i ) denotes the minimum eigenvalue of the matrix D i . Then the proposed MPC scheme (11)-(15) with S k computed as in Proposition 1 and W k chosen as in (20), renders the origin of system (10) exponentially stable. Proof: See [34]. Proposition 1 provides a periodic matrix sequence S k ensuring that problem (11) is recursively feasible. Note that in (18), the constraint S 0 = S N is implicit from Definition 1. However, Proposition 2 selects the periodic matrix sequence W k in such a way that the cost J k decreases over time, thus guaranteeing a closed-loop stability. Matrix D k in (19) is a further degree of freedom, which can be used as a tuning parameter of the control design procedure.

IV. REFERENCE TRAJECTORY GENERATION
In order to generate the reference trajectory ξ r , one has to find a solution of the ER3BP (8). Due to the time-variance of the problem, generating a periodic orbit in the ER3BP involves more complications than in the CR3BP. The distance between the primaries changes accordingly to their true anomaly, thus Lagrange points have no fixed positions. This yields an asymmetry of the problem. For this reason, halo orbit generation in the ER3BP is more involved than in CR3BP. Furthermore, in CR3BP an orbit can achieve a period which is uncorrelated to the primaries periodicity. In ER3BP this is, in general, not possible, since the right-hand side of (2) is periodic with period 2π. However, designing the reference trajectory in the ER3BP leads to significant improvements in the control performance, with respect to trajectories designed in the CR3BP, as shown in Section V.
In this work, reference periodic halo orbits in the ER3BP are generated starting from halo orbits in the CR3BP, through differential corrections and eccentricity continuation techniques, similar to [15] and [17]. Let M S be the number of the satellite revolutions around a Lagrange point and M P the number of primaries' revolutions around their barycenter. The objective is to generate an orbit with a resonance ratio M S :M P . By exploiting the mirror theorem [46], a periodic orbit with period T C = 2(M P /M S )π is generated in the CR3BP through a single-shooting algorithm, see, e.g., [47]. In particular, the latter provides the initial condition ξ r c (0) = [x 0 , 0, z 0 , 0, y ′ 0 , 0] T such that, integrating system (2) with e = 0, a perpendicular and symmetric crossing of the normal plane occurs after half of the period. At this point it is sufficient to propagate the half orbit forward for the remaining half period to have the full periodic solution ξ r c . To generate a periodic solution in the ER3BP, the rationale is similar, i.e., the objective is to find a vector of initial conditions such that, integrating the unforced version of system (2), the resulting trajectory is a closed periodic orbit. A sufficient condition for this is stated in [15]: for an orbit to be periodic in the ER3BP, it is sufficient that it has two perpendicular crossings with either the normal plane or the syzygy-axis, or both of them, when the primaries are at the apse-line.
Let χ 0 = [x 0 , z 0 , y ′ 0 ] T be the vector of the free variables at θ = θ 0 , i.e., when the Moon is either at the periapsis (θ 0 = 0) or at apoapsis (θ 0 = π). The remaining three variables are equal to zero, in order to ensure the first perpendicular crossing of the normal plane. According to the periodicity criterion, the second crossing has been imposed at θ = θ 0 + π. This means that, after the integration of (2) for θ ∈ [θ 0 , θ 0 + π], the final state must satisfy the condition T e (χ 0 ) = [y π (χ 0 ), x ′ π (χ 0 ), z ′ π (χ 0 )] T = 0. At this point it is sufficient to propagate the half trajectory forward to θ = θ 0 +2M P π to obtain the full periodic solution. In this way all obtained periodic trajectories are M S :M P resonant orbits in the ER3BP with period T = 2M P π = M S T C , i.e., an integer multiple of the primaries' periodicity.
The condition T e (χ 0 ) = 0 can be enforced by using numerical root finding techniques. However, the reliability and quality of the solution are often dependent on the initialization of the solver. In other words, using directly ξ r c as an initial guess in the ER3BP with the actual Earth-Moon eccentricity usually does not produce acceptable results. For this reason, this problem is tackled iteratively with the eccentricity used as the continuation parameter. The eccentricity is gradually increased with a fixed small step δe until the actual primaries' eccentricity e is achieved. At each iteration j, a periodic orbit is generated in the ER3BP through the described procedure with the current eccentricity e j ∈ {0, δe, 2δe . . . , e} using the orbit computed at iteration j −1 as an initial guess for the solver. The procedure is outlined in Algorithm 1.

V. NUMERICAL RESULTS
In this section, the performance of the proposed MPC scheme is evaluated through numerical simulations.

A. Simulation Model and Control Design
The proposed Sum-of-Norms MPC scheme, hereafter referred as SoN-MPC, has been tested on a high-fidelity nonlinear model, including several sources of uncertainty. The nonlinear ER3BP dynamics (2) and (3) has been corrupted by measurement noise, thruster imperfections, and fourth body perturbation. In particular, position and velocity measurements provided by the localization module are affected by additive noise with standard deviations σ r = 10 km and σ v = 0.1 m/s, respectively. These values are in the same order as those considered in [41]. In order to account for thruster imperfections, an input disturbance with standard deviation σ u = 10 −7 m/s 2 , corresponding to 1 mN, is considered. Moreover, the missions are simulated under solar gravity perturbation. Despite its long distance to the Earth-Moon system, the Sun represents the most perturbing fourth body for the three-body problem. Its huge mass yields an additional acceleration term in the ER3BP model, whose derivation is detailed in Appendix B.
As far as the reference mission is concerned, the primaries orbit features are summarized in Table I. A servicing satellite with mass m 3 = 10 000 kg is assumed to be equipped with a fixed electric thruster capable of exerting a maximum force of 1 N, resulting in a maximum acceleration of 10 −4 m/s 2 . According to the values in Table I, the maximum nondimensional acceleration, scaled as in (6), is u max = 0.0364. The primaries' orbit is sampled assuming the sampling interval θ s = 0.0491 rad, corresponding to N = 128 samples.
The SoN-MPC scheme (11)-(15) is applied with a horizon H = N , corresponding to one orbit of the primaries. This choice of the prediction horizon is motivated by the fact that low-thrust propulsion systems involve low control authority, thus problem (11) may be infeasible for shorter horizons. Moreover, halo orbits are typically characterized by highly unstable regions [4]. Therefore, it is appropriate to include the whole orbit period in the predicted dynamics, so as to prevent too aggressive control actions. Nevertheless, despite the long prediction horizon, sufficient time to solve problem (11) is available, thanks to the large sampling time (in the order of hours).
The design of the control input is performed according to Propositions 1 and 2. The gain matrix K k in (16) is chosen as the solution to the standard periodic LQR problem [45], with state weighting matrix Q T Q, with Q used in the cost (13), and input weighting equal to the identity matrix. Matrices S k+H and W k+H are chosen as in (18), (19), and (20), respectively, with the matrix D k equal to the identity matrix for all k. The LMI problem (18) is solved by using CVX [48] and the commercial solver Mosek, capable to tackle semi-definite conic programming. Its solution required about 3 s on a standard laptop. However, note that the solution of (18) can be computed offline, therefore, its computational burden is not a key issue. The solution of the MPC problem (11) is carried out by using CVX and the commercial solver Gurobi. A single MPC problem instance is composed by m(H − 1) + n H optimization variables (corresponding to 1149 variables in our case study) and the computing time for its solution is in the order of 0.5 s.
To solve the reference orbit generation problem, i.e., to implement Algorithm 1, we iteratively invoked the built-in MATLAB function fsolve, which numerically finds the solution relying on Newton's method for each eccentricity step δe = 0.001. Examples of the generated reference trajectories are reported in Figs. 2 and 7. Note that the resulting revolutions of the third body are not overlapping, opposite to what is observed in the CR3BP. This is in line with the fact that in the ER3BP, it is not possible to achieve periodicities that are not integer multiples of 2π.
Under the above conditions, two space mission scenarios are simulated.
1) Station-Keeping: The operating satellite is required to actively track a halo orbit around the collinear Lagrange point L1 with a 3:1 resonance ratio. 2) Orbital Rendezvous: The servicing satellite has to approach an unperturbed reference point, assumed to lie on a halo orbit around the L2 Lagrange point with a 2:1 resonance ratio.

B. Station-Keeping
In this first mission, the objective is to show the ability of the SoN-MPC to keep the satellite on the reference trajectory, which is the 3:1 halo orbit around the Lagrange point L1, depicted in Fig. 2. Thus, the initial state is set as ξ (0) = ξ r (0). The state weighting matrix in cost function (13) is chosen as Q = I 6×6 . In Fig. 3, the dimensional components of the position and velocity errors between the satellite and the reference are reported, while Fig. 4 shows the 2-norm of the dimensional input profile. Note that the satellite is able to accurately track the halo orbit, keeping the tracking error bound in a small range, in spite of all the considered disturbances. The oscillatory behavior of the error is mainly due to the fourth body perturbation (which rapidly leads to divergence in the absence of control, due to strong instability  of the halo reference trajectory [5]). Fig. 5 shows the trajectory of the satellite in the inertial Earth centered frame. It can be observed that the tracked halo orbit accomplishes three revolutions around the Earth.
It is worth stressing the advantages of addressing the station-keeping problem in the ER3BP setting, with respect to adopting the CR3BP one, as usually done in the literature. To this aim, the SoN-MPC approach has been applied to maintain the reference trajectory generated from the CR3BP model (i.e., the red curve in Fig. 2). The tests have been performed on the same high-fidelity simulation model, with CR3BP-specific linearized model used for prediction in (11). Interesting results are obtained in regard to the fuel consumption for the entire mission. First, it has been observed that problem (11) turns out to be infeasible after few time steps, with the maximum thrust set to 1 N. Indeed, the small control authority is not sufficient to compensate the difference between the eccentricity assumed for the trajectory generation   and the one considered in the simulation. In order to assess the control effort needed to compensate such a discrepancy, the maximum deliverable thrust has been increased to 2 N, for both scenarios. In Fig. 6, the comparison between the two resulting command profiles is shown. Significant fuel saving can be observed by tracking the trajectory generated for the ER3BP. In fact, considering k ||u(k)|| as a fuel consumption indicator, tracking the trajectory designed in the ER3BP leads to a 64% fuel consumption reduction. This corroborates the benefits of designing periodic MPC scheme using the reference orbit generated with the ER3BP model in low-thrust missions.

C. Rendezvous
In this scenario, the servicing satellite aim is to reach another satellite already orbiting on a 2:1 halo orbit near the L2 Lagrange point (shown in Fig. 7). The initial state ξ (0)  is picked from a Gaussian distribution with mean ξ r (0) + (1/d(0))[2·10 6 , −1·10 6 , 0.5·10 6 , 0, 0, 0] T , which corresponds to an initial separating distance of 2291.3 km. The covariance matrix of ξ (0) is chosen to cover a variation of 150 km on each initial position component and 1 m/s on each initial velocity component. The state weighting matrix is set to Q = blockdiag{5 · I 3×3 , I 3×3 }. As a first experiment, the SoN-MPC scheme has been tested on the ER3BP dynamics with no disturbances. The result of a typical run is shown in Fig. 8. It can be observed that the tracking error goes to zero in finite time (approximately 8 days). This is a typical feature of MPC schemes adopting sum-of-norms cost functions, as opposite to classical quadratic ones.
A Monte Carlo set of 100 simulations, in which the servicing satellite starts at different random initial conditions, has been performed on the high-fidelity simulator, including all the disturbance sources. Fig. 9 shows the distance and velocity  errors during all the simulated maneuvers. It can be seen that the rendezvous between the two satellites is always achieved after about 12 days, corresponding to less than half a lunar orbit. Moreover, the dispersion of the simulated trajectories after the initial transient is negligible, confirming the inherent robustness of the adopted MPC scheme with respect to the considered disturbances (fourth-body perturbation, thrust error, measurement noise).
In order to assess the benefits of the sum-of-norms formulation of the MPC problem, the performance of SoN-MPC is compared to that of a periodic quadratic MPC (hereafter referred to as Q-MPC). This amounts to solving problem (11) with the following cost function: where the matrix sequence k is the solution of the periodic LQR [45]. The distance and velocity errors of the two approaches for one simulated maneuver are reported in Fig. 10,  while the time history of the norm of the corresponding input acceleration vector is shown in Fig. 11. It can be observed that SoN-MPC is superior in terms of tracking errors, while the periodic Q-MPC is prone to oscillations, which are mainly due to the fourth-body perturbation. The same behavior has been observed for the 3:1 halo orbit considered in the station-keeping maneuver of Section V-B.
To further evaluate the performance of the proposed approach, both SoN-MPC and Q-MPC have been tested with a state weighting matrix Q = diag{q 1 · I 3×3 , q 2 · I 3×3 }, for different values of q 1 and q 2 . For each pair (q 1 , q 2 ), a set of 50 simulations has been performed, with the same setting as those in Fig. 9. The root mean square error (RMSE) of the steady-state distance error (from the 15th day on) is reported in Table II. It can be seen that SoN-MPC yields much more precise tracking for all the considered values of (q 1 , q 2 ).

VI. CONCLUSION
A periodic MPC scheme has been developed for halo orbit stabilization and tracking. Its ability to accurately track a periodic orbit in the ER3BP has been demonstrated through numerical simulations which included several perturbation effects. Furthermore, the same control scheme has been shown to successfully complete a simulated rendezvous mission, thereby highlighting its applicability to different low-thrust cislunar scenarios.
Possible developments of the proposed approach concern the inclusion of state constraints dictated by mission requirements and the adoption of offset-free scheme to reject persistent disturbances such as the fourth body perturbation.
Alternative techniques for generating the reference trajectory may also be considered, an example being the cooperative dual-task space framework [49]. Adaptation to system variations over time or to actuator faults are other subjects for continuing research.

A. ER3BP Model With Forcing Input
Let R = [X, Y, Z ] T be the dimensional position of the third body in the rotational frame andṘ = [Ẋ ,Ẏ ,Ż ] T its timederivative. Remember that dimensional position is related to nondimensional one as R = d(θ)r.
The coordinate system rotates at rateθ about the z-axis, so that the angular velocity vector is ω = [0, 0,θ] T and the velocity vector can be written as The position of the spacecraft with respect to the primaries can be expressed by R 1 = [X + ρd, Y, Z ] T and R 2 = [X + (ρ − 1)d, Y, Z ] T .
Denoting the kinematic energy by K, the potential by U, and the Lagrangian by L, one obtains and The forced Euler-Lagrange equations, with the components of the position vector R as the generalized coordinates, reads d dt For the sake of brevity, let us derive only (26) for the component X that yields In order to simplify (27), a transformation to the rotating and pulsating frame defined in Section II is required. Such a transformation exploits the scaling in (1), a normalization of the time by the characteristic time t ⋆ = (G(m 1 + m 2 )/d(θ )) 1/2 and then a transformation of time derivatives into derivatives with respect to true anomaly, taking into account that According to these considerations, the relationships between the dimensional, time-dependent and the nondimensional, true anomaly-dependent velocities and accelerations arė X = h p (1 + e cos θ)x ′ + e sin θ x (28)Ẍ = h 2 p 3 (1 + e cos θ ) 2 (1 + e cos θ )x ′′ + e cos θ x . (29) Then, substituting (28) and (29) in (27), one obtains x ′′ − 2y ′ − 1 1 + e cos θ x = (1 − ρ)(ρ + x) where r 1 = ((x + ρ) 2 + y 2 + z 2 ) 1/2 and r 2 = ((x + ρ − 1) 2 + y 2 + z 2 ) 1/2 . Defining the pseudo-potential as in (3) and scalingū x as in (5), (30) can be rewritten as ) which is the first equation of (2). The other equations are derived in the same way.

B. Fourth Body Perturbation
In [50], a characterization of the solar gravity acceleration in the CR3BP has been given where ω 4 is the magnitude of the nondimensional angular velocity of the Sun as viewed in the Earth-Moon rotating frame and θ 4,0 is the initial Sun angular position. The angular velocity is computed as the difference between the nondimensional mean motion of the Sun in the inertial frame centered at the Earth-Moon barycenter, i.e., n 4 ((1 + ρ 4 )/||r 4 || 3 ) 1/2 , and the nondimensional mean motion of the Earth-Moon system with respect to the same observer, that is 1. The simulation model adopted for the numerical results is the same as in (2) where the modified pseudo-potential is * = + 4 .