End-to-End Lyapunov-Based Eclipse-Feasible Low-Thrust Transfer Trajectories to NRHO

Design of low-thrust transfers to the vicinity of the Moon is of interest, with the selection of a Near-Rectilinear Halo Orbit (HALO) by NASA as a staging platform for exploration of the Moon and beyond [1]. Designing low-thrust trajectories can be quite difficult due to the combination of the very low control authority available from low-thrust propulsion systems and the highly nonlinear dynamical environment of the cislunar region. Since existing low-thrust spacecraft are powered with solar arrays, it’s essential that no solar eclipses last longer than a certain time interval. For instance, for the transfer of the Co-Manifested Vehicle (CMV) of Gateway, all eclipse durations need to be less than 90 minutes [2]. Extended operation of spacecraft within eclipses can deplete the batteries and lead to a complete loss of the vehicle.

Due to low propulsive accelerations on the order of $1.0\times 10^{-4}$ m/s², transfer trajectories can require times of flight on the order of a year or longer. Furthermore, low-thrust trajectories consist of different phases, where the primary gravitational influence shift from Earth to the Moon. Therefore, transfer trajectories are typically solved in multiple subphases. Ref. [2] gives an overview of NASA’s third and most recent Design Reference Mission (DRM) for the transfer of the CMV, which was designed in four subphases using Copernicus [3]. Indirect optimization methods are also used for designing low-thrust trajectories to quasi-periodic, near-rectilinear Halo orbits that leverages ephemeris-driven, “invariant manifold analogs” as long-duration asymptotic terminal coast arcs [4]. All discontinuous events (such as entry into and exit from Earth eclipses and throttle switches) are made smooth through the powerful and novel Composite Smooth Control (CSC) framework [5]. Ref. [6] solves a similar transfer problem in two subphases with an indirect method to optimize a powered Earth-spiral subphase that is then heuristically patched into a second ballistic subphase. Numerical continuation and homotopy methods are fundamental to the convergence of the Hamiltonian two-point boundary-value problems associated with indirect methods [7, 8, 9, 10, 11]. To consolidate the design approach, we solved a similar problem in one phase (i.e., an Earth-centered perturbed two-body model is used with perturbations due to the Moon, Sun, and Earth’s second zonal harmonic subject to Earth eclipses) with a multiple-shooting indirect method [12]. Both minimum-time and minimum-fuel solutions were achieved by starting with a high level of spacecraft acceleration and performing numerical continuation to gradually reduce the value of acceleration.

Another approach for designing low-thrust many-revolution trajectories is to use Lyapunov-based approaches. Lyapunov control (LC) is based on the Lyapunov stability theory, using which an LC law is obtained by finding the expression for control which makes the time derivative of a control-Lyapunov function (CLF) for the system negative [13, 14, 15]. The CLF is positive in terms of the states of the system and should become 0 at the equilibrium (corresponding to a desired “goal” or target state). The method can be likened to converting the second-order trajectory optimization problem into a first-order stabilization problem. We note that LC has been used extensively for low-thrust trajectory optimization [\cite[cite]{[\@@bibref{}{petropoulos_simple_2003}{}{}]}]. For instance, Ref. [16] applies Q-law [17] to solve Earth-Moon transfers in two subphases. The results are shown to serve as high-quality initial guesses for GPOPS-II [18] – a pseudospectral direct method solver. The authors leverage the computational efficiency of LC to perform an extensive trade study over potential epochs and departure orbits, allowing for an a posteriori analysis of the eclipses. Ref. [15] proposes a hybrid LC based methodology for designing Earth-Moon transfers in a full-ephemeris model. A study on the sensitivity of the LC law to missed-thrust events is also performed to demonstrate the robustness of the control law.

In this paper, we consider a single-phase design approach similar to what we considered in Ref. [12]; however, a Lyapunov control (LC) law is used instead of an indirect method to solve the trajectory optimization problem. Further, convergence of LC laws is asymptotic and depends on the rate at which the Lyapunov function value is decreased. Finite-time convergence to the goal can be achieved by parameterizing the CLF and optimizing these new parameters with respect to a cost function to obtain near-optimal solutions (for example, with respect to flight time or fuel consumption). In addition, convergence tolerances can be set that define when the propagated state is “close enough” to the goal/target state. LC laws are also straightforward to design and implement, and more importantly, are closed-loop in nature (i.e., they only depend on the current state) [19]. Thus, a motivation for this work is to rapidly solve low-thrust transfer problems. These solutions, in turn, can be used as initial guesses for other high-fidelity trajectory optimization tools that will provide more optimal solutions that precisely satisfy boundary conditions like, for instance, the indirect method in Ref. [12] or Copernicus.

An important operational constraint, for low-thrust trajectories to the Gateway, is that the duration of all eclipses has to be less than a prescribed 90-minute threshold. Mission design strategies for ensuring all solar eclipse durations are less than the prescribed time often entail generating a large set of reference trajectories for a range of departure epochs, as was done for Artemis I [20], to have a variety of options. However, many departure windows may not be feasible. Ref. [20] reports that 18% of all launch days were infeasible for Artemis I. The eclipse-duration constraint can pose more challenges for extremely low-thrust propulsion systems that require significantly longer times of flight. For low-thrust transfers departing from a Geostationary Transfer Orbit (GTO) to the Moon, Earth eclipse events highly depend on the departure epoch and GTO orientation. Adjustments can be made to these values to mitigate eclipse durations. However, even with an analyst’s extensive experience, this post-processing approach can take many iterations to identify feasible launch opportunities. Long intermittent Earth and Moon eclipses occurring in cislunar space, when the spacecraft’s relative velocity is much slower, are possible and not as easily preventable.

Incorporating eclipse-duration constraints within the trajectory optimization can could increase the number of feasible departure windows and improve optimality by allowing the trajectory optimization to be “aware” of such constraints within the optimization process. Eclipse-duration constraints are difficult to enforce due to the fact that 1) the number of eclipses is not known a priori and 2) the number and duration of eclipses can change within the trajectory optimization process. Ref. [20] outlines an effective strategy that treats eclipse durations as inequality constraints in Copernicus. The results, in the paper, indicate that it was possible to increase the number of feasible launch dates by about 20% for Artemis I. In this paper, however, we attempt to satisfy the eclipse maximum-duration constraint with a soft-penalization enforced while optimizing the parameters of the transfer problem. This eliminates an analyst-in-the-loop design approach and automates and facilitates the solution procedure.

The contributions of the paper are as follows. A LC law is derived and used to solve transfer trajectory optimization problems similar to those in Refs. [2, 6]. The LC law can only be used to transfer the spacecraft starting from a fully defined state into an orbit, i.e., in its standard formulation LC cannot be used for rendezvous type transfers unless modifications are applied to the problem formulation [21]. Therefore, always starting from a point on the NRHO, the control law is applied backward in time for a departure from a GTO to an insertion at NRHO and forward in time for a departure from NRHO to an insertion at GTO. We only consider transfer maneuvers from GTO to NRHO. The GTO departure time and true anomaly are free and the NRHO insertion time is fixed. However, we consider the NRHO insertion time to be a design parameter. Because the ephemeris-propagated and ephemeris-corrected NRHO provided in Ref. [1] is considered, the entire state on the NRHO can be defined by the time (i.e., epoch). While numerically integrating the spacecraft equations of motion, event detection is used to determine if the solutions converge, if the spacecraft intersects the surface of the Earth or Moon, and when eclipses due to the Earth and Moon occur. Coasting is enforced during eclipses and the duration of each eclipse is calculated. Particle swarm optimization [22] is used to optimize the NRHO insertion date, time of flight, and parameters of the CLF with respect to a hierarchical cost function. The cost function prioritizes 1) convergence to the target orbit, 2) satisfaction of the maximum-eclipse-duration constraint, and 3) minimization of the time of flight.

The entire transfer problem is solved in the J2000 Earth-centered inertial (ECI) frame. All accelerations are expressed with respect to this frame. The spacecraft’s motion is modeled with position, $\bm{r}=[x,\,y,\,z]^{\top}$, and velocity, $\bm{v}=[v_{x},\,v_{y},\,v_{z}]^{\top}$, vectors. The spacecraft’s state vector is $\bm{x}=[\bm{r}^{\top},\,\bm{v}^{\top}]^{\top}$ and the equations of motion are defined as,

where $t$ denotes time, $\bm{a}_{\text{kep}}$ is the two-body (Keplerian) acceleration due to the Earth, $\bm{a}_{\text{3rd}}$ is the collection of third-body gravitational perturbations, and $\bm{a}_{J_{2}}$ is the acceleration due to Earth’s $J_{2}$ gravitational perturbation. In the last acceleration term, $\hat{\bm{\alpha}}a_{\text{sc}}\delta_{\text{ecl}}$, which denotes the acceleration produced by the propulsion system, $\hat{\bm{\alpha}}$ is the thrust steering unit vector and $\delta_{\text{ecl}}\in\{0,1\}$ is the eclipse-triggered throttle factor. Since the contribution and emphasis of the work is on satisfying the maximum-eclipse-duration constraint, a constant spacecraft acceleration is assumed with its value set to $a_{\text{sc}}=1.0\times 10^{-4}$ m/s². This value is chosen to match the transfer problems in Refs. [2, 6]. Propellant-mass considerations belong to our future work. The thrust steering unit vector can freely orient in space, but is constrained to a unit vector, i.e.,

In this work, the change in spacecraft mass is not taken into account. But, our future work will investigate implementing a LC law coasting mechanism to obtain suboptimal minimum-fuel solutions, like the one that is introduced with Q-law in Ref. [23]. Earth’s two-body acceleration can be written as,

where $r=\|\bm{r}\|$ and $\mu_{\text{Earth}}$ is the gravitational parameter of the Earth. Perturbing accelerations due to the gravity of the Moon, Sun, and Jupiter are considered and written as [24, 25],

where

with $K\in\left\{\text{Moon},\,\text{Sun},\,\text{Jupiter}\right\}$, and $\bm{r}_{k}$ denotes the position of the $k$-th body with respect to the Earth. Note that this formulation avoids any numerical error due to cancellations when terms are of significantly different values [24]. The acceleration vector due to Earth’s $J_{2}$ gravitational perturbation is written as [14, 26],

where $R_{\text{Earth}}$ is the mean radius of the Earth and $J_{2}=1082.63\times 10^{-6}$.

A canonical distance unit (DU) is defined by one Earth radius, $R_{\text{Earth}}$, and a canonical time unit (TU) is defined such that the scaled value of Earth’s gravitational parameter is 1 DU³/TU². These canonical distance and time units are then used to scale all states and parameters of the dynamical model. Future work could include investigating more sophisticated scaling methods and even time regularization methods such as Ref. [27] that might make numerical integration of Eq. (1) more efficient. All planetary ephemerides and constants are obtained using NASA’s SPICE toolkit [28] and the generic kernel files¹¹1https://naif.jpl.nasa.gov/pub/naif/generic_kernels/

It is computationally inefficient to call SPICE routines during numerical integration. Thus, all ephemerides obtained positions, i.e., $\bm{r}_{\text{k}}$ appearing in Eq. (4) and in the eclipse model presented in the next section, are fitted by a spline function, which has proved to be more computationally efficient. The interpolation is performed to an accuracy on the order of 0.1 m.

In this work, eclipses due to the Earth and Moon are considered. The cylindrical eclipse model from Ref. [29] is used. The eclipse model assumes the Earth, Moon, and Sun to be perfect spheres and the spacecraft to be a point mass. Eclipse coasting is enforced during both umbra (total eclipse) and penumbra (partial eclipse). Thus, only the penumbra exits and entries are calculated, since umbra occurs inside penumbra. Figure 1 illustrates the Sun-Earth shadow geometry. Note that Figure 1 is greatly exaggerated and not drawn to scale. The same geometry is also used for modeling Moon eclipses.

Let $r_{\text{Sun}}=\|\bm{r}_{\text{Sun}}\|$ be the distance between the Earth and the Sun and let $\eta\in[0,1]$. From the geometrical proportion of the penumbral cone, we have

where the value $\eta$ can be expressed as,

Therefore, the angle that the penumbral cone makes with respect to $\bm{r}_{\text{Sun}}$ can be expressed as,

The position of the spacecraft projected onto $\bm{r}_{\text{Sun}}$ is defined as,

where $\hat{\bm{r}}_{\text{Sun}}=\bm{r}_{\text{Sun}}/r_{\text{Sun}}$. Earth shadows only occur on the side of the Earth opposite the Sun, i.e., the spacecraft can only encounter a shadow when $\bm{r}^{\top}\hat{\bm{r}}_{\text{Sun}}<0$. Let the distance between the spacecraft and the center of the penumbral cone be defined as,

and let the distance between the penumbral terminator and the center of the penumbral cone at the projected spacecraft location be defined as,

Therefore, it can be stated that the spacecraft is in the Earth shadow when $\gamma_{\text{Earth}}<\kappa_{\text{Earth}}$ and $\bm{r}^{\top}\hat{\bm{r}}_{\text{Sun}}<0$ and not in a shadow otherwise.

Let $\gamma_{\text{Moon}}$ and $\kappa_{\text{Moon}}$ denote the same definitions as Eqs. (11) and (12), respectively, but for the Moon-Sun eclipse model. Also, let $\bm{r}_{\text{sc/Moon}}=\bm{r}-\bm{r}_{\text{Moon}}$ be the position of the spacecraft with respect to the Moon and $\hat{\bm{r}}_{\text{Sun/Moon}}=\left(\bm{r}_{\text{Sun}}-\bm{r}_{\text{Moon}}\right)/\left\|\bm{r}_{\text{Sun}}-\bm{r}_{\text{Moon}}\right\|$ be the unit vector pointing towards the Sun with respect to the Moon. The following switching functions can then be defined

The eclipse-triggered throttle factor in Eq. (1), $\delta_{\text{ecl}}$, can be defined as a multiplication of two factors as,

where

and

During numerical integration of Eq. (1), an event-detection algorithm is used to determine the exact time of each eclipse entry, $t_{\text{ecl},n}^{-}$, and exit, $t_{\text{ecl},n}^{+}$, for all $n$ eclipses with $n=1,\dots,N$. The duration of each eclipse, denoted as $d_{\text{ecl},n}=t_{\text{ecl},n}^{+}-t_{\text{ecl},n}^{-}$, is also calculated. Note that the total number of Earth and Moon eclipses, $N$, is not known in advance. The eclipse constraint, that all eclipses be less than 90 minutes [2], is considered in this work and can be expressed formally as,

Unlike the eclipse model used by the authors in Ref. [12], the domain for the eclipse model from Ref. [29] is defined interior to the occulting bodies. This was the main reason that we adopted this eclipse model. While optimizing all the parameters of the transfer problem, event detection is used to stop integration when the spacecraft intersects the Earth or Moon. This logic works most of the time, however, sometimes the dynamics become significantly nonlinear around the Moon and the event-detection algorithm misses the intersection event. When the model from Ref. [30] is used, NaN’s are returned since the model is undefined for the domain inside the occulting body, which breaks the optimization routine. A similar problem is reported with the eclipse model used in [20] along with a strategy for overcoming it. The previously presented eclipse model circumvents this problem altogether.

The transfer problem, departing from a GTO at a date $t_{0}$ and inserting into NRHO at a later date $t_{f}$, is solved backwards in time over the time horizon

The GTO orbital parameters are taken from Ref. [2] with apogee and perigee altitudes of 33,900 km and 350 km, respectively, and an inclination of 28.5$\degree$. Because it is stated in Ref. [2] that the right ascension of the ascending node (RAAN) and argument of perigee (ARGP) of the GTO are unrestricted for the initial analysis, only the specific angular momentum, eccentricity, and inclination are considered in the boundary conditions for the GTO, leaving the RAAN and ARGP as free parameters. True anomaly is also free, however, this fact is inherent to the LC law that will be used to solve the trajectory optimization problem.

The boundary conditions on the GTO are therefore defined as

where subscript ‘GTO’ denotes values of the GTO and the specific angular momentum, $h$, eccentricity, $e$, and inclination, $i$, are defined as [31],

where $v=\|\bm{v}\|$ and $h_{\bm{z}}$ denotes the $z$ component of the specific angular momentum vector in the J2000 ECI frame.

It is expected for the spacecraft’s osculating orbit, with respect to the Earth, to become hyperbolic close to the Moon. Thus, the angular momentum was selected as opposed to, for example, the semimajor axis. The semiparameter, defined as $p=h^{2}/\mu_{\text{Earth}}$ [31], could also be an acceptable element to target. Ultimately, the goal is to target the size, shape, and inclination of the GTO only. A variety of other boundary conditions could be formulated of which some could lead to a better CLF and, subsequently, a better control law and therefore will be investigated in our future work.

The NRHO ephemeris is obtained from the kernel file²²2https://naif.jpl.nasa.gov/pub/naif/misc/MORE_PROJECTS/DSG/

described in Ref. [1]. The boundary condition at $t_{f}$ is defined as

The minimum-time constant-acceleration transfer trajectory optimization problem can be stated as,

A parameterized LC law based on the goal defined by Eq. (20) will be derived and substituted into Eq. (1). Eq. (23) can then be solved as a parameter optimization problem using a heuristic algorithm in which Eq. (1) is integrated over the time horizon given by Eq. (19) with the initial condition given by Eq. (22).

We note that Eq. (18) is quite challenging to enforce directly since the number of eclipses, $N$, is not known a priori and can also change during the iterations of the optimization process. Instead, Eq. (18) is enforced as a soft penalty along with a penalty to further promote satisfaction of Eq. (20) since it is not guaranteed. These two penalties and the time of flight are encoded into a single cost function that the heuristic algorithm minimizes. This cost function will be explained in detail after the control law is derived.

The control-Lyapunov function (CLF) is defined as,

where the constraint vector, $\bm{w}(\bm{x})$, is defined as,

Note that no scaling is performed on Eq. (25) because the canonical scaling method ensures that $h$ has the same order of magnitude of $e$ and $i$. Instead of using a diagonal parameter matrix $\bm{K}$, a full parameter matrix is used to consider a larger family of CLFs as it is proposed in [32]. The $3\times 3$ positive-definite matrix $\bm{K}$ is defined through a novel eigendecomposition method as,

where the column vectors of $\bm{Q}$ make up the eigenvectors of $\bm{K}$ and $\bm{\Lambda}$ is a diagonal matrix of the eigenvalues of $\bm{K}$. This parameterization is based on Ref. [32] and allows an efficient (i.e., minimal number of parameters) representation of a full matrix that is guaranteed to be positive-definite subject only to bounds on its parameters. The eigenvalue matrix, $\bm{\Lambda}$, is simple to construct, i.e.,

where the parameters $k_{1}$, $k_{2}$, and $k_{3}$ are constrained to being real and positive. The matrix $\bm{Q}$ can be generated as a rotation matrix and Ref. [32] outlines a generalized way to generate rotation matrices in $n$-dimensions. Because $\bm{Q}$ is $3\times 3$, in this work, the method in Ref. [32] reduces to any standard Euclidean rotation matrix parameterized by 3 angle-like parameters, $k_{4}$, $k_{5}$, and $k_{6}$.

A control law is derived by making the time-derivative of Eq. (24), $dV/dt=\dot{V}$, as negative as possible subject to Eq. (2), i.e., we have

In Eq. (28), $\dot{V}$ can be found through the chain rule as,

Since $V$ does not explicitly depend on time, Eq. (29) reduces to

It can be shown that Eq. (28) is pointwise satisfied with the selection of thrust steering vector as,

Because the transfer problem is being solved backwards in time, the sign of Eq. (31) should be reversed to ensure $V$ approaches 0 at the GTO departure time $t_{0}$. The resulting control law is

Eq. (32) is calculated in CasADi [33], a symbolic framework that uses automatic differentiation. Note that if the transfer problem that departs at NRHO and arrives at GTO were solved, then Eq. (31) would be used.

Due to the very low control authority available from the low-thrust propulsion system (compared to the highly-perturbed dynamical model), the CLF time derivative, Eq. (30), may become positive over one (or more) finite intervals, thereby not guaranteeing the system to converge as stated by the Lyapunov stability theory. However, this doesn’t guarantee nonconvergence either. The results of Ref. [15] show that converged solutions can still be found as long as the CLF time derivative is negative almost everywhere except on a finite number of small intervals. This property is also observed to be held in our numerical results.

After deriving the LC law, the trajectory optimization problem in Eq. (23) can now be stated as a parameter optimization problem (POP). The parameters considered are the NRHO insertion date, $t_{f}$, bounded by

the time of flight, $\Delta t$, bounded by

and finally the 6 CLF parameters in Eq. (26), bounded by

These parameters are optimized using MATLAB’s particle swarm optimization (PSO), a stochastic optimization algorithm that optimizes a scalar cost function subject only to bounds on the design variables. Under this parametrization, the time horizon in Eq. (19) can be expressed as,

An important step in the resulting POP is to accurately solve for the initial value problem (IVP) given by the set of ordinary differential equation (ODEs) in Eq. (1) with the control law in Eq. (32) and boundary condition given by Eq. (22) over the time horizon in Eq. (36), i.e.,

MATLAB’s variable-step variable-order nonstiff ODE integrator ode113 is used with an absolute and relative tolerance of $1.0\times 10^{-10}$. Our extensive numerical studies indicate that this integrator performed most efficiently with the prescribed tolerances against the rest of MATLAB’s ODE integrators.

The event-detection capability of ode113 is used extensively while solving Eq. (37). The method works by monitoring the sign of an $M$ number of scalar event functions, $e_{m}$ for $m=1,\dots,M$. When the $m$-th function changes sign, a regula falsi algorithm is used to find the precise location of the zero of $e_{m}$, and, if all the corresponding termination conditions are met, integration stops. In our paper, there are $M=7$ event functions defined as follows,

Eqs. (38a) and (38b) monitor if the spacecraft has intersected 200 km above the surface of the Earth and Moon. If either of their signs change, then integration stops and the cost function is appropriately updated and returned to PSO. Eqs. (38c), (38d), and (38e) monitor if the solution has converged or not (i.e., orbit insertion has been achieved or not). If any one of their signs become negative while the other two are also negative, then integration stops and the cost function is appropriately updated and returned to PSO. The tolerance $\epsilon=1.0\times 10^{-3}$ was chosen as it provides a balance between convergence speed and accuracy; however, future work should investigate using different convergence criteria.

Eqs. (38f) and (38g) are from Eqs. (13) and (14), respectively, and determine if the spacecraft is in an eclipse or not. If Eq. (13) (resp. Eq. (14)) changes sign while $S_{\text{ecl,Earth,2}}$ (resp. $S_{\text{ecl,Moon,2}}$) is negative, then integration is terminated. However, integration is then restarted from the same time and state. This logic is followed so that the discrete function in Eq. (15) is modeled as accurately as possible.

In formulating the cost function, let $t_{\text{end}}$ denote the final time returned by Eq. (37) under the event-detection logic, i.e., $t_{\text{end}}$ always satisfies $t_{f}-\Delta t\leq t_{\text{end}}\leq t_{f}$. The first priority of the cost function is to ensure the solution converges. If a solution to Eq. (37) does not satisfy the constraint below,

then, the value of cost, $J_{1}$, defined as,

is returned. Because of the highly nonlinear dynamics in the vicinity of the Moon, it was found that solutions commonly intersect the surface of the Moon. Thus, if integration was stopped due to Eq. (38b) becoming negative, then,

is returned as the cost function to PSO as a penalization.

If Eq. (39) is satisfied for a solution, but Eq. (18) is not, then

is returned as the cost function to PSO. Note that it appears there’s a possibility for a division by zero in Eq. (42), however, the expression in Eq. (42) is not evaluated if Eq. (18) is satisfied, which eliminates this possibility. Also, Eq. (42) is made negative to differentiate it from Eqs. (40) and (41), but inverted so that the violated eclipse durations are still minimized.

Finally, if Eqs. (39) and (18) are satisfied for a solution, then

is returned as the cost function to PSO. This cost is also made negative to differentiate from Eqs. (40) and (41) and inverted so that time of flight is minimized. However, to ensure it is differentiated from Eq. (42), it has units of $1/[\text{years}]$ while Eq. (42) has units of $1/[\text{seconds}]$ so that they are on different orders of magnitude. Further, to ensure that minuscule eclipse violations aren’t interpreted as extremely low times of flights, if $J_{3}<J_{4}$ occurs for a converged solution, then Eq. (43) is returned instead of Eq. (42). While this allows for solutions with eclipses longer than 90 minutes to be deemed feasible, these solutions will only be infeasible by a duration on the order of a second. Note that the cost function used in this work is not continuous due to the logic involved, however, stochastic optimization algorithms, such as PSO, can deal with discontinuous cost functions.

Because LC laws are prone to extreme oscillations/chattering at the end of the maneuver, the step size of variable step integrators can become minuscule and halt progress [13]. To overcome this issue, integration of ode113 is stopped when a certain number of function evaluations is reached. A simple logic is implemented inside ode113 and, if triggered, then Eq. (40) is simply returned as the cost function to PSO. In this paper, $1.0\times 10^{5}$ function evaluations were arbitrarily chosen and found to provide acceptable results; however, different values may further benefit the algorithm given that the number of iterations is a problem-dependent number. The event-detection logic and cost function values are summarized in Algorithm 1.

The results are presented for a transfer problem with a fixed NRHO insertion date and then with a variable NRHO insertion date. These two cases are considered to demonstrate the impact of the insertion date and the types of solutions that can be achieved. The fixed NRHO insertion date was arbitrarily selected as $t_{f}=$ 2026 DEC 06 00:00:00 UTC. This coincides with a state that is roughly at apolune on the NRHO. When $t_{f}$ is allowed to vary, it is bounded between $t_{f,\text{lb}}=$ 2026 NOV 06 00:00:00 UTC and $t_{f,\text{ub}}=$ 2027 JAN 05 00:00:00 UTC, or, $t_{f,\text{lb}}=t_{f}-30$ days and $t_{f,\text{ub}}=t_{f}+30$ days. The value of 30 days was selected because it is approximately equal to 1 period of the Moon’s orbit around the Earth, giving a variety of phasing possibilities to the solution space. The time of flight for all cases was bounded by $\Delta t_{\text{lb}}=200$ days and $\Delta t_{\text{ub}}=400$ days. Simulations were performed on a 2023 MacBook Pro with the Apple M2 Pro chip, which allows for 10 “workers” in MATLAB to run PSO in parallel.

We presented a methodology for efficiently finding low-thrust spacecraft transfer trajectories under constant acceleration from a geostationary transfer orbit (GTO) to the near-rectilinear halo orbit (NRHO) earmarked for NASA’s Gateway. The method is based on a closed-loop control law derived from a novel parameterization of quadratic control-Lyapunov functions. This control law is applied in a backward-in-time sense to generate solutions departing from the GTO and inserting into the NRHO. Solutions may also be obtained that depart from the NRHO and insert into the GTO.

To solve the resulting trajectory optimization problems, the parameters of the control law, time of flight, and NRHO insertion date are all optimized simultaneously with a stochastic optimization algorithm – particle swarm optimization (PSO). The cost function is designed to prioritize 1) convergence to the target orbit, 2) satisfaction of the constraint that all eclipse durations be less than 90 minutes, and 3) minimization of the time of flight. Results indicate that eclipse-feasible solutions can be obtained on the order of 10 minutes with the processing power of a personal laptop computer. Solutions obtained can serve as high-quality initial guesses to NASA’s high-fidelity trajectory optimization tools such as Copernicus and GMAT.

We would like to acknowledge Saeid Tafazzol for his useful suggestions and discussions on efficiently parameterizing $n$-dimensional positive-definite matrices.