# Second-order Averaging of Low-Thrust Transfers in Fixed Time

Lamberto Dell'Elce, Alesia Herasimenka, Aaron J Rosengren, Nicola Baresi

► To cite this version:

Lamberto Dell'Elce, Alesia Herasimenka, Aaron J Rosengren, Nicola Baresi. Second-order Averaging of Low-Thrust Transfers in Fixed Time. Journal of Guidance, Control, and Dynamics, In press. hal-04387944

HAL Id: hal-04387944

https://inria.hal.science/hal-04387944v1

Submitted on 11 Jan 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Second-order Averaging of Low-Thrust Transfers in Fixed Time*

Lamberto Dell'Elce† and Alesia Herasimenka‡

Inria center of University Côte d'Azur, Sophia Antipolis, 06902, France

Aaron J. Rosengren§

University of California San Diego, La Jolla, CA 92093, USA.

Nicola Baresi¶

University of Surrey, Guildford, GU2 7XH, UK

An algorithm for the second-order approximation of fixed-time, low-thrust orbital transfers and generation of Bacon plots is proposed. After averaging the extremal flow of the optimal-control Hamiltonian, a one-parameter family of solutions of a reduced-order, two-point boundary value problems is computed by means of a differential-continuation scheme. Sensitivities of the shooting function are then used in conjunction with an *ad hoc* near-identity transformation between averaged and osculating variables to achieve an accurate solution for all longitudes of the departure and arrival orbits. Hence, a single simplified shooting problem has to be solved to approximate the solution for any combination of departure and arrival dates, and to draw a Bacon plot for an arbitrary launch window. Both the averaged flow and the near-identity transformation are efficiently evaluated via the fast-Fourier transform algorithm, yielding a fully-numerical procedure.

## Nomenclature

* $\epsilon$ Specific thrust magnitude [$\frac{N}{kg}$]

* $m$ Residual mass ratio [–]

* $q, Q$ Osculating and averaged sets of Keplerian integrals of motion

* $\phi, \Phi$ Osculating and averaged mean longitude [deg]

* $t_0, t_f$ Departure and arrival dates

* $\beta$ Inverse of the effective exhaust velocity [$\frac{s}{m}$]

* $u$ Control variable

* $f_j, g_j$ $j$-th vector fields of Gauss variational equations

<footnote>*Presented as Paper 23-254 at the 2023 AAS/AIAA Spaceflight Mechanics Meeting, Austin, TX, January 15-19, 2023.</footnote>

<footnote>†Researcher, McTAO team & Jean Alexandre Dieudonné Laboratory (LJAD); lamberto.dell-elce@inria.fr.</footnote>

<footnote>‡PhD student, McTAO team & Jean Alexandre Dieudonné Laboratory (LJAD).</footnote>

<footnote>§Assistant professor, Department of Mechanical and Aerospace Engineering.</footnote>

<footnote>¶Lecturer in Orbital Mechanics, Surrey Space Centre.</footnote>
$\omega$ Mean motion [$\frac{\text{rad}}{\text{s}}$]

$\mathcal{H}$ Optimal control Hamiltonian

$\mathbf{p}, \mathbf{P}$ Costate variables of $\mathbf{q}$ and $\mathbf{Q}$

$\lambda, \Lambda$ Costate variables of $\varphi$ and $\Phi$

$\gamma, \Gamma$ Osculating and averaged costate variables of $m$

$S$ Shooting function

$\nu(\cdot)$ Short-periodic variation of $(\cdot)$

$\delta(\cdot)$ Second-order correction of $(\cdot)$

$\tau$ Re-scaled time [$\frac{m}{s}$]

$i$ Imaginary unit

# I. Introduction

Lambert's problem, in its classical form, consists in finding a Keplerian orbit joining two position vectors in a given transfer time. Solutions of this problem are extensively used for preliminary mission design since they offer the identification of launch opportunities and a rough evaluation of their fuel cost by assuming impulsive maneuvers at the two boundary points. Specifically, the concept of "launch windows" and the use of graphical representations to find feasible trajectories for planetary missions were introduced in the early space age [1]. These graphical representations, which include what are now called pork-chop plots, help mission planners to visualize and analyze the trade-offs between departure and arrival dates, taking into account the positions and velocities of the planets involved. Over the years, these plots were refined and adapted for specific missions, and are now standard tools in the field of space exploration. Pork-chop plots are often generated to illustrate total $\Delta V$ cost as a function of the departure and arrival dates by recursively solving Lambert's problems. When considering the low-thrust counterpart of this transfer problem in fixed time, the entire state vector is imposed at the two boundaries, i.e., both position and velocity of the satellite have to match the ones of departure and arrival bodies, because impulsive maneuvers are not allowed. In contrast to the original problem, no exact closed-form solution is available, so that assumptions on the shape of the orbit and thrust direction are introduced to achieve efficient approximations [2–4]. A similar, shape-based approach was used in [5] to draw the equivalence of pork-chop charts for low-thrust problems, which are referred to as Bacon plots [6]. Multiple-shooting techniques where used in [7] to tackle optimization of low-thrust transfers; however, this approach allowed the solution of only a single point of a Bacon plot. Averaging techniques were applied to the solution of low-thrust orbital transfers in [8, Chapter 5] and [9], but the costate of the longitude was set to zero all over the trajectory, yielding poor accuracy on the departure and arrival longitudes of the maneuver, as discussed in [10]. Averaging techniques for problems with orbital perturbations, shadowing effects, and constraint on the thrust direction were used in [11], but, again, the role
of the costate of the fast variable was ignored. This simplification is not introduced in the present paper, and we will
show that it entails major benefits on the approximation of transfers occurring in “few” revolutions (about ten), which
are maneuvers essentially ignored in the literature on averaging. Direct techniques were also used in conjunction with
averaging in [12]. In this case, a harmonic parametrization of the control variable is introduced first, and optimization
of the coefficients is carried out without introducing costate variables.

In our view, two problems need to be addressed when tackling low-thrust transfers for a fixed maneuvering time.
First, the minimum thrust magnitude necessary to carry out the maneuver has to be identified. Solutions to this problem
are characterized by the absence of coasting arcs and the exploitation of the maximum control force throughout the
entire trajectory. Second, once a sufficiently-large thrust is chosen, minimum-energy maneuvers can be found [13–15].
This work focuses on the first problem, which is of interest because it offers a lower bound on the thrust that cannot
be violated when optimizing other cost functions and it may serve as initial guess for the minimum fuel problem. A
numerical methodology based on the averaging of the extremal flow of the optimal-control system [10] is proposed.
First, a reduced-order solution of the averaged two-point boundary value problem (TPBVP) parametrized by the costate
of the fast variable is solved. This step requires the solution of a single shooting problem followed by a numerical
continuation procedure. This problem is independent of the thrust magnitude. Second, the sensitivity of the shooting
function are computed. These sensitivities are then used to compute perturbations of the averaged TPBVP associated to
short-periodic variations and second-order terms, which are hereby retained to obtain a first-order approximation of the
fast variable from the averaged solution. This information is finally used to find the minimum thrust required for the
transfer. Steps 1 and 2 are independent of the boundary conditions on the longitude, which is a major advantage of the
methodology since these are the only two computationally demanding steps. In this way, a single one-parameter family
of solutions of the averaged problem is sufficient to approximate the solution of the original problem for any departure
and arrival windows and any resolution of the grid.

The orbital-transfer problem that we tackle is formulated in Section II. Its averaged counterpart, the equations of
motion of first and second-order terms, and formal solutions of short-periodic variations are deduced in Section III.
These ingredients are used to develop an algorithm for the generation of Bacon plots, which is detailed in Section IV.
Finally, the methodology is applied to a geostationary transfer orbit (GTO) to geostationary orbit (GEO) transfer, and an
Earth-Venus transfer in Section V. Discontinuities of the minimum thrust magnitude map are discussed and related to the
cut locus of the problem, which is an aspect that was not covered in the literature on multi-revolution orbital transfers.
## II. Low-thrust orbital transfer in fixed time

Denote by $\varphi \in \mathbb{S}^1$ the mean longitude and by $\mathbf{q} \in \mathcal{M} \subset \mathbb{R}^n$ a set of integrals of motions of the two-body problem. Specifically, we use equinoctial elements, namely,

$$ \mathbf{q} = \left[ a, e \cos \theta, e \sin \theta, \tan \frac{i}{2} \cos \Omega, \tan \frac{i}{2} \sin \Omega \right]^T, \quad (1) $$

where $a, e, \theta, i, \Omega$ denote the semi-major axis, eccentricity, longitude of periapsis, inclination, and longitude of ascending node of the orbit, respectively. Hereafter, $\mathbf{q}$ and $\varphi$ are referred to as slow and fast variables, respectively. Let $t_0$ and $t_f$ be the desired departure and arrival dates. We are interested in finding the minimum thrust magnitude necessary to carry out the transfer between two Keplerian orbits with elements $(\mathbf{q}_0, \varphi_0(t_0))$ and $(\mathbf{q}_f, \varphi_f(t_f))$ in a time $t_f - t_0$. No orbital perturbation is considered. The notation $\varphi_0(t_0)$ denotes the value of the mean longitude of the departure orbit at epoch $t_0$. For example, if an interplanetary transfer is envisaged, $\varphi_0(t)$ is given by the ephemeris of the departure planet at time $t$. Analogous notation holds for $\varphi_f(t_f)$. The optimal-control problem that we tackle is thus

min $\varepsilon$ subject to:

$$ \begin{align} \frac{d\mathbf{q}}{dt} &= \frac{\varepsilon}{m} \sum_{j=1}^{3} \mathbf{f}_j(\mathbf{q}, \varphi) u_j, \notag \\ \frac{d\varphi}{dt} &= \omega(\mathbf{q}) + \frac{\varepsilon}{m} \sum_{j=1}^{3} g_j(\mathbf{q}, \varphi) u_j, \notag \\ \frac{dm}{dt} &= -\varepsilon\beta \|\mathbf{u}\|, \tag{2} \\ \|\mathbf{u}\| \le 1 \qquad \forall t \in [t_0, t_f], \notag \end{align} $$

$$ \mathbf{q}(t_0) = \mathbf{q}_0, \quad \varphi(t_0) = \varphi_0(t_0), \quad m(t_0) = 1 $$

$$ \mathbf{q}(t_f) = \mathbf{q}_f, \quad \varphi(t_f) = \varphi_f(t_f); $$

where $\mathbf{u}$ is the control variable taking values in the unit sphere of $\mathbb{R}^3$, $\varepsilon$ is the specific thrust, which is obtained by dividing the thrust magnitude by the initial mass of the satellite, $m$ is the instantaneous to initial mass ratio of the satellite, and $\beta$ is the inverse of the effective exhaust velocity*. For low-thrust propulsion, $\beta \ll 1$. Vector fields $\mathbf{f}_j$ and $g_j$ are periodic with respect to the fast variable, and they stem from the Gauss variational equations (GVE) of the chosen set of orbital elements.

<footnote>*Namely, $\beta = \frac{1}{g_0 I_{sp}}$, where $g_0 = 9.81 \frac{m}{s^2}$ and $I_{sp}$ denote the Earth's gravitational acceleration at sea level and the specific impulse of the thruster, respectively.</footnote>
The pre-Hamiltonian of Problem (2) is defined as

$$ \mathcal{H}' = \lambda \omega(\mathbf{q}) + \frac{\varepsilon}{m} \sum_{j=1}^{3} [\mathbf{p} \cdot \mathbf{f}_j(\mathbf{q}, \varphi) + \lambda g_j(\mathbf{q}, \varphi)] u_j - \varepsilon \beta \gamma \|\mathbf{u}\|, \quad (3) $$

where $\mathbf{p}$, $\lambda$, and $\gamma$ denote adjoints of the slow and fast variables, and of the mass, respectively. Application of the Pontryagin maximum principle (PMP) [16] yields the control action of the extremal flow $\mathbf{u}(\mathbf{q}, \varphi, \mathbf{p}, \lambda)$, namely

$$ \mathbf{u}^* = \arg \max_{\|\mathbf{u}\| \le 1} \mathcal{H}'(\mathbf{q}, \mathbf{p}, \varphi, \lambda, m, \gamma, \varepsilon, \beta, \mathbf{u}) = \begin{cases} \frac{\mathbf{K}}{K} & \text{if } K - \beta \gamma m \ge 0, \\ \mathbf{0} & \text{otherwise,} \end{cases} \quad (4) $$

where the j-th component of $\mathbf{K}$ is $K_j = \mathbf{p} \cdot \mathbf{f}_j(\mathbf{q}, \varphi) + \lambda g_j(\mathbf{q}, \varphi)$, and $K := \|\mathbf{K}\|$. The transversality condition on $\gamma(t_f) = 0$ (stemming from not constraining $m(t_f)$) and the non-negativity of the time derivative of $\gamma$, namely

$$ \frac{d\gamma}{dt} = \begin{cases} \varepsilon\beta\frac{K}{m^2} \ge 0 & \text{if } K - \beta\gamma m \ge 0, \\ 0 & \text{otherwise,} \end{cases} \quad (5) $$

imply that $\gamma(t) \le 0 \forall t \in [t_0, t_f]$ and, consequently, that the switching condition $K - \beta\gamma m = 0$ is never verified all over the transfer. Hence, replacing $\mathbf{u}^*$ into Eq. (3) yields the Hamiltonian of the extremal flow:

$$ \mathcal{H} = \lambda \omega(\mathbf{q}) + \varepsilon \left( \frac{K(\mathbf{q}, \mathbf{p}, \varphi, \lambda)}{m} - \beta\gamma \right). \quad (6) $$

Transversality conditions are such that adjoints are free at both $t_0$ and $t_f$ [17], so that necessary conditions for the optimality of Problem (2) are satisfied by finding a zero of the shooting function $S(\mathbf{p}(t_0), \lambda(t_0), \varepsilon)$ defined as

$$ S(\mathbf{p}(t_0), \lambda(t_0), \varepsilon) = \begin{bmatrix} \mathbf{q}^*(t_f) - \mathbf{q}_f \\ \varphi^*(t_f) - \varphi_f(t_f) \\ \|\mathbf{p}(0)\| - 1 \end{bmatrix}, \quad (7) $$

where $\mathbf{q}^*(t)$ and $\varphi^*(t)$ denote trajectories of $\mathbf{q}$ and $\varphi$, respectively, emanated by $\mathcal{H}$ with initial conditions $\mathbf{q}_0$, $\varphi_0$, $\mathbf{p}(t_0)$, $\lambda(t_0)$, and small parameter $\varepsilon$. The last equation of $S$ is an arbitrary normalizing condition due to the homogeneity of $\mathcal{H}$ with respect to the adjoints [18]. The condition $\gamma(t_f) = 0$ and the related shooting variable $\gamma(t_0)$ are not included in Eq. (7), since the value of $\gamma$ has no impact on the equations of motion of the other variables. If desired, $\gamma(t_0)$ can be inferred a posteriori by backward integration of Eq. (5). Finally, we note that the inequality $m(t_f) > 0$ should also be verified. This can be automatically achieved by means of the change of time variable detailed in Section III.A.
### III. Averaged problem

States and adjoints are decomposed as

$$ 
\begin{aligned}
q &= Q + \varepsilon [\delta Q + v_Q(Q, P, \hat{\varphi})], \\
p &= P + \varepsilon [\delta P + v_P(Q, P, \hat{\varphi})], \\
\varphi &= \varphi_0 + \frac{\Phi}{\varepsilon} + \delta\Phi + \varepsilon v_\Phi(Q, P, \hat{\varphi}), \\
\lambda &= \varepsilon (\Lambda + v_\Lambda(Q, P, \hat{\varphi})) + \varepsilon^2 \delta\Lambda (Q, P, \hat{\varphi}, \Lambda, \delta Q, \delta P, \delta\Gamma), \\
\gamma &= \Gamma + \varepsilon [\delta\Gamma + v_\Gamma(Q, P, \hat{\varphi})].
\end{aligned}
\quad (8) $$

Hereafter, $Q$ and $\Phi$ are referred to as the slow and fast averaged states, respectively, $P$ and $\Lambda$ are their adjoints, and $v_\cdot$ denotes short-periodic variations. No decomposition of the mass variable is introduced since its trajectories are linear functions of time. Second-order terms,† namely $\delta Q$, $\delta P$, $\delta\Phi$, $\delta\Lambda$, and $\delta\Gamma$, are necessary to have a precise enough approximation of the fast variable; that is such that its error is $O(\varepsilon)$ [19]. The power of $\varepsilon$ that multiplies each term of the decomposition (i.e., $\varepsilon^{-1}$, $\varepsilon^0$, $\varepsilon^1$, $\varepsilon^2$) is used to make explicit each respective formal order of magnitude, so that $Q, P, \Phi, \Lambda, \Gamma, \delta Q, \delta P, \delta\Phi, \delta\Lambda, \delta\Gamma$, and any $v_\cdot$ (once stripped out of their power of $\varepsilon$) are all $O(1)$ throughout the entire trajectory. Specifically, the gross behavior of the fast variable is captured by $\Phi$, so that $\varepsilon\varphi(t) = \Phi(t) + O(\varepsilon) = O(1)$ when $t = O(\varepsilon^{-1})$. Similarly, the adjoint $\Lambda$ is well defined when $\varepsilon$ approaches zero as discussed in [10], where $\lambda$ is shown to be of order $\varepsilon$ throughout the entire trajectory. Short-periodic variations and $\delta\Lambda$ are evaluated at the angle:

$$ \hat{\varphi} = \varphi_0 + \frac{\Phi}{\varepsilon} + \delta\Phi. \quad (9) $$

The decomposition in Eq. (8) is not standard since second-order averaging is generally carried out by including second-order terms in the equations of motion of the averaged state without introducing new variables such as $\delta Q$ and $\delta P$. The main advantage of the proposed decomposition is that both the averaged counterpart of Eq. (7) and the computation of second-order terms are $\varepsilon$-independent. We also note that second-order short-periodic variations are neglected since they do not affect the order of magnitude of the error in $\varphi$, so that Eq. (8) is not a full second-order approximation.

The averaging operator acting on a generic function $f(\varphi)$, periodic with respect to $\varphi$, is denoted with an overline, and it is defined as

$$ \bar{f} := \frac{1}{2\pi} \int_{0}^{2\pi} f(\varphi) \, d\varphi. \quad (10) $$

When the operator is applied to a function with other arguments than $\varphi$; these arguments are kept constant in the integral.

<footnote>†The reference "second order" may stride with the order of magnitude of these terms in the expansion of Eq. (8). Instead, they are referred to as second-order terms because time derivatives of $\varepsilon\delta Q$ and $\varepsilon\delta P$ are $O(\varepsilon^2)$. For the sake of clarity and without arguing further on this notation, $\delta\Phi$ and $\delta\Lambda$ are also referred to as "second-order terms" for analogy.</footnote>
The equations of motion of averaged, short-periodic, and second-order terms are deduced by first substituting Eq. (8) into the Hamilton equations associated to $\mathcal{H}(q, p, \varphi, \lambda, m, \gamma, \varepsilon, \beta)$, and then by expanding these equations up to the second-order in $\varepsilon$. This yields:

$$ \mathcal{H} = \varepsilon \left[ (\Lambda + \nu_{\Lambda}) \omega(\mathbf{Q}) + \frac{K(\mathbf{Q}, \mathbf{P}, \hat{\varphi}, 0)}{m} - \beta\Gamma \right] \\ + \varepsilon^2 \left[ \delta\Lambda \omega(\mathbf{Q}) + (\Lambda + \nu_{\Lambda}) \frac{\partial \omega}{\partial \mathbf{Q}} \cdot (\delta\mathbf{Q} + \nu_{\mathbf{Q}}) - \beta (\delta\Gamma + \nu_{\Gamma}) \right. \\ \left. + \frac{1}{m} \left( \frac{\partial K}{\partial \mathbf{Q}} \cdot (\delta\mathbf{Q} + \nu_{\mathbf{Q}}) + \frac{\partial K}{\partial \mathbf{P}} \cdot (\delta\mathbf{P} + \nu_{\mathbf{P}}) + \frac{\partial K}{\partial \lambda} (\Lambda + \nu_{\Lambda}) + \frac{\partial K}{\partial \hat{\varphi}} \nu_{\varphi} \right) \right]. \quad (11) $$

Time derivatives of secular and second-order terms are finally achieved by averiging the resulting equations and collecting terms of the appropriate order (hence, they are independent of the fast variable by construction), whereas the formal solution of short-periodic variations is obtained by compensating first-order oscillatory dynamics. The outcome of this process is detailed in the reminder of the section.

To simplify the notation, omitted arguments of functions and their derivatives are understood as the averaged value of slow variables, namely $\mathbf{Q}, \mathbf{P}, \varepsilon\Lambda = 0$, and fast variable $\hat{\varphi}$. For example $K = K(\mathbf{Q}, \mathbf{P}, \hat{\varphi}, 0)$ and $\frac{\partial K}{\partial \mathbf{Q}} = \frac{\partial K}{\partial \mathbf{q}}|_{\mathbf{q}=\mathbf{Q}, \mathbf{p}=\mathbf{P}, \hat{\varphi}=\hat{\varphi}, \lambda=0}$.

## A. First-order terms and the averaged-shooting problem

Targeting $\varepsilon$-independent equations of motion, we introduce the re-scaling of time variable

$$ m \, d\tau = \varepsilon \, dt. \quad (12) $$

Integration of Eq. (12),

$$ \int_{0}^{\tau_f} m \, d\tau = \varepsilon (t_f - t_0), \quad (13) $$

suggests that the original minimum-$\varepsilon$ problem with fixed transfer time can be replaced by an $\varepsilon$-independent, averaged problem with Lagrangian cost $m$ and Hamiltonian

$$ \bar{\mathcal{H}} = m [\omega(\mathbf{Q})\Lambda - \beta\Gamma - p_0] + \bar{K}(\mathbf{Q}, \mathbf{P}, 0), \quad (14) $$
where $p_0 \ge 0$ denotes the adjoint to the transformed Lagrangian-cost function. The variable $\Phi$ is cyclic and its momentum $\Lambda$ is thus a constant of motion. Hence, the time derivatives of averaged states and adjoints are

$$ 
\begin{aligned}
\frac{d\mathbf{Q}}{d\tau} &= \frac{\partial \bar{\mathbf{K}}}{\partial \mathbf{P}}, & \frac{d\Phi}{d\tau} &= m\omega(\mathbf{Q}), & \frac{dm}{d\tau} &= -\beta m, \\
\frac{d\mathbf{P}}{d\tau} &= -m\frac{\partial \omega}{\partial \mathbf{Q}}\Lambda - \frac{\partial \bar{\mathbf{K}}}{\partial \mathbf{Q}}, & \frac{d\Lambda}{d\tau} &= 0, & \frac{d\Gamma}{d\tau} &= -\omega(\mathbf{Q})\Lambda + \beta\Gamma + p_0.
\end{aligned}
\quad (15) $$

Because $\Gamma$ has no impact on the motion of the other variables, both its initial value and $p_0$ can be excluded from the shooting variables, so that the averaged counterpart of Eq. (7) is

$$ \bar{S}(\mathbf{P}(0), \tau_f | \mathbf{q}_0, \mathbf{q}_f, \beta, \Lambda) = \begin{bmatrix} \mathbf{Q}^*(\tau_f) - \mathbf{q}_f \\ \|\mathbf{P}(0)\| - 1 \end{bmatrix}, \quad (16) $$

where, similarly to Eq (7), $\mathbf{Q}^*(\tau)$ denotes the trajectory emanated by $\bar{\mathcal{H}}$ at time $\tau$ with initial condition $\mathbf{q}_0$. We stress again that $\bar{S}$ is independent of both $\varepsilon$ and $\Phi$, and it is a minimum-time problem. Conversely, the maneuvering time is fixed in both Problem (2) and its shooting function, Eq. (7). This dualism stems from the relation between minimum-time flow and the boundary of the reachable set [20, Chap 10, Theorem 1]. The notation $\bar{S}(\mathbf{P}(0), \tau_f | \mathbf{q}_0, \mathbf{q}_f, \beta, \Lambda)$ is introduced to emphasize that $\mathbf{P}(0)$ and $\tau_f$ are shooting variables, whereas $\mathbf{q}_0, \mathbf{q}_f, \beta$, and $\Lambda$ are parameters. The original and averaged shooting functions (Eqs (7) and (16), respectively) have the same boundary conditions of slow variables, namely $\mathbf{q}_0$ and $\mathbf{q}_f$. Corrections due to short-periodic variations and second-order terms are accounted for by computing sensitivities of $\bar{S}$, as detailed in the Sections III.B and III.C, respectively.

## B. Short-periodic variations

The time derivative of short-periodic variations is

$$ \frac{d\varepsilon\mathbf{v}(\cdot)}{d\tau} = \frac{d\Phi}{d\tau} \frac{\partial \mathbf{v}(\cdot)}{\partial \hat{\varphi}} + \varepsilon \left[ \frac{d\delta\Phi}{d\tau} \frac{\partial \mathbf{v}(\cdot)}{\partial \hat{\varphi}} + \frac{d\mathbf{Q}}{d\tau} \cdot \frac{\partial \mathbf{v}(\cdot)}{\partial \mathbf{Q}} + \frac{d\mathbf{P}}{d\tau} \cdot \frac{\partial \mathbf{v}(\cdot)}{\partial \mathbf{P}} \right] = m\omega \frac{\partial \mathbf{v}(\cdot)}{\partial \hat{\varphi}} + O(\varepsilon). \quad (17) $$

The average of Eq. (17) is zero because the derivatives of the averaged and second-order variables do not depend on the fast variable, and $\overline{\mathbf{v}(\cdot)} = 0$ by construction.
Short-periodic variations are obtained by solving

$$
\begin{align}
m \omega \frac{\partial \mathbf{v}_Q}{\partial \hat{\varphi}} &= \frac{\partial K}{\partial \mathbf{P}} - \frac{\partial \bar{K}}{\partial \mathbf{P}}, \nonumber \\
m \omega \frac{\partial \mathbf{v}_P}{\partial \hat{\varphi}} &= -\frac{\partial K}{\partial \mathbf{Q}} + \frac{\partial \bar{K}}{\partial \mathbf{Q}} - m \frac{\partial \omega}{\partial \mathbf{Q}} v_{\Lambda}, \nonumber \\
m \omega \frac{\partial v_{\Phi}}{\partial \hat{\varphi}} &= \frac{\partial K}{\partial \lambda} - \frac{\partial \bar{K}}{\partial \lambda} + m \frac{\partial \omega}{\partial \mathbf{Q}} \cdot \mathbf{v}_Q, \tag{18} \\
m \omega \frac{\partial v_{\Lambda}}{\partial \hat{\varphi}} &= -\frac{\partial K}{\partial \hat{\varphi}}, \nonumber \\
m \omega \frac{\partial v_{\Gamma}}{\partial \hat{\varphi}} &= -m \omega v_{\Lambda}; \nonumber
\end{align}
$$

with the conditions to resolve the integration constants:

$$
\overline{v_Q} = 0, \quad \overline{v_P} = 0, \quad \overline{v_\Phi} = 0, \quad \overline{v_\Lambda} = 0, \quad \overline{v_\Gamma} = 0. \tag{19}
$$

Denote by $f^{(k)}$ the $k$-th coefficient of the Fourier series of a smooth periodic function $f(\varphi)$, namely $f(\varphi) = \sum_{k \in \mathbb{Z}} f^{(k)} e^{ik\varphi}$, and by $i$ the imaginary unit. The formal solution of Eqs. (18) and (19) is

$$
\begin{align}
\mathbf{v}_Q &= -\frac{i}{m\omega} \sum_{k \in \mathbb{Z} \setminus \{0\}} \frac{1}{k} \frac{\partial K^{(k)}}{\partial \mathbf{P}} e^{ik\hat{\varphi}}, \\
\mathbf{v}_P &= \frac{i}{m\omega} \sum_{k \in \mathbb{Z} \setminus \{0\}} \frac{1}{k} \left( \frac{\partial K^{(k)}}{\partial \mathbf{Q}} - \frac{K^{(k)}}{\omega} \frac{\partial \omega}{\partial \mathbf{Q}} \right) e^{ik\hat{\varphi}}, \\
\mathbf{v}_{\Phi} &= -\frac{i}{m\omega} \sum_{k \in \mathbb{Z} \setminus \{0\}} \frac{1}{k} \left( \frac{\partial K^{(k)}}{\partial \lambda} - \frac{i}{k\omega} \frac{\partial \omega}{\partial \mathbf{Q}} \cdot \frac{\partial K^{(k)}}{\partial \mathbf{P}} \right) e^{ik\hat{\varphi}}, \tag{20} \\
\mathbf{v}_{\Lambda} &= -\frac{K(\mathbf{Q}, \mathbf{P}, \hat{\varphi}, 0) - \bar{K}(\mathbf{Q}, \mathbf{P}, 0)}{m\omega(\mathbf{Q})}, \\
\mathbf{v}_{\Gamma} &= -\frac{i}{m\omega} \sum_{k \in \mathbb{Z} \setminus \{0\}} \frac{K^{(k)}}{k} e^{ik\hat{\varphi}},
\end{align}
$$

which can be verified by substituting it into Eq. (18). We note that $v_\Lambda$ establishes a first-order equivalence between the original and averaged Hamiltonians. Short-periodic variations can be efficiently evaluated by means of the fast-Fourier transform (FFT) algorithm [21].
### C. Second-order terms

The equations of motion of $\delta Q$, $\delta P$, and $\delta \Phi$ are obtained by expanding Eq. (6) up to the second order in $\epsilon$, averaging, and discarding first-order terms (which where already used in the previous sub-sections). This yields‡

$$ \begin{aligned} \frac{d\delta Q}{d\tau} &= \overline{\frac{\partial^2 K}{\partial Q \partial P} (\delta Q + \nu_Q)} + \overline{\frac{\partial^2 K}{\partial P^2} (\delta P + \nu_P)} + \overline{\frac{\partial^2 K}{\partial \varphi \partial P} \nu_\Phi} + \overline{\frac{\partial^2 K}{\partial \lambda \partial P} (\Lambda + \nu_\Lambda)}, \\ \frac{d\delta P}{d\tau} &= -\overline{\frac{\partial^2 K}{\partial Q^2} (\delta Q + \nu_Q)} - \overline{\frac{\partial^2 K}{\partial P \partial Q} (\delta P + \nu_P)} - \overline{\frac{\partial^2 K}{\partial \varphi \partial Q} \nu_\Phi} - \overline{\frac{\partial^2 K}{\partial \lambda \partial Q} (\Lambda + \nu_\Lambda)} \\ &\quad - m \overline{\frac{\partial^2 \omega}{\partial Q} (\delta Q + \nu_Q) (\Lambda + \nu_\Lambda)} - m \overline{\frac{\partial \omega}{\partial Q} \delta \Lambda}, \end{aligned} \quad (21) $$

$$ \frac{d\delta\Phi}{d\tau} = m \frac{\partial \omega}{\partial Q} \cdot \delta Q + \frac{\partial \bar{K}}{\partial \lambda}, $$

$$ \frac{d\delta\Gamma}{d\tau} = -\frac{\partial \omega}{\partial Q} \cdot \overline{(\delta Q + \nu_Q)(\Lambda + \nu_\Lambda)} + \beta \delta\Gamma. $$

Higher-order terms of Eq. (17) do not contribute to Eq. (21) because their average is zero, as previously mentioned. If the determination of second-order, short-periodic variations were envisaged, they could not be neglected.

The expansion of the Hamiltonian up to the second order in $\epsilon$ enables an explicit expression of $\delta\Lambda$ as a function of the other variables, namely

$$ \delta\Lambda = -\frac{1}{m\omega} \left[ \left( \frac{\partial K}{\partial Q} + m(\Lambda + \nu_\Lambda) \frac{\partial \omega}{\partial Q} \right) \cdot (\delta Q + \nu_Q) \right. \\ \left. + \frac{\partial K}{\partial P} \cdot (\delta P + \nu_P) + \frac{\partial K}{\partial \Phi} \nu_\Phi + \frac{\partial K}{\partial \lambda} (\Lambda + \nu_\Lambda) - m\beta (\delta\Gamma + \nu_\Gamma) \right]. \quad (22) $$

Denoting

$$ \begin{aligned} \nu_Q^{(k)} &:= -\frac{i}{m \omega k} \frac{\partial K^{(k)}}{\partial P}, \\ \nu_P^{(k)} &:= \frac{i}{m \omega k} \left( \frac{\partial K^{(k)}}{\partial Q} - \frac{K^{(k)}}{\omega} \frac{\partial \omega}{\partial Q} \right), \end{aligned} \quad (23) $$

$$ \nu'_{\Phi}^{(k)} := \frac{1}{m\omega} \left( m \frac{\partial \omega}{\partial Q} \cdot \nu_Q^{(k)} + \frac{\partial K^{(k)}}{\partial \lambda} \right), $$

the k-th coefficients of the series in Eq. (20), the averaged value of $\delta\Lambda$ (which appears in the time derivative of $\delta P$) is given by

$$ \begin{aligned} \overline{\delta\Lambda} = &-\frac{1}{m\omega} \left[ \left( \frac{\partial K^{(0)}}{\partial Q} + m\Lambda \frac{\partial \omega}{\partial Q} \right) \cdot \delta Q + \frac{\partial K^{(0)}}{\partial P} \cdot \delta P + \frac{\partial K^{(0)}}{\partial \lambda} \Lambda - m\beta \delta\Gamma \right] \\ &-\frac{1}{m\omega} \sum_{\mathcal{Z}\setminus\{0\}} \left[ \left( \frac{\partial K^{(-k)}}{\partial Q} + m\nu_{\Lambda}^{(-k)} \frac{\partial \omega}{\partial Q} \right) \cdot \nu_Q^{(k)} + \frac{\partial K^{(-k)}}{\partial P} \cdot \nu_P^{(k)} - \left( \nu'_{\Phi}^{(-k)} + \frac{1}{m\omega} \frac{\partial K^{(-k)}}{\partial \lambda} \right) K^{(k)} \right], \end{aligned} \quad (24) $$

<footnote>‡The convention $\left[ \frac{\partial^2 K}{\partial Q \partial P} \right]_{[l,m]} = \frac{\partial^2 K}{\partial X_m \partial P_l}$ is used for second derivatives. Indeed, $\frac{\partial^2 K}{\partial Q \partial P} = \left( \frac{\partial^2 K}{\partial P \partial Q} \right)^T$.</footnote>
where we extensively used the integration rule of the average product of two periodic functions $f(\varphi)$ and $g(\varphi)$ as

$$ \overline{f(\varphi) g(\varphi)} = \sum_{k \in \mathbb{Z}} f^{(k)} g^{(-k)}. \qquad (25) $$

Substituting Eq.(20) into (21), integrating by parts terms that include derivatives of $\varphi$, and again using the rule of Eq. (25), yields

$$ \begin{aligned} \frac{d}{d\tau} \delta \mathbf{Q} = & \frac{\partial^2 K^{(0)}}{\partial \mathbf{Q} \partial \mathbf{P}} \delta \mathbf{Q} + \frac{\partial^2 K^{(0)}}{\partial \mathbf{P}^2} \delta \mathbf{P} + \frac{\partial^2 K^{(0)}}{\partial \lambda \partial \mathbf{P}} \Lambda \\ & + \sum_{k \in \mathbb{Z} \setminus \{0\}} \left( \frac{\partial^2 K^{(-k)}}{\partial \mathbf{Q} \partial \mathbf{P}} \mathbf{v}_{\mathbf{Q}}^{(k)} + \frac{\partial^2 K^{(-k)}}{\partial \mathbf{P}^2} \mathbf{v}_{\mathbf{P}}^{(k)} - \frac{\partial K^{(-k)}}{\partial \mathbf{P}} \mathbf{v}_{\Phi}^{(k)} + \frac{\partial^2 K^{(-k)}}{\partial \lambda \partial \mathbf{P}} \mathbf{v}_{\Lambda}^{(k)} \right), \\ \frac{d}{d\tau} \delta \mathbf{P} = & - \left( \frac{\partial^2 K^{(0)}}{\partial \mathbf{Q}^2} + m \Lambda \frac{\partial^2 \omega}{\partial \mathbf{Q}^2} \right) \delta \mathbf{Q} - \frac{\partial^2 K^{(0)}}{\partial \mathbf{P} \partial \mathbf{Q}} \delta \mathbf{P} - \frac{\partial^2 K^{(0)}}{\partial \lambda \partial \mathbf{Q}} \Lambda - m \frac{\partial \omega}{\partial \mathbf{Q}} \delta \Lambda \\ & - \sum_{k \in \mathbb{Z} \setminus \{0\}} \left[ \left( \frac{\partial^2 K^{(-k)}}{\partial \mathbf{Q}^2} + m \mathbf{v}_{\Lambda}^{(-k)} \frac{\partial^2 \omega}{\partial \mathbf{Q}^2} \right) \mathbf{v}_{\mathbf{Q}}^{(k)} + \frac{\partial^2 K^{(-k)}}{\partial \mathbf{P} \partial \mathbf{Q}} \mathbf{v}_{\mathbf{P}}^{(k)} - \frac{\partial K^{(-k)}}{\partial \mathbf{Q}} \mathbf{v}_{\Phi}^{(k)} + \frac{\partial^2 K^{(-k)}}{\partial \lambda \partial \mathbf{Q}} \mathbf{v}_{\Lambda}^{(k)} \right], \end{aligned} \qquad (26) $$

$$ \frac{d}{d\tau} \delta\Phi = m \frac{\partial \omega}{\partial \mathbf{Q}} \cdot \delta\mathbf{Q} + \frac{\partial K^{(0)}}{\partial \lambda}, $$

$$ \frac{d}{d\tau} \delta\Gamma = - \frac{\partial \omega}{\partial \mathbf{Q}} \cdot \left( \delta\mathbf{Q} \Lambda + \sum_{k \in \mathbb{Z} \setminus \{0\}} \mathbf{v}_{\mathbf{Q}}^{(-k)} \mathbf{v}_{\Lambda}^{(k)} \right) + \beta \delta\Gamma. $$

Equation (26) is a linear system in $\delta \mathbf{Q}$, $\delta \mathbf{P}$, $\delta \Phi$, and $\delta \Gamma$. Although these expressions appear cumbersome, they are obtained by elementary combinations of Fourier coefficients of $K$ and their derivatives. As such, this formulation lends itself to a practical implementation by leveraging on automatic differentiation (AD) of the Hamiltonian [13] and FFT.

## IV. Algorithm for the generation of Bacon plots

We propose a four-step algorithm for the generation of approximated solutions of Problem (2). A valuable feature of the algorithm is that it relies on a one-parameter family of solutions of the averaged TPBVP that is independent of $\varepsilon$, $\beta$, $t_0$, $t_f$, and boundary longitudes. Almost the entire computational cost of the algorithm is on Steps 1 and 2, which are independent of these parameters. Accurate solutions for arbitrary combinations of these variables can be efficiently evaluated (Steps 3 and 4) with almost no additional computational burden.

The algorithm is articulated as follows:

1) Compute a family of solutions of the averaged problem parametrized by $\Lambda$, namely

$$ \begin{gathered} \forall \Lambda \in [\Lambda^{\min}, \Lambda^{\max}] \text{ find } \mathbf{P}_0^{\Lambda}, \tau_f^{\Lambda} \text{ s.t.:} \\ \bar{S}(\mathbf{P}_0^{\Lambda}, \tau_f^{\Lambda} | \mathbf{q}_0, \mathbf{q}_f, 0, \Lambda) = 0. \end{gathered} \qquad (27) $$

After finding the solution for $\Lambda = 0$ (this is the only shooting problem that has to be solved without a priori
knowledge), a numerical continuation scheme is used to obtain solutions in the range $[\Lambda^{min}, \Lambda^{max}]$ ($\Lambda$ serves as the continuation parameter). Although fold bifurcations may occur and multiple solutions may be found for the same value of $\Lambda$, we prefer to formally use $\Lambda$ as continuation parameter, instead of introducing a pseudo-arc length variable, to streamline the notation.

We note that the same initial conditions of the osculating problem, $q_0$ and $q_f$, are used here, i.e., short-periodic variations are neglected at this step and, as such, these solutions are independent of the boundary phases, $\varphi_0(t_0)$ and $\varphi_f(t_f)$. In addition, $\beta$ is set to 0 at this stage, i.e., mass is constant.

2) Evaluate sensitivities of the averaged-shooting function with respect to both its variables and parameters and second-order variations $\delta Q_2^\Lambda$ and $\delta \Phi_2^\Lambda$, which are obtained by integrating Eq. (21) with initial conditions $\delta Q = \delta P = 0$ and $\delta \Phi = \delta \Gamma = 0$. This yields

$$ \begin{array}{l} \forall \Lambda \in [\Lambda^{min}, \Lambda^{max}] \text{ evaluate} \\ \frac{\partial \bar{S}}{\partial \tau_f}, \frac{\partial \bar{S}}{\partial \mathbf{P}_0}, \frac{\partial \bar{S}}{\partial q_0}, \frac{\partial \bar{S}}{\partial q_f}, \frac{\partial \bar{S}}{\partial \beta}, \frac{\partial \Phi_f}{\partial \tau_f}, \frac{\partial \Phi_f}{\partial \mathbf{P}_0}, \frac{\partial \Phi_f}{\partial q_0}, \frac{\partial \Phi_f}{\partial \beta}, \\ \delta Q_2^\Lambda = \delta Q (\tau_f^\Lambda | q_0, \mathbf{P}_0^\Lambda, 0, \Lambda), \\ \delta \Phi_2^\Lambda = \delta \Phi (\tau_f^\Lambda | q_0, \mathbf{P}_0^\Lambda, 0, \Lambda), \end{array} \qquad (28) $$

where $\Phi_f = \Phi (\tau_f^\Lambda | q_0, \mathbf{P}_0^\Lambda, 0, \Lambda)$ is the final re-scaled fast variable of the averaged system. When AD is used to achieve exact derivatives of Eq. (15), this step requires the integration of the linearized system for each solution $(\mathbf{P}^\Lambda, \tau_f^\Lambda)$ and of Eq. (26).

3) Compute variations of $(\mathbf{P}^\Lambda, \tau_f^\Lambda)$ due to the desired value of $\beta$ and of the boundary phases, $\varphi_0(t_0)$ and $\varphi_f(t_f)$. Because low-thrust propulsion is characterized by $\beta \ll 1$, solutions for arbitrary (small) $\beta$ are achieved by computing sensitivities of the shooting function to this parameter$^8$. We also note that the mass-dependent term of the time derivative of $\mathbf{P}$ in Eq. (15) vanishes when $\Lambda = 0$, which implies that solutions of the averaged problem are independent of $\beta$ in this case, and suggests a moderate sensitivity to the $\beta$ parameter for small $\Lambda$. Denote by $\delta Q_0$ and $\delta Q_f$ the perturbations of the boundary conditions due to both short-periodic variations and second-order terms, i.e.,

$$ \begin{array}{l} \delta Q_0 = -\nu_Q (q_0, \mathbf{P}_0^\Lambda, \varphi_0(t_0)), \\ \delta Q_f = -\nu_Q (q_f, \mathbf{P} (\tau_f^\Lambda), \varphi_f(t_f)) - \delta Q_2^\Lambda, \end{array} \qquad (29) $$

<footnote>$^8$If solutions for large $\beta$ are envisaged, the computation of the one-parameter family $(\mathbf{P}_0^\Lambda, \tau_f^\Lambda)$ at Step 1 can be carried out by using the desired value of $\beta$. The starting point of the continuation, $(\mathbf{P}_0^0, \tau_f^0)$, remains unchanged because trajectories of $\mathbf{Q}$ and $\mathbf{P}$ are mass-independent when $\Lambda = 0$.</footnote>
Hence, this step consists in solving:

$$
\begin{align}
\text{Given } \beta, \varphi_0(t_0), \varphi_f(t_f) \text{ find } \delta\tau_f^\Lambda, \delta\mathbf{P}_0^\Lambda, \delta\Phi_f^\Lambda \text{ s.t.:} \nonumber \\
\frac{\partial \bar{S}}{\partial \tau_f} \delta\tau_f^\Lambda + \frac{\partial \bar{S}}{\partial \mathbf{P}_0} \delta\mathbf{P}_0^\Lambda &= -\frac{\partial \bar{S}}{\partial \mathbf{q}_0} \delta\mathbf{Q}_0 - \frac{\partial \bar{S}}{\partial \mathbf{q}_f} \delta\mathbf{Q}_f - \frac{\partial \bar{S}}{\partial \beta} \beta, \tag{30} \\
\delta\Phi_f^\Lambda &= \frac{\partial \Phi_f}{\partial \tau_f} \delta\tau_f + \frac{\partial \Phi_f}{\partial \mathbf{P}_0} \delta\mathbf{P}_0 + \frac{\partial \Phi_f}{\partial \mathbf{q}_0} \delta\mathbf{Q}_0 + \frac{\partial \Phi_f}{\partial \beta} \beta + \delta\Phi_2^\Lambda, \nonumber
\end{align}
$$

whose computational complexity is as modest as the solution of a six-dimensional linear system.

4) Identify minimum ε and the good member of the family of averaged solutions, namely

$$
\text{Find } \Lambda \text{ s.t.:} \qquad \varphi_0(t_0) + \frac{\Phi_f^\Lambda}{\varepsilon^\Lambda} + \delta\Phi_f^\Lambda - \varphi_f(t_f) = 0 \quad \mod 2\pi, \qquad (31)
$$

where $\varepsilon^\Lambda$ is chosen such that the final masses of the averaged and osculating trajectories, i.e.,

$$
m(\tau_f^\Lambda + \varepsilon^\Lambda \delta\tau_f^\Lambda) = e^{-\beta(\tau_f^\Lambda + \varepsilon^\Lambda \delta\tau_f^\Lambda)} \approx (1 - \varepsilon^\Lambda \beta \delta\tau_f^\Lambda) e^{-\beta \tau_f} \quad \text{and} \quad m(t_f - t_0) = 1 - \varepsilon^\Lambda \beta (t_f - t_0), \quad (32)
$$

respectively, equal at first order, yielding

$$
\varepsilon^{\Lambda} = \begin{cases} \frac{\tau_f^{\Lambda}}{t_f - t_0 - \delta\tau_f^{\Lambda}} & \text{if } \beta = 0, \\ \frac{1}{\beta} \frac{1 - e^{-\beta\tau_f^{\Lambda}}}{t_f - t_0 - e^{-\beta\tau_f^{\Lambda}}\delta\tau_f^{\Lambda}} & \text{otherwise.} \end{cases} \tag{33}
$$

Equation (31) can be interpreted as a matching between the angular length estimated by the averaged system and the desired one given by the position on the orbit of the departure and arrival points. Second-order variations are necessary to have an approximation of the longitude that converges as $\varepsilon$ approaches zero.

The software Hampath<sup>¶</sup> described in [13] was used in our numerical implementation to tackle Steps 1 and 2. Specifically, Hampath leverages on AD to compute derivatives of both $\overline{\mathcal{H}}$ and $\overline{S}$. Hampath only requires for input the maximized Hamiltonian and the shooting function, and it also offers a continuation algorithm based on differential homotopy that can be used to compute the one-parameter family of solutions $(\mathbf{P}_0^\Lambda, \tau_f^\Lambda)$.

We stress that all steps only use the averaged system and the near-identity transformation. Departure and arrival dates only appear in Steps 3 and 4, so that outcomes of Steps 1 and 2, which involve the only computationally demanding tasks, need to be evaluated only once for any combination of $\beta, t_0$, and $t_f$.

<footnote>¶Website: https://www.hampath.org/. New developments of Hampath are available in the control-toolbox library: https://ct.gitlabpages.inria.fr/gallery/</footnote>
<table><caption>**Table 1** GTO-GEO transfer. Elements of departure and arrival orbits.</caption><thead><tr><td>Variable</td><td>GTO</td><td>GEO</td></tr></thead><tbody><tr><td>Semi-major axis [km]</td><td>12000</td><td>42165</td></tr><tr><td>Eccentricity [-]</td><td>0.75</td><td>0</td></tr><tr><td>Inclination [deg]</td><td>4</td><td>0</td></tr><tr><td>Argument of perigee [deg]</td><td>0</td><td>0</td></tr><tr><td>Right ascension of ascending node [deg]</td><td>30</td><td>0</td></tr></tbody></table>

## V. Examples

Two examples are illustrated. The first one consists of a GTO to GEO transfer, and it is characterized by a large number of revolutions. The second one considers an Earth-Venus transfer, which is performed in just a few orbits.

### A. GTO to GEO transfer

Table 1 lists simulation parameters of this case study. The one-parameter family of averaged solutions computed at Step 1 is illustrated in Fig. 1. Because this transfer occurs in several hundreds of orbits, the continuation is stopped for a small value of $Λ$, since several local solutions of the original problem are expected to be found in proximity of $Λ = 0$ when the final angular length (approximated by $Φ/ε$) is large in Eq.(31) (note the presence of the modulo $2π$ in this equation). In this short range of $Λ$, $τ_f^Λ$ is convex (it has a minimum for $Λ = 0$ as it could be regarded as the $Φ$-free, minimum $τ_f$ solution) and the other shooting variables are nearly linear. The conjecture on the presence of various local solutions in proximity of $Λ = 0$ is verified in Fig. 2, which depicts the solutions identified at the Step 4 for a given departure and arrival dates. Here, the solid line is the draw of Eq. (33), and black dots are solutions of Eq. (31) If no solution were found, the continuation had to be pushed further. From this point of view, the GTO-to-GEO transfer is a "simple" case study because the importance of the position on departure and arrival orbits dwindles when several revolutions are necessary to complete the maneuver. Hence, simple averaging offers a fair approximation without a real need for the inclusion of short-periodic and second-order terms.

The Bacon plots of minimum thrust and final mass ratio are illustrated in Fig. 3 and 4, respectively. Solutions along red-dashed lines have the same maneuvering time, which is larger for upper lines. Again, results are rather intuitive: level sets of the minimum thrust are nearly aligned with iso-transfer-time lines (in red) but short-periodic variations can be appreciated. The delivered mass is almost constant because $τ_f^Λ$ exhibits minor changes and $δτ_f^Λ$ is negligible with respect to $t_f - t_0$. The change of time variable carried out in Eq. (12) facilitates the integration of the trajectory when the mass approaches zero by converting the increase of the thrust-to-mass ratio into a stretching of the time variable. This numerical care also ensures that mass remains strictly positive throughout the maneuver, since $m(τ)$ is an exponential decay (whereas $m(t)$ is a linear function with negative slope).
<figure><img src="image_2.png" /><figcaption>Fig. 1 Shooting variables outcome of Step 1 for the GEO transfer.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Fig. 2 Identified solutions for the bottom-left point of the launch window in the GTO-to-GEO transfer. Red dot denotes the global optimum.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Fig. 3 Minimum thrust magnitude of the GEO transfer for various departure and arrival dates.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4 Final mass for the transfer to GEO with $\beta = 10^{-4}$ s/m.</figcaption></figure>
<table><caption>Table 2 Earth-Venus transfer. Elements of departure and arrival orbits.</caption><thead><tr><th>Variable</th><th>Earth</th><th>Venus</th></tr></thead><tbody><tr><td>Semi-major axis [AU]</td><td>1</td><td>0.723336</td></tr><tr><td>Eccentricity [-]</td><td>1.671 &middot; 10<sup>-2</sup></td><td>6.78 &middot; 10<sup>-3</sup></td></tr><tr><td>Ecliptic inclination [deg]</td><td>0</td><td>3.3947</td></tr><tr><td>Argument of perihelion [deg]</td><td>102.9</td><td>54.9</td></tr><tr><td>Longitude of ascending node [deg]</td><td>0</td><td>76.7</td></tr><tr><td>Time of periapsis [UT]</td><td>Jan 15 2007, 00:30</td><td>Oct 6 2006, 5:19</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>Fig. 5 Shooting variables outcome of Step 1 for the Earth-Venus transfer.</figcaption></figure>

## B. Earth-Venus transfer

Table 2 lists simulation parameters of this case study. Here, the transfer is achieved in few revolutions and, as such, it is necessary to compute the family of averaged solutions at Step 1 over a large range of $\Lambda$. This family is illustrated on Fig. 5. In this case, a fold bifurcation occurs when $\Lambda \approx 3$. Shooting variables exhibit large-nonlinear variations in the family, including a change of sign of $p_{I,1}^{\Lambda}$.

Figures 6a and 6b illustrate Bacon plots of $\varepsilon^*$ for $\beta = 0$ s/m and $\beta = 10^{-4}$ s/m, respectively ($\beta = 10^{-4}$ s/m corresponds to a specific impulse of about 1000 s, which is a value compatible with existing ion-thrust technology). Minimum thrust magnitude lowers when $\beta$ increases. This is expected because mass consumption facilitates the transfer by increasing the thrust-to-mass ratio at the end of the maneuver. When $\beta = 0$ s/m, the specific impulse of the thrusters approaches infinity, and mass consumption is neglected. The $\beta = 0$ s/m case is of interest as it yields a conservative (i.e., $\varepsilon$ is overestimated for any $\beta > 0$), but feasible, solution of Problem (2).
<figure><img src="image_7.png" /><figcaption>Fig. 6 Minimum thrust magnitude of the Earth-Venus transfer for various departure and arrival dates.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Fig. 7 Final mass for the Earth-Venus transfer with $\beta = 10^{-4}$ s/m.</figcaption></figure>

Fig. 7 depicts the final mass for the transfer with $\beta = 10^{-4}$ s/m. This value could be increased, for example, by solving the minimum fuel problem with values of $\varepsilon$ larger than the ones provided in Fig. 6b. Differently of the GTO-to-GEO example, variations of the final mass are appreciable, since $\tau_f^\Lambda$ spans over a large interval.

Bacon plots of $\varepsilon^*$ are not smooth, as both jump and derivative discontinuities are possible. The first ones arise when new local solutions are generated, but this does not happen in the case study at hand. Instead, discontinuities of the derivative, which appear in Fig. 6, occur when two existing local solutions have the same $\varepsilon$. This can be better appreciated by inspecting Fig. 8, which depicts the value of $\Lambda$ for the transfers with $\beta = 10^{-4}$ s/m. In fact, each point of the Bacon plot is associated to a precise member of the one-parameter family of averaged solutions. Discontinuities of $\Lambda$ (occurring at the interface of yellow and blue regions of this figure) are associated to coexistence of two local solutions with the same cost function $\varepsilon$. This mechanism is further emphasized in Fig. 9, which shows the $\varepsilon$ of local solutions along the red-dashed line of Fig. 8. Two distinct solutions are identified for some values of the departure date, but only the one with lower $\varepsilon$ is of interest. Red dot of Fig. 9 indicates a cut point, i.e., where the two local solution have the same $\varepsilon$, and one branch stops tracking the global optimum. At this point, a discontinuity occur in the charts of Fig. 6 and 8 (derivative of $\varepsilon$, and jump of $\Lambda$ and of all other shooting variables). In conclusion, discontinuities stem from crossing of the cut locus, whose immediate identification is a major asset of the proposed methodology. For the sake of completeness, we note that the achievement of global optima and identification of cut points is guaranteed under the assumption of uniqueness of the one-parameter family of averaged solutions. this is a fair assumption, considering that the role of $\Lambda$ was systematically ignored in the literature on multi-revolution transfers, and only solutions with $\Lambda = 0$ were studied.
<figure><img src="image_9.png" /><figcaption>Fig. 8 Continuation parameter, $\Lambda$, of the best solutions identified in Step 4 for the Earth-Venus transfer with $\beta = 10^{-4}$ s/m.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 9 Value of $\varepsilon$ for local solutions along the red-dashed line of Fig. 8.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Fig. 10 Identified solutions for departure and arrival dates at the point emphasized in Fig. 8. in the Earth-Venus transfer. Red dot denotes the global optimum.</figcaption></figure>

Figure 10 depicts local solutions found at the departure and arrival time located at the red dot of Fig. (8), which is close to a cut point. Differently of Fig. 2, the best-two solutions have indeed similar $ε$, but they are very distant in the family of averaged solutions, which suggests that they likely represent very different trajectories.

The relevance of second-order terms is disclosed in Fig. 11. Contributions of averaged, second-order, and short-periodic variations are pointed out in Fig. 11a, where a sample trajectory of the adjoint $p_4$ is drawn. Drift between $p_4$ and $P_4$ is of order $ε$ at the end of the maneuver, and (for this particular trajectory) including $\delta P_4$ is more important than short-periodic variations $v_{P_4}$ to reduce the error with respect to the original trajectory. In general, the contributions of the second-order terms and short-periodic variations are of the same order of magnitude. The necessity to account for second-order terms stems from the need of a first-order estimation of the fast variable. As depicted in Fig. 11b, the error between $\varphi$ and $\Phi$ converges to $\delta\Phi$ (which is of order one) as $ε$ approaches zero.

Indeed, a drawback of the proposed framework is that the quality of the approximation is tied to the small magnitude of $ε$. In fact, our approximation of slow variables is of order $ε^2$ and the one of $\varphi$ is of order $ε$ (which is an improvement with respect to simple averaging). Hence, the methodology is not adequate to approximate short transfers, which are characterized by large values of $ε$. These estimates are not quantitatively clear. Our experience suggests that the proposed algorithm offers a reliable approximation when the number of orbital revolutions is about 10, an excellent approximation for transfers occurring in more than 100 orbits (but in this case simple averaging may be sufficient for the purposes of a preliminary analysis), and it is paramount for transfers with 10 to 100 orbits, since approximations
<figure><img src="image_12.png" /><figcaption>(a) Adjoint $p_4$. Here, $\Lambda = -1$, $\beta = 10^{-4}$ s/m, and $\varepsilon = 0.05$ mN/kg.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(b) Error between true fast variable and first-order approximation, $\frac{\Phi}{\varepsilon}$ with $\Lambda = -1$ and $\beta = 10^{-4}$ s/m. Convergence to the second-order contribution $\delta\Phi$ is achieved as expected.</figcaption></figure>

<caption>Fig. 11 Importance of second-order terms.</caption>
<figure><img src="image_14.png" /><figcaption>Fig. 12 Comparison between approximated and true solutions of Problem (2).</figcaption></figure>
given by simple averaging are modest. This rule of thumbs is supported by Fig. 12, which compares the solution of the
original problem obtained with the proposed approximation (solid-red line) and the one given by simple averaging
neglecting short-periodic and second-order terms (dashed-blue line). Black-dotted line is the true solution, and the
departure and arrival dates are the ones of the red dot in Fig. 8. Our approximation outperforms simple averaging in this
case. The extremely large position error obtained with simple averaging is mostly due to the poor approximation of the
fast variable.

Finally, we note that the overall runtime of the algorithm¹ is of the order of 1 hour. This value appears large at a first view. However, the largest part of the time is spent in the computation of Step 1 ( 50 minutes) and Step 2 ( 10 minutes), which are independent of the resolution and bounds of the launch windows. In fact, the Bacon plots that we presented were computed on a grid of 360 × 360 points. If each point of the plot were evaluated in a single second (this is the order of magnitude of the time required to compute a single trajectory in [6]), computational time would be of the order of 36 hours. Drawing a full plot with a completely different launch window can be achieved at almost no additional cost. For completeness, overall runtime of the GTO-to-GEO example is about 5 minutes.

VI. Conclusion

This paper offers an efficient, second-order numerical approximation to the minimum-thrust orbital transfer problem which may serve as a tool for preliminary mission design. A major feature of the proposed methodology is that drawing a full chart of transfers (namely, for any combination of mass consumption rate, departure and arrival dates) requires the solution of a single, one-parameter family of averaged TPBVP. The proposed algorithm delivers a reliable approximation for scenarios involving a few orbital revolutions and achieves an excellent level of accuracy for transfers spanning many orbits, where simple averaging might suffice for preliminary analyses. However, the algorithm's significance becomes paramount for transfers within the intermediate range of revolutions, as approximations given by simple averaging prove to be only modest in such cases. In addition, the methodology allows the immediate identification of cut points, which enables the straightforward discarding of non-minimizing local solutions. The accuracy of the approximation deteriorates when dealing with very short transfers characterized by large values of the thrust-to-mass ratio. For these transfer, shape-based approaches should offer a more accurate approximation. In addition, orbital perturbations are neglected in this study. Their inclusion in the methodology may be possible by carrying out a two-parameter continuation of the averaged solution, where the second parameter represent the ratio between the magnitude of the thrust and of the perturbation.

The formal solution of short-periodic variations and an implementable expression of time derivatives of second-order
terms are also provided. Such novel expressions can be efficiently evaluated by means of the FFT algorithm and an AD
tool, and they can be directly applied to the more general framework of minimum-time, multi-revolution transfers.

<footnote>We used a Dell Latitude 7390 with the Matlab version of Hampath.</footnote>
References

[1] Jordan, J. F., "The application of Lambert's theorem to the solution of interplanetary transfer problems," Tech. Rep. 32-521, Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA, USA, 1964.

[2] Izzo, D., "Lambert's problem for exponential sinusoids," *Journal of guidance, control, and dynamics*, Vol. 29, No. 5, 2006, pp. 1242-1245. https://doi.org/10.2514/1.21796.

[3] Zuiani, F., Vasile, M., Palmas, A., and Avanzini, G., "Direct transcription of low-thrust trajectories with finite trajectory elements," *Acta Astronautica*, Vol. 72, 2012, p. 108-120. https://doi.org/10.1016/j.actaastro.2011.09.011, URL http://dx.doi.org/10.1016/j.actaastro.2011.09.011.

[4] Avanzini, G., Palmas, A., and Vellutini, E., "Solution of low-thrust Lambert problem with perturbative expansions of equinoctial elements," *Journal of Guidance, Control, and Dynamics*, Vol. 38, No. 9, 2015, pp. 1585-1601. https://doi.org/10.2514/1.g001018, URL https://doi.org/10.2514%2F1.g001018.

[5] Novak, D. M., and Vasile, M., "Improved Shaping Approach to the Preliminary Design of Low-Thrust Trajectories," *Journal of Guidance, Control, and Dynamics*, Vol. 34, No. 1, 2011, p. 128-147. https://doi.org/10.2514/1.50434, URL http://dx.doi.org/10.2514/1.50434.

[6] Woolley, R. C., Laipert, F., Nicholas, A. K., and Olikara, Z., "Low-thrust trajectory Bacon plots for Mars mission design," *29th AAS/AIAA Space Flight Mechanics Meeting, Ka'anapali, Hawaii*, Vol. 21, American Astronautical Society (AAS), 2019, pp. 797-814. https://doi.org/hdl.handle.net/2014/50666.

[7] Gao, Y., and Kluever, C. A., "Low-thrust interplanetary orbit transfers using hybrid trajectory optimization method with multiple shooting," *AIAA/AAS Astrodynamics Specialist Conference and Exhibit*, 2004, p. 5088. https://doi.org/10.2514/6.2004-5088.

[8] Conway, B. (ed.), *Spacecraft Trajectory Optimization*, Cambridge University Press (CUP), 2010. https://doi.org/10.1017/cbo9780511778025, URL http://dx.doi.org/10.1017/cbo9780511778025.

[9] Tarzi, Z., Speyer, J., and Wirz, R., "Fuel optimum low-thrust elliptic transfer using numerical averaging," *Acta Astronautica*, Vol. 86, 2013, pp. 95-118. https://doi.org/10.1016/j.actaastro.2013.01.003.

[10] Dell'Elce, L., Caillau, J.-B., and Pomet, J.-B., "On the convergence of time-optimal maneuvers of fast-oscillating control systems," *2021 European Control Conference (ECC)*, IEEE, 2021, pp. 2008-2013. https://doi.org/10.23919/ECC54610.2021.9655101.

[11] Geffroy, S., and Epenoy, R., "Optimal low-thrust transfers with constraints-generalization of averaging techniques," *Acta Astronautica*, Vol. 41, No. 3, 1997, p. 133-149. https://doi.org/10.1016/s0094-5765(97)00208-7, URL http://dx.doi.org/10.1016/s0094-5765(97)00208-7.

[12] Wei, P., and Changhou, L., "Direct Optimization of Low-Thrust Transfers Using Averaging Techniques Based on Fourier Series Expansion," *2011 Fourth International Joint Conference on Computational Sciences and Optimization*, IEEE, 2011. https://doi.org/10.1109/cso.2011.118, URL http://dx.doi.org/10.1109/cso.2011.118.
[13] Caillau, J.-B., Cots, O., and Gergaud, J., "Differential continuation for regular optimal control problems," *Optimization Methods and Software*, Vol. 27, No. 2, 2012, pp. 177-196. https://doi.org/10.1080/10556788.2011.593625.

[14] Chang, D. E., Chichka, D. F., and Marsden, J. E., "Lyapunov-based transfer between elliptic Keplerian orbits," *Discrete and Continuous Dynamical Systems - B*, Vol. 2, No. 1, 2002, pp. 57-67. https://doi.org/10.3934/dcdsb.2002.2.57, URL https://doi.org/10.3934%2Fdcdsb.2002.2.57.

[15] Alizadeh, I., and Villac, B. F., "Static solutions of the Hamilton-Jacobi-Bellman equation for circular orbit transfers," *Journal of Guidance, Control, and Dynamics*, Vol. 34, No. 5, 2011, pp. 1584-1588. https://doi.org/10.2514/1.53394, URL https://doi.org/10.2514%2F1.53394.

[16] Zhang, J., Xiao, Q., and Li, L., "Solution space exploration of low-thrust minimum-time trajectory optimization by combining two homotopies," *Automatica*, Vol. 148, 2023, p. 110798. https://doi.org/10.1016/j.automatica.2022.110798.

[17] Zhu, J., Trélat, E., and Cerf, M., "Geometric optimal control and applications to aerospace," *Pacific Journal of Mathematics for Industry*, Vol. 9, No. 1, 2017, pp. 1-41. https://doi.org/10.1186/s40736-017-0033-4.

[18] Jiang, F., Baoyin, H., and Li, J., "Practical Techniques for Low-Thrust Trajectory Optimization with Homotopic Approach," *Journal of Guidance, Control, and Dynamics*, Vol. 35, No. 1, 2012, p. 245-258. https://doi.org/10.2514/1.52476, URL http://dx.doi.org/10.2514/1.52476.

[19] Sanders, J. A., and Verhulst, F., *Averaging Methods in Nonlinear Dynamical Systems*, Springer New York, 1985. https://doi.org/10.1007/978-1-4757-4575-7, URL https://doi.org/10.1007%2F978-1-4757-4575-7.

[20] Jurdjevic, V., *Geometric Control Theory*, Cambridge University Press, 1996. https://doi.org/10.1017/cbo9780511530036, URL https://doi.org/10.1017%2Fcbo9780511530036.

[21] Ely, T. A., "Transforming mean and osculating elements using numerical methods," *The Journal of the Astronautical Sciences*, Vol. 62, No. 1, 2015, pp. 21-43. https://doi.org/10.1007/s40295-015-0036-2, URL https://doi.org/10.1007%2Fs40295-015-0036-2.