Suboptimal Nonlinear Model Predictive Control Strategies for Tracking Near Rectilinear Halo Orbits

1 Introduction

2 System Modelling

3 Control Design

4 Simulation Results

5 Conclusion

References

2.1 The Elliptical Restricted Three Body Problem (ER3BP)

2.2 Halo Orbits

3.1 Prediction Model

3.2 Optimal Control Problem Formulation

3.3 Controller Implementation

Abstract

Remark 1.

\PaperNumber

-XX-XX

Andrew W. Berning Jr Thanks: PhD Candidate, Aerospace Engineering, University of Michigan, 1221 Beal Avenue, Ann Arbor, MI 48109. Dominic Liao-McPherson Thanks: Research Fellow, Aerospace Engineering, University of Michigan, 1221 Beal Avenue, Ann Arbor, MI 48109. Anouck Girard Thanks: Associate Professor, Aerospace Engineering, University of Michigan, 1221 Beal Avenue, Ann Arbor, MI 48109. and Ilya Kolmanovsky Thanks: Professor, Aerospace Engineering, University of Michigan, 1221 Beal Avenue, Ann Arbor, MI 48109. Note: This research is supported by the National Science Foundation Award Number CMMI 1562209.

Near Rectilinear Halo Orbits (NRHOs), a subclass of halo orbits around the L1 and L2 Lagrange points, are promising candidates for future lunar gateways in cis-lunar space and as staging orbits for lunar missions. Closed-loop control is beneficial to compensate orbital perturbations and potential instabilities while maintaining spacecraft on NRHOs and performing relative motion maneuvers. This paper investigates the use of nonlinear model predictive control (NMPC) coupled with low-thrust actuators for station-keeping on NRHOs. It is demonstrated through numerical simulations that NMPC is able to stabilize a spacecraft to a reference orbit and handle control constraints. Further, it is shown that the computational burden of NMPC can be managed using specialized optimization routines and suboptimal approaches without jeopardizing closed-loop performance.

It has been more than 40 years since a human last traveled beyond low earth orbit (LEO) on the Apollo 17 spacecraft in 1972. Recently, there has been renewed interest in human exploration of the solar system. In particular, cis-lunar space has emerged as a focus area as evidenced by the Global Exploration Roadmap [1, 2] and the NASA Artemis program[3], including the 2024 lunar landing goal[4]. Operations in cis-lunar space will support space-based facilities for robotic and human missions to the Moon and eventually to destinations such as asteroids or Mars.

Near Rectilinear Halo Orbits (NRHOs) have emerged as promising candidates for a long term lunar gateway module and/or as staging orbits between LEO and low lunar orbit[5, 6]. NRHOs are limit cycles near the co-linear L1 and L2 Lagrange points. Certain NHROs possess several useful properties including the existence of low-energy transfer orbits [7], good stability characteristics, unobstructed views of Earth, and favorable resonance properties that allow them to avoid eclipses [8]. They were first identified in the Circular Restricted Three Body Problem (CR3BP) and are periodic natural motion trajectories in the simplified setting of the CR3BP[9].

While certain NRHOs are stable or nearly stable in restricted three body systems[5], in reality maintaining the spacecraft on them is challenging due to perturbations caused by gravitational forces from other celestial bodies, navigational errors, hardware limitations, solar radiation and magnetic forces. This provides the motivation for the development of active station-keeping algorithms that can stabilize the orbit despite these disturbances. Ideally, these algorithms should optimally balance tracking error with propellant/energy consumption. Computing actions/policies offline and uploading them reduces onboard computing requirements but may reduce robustness to disturbances. Alternatively, computing actions online can enhance the ability to react to changes and improve robustness. In particular, Model Predictive Control (MPC) is a promising methodology for online optimal control that can systematically account for constraints, nonlinearities, and both trajectory tracking and fuel minimization requirements[10, 11, 12]. However, online optimization approaches can be expensive from a computational or power consumption perspective.

There is a growing body of literature on station-keeping for NRHOs. Dynamical systems theory based techniques using Cauchy-Green Tensors and X-Axis Crossing methods have, in particular, been investigated [8, 13]. Set-invariance and analytical minimum energy-based solutions to a linearized problem that account for measurement uncertainty have been proposed[14], and linear MPC and linear quadratic regulator (LQR) approaches have been developed[15]. There is also interest in station-keeping for other halo orbits about the L1 and L2 points. The Optimal Continuation Strategies method, based on 2-point boundary value problems, was validated in-flight during NASA’s ARTEMIS mission[16] which flew longer period halo orbits around the L1 and L2 points. Discrete time sliding mode control has been applied to station-keeping on Halo and Lissajous orbits around L2 [17]. Analytical methods for station-keeping on Halo orbits in the CR3BP using the continuous-time linear quadratic regulator[18] (LQR), Floquet theory[19], and invariant manifolds[20] have also been investigated.

In this paper, we investigate the use of nonlinear model predictive control (NMPC) for station keeping on NRHOs using a low-thrust propulsion system. We show that NMPC is able to successfully maintain spacecraft flight along a halo orbit, satisfy thrust constraints, and demonstrates a high degree of robustness to disturbances and model mismatch. Moreover, while low-thrust actuators are well suited for long duration station-keeping, they require more frequent control updates, leading to more demanding onboard computing requirements. As such, we leverage new optimization algorithms[21] and suboptimal MPC[22], to demonstrate that it is possible to meet closed-loop performance requirements using little computational power. Our approach contrasts with existing literature which often assumes impulsive thrusters and uses control methodologies based on linearized models. By demonstrating the computational feasibility of NMPC for this problem, we open the door to future research exploiting the flexibility, i.e., nonlinear cost functions and constraints, of NMPC to optimize high-level objectives such as minimizing propellant consumption.

When considering spaceflight in the cis-lunar flight regime, it is most natural to consider the restricted three body problem in which the third body is of negligible mass compared to the two primary bodies (as is the case for a spacecraft, the moon, and the earth)[23]. The two restricted three body problems considered in this work are the circular restricted three body problem (for reference trajectory generation and control prediction model), in which the primary and secondary bodies are assumed to travel in circular orbits about their barycenter, and the elliptical restricted three body problem (for simulation and validation), in which the orbits of the primary and secondary bodies are assumed to have nonzero eccentricities.

For both systems of equations, described below, the primary and secondary bodies lie along the $x$ axis with the barycenter at the origin. The $z$ axis points in the direction of angular momentum of the primary-secondary system, and the $y$ axis completes the orthogonal frame.

The equations of motion of the spacecraft in the pulsating frame of ER3BP and in non-dimensionalized distance and time units are given by:


$\displaystyle x^{\prime\prime}-2y^{\prime}$	$\displaystyle=\frac{1}{1+e\cos(\theta)}\frac{\partial U}{\partial x}+u_{x}$	(1a)
$\displaystyle y^{\prime\prime}+2x^{\prime}$	$\displaystyle=\frac{1}{1+e\cos(\theta)}\frac{\partial U}{\partial y}+u_{y}$	(1b)
$\displaystyle z^{\prime\prime}+z$	$\displaystyle=\frac{1}{1+e\cos(\theta)}\frac{\partial U}{\partial z}+u_{z}$	(1c)

where $(\cdot)^{\prime}$ denotes differentiation with respect to the true anomaly of the primaries $\theta$ , $u=(u_{x},u_{y},u_{z})$ are the control accelerations provided by thrusters and

U(x,y,z)=\frac{1}{2}(x^{2}+y^{2})+\frac{1-\mu}{\|(x+\mu,y,z)\|_{2}}+\frac{\mu}{\|(x-1+\mu,y,z)\|_{2}},

(2)

is the pseudo-potential function, shown in Figure 1. Note that the independent variable in (1) is the true anomaly which is related to time, $t$ by

t^{\prime}=\frac{1}{(1+e\cos(\theta))^{2}}.\\

(3)

Refer to caption — Figure 1: Contour plot of the pseudo-potential (2) with Lagrangian points L1 and L2 marked.

In the absence of eccentricity, i.e., with $e=0$ , the ER3BP reduces to the circular restricted three body problem (CR3BP):


$\displaystyle x^{\prime\prime}-2y^{\prime}=$	$\displaystyle~\frac{\partial U}{\partial x}+u_{x}$	(4a)
$\displaystyle y^{\prime\prime}+2x^{\prime}=$	$\displaystyle\frac{\partial U}{\partial y}+u_{y},$	(4b)
$\displaystyle z^{\prime\prime}+z=$	$\displaystyle\frac{\partial U}{\partial z}+u_{z}.$	(4c)

Note that in the circular case, the rate of change of time with respect to true anomaly is unity and so, barring wrap-around issues ( $\theta\in[0,2\pi]$ ) , $t=\theta$ .

The equations of motion can be simplified by introducing the velocity $v=(x^{\prime},y^{\prime},z^{\prime})$ and position $r=(x,y,z)$ vectors. Then (1) becomes

	$\begin{bmatrix}r^{\prime}\\ v^{\prime}\end{bmatrix}=\begin{bmatrix}0&I\\ A_{21}&A_{22}\end{bmatrix}\begin{bmatrix}r\\ v\end{bmatrix}+\begin{bmatrix}0\\ I\end{bmatrix}\left(\frac{\nabla U(r)}{1+e\cos(\theta)}+u\right)$	(5a)
where
	$A_{21}=\begin{bmatrix}0&0&0\\ 0&0&0\\ 0&0&1\end{bmatrix},\quad A_{22}=\begin{bmatrix}0&2&0\\ -2&0&0\\ 0&0&0\end{bmatrix}.$	(5b)

Finally, introducing the state vector $\xi=(r,v)$ , (5) can be written compactly as

\xi^{\prime}=f_{c}(\theta,\xi,u,e).

(6)

Halo orbits are families of periodic orbits near the collinear Lagrange points in the CR3BP[9]. Members of these orbital families that occur closer to the secondary body begin to exhibit properties of NRHOs, traveling very nearly in a plane normal to the orbital plane of the primaries.

When considering the construction of these orbits in CR3BP we exploit the following property: the velocities $x^{\prime}$ and $z^{\prime}$ are equal to zero when the orbit crosses the $x\mbox{-}z$ plane ( $y=0$ ) at apoapsis and periapsis. Additionally, orbits are symmetric across the $x\mbox{-}z$ plane, so if we assume initial conditions of $X_{0}=[x_{0},0,z_{0},0,y^{\prime}_{0},0]^{\rm T}$ and integrate until the trajectory again crosses the $x\mbox{-}z$ plane at $X_{f}=[x_{f},0,z_{f},x^{\prime}_{f},y^{\prime}_{f},z^{\prime}_{f}]^{\rm T}$ , we only need to enforce $x^{\prime}_{f}=z^{\prime}_{f}=0$ to ensure a periodic orbit. In this work, a single-shooting method is used to solve this problem, iterating on the initial conditions $z_{0}$ and $y^{\prime}_{0}$ for a given $x_{0}$ until the condition $x^{\prime}_{f}=z^{\prime}_{f}=0$ is met to within some specified tolerance. Finally, in an attempt to reduce the discontinuity when a controller is tracking the reference trajectory over multiple orbital periods, the state at each time instant is shifted linearly such that the shift at $X_{0}$ is zero and the shift at $X_{f}$ is such that $X_{f}=X_{0}$ . This results in a reference trajectory that is no longer a natural motion trajectory, but it eliminates the aforementioned discontinuity. A family of periodic orbits about L2 is shown in Figure 2.

The single-shooting method utilized in this work is limited to finding periodic orbits in system (1) only for $e=0$ . Future work will include the use of a more sophisticated multiple-shooting and homotopy methods to enable finding periodic reference trajectories in ER3BP for $e\neq 0$ .

This section describes the proposed controller. An MPC controller has three constituent components, a prediction model used to evaluate the impact of control actions, an optimal control problem (OCP) formulation encapsulating the control problem, and an implementation strategy for solving the OCPs online.

MPC uses a prediction model to estimate the response of the system to control actions. Such a model is typically of the form

\xi_{k+1}=f(\theta_{k},\xi_{k},u_{k}),

(7)

where $\xi_{k}=\xi(\theta_{k})$ and $u_{k}=u(\theta_{k})$ . We derive a suitable prediction model in this form by discretizing the CR3BP, i.e., (5) with $e=0$ . The control signal is held constant over each interval $[\theta_{k},\theta_{k+1}]$ so the exact model is given by the state transition equation

\xi_{k+1}=\phi(\theta_{k},\xi_{k},u_{k})=\xi_{k}+\int_{\theta_{k}}^{\theta_{k+1}}f_{c}(\theta,\xi(\theta),u_{k},0)~d\theta.

(8)

We approximate (8) numerically using a standard 4^th order Runge-Kutta (RK4) method:

	$\displaystyle\xi_{k+1}=f(\theta_{k},\xi_{k},u_{k})=\xi_{k}+\frac{\Delta\theta}{\\|b\\|_{1}}\sum_{i=1}^{4}b_{i}~k_{i}(\theta_{k},\xi_{k},u_{k}),$		(9)
	$\displaystyle k_{i}=f_{c}(\theta_{k}+\Delta\theta~c_{i},\xi_{k}+\Delta\theta~a_{i}~k_{i-1},u_{k},0),~k_{0}=0,$

with a uniform step size $\Delta\theta>0$ and where $a=(0,0.5,0.5,1)$ , $b=(1,2,2,1)$ , and $c=(0,0.5,0.5,1)$ are the coefficients of the method.

Throughout this paper, the CR3BP is used as the prediction model.

In MPC the feedback law is defined implicitly through the solution of a receding horizon optimal control problem (OCP). Let $N>0$ be the length of the prediction horizon, $i$ be the index along the prediction horizon, and $k$ be the discrete true anomaly index. We use the notation $\xi_{i|k}$ to denote the predicted state $i$ steps into the prediction horizon at the sampling instant $k$ and denote the planned control actions $u_{i|k}$ analogously. The OCP formulation is then

where $\theta_{i|k}=\theta_{k}+i\Delta\theta$ , $\tilde{\xi}=(\xi_{0|k},\ldots,\xi_{N|k})$ and $\tilde{u}=(u_{0|k},\ldots,u_{N-1|k})$ are the optimization variables, $\bar{\xi}_{i|k}$ is the desired state, $u_{max}>0$ is an upper bound on the input acceleration, $f$ is defined in (9), and $Q=Q^{T}\succ 0$ and $R=R^{T}\succ 0$ are weighting matrices. The control input is then $u_{k}=u^{*}_{0|k}$ where $(\cdot)^{*}$ denotes a minimizer of (10).

An efficient method for solving (10) is essential for implementing NMPC. In this paper, we use a suboptimal variant of Sequential Quadratic Programming (SQP) algorithm that exploits the time sequential structure of MPC. The OCP (10) can be written compactly as the nonlinear programming problem

where $w=(\tilde{\xi},\tilde{u})$ . The NLP is solved using the sequential quadratic programming iteration

w_{i+1|k}=w_{i|k}+d^{*},\quad\lambda_{i+1|k}=\lambda^{*},\quad v_{i+1|k}=v^{*},

(12)

where $\lambda$ and $v$ are dual variables associated with the equality and inequality constraints, respectively, $z_{i|k}=(w_{i|k},\lambda_{i|k},v_{i|k})$ is the solution estimate at the sampling instant $\theta_{k}$ after $i$ SQP iterations and $(d^{*},\lambda^{*},v^{*})$ is the primal-dual solution to the following quadratic program (QP)

In time-distributed optimization (TDO), the SQP algorithm is limited to $\ell>0$ iterations and the final solution guess from previous timestep is used to warmstart the SQP algorithm, i.e., $z_{0|k+1}=z_{\ell|k}$ . This leads to a coupled plant-optimizer system as shown in Figure 3. Under appropriate assumptions it is possible to show that TDO based MPC recovers the stability and robustness properties of optimal MPC using a finite number of iterations[22].

Figure 3: A comparison between optimal MPC, which is a static feedback law, and time-distributed MPC, which is a dynamic compensator. The operator

\mathcal{T}_{M}

represents a fixed number of SQP iterations and

\Xi

selects the control input from the full solution.

Analytic derivatives, potentially coupled with symbolic optimization[24], are important for accelerating the SQP routine. The cost function and the inequality constraints are quadratic and linear functions so their derivatives are readily available. The derivative of the equality constraints, i.e., of the prediction model, are more complicated due to the RK4 integration scheme used in (9). Algorithm 1 evaluates the sensitivities $\nabla_{\xi}f$ and $\nabla_{u}f$ . Finally, the quadratic programming algorithm used to solve (13) significantly influences the overall computation time. We use the FBstab method[21] which can exploit the sparsity structure of (13) and be easily warmstarted between SQP iterations.

Algorithm 1 Sensitivities of RK4 Integration

\Delta\theta

\xi

u

\theta

f(\theta,\xi,u)

\nabla_{\xi}f(\theta,\xi,u)

\nabla_{u}f(\theta,\xi,u)

k_{0}=0

A_{0}=0

B_{0}=0

4: for

i=1,\ldots 4

\theta_{i}=\theta+c_{i}\Delta\theta

\xi_{i}=\xi+a_{i}\Delta\theta k_{i-1}

k_{i}=f_{c}(\theta_{i},\xi_{i},u)

A_{i}=\nabla_{\xi}f_{c}(\theta_{i},\xi_{i},u)\left[I+a_{i}\Delta\theta A_{i-1}\right]

B_{i}=\nabla_{\xi}f_{c}(\theta_{i},\xi_{i},u)\left[a_{i}\Delta\theta B_{i-1}\right]+\nabla_{u}f_{c}(\theta_{i},\xi_{i},u)

10: end for

11:

f(\theta,\xi,u)=\xi+\frac{\Delta\theta}{\|b\|_{1}}\sum_{i=0}^{4}b_{i}k_{i}

12:

\nabla_{\xi}f(\theta,\xi,u)=I+\frac{\Delta\theta}{\|b\|_{1}}\sum_{i=0}^{4}b_{i}A_{i}

13:

\nabla_{u}f(\theta,\xi,u)=\frac{\Delta\theta}{\|b\|_{1}}\sum_{i=0}^{4}b_{i}B_{i}

The plant and controller are both implemented in Simulink 2019b MATLAB function blocks and compiled into C code. We use a MATLAB implementation of FBstab⁰⁰ 0 https://github.com/dliaomcp/fbstab-matlab that is compatible with automatic code generation.

The specific parameters used for simulations in this work are presented in Table 1. Simulations are performed under white Gaussian process noise with zero mean and variance $\tau$ . The thrust limits and mass used in simulation are in line with expected capabilities of the Advanced Electric Propulsion System [25].

$m$	$10,000~kg$	Spacecraft mass
$N$	$35$	OCP horizon length
$M$	$3$	Maximum number of SQP iterations
$\Delta\theta$	$0.01~rad$	Discrete step size
$\mu$	$0.012$	Mass ratio
$e$	$0.055$	Eccentricity
$\tau$	$10^{-10}~\frac{km}{s^{2}}$	Process noise variance
$Q$	$10^{3}\begin{bmatrix}10\mathbb{I}_{3}&0_{3}\\ 0_{3}&\mathbb{I}_{3}\end{bmatrix}$	Cost function weighting
$R$	$\mathbb{I}_{3}$	Cost function weighting
$u_{max}$	$2000~mN$	Max control constraint

Table 1: Simulation parameters

Figures 4–7 show the results of simulating the CR3BP, i.e., (1) with $e=0$ , in closed-loop with the NMPC controller described in the Control Design Section. In this case the reference trajectory is a periodic natural motion trajectory of the CR3BP so the spacecraft is able to closely track the reference using little control input.

Figure 4: State trajectories for spacecraft in CR3BP tracking halo orbit.

Figure 6: State error trajectories for spacecraft in CR3BP tracking halo orbit.

Figure 7: Control history and OCP residual

||F||

for spacecraft in CR3BP tracking halo orbit.

To test the robustness of the NMPC control strategy we performed closed-loop simulations using the CR3BP for a variety of initial conditions. To be specific, we sampled 10 initial conditions from the uniform distribution over the hypercube

\mathcal{H}=\left\{x~|~-\begin{bmatrix}1_{3}\Delta r_{max}\\ 1_{3}\Delta v_{max}\par\end{bmatrix}\leq x-x_{0}\leq\begin{bmatrix}1_{3}\Delta r_{max}\\ 1_{3}\Delta v_{max}\par\end{bmatrix}\right\}

(14)

where $x_{0}\approx(0.9878,0,0.0290,0,0.8763,0)$ is the initial condition for the Halo orbit, $\Delta r_{max}=500~km$ , and $\Delta v_{max}=0.01~km/s$ . The results are shown in Figures 8 – 10, the spacecraft is able to reject the impact of the initial disturbance and successfully converge to the reference trajectory. As expected, this requires significant control effort.

Figure 8: Spacecraft state error trajectories in the CR3BP for a variety of initial conditions.

Figure 9: Spacecraft control inputs in the CR3BP for a variety of initial conditions.

Finally, we demonstrate the robustness of the controller to model mismatch. Figures 11 – 13 illustrate the closed-loop response of the controller when the simulation model eccentricity is set to $e=0.055$ . Tracking performance is significantly degraded and control effort is higher, however the controller is still able to stabilize the orbit. A comparison of Figures 6 and 12 shows the difference in tracking error between the two simulations, and a comparison of Figures 7 and 13 illustrates the difference in control effort and constraint handling. In the elliptic case, the spacecraft uses approximately an order of magnitude higher control effort ( $2N$ vs roughly $200mN$ in the $e=0$ case) to stabilize the system with a nonzero eccentricity. This robustness to large model mismatch and the ability to seamlessly accommodate changes in control constraints is an advantage of the MPC methodology. In the future, tracking error could be reduced by incorporating non-zero eccentricity values into the prediction model (8) and control utilization can be reduced through a better reference trajectory.

Figure 12: State error trajectories for spacecraft in ER3BP tracking halo orbit.

Figure 13: Control history and OCP residual

||F||

for spacecraft in ER3BP tracking halo orbit.

The controller has been implemented in MATLAB/Simulink using the FBstab [21] quadratic programming solver and translated into C code by Simulink 2019a. The mean execution time of the controller found to be $8.6~ms$ on a 2019 Macbook Pro with 2.4 GHz i9 processor. A preliminary clock scaling analysis suggests a controller execution time of approximately $100~ms$ on a 200 MHz RAD750 radiation hardened computer. Given a true anomaly sampling period of $\Delta\theta=0.01$ , which implies a sampling period of approximately $1$ hour, this indicates that the computational burden of the NMPC strategy is likely to be manageable.

In order to compare the spacecraft’s performance for different values of $\ell$ , we introduce a cost function $J$ , analogous to (10a),

J=\sum_{j}||\xi_{j}-\bar{\xi}_{j}||_{Q}^{2}+||u_{j}||_{R}^{2}

(15)

for all discrete $\theta$ instants from $\theta=0$ to the end of the simulation (here, 5 full revolutions of the halo orbit). The results of this sensitivity study are shown in Figure 14, showing minimal optimality benefits to running the SQP scheme for more than three iterations. This justifies the use of $\ell=3$ and illustrates the computational efficiency benefits of the TDO scheme.

Figure 14: Sensitivity study between max SQP iterations

\ell

and relative optimality of spacecraft trajectory tracking

J

In this work we have presented a nonlinear MPC-based approach to the station-keeping of halo orbits in the Earth-Moon system. A method is described for generating periodic halo orbits in the CR3BP and the same equations of motion, in combination with an RK4 integration scheme, are used in the prediction model for the controller. This nonlinear optimization problem is solved with time-distributed SQP techniques utilizing the FBstab quadratic programming algorithm. Finally, the controller is validated in numerical simulations of the CR3BP and ER3BP with process noise across several full revolutions of the halo orbit.

The results showed that even for large disturbances and model mismatch, the NMPC scheme is capable of stabilizing the system about the reference trajectory while satisfying control constraints. Additionally, it is shown that a TDO approach to the NMPC problem allows for the use of fewer computational resources without significant reduction in performance. This work illustrates the computation feasibility of NMPC for this application which opens the door to more advanced formulations that directly incorporate high level objectives such as rendezvous or minimal propellant consumption.

$\displaystyle\min_{\tilde{\xi},\tilde{u}}~~~$

$\displaystyle||\xi_{i|k}-\bar{\xi}_{i|k}||_{Q}^{2}+\frac{1}{2}\sum_{i=0}^{N-1}||\xi_{i|k}-\bar{\xi}_{i|k}||_{Q}^{2}+||u_{i|k}||_{R}^{2}$

$\displaystyle s.t.~~$

$\displaystyle\xi_{0|k}=\xi_{k},$

$\displaystyle\xi_{i+1|k}=f(\theta_{i|k},\xi_{i|k},u_{i|k}),~~i=0,\dots,N{-}1,$

$\displaystyle\|u_{i|k}\|_{\infty}\leq u_{max},~~i=0,\dots,N,$

[1] B. Hufenbach, K. Laurini, N. Satoh, C. Lange, R. Martinez, J. Hill, M. Landgraf, and A. Bergamasco, “International missions to lunar vicinity and surface-near-term mission scenario of the Global Space Exploration Roadmap,” IAF 66th International Astronautical Congress, 2015.
[2] K. C. Laurini, B. Hufenbach, J. Hill, and A. Ouellet, “The global exploration roadmap and expanding human/robotic exploration mission collaboration opportunities,” IAF 66th International Astronautical Congress, 2015.
[3] M. Woodard, D. Folta, and D. Woodfork, “ARTEMIS: the first mission to the lunar libration orbits,” 21st International Symposium on Space Flight Dynamics, Toulouse, France, 2009.
[4] N. Aeronautics and S. A. (NASA), “What is Artemis?,” https://www.nasa.gov/what-is-artemis. Accessed: 2020-05-15.
[5] E. M. Zimovan, K. C. Howell, and D. C. Davis, “Near rectilinear halo orbits and their application in cis-lunar space,” 3rd IAA Conference on Dynamics and Control of Space Systems, Moscow, Russia, 2017, p. 20.
[6] R. Whitley and R. Martinez, “Options for staging orbits in cislunar space,” 2016 IEEE Aerospace Conference, IEEE, 2016, pp. 1–9.
[7] L. Bury and J. W. McMahon, “Landing Trajectories to Moons from the Unstable Invariant Manifolds of Periodic Libration Point Orbits,” AIAA Scitech 2020 Forum, 2020, p. 2181.
[8] D. Guzzetti, E. M. Zimovan, K. C. Howell, and D. C. Davis, “Stationkeeping analysis for spacecraft in lunar near rectilinear halo orbits,” 27th AAS/AIAA Space Flight Mechanics Meeting, 2017, pp. 1–20.
[9] K. Howell and J. Breakwell, “Almost rectilinear halo orbits,” Celestial mechanics, Vol. 32, No. 1, 1984, pp. 29–52.
[10] J. B. Rawlings and D. Q. Mayne, Model predictive control: Theory and design. Nob Hill Pub., 2009.
[11] L. Grüne and J. Pannek, “Nonlinear model predictive control,” Nonlinear Model Predictive Control, pp. 45–69, Springer, 2017.
[12] M. Ellis, H. Durand, and P. D. Christofides, “A tutorial review of economic model predictive control methods,” Journal of Process Control, Vol. 24, No. 8, 2014, pp. 1156–1178.
[13] D. Davis, S. Bhatt, K. Howell, J.-W. Jang, R. Whitley, F. Clark, D. Guzzetti, E. Zimovan, and G. Barton, “Orbit maintenance and navigation of human spacecraft at cislunar near rectilinear halo orbits,” 2017.
[14] V. Muralidharan, A. Weiss, and U. V. Kalabic, “Control Strategy for Long-Term Station-Keeping on Near-Rectilinear Halo Orbits,” AIAA Scitech 2020 Forum, 2020, p. 1459.
[15] U. Kalabic, A. Weiss, S. Di Cairano, and I. Kolmanovsky, “Station-keeping and momentum-management on halo orbits around L2: Linear-quadratic feedback and model predictive control approaches,” Proc. AAS Space Flight Mechanics Meeting, 2015, pp. 15–307.
[16] D. C. Folta, T. A. Pavlak, A. F. Haapala, K. C. Howell, and M. A. Woodard, “Earth–Moon libration point orbit stationkeeping: theory, modeling, and operations,” Acta Astronautica, Vol. 94, No. 1, 2014, pp. 421–433.
[17] Y. Lian, G. Gómez, J. J. Masdemont, and G. Tang, “Station-keeping of real Earth–Moon libration point orbits using discrete-time sliding mode control,” Communications in nonlinear science and numerical simulation, Vol. 19, No. 10, 2014, pp. 3792–3807.
[18] J. V. Breakwell, A. A. Kamel, and M. J. Ratner, “Station-keeping for a translunar communication station,” Celestial Mechanics, Vol. 10, No. 3, 1974, pp. 357–373.
[19] C. Simó, G. Gómez, J. Llibre, R. Martinez, and J. Rodriguez, “On the optimal station keeping control of halo orbits,” Acta Astronautica, Vol. 15, No. 6-7, 1987, pp. 391–397.
[20] G. Gómez, J. Llibre, R. Martínez, and C. Simó, “Station keeping of a quasiperiodic halo orbit using invariant manifolds,” Proceed. 2nd Internat. Symp. on spacecraft flight dynamics, Darmstadt, 1986, pp. 65–70.
[21] D. Liao-McPherson and I. Kolmanovsky, “FBstab: A proximally stabilized semismooth algorithm for convex quadratic programming,” Automatica, Vol. 113, 2020, p. 108801, https://doi.org/10.1016/j.automatica.2019.108801.
[22] D. Liao-McPherson, M. Nicotra, and I. Kolmanovsky, “Time-distributed optimization for real-time model predictive control: Stability, robustness, and constraint satisfaction,” Automatica, Vol. 117, 2020, p. 108973, https://doi.org/10.1016/j.automatica.2020.108973.
[23] R. Broucke, “Stability of periodic orbits in the elliptic, restricted three-body problem.,” AIAA journal, Vol. 7, No. 6, 1969, pp. 1003–1009.
[24] K. Walker, B. Samadi, M. Huang, J. Gerhard, K. Butts, and I. Kolmanovsky, “Design environment for nonlinear model predictive control,” tech. rep., SAE Technical Paper, 2016.
[25] D. A. Herman, T. A. Tofil, W. Santiago, H. Kamhawi, J. E. Polk, J. S. Snyder, R. R. Hofer, F. Q. Picha, J. Jackson, and M. Allen, “Overview of the Development and Mission Application of the Advanced Electric Propulsion System (AEPS),” 2018.