Optimal Impulsive Orbit Transfers from Gateway to Low Lunar Orbit

: Gateway represents a key element of the Artemis program for the upcoming lunar exploration aimed at establishing a sustainable presence by the mid-2030s. This paper investigates minimum-fuel bi-impulsive orbit transfers from Gateway to low lunar orbits (LLOs) with a maximum time of ﬂight of 48 h. Two distinct scenarios are analyzed: (i) target orbits with free right ascension of the ascending node (RAAN), and (ii) target orbits with speciﬁed RAAN. For case (i), a global optimization technique based on a heuristic algorithm is exploited to obtain the minimum-fuel transfer. Several inclinations of the target orbit are considered. For case (ii), two distinct techniques are proposed: (a) a purely heuristic approach, and (b) a semi-analytical method based on local reﬁnement of a Lambert-based solution. Numerical propagations are conducted in all scenarios in a high-ﬁdelity framework that includes all relevant perturbations. A comparison between the different strategies and the related numerical results is provided.


Introduction
NASA's Artemis program represents a paradigm shift in human space exploration, aiming at returning men to the lunar surface and establishing a sustainable presence by the mid-2030s [1].Gateway is a pivotal component of the Artemis program, acting as a permanent orbital outpost around the Moon.This innovative space station will travel a nearrectilinear halo orbit (NRHO), whose stability properties were verified by the CAPSTONE mission [2].Its purpose is to serve as a staging point for crewed missions to the lunar surface, as well as a hub for scientific research and a gateway for future missions beyond the Moon.
The design of transfer trajectories in cislunar space is crucial for successful mission planning.To this aim, different strategies have been investigated for either stable or unstable orbits.The underlying manifold structures can be leveraged for advantageous departure from or approaching toward unstable orbits.Howell et al. [3] proposed low-cost transfers between the Earth-Moon and Sun-Earth systems by means of transit orbits between the respective manifolds.Alessi et al. [4] advanced a two-impulse transfer strategy between low earth orbits (LEOs) and Lissajous orbits around either L 1 or L 2 Lagrangian points in the circular restricted 3-body problem (CR3BP) dynamical framework, exploiting the stable invariant manifold of the target Lissajous orbits.Pontani et al. introduced a polyhedral representation technique for invariant manifolds, useful to classify the trajectories near L 1 and to identify optimal two-impulse and low-thrust transfers from LEOs to Lyapunov orbits [5].Further contributions investigated orbit transfers that connect a variety of Earth orbits to libration point orbits, in particular near-rectilinear halo orbits [6][7][8].In a recent publication, Muralidharan and Howell focused on leveraging stretching directions as a methodology for departure and trajectory design applications [9].
Furthermore, transfer trajectories have been identified for stable orbits, which do not have manifold structures, by means of either low-thrust or impulsive maneuvers.Parish et al. [10] computed low-thrust transfers between distant retrograde orbits (DROs) and halo orbits in the Earth-Moon system.Pino et al. [11] exploited energy surfaces to identify intermediate orbits for low-thrust orbit transfers.Pritchett et al. [12] employed a collocation algorithm for low-thrust transfers.Das et al. [13] illustrated trajectory design techniques using machine learning applied to dynamical structures, combined with numerical optimization.McCarty et al. [14] leveraged the Monotonic Basin Hopping (MBH) optimization algorithm to identify low-thrust transfers from a Gateway-like NRHO to a DRO.Several transfer options have been presented for stable orbits by means of impulsive maneuvers as well.Vutukuri [15] examined intermediate resonant orbits and their manifolds to achieve transfers between stable orbits.Intermediate dynamical structures were employed also by Zimovan et al. [16], using higher-period orbits.
Low lunar orbits (LLOs) provide close-up views of the lunar surface, and represent the most convenient parking orbits for lunar descent and landing [17].Trofimov et al. [18] investigated a variety of two-impulse transfers from Gateway to an elliptic parking orbit, as well as a strategy for direct lunar landing from Gateway, in a high-fidelity dynamical framework.Moreover, high-thrust transfers between Gateway and LLOs are of great interest as they need a limited time of flight (ToF), suitable for crewed missions.However, few studies have investigated this specific topic.Rozek et al. [19] identified two-impulse transfers by means of an elitist non-dominated sorting genetic algorithm and improved solutions using gradient-based optimization.Lu et al. [20] combined local gradient optimization and a numerical continuation strategy to generate optimal direct transfer trajectories, setting the maximum ToF to six days.Bucchioni et al. [21] compared Lambert/differential corrections, Hohmann/differential corrections, and direct numerical optimization to design similar orbit transfers.Zeng et al. [22] focused on impulsive orbit transfers toward a circular LLO in the context of one-way and roundtrip missions from a variety of NRHOs.The great majority of the previously mentioned works (with the remarkable exceptions of Refs.[9,18]) were developed in the dynamical framework of CR3BP, introducing inevitable approximations.
This paper addresses minimum-fuel direct two-impulse transfers from Gateway to circular LLOs at an altitude of 200 km, with a maximum duration of 2 days.These trajectories are especially suitable for crewed missions to the Moon.Unlike Ref. [22], this work employs a high-fidelity dynamical framework based on ephemeris of Gateway and celestial bodies, with the inclusion of all relevant orbit perturbations that affect orbit dynamics.In this work, the term "direct" refers to those transfers that avoid multiple revolutions around the Moon before reaching the target orbit.Specifically, this study has the following objectives: • develop a high-fidelity dynamical framework in cislunar space, to investigate orbit transfers, using modified equinoctial elements for orbit dynamics propagations; • determine optimal two-impulse orbit transfers from Gateway to a set of LLOs with different inclinations, using a new heuristic optimization approach; • determine optimal two-impulse orbit transfers from Gateway to polar LLOs with specified RAAN, using two distinct numerical techniques, i.e., (i) Lambert-based guess generation, and (ii) heuristic guess generation, with subsequent refinement (for both (i) and (ii)); • identify the existing options for optimal orbit transfers, while comparing the related propellant consumption in relation to the final orbit.
The paper is organized as follows.Section 2 describes the reference frames and the orbit dynamics using modified equinoctial elements, with the inclusion of all the relevant orbit perturbations.The orbital motion of Gateway is illustrated in Section 3. Section 4 provides a detailed description of the minimum-fuel direct transfer design, considering RAAN as a free parameter.Section 5 deals with minimum-fuel direct transfers with specified final RAAN.Final remarks are reported in Section 6.

Orbit Dynamics
Section 2 is focused on orbit dynamics.The spacecraft is modeled as a point mass and its motion around the Moon is subject to several perturbations.

Reference Frames
As a preliminary step, a set of reference frames is introduced in order to properly describe the dynamics of the spacecraft.In the following, the notation R j pαqpj " 1, 2, 3q denotes an elementary rotation about axis j by angle α:

•
The local-vertical local-horizontal (LVLH) pr, θ, ĥq frame is attached to the center of mass of the spacecraft.The unit vector r is directed from the Moon center to the instantaneous position of the vehicle; ĥ is aligned with the spacecraft's orbit angular momentum, while θ completes the right-handed set.If Ω denotes the RAAN, i the inclination of the orbital plane, ω the argument of periselenium, and θ the true anomaly, then the following rotations link MCI to the LVLH-frame: • The local-horizontal (LH) frame pr, Ê, Nq is attached to the center of mass of the spacecraft.The unit vectors Ê and N are aligned with the local (lunar) east and north directions, respectively.The LH and LVLH frames can be linked through the following relation: where ζ is the heading angle.
Figure 1 shows the reference frames and the related angles.

Equations of Motion
A set of modified equinoctial elements x i (i " 1, . . ., 6) is introduced in order to describe the motion of the spacecraft, which is modeled as a point mass: x 2 " e cos pΩ `ωq, x 3 " e sin pΩ `ωq, (5) x 5 " tan i 2 sin Ω, where a is the semi-major axis and e the eccentricity.Using the modified equinoctial elements, the dynamical equations are written as where z " " x 1 x 2 x 3 x 4 x 5 ‰ T , a " " a r a θ a h ‰ T is the vector of non-Keplerian accelerations expressed in the LVLH frame, η " 1 `x2 cos x 6 `x3 sin x 6 , while with µ K denoting the lunar gravitational parameter.The time derivative for the orbital element x 6 is instead given by The spacecraft is equipped with a chemical high-thrust propulsion system.As the thrusting arcs have very limited duration with respect to the total transfer time, all the maneuvers are modeled using the impulsive thrust approximation.

Orbit Perturbations
This study addresses the presence of perturbing effects.Other than the lunar gravitational field, which includes several harmonics of the selenopotential, the motion of the spacecraft is affected by the gravitational pull of the Earth and Sun.

Third-Body Perturbations
The gravitational pull of the Earth and Sun can be modeled as third-body perturbations.If r denotes the instantaneous position vector of the spacecraft with respect to the Moon, while r B denotes the vector position of the third body (B) with respect to the Moon, then the perturbing acceleration a 3B due to a third body can be expressed in the MCI frame, following the Battin-Giorgi approach [24,25]: where µ B is the third-body gravitational parameter and If a 3E and a 3S , respectively, denote the third-body accelerations due to the Earth and Sun, then the total perturbing acceleration due to the gravitational pulls of third bodies can be computed as a 3 " a 3E ` a 3S .Finally, by means of Equation ( 1), the acceleration is expressed in the LVLH frame.As a matter of fact, if pa

Harmonics of the Selenopotential
This research employs the Lunar Prospector LP100K model for the Moon's gravitational field, which provides the coefficients J lm and λ lm of zonal, tesseral, and sectoral harmonics up to degree 100 [26].However, only the relevant harmonics are retained in the dynamical model, that is, those associated with coefficients |J lm | ą 5 ˆ10 ´6 [27].These coefficients allow expanding the selenopotential U in terms of Legendre polynomials P l,m : where R K is the lunar equatorial radius, φ the latitude, and λ g the geographical longitude, taken from the lunar reference meridian.This intersects the line that connects the Earth and Moon at all times and is directed toward the far side of the Moon.The gravitational acceleration can then be obtained in the LH frame: The perturbing acceleration a H is Finally, by means of Equation ( 2), it is straightforward to obtain the components of the perturbing acceleration in the LVLH frame.

Orbital Motion of Gateway
Gateway's reference orbit is a southern L 2 NRHO.Several criteria were considered in selecting this orbit [28]: (i) an L 2 -family NRHO provides enhanced visibility of the far side of the Moon for communications compared to the L 1 family [17]; (ii) the L 2 NRHOs enjoy more favorable stability properties compared to the L 1 NRHOs, resulting in lower orbit maintenance requirements [29,30]; (iii) a southern NRHO offers lower delta-v and propellant requirements for spacecraft return to Earth's northern hemisphere; (iv) the southern family also allows very good communications coverage of the lunar south pole, which is a region of high interest for the presence of large water ice deposits in the shadows of craters [31]; and (v) the 9:2 lunar synodic resonance (LSR) with the orbital period of the Moon results in the avoidance of eclipses by the Earth [32].As a matter of fact, since Earth eclipses can last well beyond current hardware limitations, avoiding them was a primary design goal for this reference orbit.While lunar eclipses are relatively common, with several occurring each year, they always last less than 80 min and pose no threat.
The resulting orbit has a periselenium radius that varies from 3196 to 3557 km, with an average value of 3366 km, and an average orbital period of 6.562 days.

Gateway Orbit Dynamics from Ephemeris
The ephemeris of Gateway is contained in an SPK-type file, produced by NASA's Navigation and Ancillary Information Facility (NAIF) [33].The time span of the ephemeris extends from 2 January 2020 at 08:09:36 UTC to 11 February 2035 at 03:59:59 UTC.The numerical model includes the perturbing effects of the gravitational fields of the Earth, Sun, and Jupiter, as well as the harmonics of the selenopotential up to order and degree 8.This orbit considers small aposelenium maneuvers of a few millimeters per second per revolution, which are needed for stationkeeping purposes.As a matter of fact, the orbit diverges after more than 9 revolutions without maintenance [34].Figure 2 shows the reference trajectory in the Earth-Moon synodic reference frame.Its inspection reveals a significant impact of the perturbations: the trajectory, which is periodic in the CR3BP, becomes quasi-periodic in the accurate dynamical framework.Therefore, stationkeeping is necessary to maintain the orbit close to the reference one, obtained in CR3BP.The existence of a periodic solution in the CR3BP remains as a quasi-periodic solution even when perturbations are included, in agreement with the Kolmogorov-Arnold-Moser (KAM) theorem [35][36][37].Indeed, this theorem guarantees the persistence of quasi-periodic motions under small perturbations.

Gateway Orbit Dynamics from High-Fidelity Model
Gateway's trajectory can also be reproduced using the high-fidelity model, i.e., by integrating Equation ( 9).This allows an evaluation of the accuracy of the model and a deeper analysis of the orbital motion of the space station.
The integration is performed using MATLAB 2023b ode113, with both absolute and relative tolerances set to 10 ´13 .The time frame selected for the integration is slightly longer than the period of the station's orbit around the Moon, starting on 17 May 2025 at 10:00:00 UTC and ending on 24 May 2025 at 10:00:00 UTC.The high-fidelity model provides insight into the total perturbing acceleration a p , which is shown in Figure 4. Its maximum value corresponds to the aposelenium, which is located at about 3.5 days.Furthermore, the average magnitude (3.896 ˆ10 ´4 m/s 2 ) is close to the thrust acceleration produced by a low-thrust engine.This confirms that cislunar space is a highly perturbed environment.Moreover, it is worth pointing out that the osculating RAAN varies by 360 degrees in about 28 days.This behavior is shown in Figure 5.

Minimum-Fuel Orbit Transfers (Free RAAN)
This section addresses minimum-fuel two-impulse transfers from Gateway to circular LLOs, thereby minimizing the total change in velocity.Several inclinations of the target LLO are considered, with the RAAN left as a free parameter, in order to identify the most convenient orientation of the orbital plane.

Formulation of the Problem
The selected target LLO is a circular orbit, with an altitude of h ˚" 200 km, i.e., with radius r c :" R K `h˚.Ten equally spaced inclinations are considered in the range r0, 90s degrees.Two three-dimensional velocity changes, denoted as ∆ V 1 and ∆ V 2 , are applied according to the following transfer strategy: • ∆ V 1 is applied at the initial time t 0 along the orbit of Gateway to inject the spacecraft into a transfer trajectory.Consequently, the spacecraft's dynamical state after the maneuver is given by where ( X 0 , V 0 ) and ( X 0 , V 0 ), respectively, indicate the position and velocity vectors before and after the instantaneous application of the velocity change at t 0 ; • ∆ V 2 is applied at the end of the transfer time to insert the spacecraft into the desired LLO.
The transfer design aims to determine the optimal components of ∆ V 1 and the optimal initial time t 0 to start the transfer.The optimization problem is formulated as Hence forward, ∆V 1 and ∆V 2 will indicate the magnitudes of the two velocity changes.The method of solution is presented in the following paragraph.

Method of Solution
As a first step, the desired inclination i d of the final LLO is selected.The problem is solved through numerical optimization, which iteratively performs the following steps: 1.
After applying ∆ V 1 at t 0 , the spacecraft is inserted into the transfer.Its state is propagated forward in time according to the equations of motion of the high-fidelity model (cfr.Section 2.2), with a maximum ToF of 48 h.The integration ends when the target sphere of radius r c is reached.

2.
The value of Ω f for the target orbit can be computed from where r f is the spacecraft position vector at the end of the transfer, whose components are expressed in the MCI frame, while ĥ " identifies the direction of the angular momentum associated with the orbit containing the final spacecraft position.

3.
Equation ( 21) yields two solutions for the RAAN, denoted by (Ω 1 , Ω 2 ).The corresponding values of the argument of latitude (θ t,1 , θ t,2 ) are then calculated.For each solution (Ω k , θ t,k ), the comparison point P k along the LLO is determined via a single-parameter search, aimed at minimizing the cost function where ∆r is the distance of P k from the final position (identified with r f q, while k r and k 2 are weighting coefficients with proper units.The single parameter of interest is the argument of latitude, identifying P k along the desired final LLO.

4.
The primary optimization is conducted using a cost function defined as where ∆h represents the difference between the altitude reached at the end of the transfer and the desired altitude h ˚, while c 1 , c 2 , c h , and c r are weighting coefficients with proper units.The term ∆h guides the trajectory toward the target sphere.When i d ă 90 deg, including ∆r in the cost function yields convergence when the intersection between the transfer and target spheres occurs at latitudes higher than the inclination of the target orbit.

5.
The cost function J is evaluated for each pair pΩ k , θ t,k q, k " 1, 2, and the values that minimize J are selected.
The entire optimization process consists of two consecutive steps: (i) A first global optimization of J is performed by means of the Stochastic Fractal Search (SFS) algorithm.Its stochastic approach enhances the exploration phase and reduces the possibility of becoming stuck in local minima.A detailed description of the SFS algorithm is reported in Appendix A. (ii) Then, the MATLAB built-in routine fmincon is employed for solution refinement, with an equality constraint on the final altitude, and an inequality constraint on the latitude of r f , i.e., |φpt f q| ď i d .The latter inequality constraint ensures that an LLO with inclination i d passes through r f .
For the numerical simulations, the following canonical units are employed: • The distance unit (DU) is equal to the lunar mean equatorial radius, DU " 1738 km; • The time unit (TU) is chosen to yield a unit gravitational parameter, TU " with µ K " 4902.8 km 3 /s 2 .

Numerical Results
Figure 6 shows ∆V tot " ∆V 1 `∆V 2 , with respect to the inclination of the LLO.It is clear that ∆V tot increases as the inclination decreases.This behavior is expected since the osculating inclination remains relatively close to 90 deg.Therefore, transfers to polar or near-polar LLOs require a smaller change in orbital plane and reduced fuel consumption as a result.Table 1 collects the relevant data for the optimal transfers.In particular, the overall velocity change for the final polar LLO is very close to the values found by Bucchioni and Innocenti [21] and Lu et al. [20], who found the values 661 m/s and 650 m/s, respectively.However, the latter works assumed a target LLO similar though not identical to that considered in this work, while using the CR3BP model.Figure 7 shows the transfer trajectories in the MCI frame for the two limiting cases, i.e., polar and equatorial LLOs.In both cases, the orbital plane of the transfer is close to the orbital plane of the target LLO.

Minimum-Fuel Orbit Transfers (Specified RAAN)
The aim of this section is to identify minimum-fuel transfers from Gateway to polar circular LLOs with a maximum ToF of 48 h.Unlike the previous analysis, the RAAN value is no longer a free parameter.Eight scenarios are analyzed, corresponding to different specified RAAN values, evenly distributed in the range r0, 360s degrees.Because the osculating RAAN varies by 360 degrees in about 28 days (cf. Figure 5), the transfer starting time is sought in the same interval.Polar orbits are selected as target LLOs because (i) polar inclination requires a lower ∆V tot (cf.Table 1), and (ii) they guarantee the coverage of lunar poles, especially the South pole, which harbors vital water ice deposits [38].The problem is solved through two different methods, which are described below.

Formulation of the Problem
The goal is to minimize the total ∆V for a bi-impulsive transfer from Gateway to target LLO.The two velocity changes ∆ V 1 and ∆ V 2 are applied as with the previous strategy (cf.Section 4.1).The final orbit must not only be circular but also polar, with the specified RAAN value Ω d .The argument of latitude along the target LLO remains a free parameter.The optimization problem can be formulated as min For this purpose, a metaheuristic method could be directly employed, as in Section 4.However, an alternative two-step approach is here proposed: (i) identifying a preliminary solution through a global optimizer using the Lambert algorithm assuming Keplerian motion (cf.Section 5.2); and (ii) refining the solution through a local optimizer in a high-fidelity dynamical framework (cf.Section 5.2).This method is described in the following paragraph.

Method 1: Lambert-Based Guess Generation and High-Fidelity Refinement
The initial guess solution is sought through a Lambert problem solving algorithm [39], which allows the determination of the conic connecting two points in three-dimensional space for a given time of flight ∆t.Given the initial and final positions along with a time span, the algorithm efficiently computes the Keplerian transfer trajectory.
In our problem, the initial position and velocity of the spacecraft along the orbit of Gateway are uniquely determined by the designated start time t 0 .The final state at time t f is specified by the argument of latitude θ f of the LLO, whose orientation is fixed by inclination and RAAN.
Therefore, the algorithm requires the following inputs: • The initial state p X 0 , V 0 q, identified by t 0 ; • The final state p X f , V f q, identified by θ f ; • The transfer time, ∆t " t f ´t0 .
Both state vectors are referenced in the MCI frame.The algorithm outputs the orbit elements of the transfer trajectory, allowing the determination of ∆ V 1 and ∆ V 2 for Gatewayto-transfer orbit insertion at t 0 and for transfer orbit-to-LLO insertion at t f .
The optimization problem, as formulated in Equation ( 27), is solved by means of the SFS, with the two velocity changes determined by the Lambert algorithm.Its introduction is motivated by the immediate determination of the transfer trajectory within a Keplerian framework.Consequently, propagation becomes unnecessary when SFS is used for optimization, making the search for a solution more efficient.
The initial guess, which solves the problem ( 27) in a Keplerian model, is employed to refine the solution within the high-fidelity framework (cf.Section 2.3).This refinement involves solving the optimization problem (20) through MATLAB fmincon.
While the guess value t 0 is directly obtained from the global optimization procedure, the three components of ∆ V 1 in the MCI frame are computed through the Lambert algo-rithm, starting from the previous optimal solution.The optimization process advances by propagating the perturbed motion of the spacecraft using Equation ( 9) after the impulsive application of ∆ V 1 .Due to the non-Keplerian dynamical environment, reaching the final position at θ f within ∆t is unlikely.Consequently, the propagation is stopped when an altitude of 200 km is reached or at the periselenium of the transfer orbit.Finally, ∆ V 2 is applied.This method ensures targeting a circular and polar orbit, although it may not meet the prescribed RAAN condition.Therefore, within fmincon, the difference between the final and the desired RAAN is enforced as a nonlinear constraint.

Numerical Results Using Method 1
As a first step, the SFS algorithm is employed to obtain the initial guess as a solution of the Keplerian Lambert problem.The tuning parameters are set as follows:  From inspection of Table 2, it is evident that most trajectories yield the maximum ToF (48 h), corresponding to the upper bound of the search interval, as a result of the optimization process.
Secondly, the fmincon routine is employed to refine the solution in the high-fidelity model, while enforcing the nonlinear constraint on the final RAAN: ˇˇ1 ´cos pΩ d ´Ω f q ˇˇď 1.5 ¨10 ´6.(28) This tolerance yields an error on the final RAAN less than 0.1 deg.The search space is identified by lower and upper bounds on the variation of the initial guess time t pGq 0 (rt pGq 0 ´δt 0 , t pGq 0 `δt 0 s), as well as on the ∆ V 1 components along the LVLH frame (respectively, adding ˘δV r , ˘δVθ, ˘δV h ), where δt 0 " 1 h, δV r " δV θ " δV h " 50 m/s.However, a further reduction in the search space was necessary in order to guarantee convergence for the cases associated with Ω d " 270, 315 deg.Table 3 collects the results of the solution refinement.The algorithm is capable of generating feasible transfers that satisfy the requirements on the final orbit.
Figure 8 shows the stream of trajectories associated with different RAAN values of the target LLO.Interestingly, the trajectories associated with the lowest ∆V tot , e.g., Ω d " 90 or 270 deg (cf.Table 3) exhibit a reduced change in direction after the first velocity change; instead, the greatest values of ∆V tot are associated with a significant change in orbital plane.Moreover, orbit injection always occurs in the north hemisphere.

Method 2: High-Fidelity Heuristic Guess and Local Refinement
The heuristic algorithm (cf.Section 4) is here exploited to obtain the global optimal parameters for the sake of comparison.However, as Ω f is fixed, there is no need to solve Equation (22).Thus, a single value is determined for θ t and for the cost function.The SFS algorithm is initially set to the same parameters as in Section 5.3.However, as the cost function exhibits several local minima, some scenarios required several attempts and an appropriate tuning of the parameters in order to increase the exploration capabilities of the algorithm.The weighting coefficients k 1 , k 2 and c 1 , c 2 , c r , c h (cf.Section 4.2) are set to unit values.As illustrated in Section 4.2, the solution of global optimization is refined by means of fmincon.

Numerical Results Using Method 2
Table 4 collects the outcomes of the optimization process based on high-fidelity heuristic guess generation and subsequent local refinement.

Comparison of Results
Figure 9 shows a comparison between the two strategies for each scenario.Remarkably, both approaches yield very close results in terms of ∆V tot .A close inspection of Tables 3 and 4 highlights that, even though the starting times may be different, they provide similar ∆V tot .In particular, the Lambert algorithm is capable of providing a good initial guess solution, though in a Keplerian dynamical framework.As a matter of fact, the semi-analytical approach is proven to guarantee a faster convergence with respect to the heuristic method, as the global exploration does not require the integration of the high-fidelity equations of motion.As an example, a single iteration of the SFS algorithm on a 48-core computer requires on average 1.58 s for the Lambert-based guess generation and 67.17 s for the heuristic optimization.On the other hand, the accuracy of the Lambert-based guess solution is expected to decrease when trajectories associated with longer transfer times are sought, because of the action of perturbations.

Concluding Remarks
This research addresses minimum-fuel orbit transfers from Gateway to LLOs associated with a maximum duration of 48 h.The spacecraft is equipped with a chemical engine.The dynamical framework includes all relevant gravitational perturbations, including third-body effects of the Earth and Sun.Two different cases are explored: first, a global optimization technique based on a heuristic algorithm is developed to obtain the optimal transfer toward an LLO with free RAAN.Several scenarios are investigated, associated with different inclinations of the target orbit; the numerical results prove that the algorithm is capable of identifying the optimal transfer leading to the desired orbit.Unsurprisingly, polar target orbits are reached with lower fuel expenditure.Then, a semi-analytical ap-proach is developed, aimed at achieving the optimal minimum-fuel transfer toward a polar target orbit with specified RAAN.This strategy relies on a numerical solver based on the Lambert theorem, which provides an initial guess.This solution is then refined by local numerical optimization in the high-fidelity dynamical framework.The results attained with this semi-analytical approach are shown to be very similar (in most cases) to those found with high-fidelity heuristic guess generation.The latter approach, however, turns out to be much more computationally expensive.one that minimizes the fitness function) is selected.The central idea of SFS consists in generating individuals for the next generation through: • An initial diffusing phase, in which the existing individuals are "diffused" around their current position.To generate the branch-like fractal shape, the diffusion-limited aggregation (DLA) [41] model is employed, by exploiting the Gaussian walk statistical technique.Only the best individual generated from this diffusion is considered, while the others are discarded in order to avoid a dramatic increase in the number of individuals.• A subsequent updating phase, in which the less well performing candidates are replaced by other individuals, randomly generated.This updating step aims at introducing a way for the individuals to exchange information, to speed up convergence toward the optimal solution.
Unconstrained parameter optimization problems can be stated as follows: determine the optimal values of the n unknown parameters χ 1 , . . ., χ n such that the fitness function J is minimized.The time evolution of the dynamical system under consideration depends on χ 1 , . . ., χ n , which are constrained to their respective ranges, The initial population of SFS is randomly generated by introducing N individuals, whose parameters are (stochastically) uniformly distributed in the respective search spaces, defined by Equation (A1).The steps that follow comprise the generic iteration (denoted with index g), which leads to generating a new populations of individuals.

2.
The best individual, which corresponds to the current minimum value of the fitness function, is identified and is associated with index i opt .3.
Diffusion phase.For i " 1, . . ., N: (a) The "diffusion limited aggregation" growth process employs the Gaussian distribution to generate new individuals; N D is the number of individuals generated by diffusion of individual i.For j " 1, . . ., N D , Ypi, jq " # gpµ opt , σq `pε χ pgq pi opt q ´ε1 χ pgq piqq if ε 2 ă P GW where i opt " arg min i"1, ..., N Jpgq i , P GW denotes a threshold value, whereas ε, ε 1 , and ε 2 are three random variables with uniform distribution in r0, 1s.Symbol gpµ, σq denotes a vector with random components g k , subject to Gaussian distribution, with mean value µ k (component of µ) and standard deviation σ " ˇˇˇl og g g ´χpgq piq ´χpgq pi opt q ¯ˇˇˇ, (A4) The term log g{g aims at decreasing the size of the Gaussian jumps, to promote more localized search as the generation index g increases.
respectively, denote the components of a 3 in the MCI and LVLH frames, then

Figure 2 .
Figure 2. Gateway trajectory in the Earth-Moon synodic reference frame, for ten periods.

Figure 3
Figure3displays the relative distance ∆r between the Gateway position obtained by numerical integration in the high-fidelity model and the ephemeris data.The mismatch is mainly due to a different model for the lunar gravity field, which has a large effect around periselenium passages.In addition, the high-fidelity model does not include any aposelenium maneuvers, unlike the ephemeris data.

Figure 3 .
Figure 3. Norm of Gateway position difference between high-fidelity model and ephemeris data for seven days.

Figure 4 .
Figure 4. Norm of the total perturbing acceleration acting on Gateway for seven days.

Figure 5 .
Figure 5.Time history of the RAAN in 28 days.

Figure 6 .
Figure 6.∆V tot as a function of LLO inclination.

Figure 7 .
Figure 7. Transfers from Gateway NRHO toward an LLO, in the MCI frame.

Figure 8 .
Figure 8. Stream of trajectories associated with different RAAN of the target orbit (with zoom on the injection positions in the inset).The trajectory of Gateway is depicted with a dashed black line.

Figure 9 .
Figure 9.Comparison between ∆V tot for different values of RAAN.

Table 1 .
Values of ∆V tot , starting epoch, and ToF according to the inclination of the LLO.

Table 2
collects the results of the guess generation based on the Lambert algorithm.

Table 2 .
Initial guess solutions obtained through Lambert algorithm.

Table 3 .
Results of solution refinement using fmincon.

Table 4 .
Refined results of heuristic global optimization.