End-to-End Lyapunov-Based Eclipse-Feasible Low-Thrust Transfer Trajectories to NRHO
Nicholas P. Nurre Ehsan Taheri








End-to-End Lyapunov-Based Eclipse-Feasible Low-Thrust Transfer Trajectories to NRHO
    Nicholas P. NurreGraduate Student, Department of Aerospace Engineering, Auburn University, 141 Engineering Dr Auburn, AL 36849.  and  
Ehsan TaheriAssistant Professor, Department of Aerospace Engineering, Auburn University, 141 Engineering Dr Auburn, AL 36849.
    September 5, 2024
========================================================================================================================================================================================================================================================================

 		

§ ABSTRACT 
Generating low-thrust transfer trajectories between Earth and the Near Rectilinear Halo Orbit (NRHO), that is selected for NASA's Gateway, can be challenging due to the low control authority available from the propulsion system and the important operational constraint that the duration of all eclipses has to be less than a prescribed 90-minute threshold. We present a method for generating eclipse-feasible, minimum-time solutions to the aforementioned trajectory design problem using a Lyapunov control law. Coasting is enforced during solar eclipses due to both the Earth and Moon. We used particle swarm optimization to optimize the NRHO insertion date, time of flight, and control law parameters according to a cost function that prioritizes 1) convergence to the target orbit, 2) satisfaction of eclipse-duration constraints, and 3) minimization of time of flight. Trajectories can serve as initial guesses for NASA's high-fidelity trajectory design tools such as Copernicus and GMAT.




§ INTRODUCTION

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 <cit.>. 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 <cit.>. 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×10^-4 m/s^2, 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. <cit.> 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 <cit.>. 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 <cit.>. 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 <cit.>. Ref. <cit.> 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 <cit.>. 

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  <cit.>. 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 <cit.>. 
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 <cit.>. For instance, Ref. <cit.> applies Q-law <cit.> to solve Earth-Moon transfers in two subphases. The results are shown to serve as high-quality initial guesses for GPOPS-II <cit.> – 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. <cit.> 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. <cit.>; 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) <cit.>. 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. <cit.> 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 <cit.>, to have a variety of options. However, many departure windows may not be feasible. Ref. <cit.> 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. <cit.> 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. <cit.>. 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 <cit.>. 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. <cit.> 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 <cit.> 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.





§ DYNAMICAL MODEL

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, r=[x, y, z]^⊤, and velocity, v=[v_x, v_y, v_z]^⊤, vectors. The spacecraft's state vector is x=[r^⊤, v^⊤]^⊤ and the equations of motion are defined as,

    ẋ = f(t,x,α̂) = [                                    v; a_kep + a_3rd + a_J_2 + α̂ a_scδ_ecl ],

where t denotes time, a_kep is the two-body (Keplerian) acceleration due to the Earth, a_3rd is the collection of third-body gravitational perturbations, and a_J_2 is the acceleration due to Earth's J_2 gravitational perturbation. In the last acceleration term, α̂ a_scδ_ecl, which denotes the acceleration produced by the propulsion system, α̂ is the thrust steering unit vector and δ_ecl∈{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_sc=1.0×10^-4 m/s^2. This value is chosen to match the transfer problems in Refs. <cit.>. 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.,

    α̂^⊤α̂ = 1.


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. <cit.>. Earth's two-body acceleration can be written as,

    a_kep = -μ_Earth/r^3r,

where r=r and μ_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 <cit.>,

    a_3rd = -∑_k∈ Kμ_kr + F(q_k)r_k/r - r_k^3,

where

    F(q_k)     = q_k(3 + 3q_k + q_k^2/1+(1+q_k)^3/2),     q_k     = r^⊤(r-2r_k)/r_k^⊤r_k,

with K∈{Moon, Sun, Jupiter}, and 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 <cit.>. The acceleration vector due to Earth's J_2 gravitational perturbation is written as <cit.>,

    a_J_2 = 3J_2μ_Earth R_Earth^2/2r^4[ x/r(5x^2/r^2-1), y/r(5y^2/r^2-1),  z/r(5z^2/r^2-3) ]^⊤,

where R_Earth is the mean radius of the Earth and J_2=1082.63×10^-6. 

A canonical distance unit (DU) is defined by one Earth radius, R_Earth, and a canonical time unit (TU) is defined such that the scaled value of Earth's gravitational parameter is 1 DU^3/TU^2. 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. <cit.> that might make numerical integration of Eq. (<ref>) more efficient. All planetary ephemerides and constants are obtained using NASA's SPICE toolkit <cit.> and the generic kernel files[<https://naif.jpl.nasa.gov/pub/naif/generic_kernels/>] 
    , , , and .
 

It is computationally inefficient to call SPICE routines during numerical integration. Thus, all ephemerides obtained positions, i.e., r_k appearing in Eq. (<ref>) 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.



§ ECLIPSE MODEL

In this work, eclipses due to the Earth and Moon are considered. The cylindrical eclipse model from Ref. <cit.> 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 <ref> illustrates the Sun-Earth shadow geometry. Note that Figure <ref> is greatly exaggerated and not drawn to scale. The same geometry is also used for modeling Moon eclipses.



Let r_Sun=r_Sun be the distance between the Earth and the Sun and let η∈ [0,1]. From the geometrical proportion of the penumbral cone, we have

    (1-η)r_Sun/2R_Earth=η r_Sun/2 R_Sun,

where the value η can be expressed as,

    η = R_Sun/R_Earth+R_Sun.












Therefore, the angle that the penumbral cone makes with respect to r_Sun can be expressed as,

    θ_p = sin^-1(R_Earth/(1-η)r_Sun) = sin^-1(R_Earth + R_Sun/r_Sun).










The position of the spacecraft projected onto r_Sun is defined as,

    r̅ = (r^⊤r̂_Sun)r̂_Sun,

where r̂_Sun=r_Sun/r_Sun. Earth shadows only occur on the side of the Earth opposite the Sun, i.e., the spacecraft can only encounter a shadow when r^⊤r̂_Sun<0. Let the distance between the spacecraft and the center of the penumbral cone be defined as,

    γ_Earth = r - r̅,

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

    κ_Earth = ((1-η)r_Sun + r̅)tan(θ_p).


Therefore, it can be stated that the spacecraft is in the Earth shadow when γ_Earth < κ_Earth and r^⊤r̂_Sun<0 and not in a shadow otherwise. 

Let γ_Moon and κ_Moon denote the same definitions as Eqs. (<ref>) and (<ref>), respectively, but for the Moon-Sun eclipse model. Also, let r_sc/Moon=r - r_Moon be the position of the spacecraft with respect to the Moon and r̂_Sun/Moon = (r_Sun - r_Moon)/r_Sun - r_Moon be the unit vector pointing towards the Sun with respect to the Moon. The following switching functions can then be defined

    S_ecl,Earth,1    = γ_Earth - κ_Earth ,     S_ecl,Earth,2    = r^⊤r̂_Sun, 
    
        S_ecl,Moon,1    = γ_Moon - κ_Moon ,     S_ecl,Moon,2    = r_sc/Moon^⊤r̂_Sun/Moon.


The eclipse-triggered throttle factor in Eq. (<ref>), δ_ecl, can be defined as a multiplication of two factors as,

    δ_ecl = δ_ecl,Earth×δ_ecl,Moon,

where

    δ_ecl,Earth =  0,        S_ecl,Earth,1<0 and S_ecl,Earth,2<0, 
     1,    else,

and 

    δ_ecl,Moon =  0,        S_ecl,Moon,1<0 and S_ecl,Moon,2<0, 
     1,    else.


During numerical integration of Eq. (<ref>), an event-detection algorithm is used to determine the exact time of each eclipse entry, t_ecl,n^-, and exit, t_ecl,n^+, for all n eclipses with n=1,…,N. The duration of each eclipse, denoted as d_ecl,n = t_ecl,n^+ - t_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 <cit.>, is considered in this work and can be expressed formally as,

    d_ecl,n≤ 90 minutes, ∀ n = 1,…,N.


Unlike the eclipse model used by the authors in Ref. <cit.>, the domain for the eclipse model from Ref. <cit.> 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. <cit.> is used, '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 <cit.> along with a strategy for overcoming it. The previously presented eclipse model circumvents this problem altogether.



§ TRAJECTORY OPTIMIZATION PROBLEM

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 

    t∈[t_f,t_0],        t_f > t_0.


The GTO orbital parameters are taken from Ref. <cit.> with apogee and perigee altitudes of 33,900 km and 350 km, respectively, and an inclination of 28.5. Because it is stated in Ref. <cit.> 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

    h(x(t_0))     = h_GTO,     e(x(t_0))     = e_GTO,     i(x(t_0))     = i_GTO,

where subscript `GTO' denotes values of the GTO and the specific angular momentum, h, eccentricity, e, and inclination, i, are defined as <cit.>,

    h     = r×v,     e     = (v^2-μ_Earth/r)r - (r^⊤v)v/μ_Earth,     i     = cos^-1(h_ẑ/h),

where v=v and h_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/μ_Earth <cit.>, 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[<https://naif.jpl.nasa.gov/pub/naif/misc/MORE_PROJECTS/DSG/>] 
    
 described in Ref. <cit.>. The boundary condition at t_f is defined as

    x(t_f) = x_NRHO(t_f).


The minimum-time constant-acceleration transfer trajectory optimization problem can be stated as,
 
    
    min_α̂,t_0,t_f      J = t_f-t_0, 
    s.t.,     Eqs. (<ref>),(<ref>),(<ref>),(<ref>),(<ref>),(<ref>).



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

We note that Eq. (<ref>) 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. (<ref>) is enforced as a soft penalty along with a penalty to further promote satisfaction of Eq. (<ref>) 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.



§ LYAPUNOV CONTROL LAW

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

    V(x) = 1/2w(x)^⊤Kw(x),

where the constraint vector, w(x), is defined as,

    w(x) = [ h(x) - h_GTO, e(x) - e_GTO,  i(x) - i_GTO ]^⊤.


Note that no scaling is performed on Eq. (<ref>) 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 K, a full parameter matrix is used to consider a larger family of CLFs as it is proposed in <cit.>. The 3×3 positive-definite matrix K is defined through a novel eigendecomposition method as,

    K= QΛQ^⊤,

where the column vectors of Q make up the eigenvectors of K and Λ is a diagonal matrix of the eigenvalues of K.  This parameterization is based on Ref. <cit.> 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, Λ, is simple to construct, i.e.,

    Λ = [ k_1   0   0;   0 k_2   0;   0   0 k_3 ],

where the parameters k_1, k_2, and k_3 are constrained to being real and positive. The matrix Q can be generated as a rotation matrix and Ref. <cit.> outlines a generalized way to generate rotation matrices in n-dimensions. Because Q is 3×3, in this work, the method in Ref. <cit.> 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. (<ref>), dV/dt=V̇, as negative as possible subject to Eq. (<ref>), i.e., we have

    α̂^* = _α̂=1V̇.




In Eq. (<ref>), V̇ can be found through the chain rule as,

    V̇ = ∂ V/∂x∂x/∂ t + ∂ V/∂ t.


Since V does not explicitly depend on time, Eq. (<ref>) reduces to 

    V̇ = ∂ V/∂rv + ∂ V/∂v(a_kep + a_3rd + a_J_2 + α̂a_scδ_ecl).


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

    α̂^* = -(∂ V/∂v/∂ V/∂v)^⊤.


Because the transfer problem is being solved backwards in time, the sign of Eq. (<ref>) should be reversed to ensure V approaches 0 at the GTO departure time t_0. The resulting control law is

    α̂^* = (∂ V/∂v/∂ V/∂v)^⊤.


Eq. (<ref>) is calculated in CasADi <cit.>, 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. (<ref>) 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. (<ref>), 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. <cit.> 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.



§ PARAMETER OPTIMIZATION PROBLEM

After deriving the LC law, the trajectory optimization problem in Eq. (<ref>) can now be stated as a parameter optimization problem (POP). The parameters considered are the NRHO insertion date, t_f, bounded by

    t_f,lb≤ t_f ≤ t_f,ub,

the time of flight, Δ t, bounded by 

    Δ t_lb≤Δ t ≤Δ t_ub,

and finally the 6 CLF parameters in Eq. (<ref>), bounded by

    0 < k_i ≤ k_ub, for i=1,2,3,        0 ≤ k_4 ≤π,        0 ≤ k_j ≤ 2π, for j=5,6.


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. (<ref>) can be expressed as,

    t ∈ [t_f,t_f - Δ t].


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. (<ref>) with the control law in Eq. (<ref>) and boundary condition given by Eq. (<ref>) over the time horizon in Eq. (<ref>), i.e.,

    ẋ    = f(t,x,α̂^*;K),       x(t_f) = x_NRHO(t_f),     t    ∈[t_f,t_f-Δ t].


MATLAB's variable-step variable-order nonstiff ODE integrator  is used with an absolute and relative tolerance of 1.0×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  is used extensively while solving Eq. (<ref>). The method works by monitoring the sign of an M number of scalar event functions, e_m for m=1,…,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,
 
    
    e_1     = r - R_Earth - 200 km, 
    
            e_2     = r - R_Moon - 200 km, 
    
            e_3     = |h - h_GTO| - ϵ, 
    
            e_4     = |e - e_GTO| - ϵ, 
    
            e_5     = |i - i_GTO| - ϵ, 
    
            e_6     = S_ecl,Earth,1, 
    
            e_7     = S_ecl,Moon,1.



Eqs. (<ref>) and (<ref>) 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. (<ref>), (<ref>), and (<ref>) 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 ϵ=1.0×10^-3 was  chosen as it provides a balance between convergence speed and accuracy; however, future work should investigate using different convergence criteria. 

Eqs. (<ref>) and (<ref>) are from Eqs. (<ref>) and (<ref>), respectively, and determine if the spacecraft is in an eclipse or not. If Eq. (<ref>) (resp. Eq. (<ref>)) changes sign while S_ecl,Earth,2 (resp. S_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. (<ref>) is modeled as accurately as possible.

In formulating the cost function, let t_end denote the final time returned by Eq. (<ref>) under the event-detection logic, i.e., t_end always satisfies t_f-Δ t≤ t_end≤ t_f. The first priority of the cost function is to ensure the solution converges. If a solution to Eq. (<ref>) does not satisfy the constraint below,

    |w(x(t_end))|<ϵ,

then, the value of cost, J_1, defined as,

    J_1=w(x(t_end)),

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. (<ref>) becoming negative, then, 

    J_2 = 1000w(x(t_end)),

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

If Eq. (<ref>) is satisfied for a solution, but Eq. (<ref>) is not, then 

    J_3 = -(∑_n = 1^N max(d_ecl,n-90 minutes,0))^-1,

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

Finally, if Eqs. (<ref>) and (<ref>) are satisfied for a solution, then 

    J_4 = -(t_f-t_end)^-1,

is returned as the cost function to PSO. This cost is also made negative to differentiate from Eqs. (<ref>) and (<ref>) and inverted so that time of flight is minimized. However, to ensure it is differentiated from Eq. (<ref>), it has units of 1/[years] while Eq. (<ref>) has units of 1/[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. (<ref>) is returned instead of Eq. (<ref>). 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 <cit.>. To overcome this issue, integration of  is stopped when a certain number of function evaluations is reached. A simple logic is implemented inside  and, if triggered, then Eq. (<ref>) is simply returned as the cost function to PSO. In this paper, 1.0×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 <ref>.







§ RESULTS

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,lb= 2026 NOV 06 00:00:00 UTC and t_f,ub= 2027 JAN 05 00:00:00 UTC, or, t_f,lb=t_f - 30 days and t_f,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 Δ t_lb=200 days and Δ t_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.



 §.§ Fixed NRHO Insertion Date

For this transfer problem, PSO was run 5 times with a swarm size of 500 and a maximum time limit of 1 hour. The best run provided an eclipse-feasible solution, with a time of flight of Δ t= 321.17 days. The trajectory for this solution, in the J2000 ECI frame, is shown in Figure <ref>. The trajectory is also shown in the Earth-centered Sun-Earth rotating frame in Figure <ref> and in the Moon-centered Earth-Moon rotating frame in Figure <ref>. These frames are denoted by the unit vectors {x̂_ECR, ŷ_ECR, ẑ_ECR} and {x̂_MCR, ŷ_MCR, ẑ_MCR}, respectively, where the subscript `ECR' denotes Earth-centered rotating and `MCR' denotes Moon-centered rotating. Note that the legend in Figure <ref> also applies to Figure <ref> and Figure <ref>.





Figures <ref>, <ref>, and <ref> show the time histories of the specific angular momentum, eccentricity, and inclination, respectively. The time histories are plotted in the “forward sense of time,” i.e., the x-axis is 0 when the spacecraft is at the GTO and is Δ t when the spacecraft is at the NRHO. Figures <ref> and <ref> show the CLF value (Eq. (<ref>)) and the CLF time derivative value (Eq. (<ref>)), respectively, and are also plotted in the forward sense of time. Note that because this problem is being solved backwards in time, the sign of the control law was reversed, so, ideally, the function in Fig. <ref> is positive definite. It can be observed though that the CLF time derivative, in fact, becomes negative over multiple intervals. However, the trajectory still converges to the target orbit.








































Figures <ref> and <ref> show the eclipse-triggered throttle factor and eclipse switching functions for the solution. Figure <ref> shows the duration of each eclipse. While included in the model, no eclipses due to the Moon occur in this solution. The first feasible solution from PSO in the same run was obtained in about 15 minutes and had a time of flight of Δ t= 332.72 days. Figure <ref> shows the duration of each of its eclipses. Comparing Figure <ref> and Figure <ref>, it can be interpreted for this particular case that PSO improves the time of flight by increasing the eclipse durations as much as possible. 




































 §.§ Variable NRHO Insertion Date

For this transfer problem, PSO was ran once with an increased swarm size of 1000 to account for the extra parameter, t_f, being optimized. The first eclipse-feasible solution was obtained after about 1 hour and 45 minutes. The NRHO insertion date is t_f= 2026 NOV 06 00:01:35 UTC with a time of flight of Δ t= 314.03 days. The best eclipse-feasible solution was obtained after about 3 hours and 30 minutes and had an NRHO insertion date of t_f= 2026 NOV 06 00:00:00 UTC and a time of flight of Δ t= 303.23 days. The NRHO insertion date of many of the eclipse-feasible solutions tends towards t_f,lb, suggesting that shifting t_f,lb and t_f,ub may provide better solutions. Figure <ref> shows the trajectory for the best solution in the J2000 ECI frame. Figure <ref> shows the duration of each of the eclipses for the solution. An interesting aspect of the solution is that Moon eclipses occur in this transfer solution, with one Moon eclipse  occurring very close to the NRHO insertion being the limiting eclipse duration.



It is hypothesized that it took longer to reach convergence because of the sensitivity introduced by making t_f a decision variable. Small changes in t_f can potentially cause large changes in the initial conditions, x(t_f), depending on how close x(t_f) is to perilune. This means that changes to t_f itself would require changes to the other parameters for the solution to converge. However, in the current formulation, all parameters are being optimized at the same level. It is hard to fully characterize this hypothesis without performing many runs of PSO, since PSO is a stochastic optimization algorithm. One future work is to investigate optimizing t_f using a bi-level optimization algorithm. Nonetheless, considering t_f as a design variable shows that it is possible to obtain solutions with a better time of flight within a new departure window. A potential use of our proposed formulation is to make the bounds t_f,lb and t_f,ub large to find a departure window over a wide range of time, or the bounds can be made small to improve optimality once an initial solution is obtained.




§ CONCLUSION

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.










§ ACKNOWLEDGMENT

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












AAS_publication   
   

