# Efficient Numerical Solution of the Low-Thrust Lambert Problem

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

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

To cite this version:

Lamberto Dell'Elce, Alesia Herasimenka, Aaron J Rosengren, Nicola Baresi. Efficient Numerical Solution of the Low-Thrust Lambert Problem. SFMM 2023 / 33rd AAS/AIAA Space Flight Mechanics Meeting, AAS, Jan 2023, Austin, United States. hal-04050277

HAL Id: hal-04050277

https://inria.hal.science/hal-04050277

Submitted on 29 Mar 2023

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.
EFFICIENT NUMERICAL SOLUTION OF THE LOW-THRUST
LAMBERT PROBLEM

Lamberto Dell'Elce*, Alesia Herasimenka,
Aaron J. Rosengren‡ and Nicola Baresi§

An algorithm for the numerical solution of low-thrust Lambert's problem is pro-
posed. 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 achieved 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 (*i.e.*, to draw a so-called "pork-chop
chart"). Both the averaged flow and the near identity transformation are efficiently
evaluated via the fast Fourier transform algorithm, yielding a fully numerical pro-
cedure.

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. Pork-chop plots are
often generated to illustrate total ΔV cost as a function of departure and arrival dates by recursively
solving Lambert's problems. Recent interest was brought to the low-thrust counterpart of the prob-
lem. Because impulsive maneuvers are not allowed, 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. Differently of the original problem, no exact closed form solution is available, so that
assumptions on the shape and thrust direction are introduced to achieve efficient approximations¹.²
Equivalent of pork-chop charts for low-thrust problems are referred to as bacon plots.³

In our view, two problems need to be addressed when tackling low-thrust transfers on a fixed
maneuvering time. First, minimum thrust magnitude necessary to carry out the maneuver has to
be identified, which, in turn, is a reachability problem. Second, once a sufficiently-large thrust is
chosen, minimum energy maneuvers can be found. This work focuses on the first problem. A
numerical methodology based on the averaging of the extremal flow of the optimal control system⁴
is proposed. First, a reduced-order solution of the averaged two-point boundary value problem

<footnote>*Researcher, Inria center at Université Côte d'Azur, 06902 Valbonne, France</footnote>

<footnote>†PhD researcher, Université Côte d'Azur, CNRS, Inria, LJAD, 06100 Nice, France; ESA contract no 4000134950/21/NL/GLC/my</footnote>

<footnote>‡Professor, University of California San Diego, La Jolla, CA 92093, USA</footnote>

<footnote>§Lecturer in Orbital Mechanics, Surrey Space Centre, University of Surrey, GU2 7XH, Guildford UK</footnote>
(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, sensitivity of the shooting functions are computed. These sensitivities are then used to compute perturbations of the averaged TPBVP associated to short-periodic variations and second-order terms. Finally, the thrust magnitude is identified by reconstructing a first-order approximation of the fast variables from the averaged solution.

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 dates. The methodology is finally applied to an Earth-Venus transfer. Discontinuities of the minimum thrust magnitude map are discussed and related to the cut locus of the problem. We ignore mass consumption and variable thrust profile in this case study because it is problem dependent. However, their inclusion is straightforward by appending mass to the vector of slow variables.

## LOW-THRUST LAMBERT'S PROBLEM

Denote by $x \in \mathcal{M} \subset \mathbb{R}^n$ a set of integrals of motions of the two body problem (specifically, we use equinoctial elements) and by $\varphi \in \mathbb{S}^1$ the mean longitude. Hereafter, $x$ 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, $\varepsilon$, necessary to carry out the transfer between two Keplerian orbits with elements $(x_0, \varphi_0(t_0))$ and $(x_f, \varphi_f(t_f))$ in a time $t_f - t_0$. 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

$$ \begin{array}{ll} \min \quad & \varepsilon \quad \text{subject to:} \\ & \displaystyle \frac{d\mathbf{x}}{dt} = \varepsilon \sum_{i=1}^{m} f_i(\mathbf{x}, \varphi) u_i \\ & \displaystyle \frac{d\varphi}{dt} = \omega(\mathbf{x}) + \varepsilon \sum_{i=1}^{m} g_i(\mathbf{x}, \varphi) u_i \qquad (1) \\ & \mathbf{x}(t_0) = \mathbf{x}_0, \qquad \varphi(t_0) = \varphi_0(t_0) \\ & \mathbf{x}(t_f) = \mathbf{x}_f, \qquad \varphi(t_f) = \varphi_f(t_f) \\ & \|\mathbf{u}\| \le 1 \qquad \forall t \in [t_0, t_f] \end{array} $$

where $m = 3$, and $\mathbf{u}$ is the control variable taking values in the unit sphere of $\mathbb{R}^m$. Vector fields $f_i$ and $g_i$ are periodic with respect to the fast variable, and they stem from Gauss variational equations (GVE) of the chosen set of orbital elements.

The pre-Hamiltonian of Problem (1) is defined as

$$ \mathcal{H}' = \lambda \omega(I) + \varepsilon \sum_{i=1}^{m} [\mathbf{p} \cdot \mathbf{f}_i(\mathbf{x}, \varphi) + \lambda g_i(\mathbf{x}, \varphi)] u_i \quad (2) $$

where $\mathbf{p}$ and $\lambda$ denote adjoints of slow and fast variables, respectively. Application of Pontryagin
maximum principle (PMP) yields the control action of the extremal flow $u^*(x, \varphi, p, \lambda)$, namely

$$u^* = \arg \max_{\|u\| \le 1} \mathcal{H}'(x, p, \varphi, \lambda, \varepsilon, u) = \frac{K}{\|K\|}$$

where the i-th component of $K$ is $K_i = p \cdot f_i(x, \varphi) + \lambda g_i(x, \varphi)$. Replacing $u^*$ into Eq. (2) yields the Hamiltonian of the extremal flow:

$$\mathcal{H} = \lambda \omega(x) + \varepsilon K(x, p, \varphi, \lambda) \quad (3)$$

where $K := \|K\|$. Transversality conditions are such that adjoints are free at both $t_0$ and $t_f$, so that necessary conditions for the optimality of Problem (1) are satisfied by finding a zero of the shooting function $S(p(t_0), \lambda(t_0), \varepsilon)$ defined as

$$S(p(t_0), \lambda(t_0), \varepsilon) = \begin{bmatrix} x^*(t_f) - x_f \\ \varphi^*(t_f) - \varphi_f(t_f) \\ \|p(0)\| - 1 \end{bmatrix} \quad (4)$$

where $x^*(t)$ and $\varphi^*(t)$ denote trajectories of $x$ and $\varphi$ emanated by $\mathcal{H}$ with initial conditions $x_0, \varphi_0, 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 adjoints.

## AVERAGED PROBLEM

States and adjoints are decomposed as

$$\begin{align}
x &= X + \varepsilon [\delta X + \nu_X(X, P, \hat{\varphi})] \nonumber \\
p &= P + \varepsilon [\delta P + \nu_P(X, P, \hat{\varphi})] \nonumber \\
\varphi &= \varphi_0 + \frac{\Phi}{\varepsilon} + \delta\Phi + \varepsilon\nu_\Phi(X, P, \hat{\varphi}) \tag{5} \\
\lambda &= \varepsilon (\Lambda + \nu_\Lambda(X, P, \hat{\varphi})) + \varepsilon^2\delta\Lambda (X, P, \hat{\varphi}, \Lambda, \delta X, \delta P) \nonumber
\end{align}$$

Hereafter, $X$ and $P$ are referred to as averaged state and adjoint, respectively, and $\nu_\cdot$ denotes short-periodic variations. Second-order terms*, namely $\delta X, \delta P, \delta\Phi$, and $\delta\Lambda$, are necessary to have a precise enough approximation of the fast variable. 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 $X, P, \Phi, \Lambda, \delta X, \delta P, \delta\Phi, \delta\Lambda$, and any $\nu_\cdot$ (once stripped out of their power of $\varepsilon$) are all $O(1)$ throughout the entire trajectory. Specifically, 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, adjoint $\Lambda$ is well defined when $\varepsilon$ approaches zero as discussed in.⁴ Short-periodic variations and $\delta\Lambda$ are evaluated at the angle:

$$\hat{\varphi} = \varphi_0 + \frac{\Phi}{\varepsilon} + \delta\Phi$$

The decomposition in Eq. (5) 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

<footnote>*The reference "second order" may stride with the order of magnitude these terms in the expansion of Eq. (5). Instead, they are referred to as second order terms because time derivatives of $\varepsilon\delta X$ and $\varepsilon\delta P$ are $O(\varepsilon^2)$. For the sake of clarity and without arguing further on this notation, also $\delta\Phi$ and $\delta\Lambda$ are referred to as "second order terms" for analogy.</footnote>
new variables such as $\delta X$ and $\delta P$. The main advantage of the proposed decomposition is that both the averaged counterpart of Problem (4) and the computation of second-order terms are $\epsilon$-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. (5) is not a full second-order approximation.

The average 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 $$

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

Equations of motion of averaged, short-periodic, and second-order terms are deduced by first substituting Eq. (5) into the Hamilton equations associated to $\mathcal{H}(x, p, \varphi, \lambda, \epsilon)$, and then by expanding these equations up to the second-order in $\epsilon$. Time derivatives of secular and second-order terms are finally achieved by averaging the resulting equations and collecting terms of the appropriate order (hence, they are independent of the fast variable by construction), whereas formal solution of short periodic variations is obtained by compensating first-order oscillatory dynamics. Outcome of this process is detailed in the remainder of the section.

To simplify the notation, omitted arguments of functions and their derivatives are understood as the average value of slow variables, namely $X, P, \epsilon\Lambda = 0$, and fast variable $\hat{\varphi}$. For example $K = K(X, P, \hat{\varphi}, 0)$, $\bar{K} = \bar{K}(X, P, 0)$, $\frac{\partial K}{\partial \bar{X}} = \frac{\partial K}{\partial x}\bigg|_{x=X,p=P,\varphi=\hat{\varphi},\lambda=0}$, $\omega = \omega(X)$, and $\frac{\partial \omega}{\partial \bar{X}} = \frac{\partial \omega}{\partial x}\bigg|_{x=X}$.

## First-order terms and averaged shooting problem

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

$$ d\tau = \epsilon dt $$

Motion of averaged states and adjoints is given by

$$ \frac{d\mathbf{X}}{d\tau} = \frac{\partial \bar{K}}{\partial \mathbf{P}}, \quad \frac{d\mathbf{P}}{d\tau} = -\frac{\partial \omega}{\partial \mathbf{X}}\Lambda - \frac{\partial \bar{K}}{\partial \mathbf{X}}, \quad \frac{d\Phi}{d\tau} = \omega(\mathbf{X}), \quad \frac{d\Lambda}{d\tau} = 0 \quad (6) $$

Equations (6) are hamiltonian with the first-order averaged Hamiltonian defined as

$$ \bar{\mathcal{H}} = \omega(\mathbf{X})\Lambda + \bar{K}(\mathbf{X}, \mathbf{P}, 0) $$

Variable $\Phi$ is cyclic and its momentum $\Lambda$ is a constant of motion.

The averaged shooting function, counterpart of Eq. (4), is

$$ \bar{S}(\mathbf{P}(0), \tau_f | \mathbf{x}_0, \mathbf{x}_f, \Lambda) = \begin{bmatrix} \mathbf{X}^*(\tau_f) - \mathbf{x}_f \\ \|\mathbf{P}(0)\| - 1 \end{bmatrix} \quad (7) $$

We note that Problem (7) is independent of both $\epsilon$ and $\Phi$, and it is a minimum time problem. Conversely, maneuvering time is fixed in both Problem (1) and its shooting function, Eq. (4). This
dualism stems from the relation between minimum-time flow and boundary of the reachable set [5, Chap 10, Theorem 1]. The notation $\bar{S}(P(0), \tau_f|x_0, x_f, \Lambda)$ is introduced to emphasize that $P(0)$ and $\tau_f$ are shooting variables, whereas $x_0, x_f$, and $\Lambda$ are parameters. Problems (4) and (7) have the same boundary conditions of slow variables, namely $x_0$ and $x_f$. Corrections due to short-periodic variations and second-order terms are accounted for by computing sensitivities of $\bar{S}$, as detailed in the next Section.

## Short-periodic variations

Time derivative of short periodic variations is

$$ \frac{d \nu_{(\cdot)}}{d \tau} = \frac{d \Phi}{d \tau} \frac{\partial \nu_{(\cdot)}}{\partial \hat{\varphi}} + \varepsilon \left[ \frac{d \delta \Phi}{d \tau} \frac{\partial \nu_{(\cdot)}}{\partial \hat{\varphi}} + \frac{d \mathbf{X}}{d \tau} \cdot \frac{\partial \nu_{(\cdot)}}{\partial \mathbf{X}} + \frac{d \mathbf{P}}{d \tau} \cdot \frac{\partial \nu_{(\cdot)}}{\partial \mathbf{P}} \right] = \omega \frac{\partial \nu_{(\cdot)}}{\partial \hat{\varphi}} + \mathcal{O}(\varepsilon) \quad (8) $$

The average of Eq. (8) is zero because derivatives of averaged and second-order variables do not depend on the fast variable and $\bar{\nu}_{(\cdot)} = 0$ by construction.

Short-periodic variations are obtained by solving

$$ \begin{aligned} \omega \frac{\partial \nu_X}{\partial \hat{\varphi}} &= \frac{\partial K}{\partial \mathbf{P}} - \frac{\partial \bar{K}}{\partial \bar{\mathbf{P}}} \\ \omega \frac{\partial \nu_P}{\partial \hat{\varphi}} &= -\frac{\partial K}{\partial \mathbf{X}} + \frac{\partial \bar{K}}{\partial \bar{\mathbf{X}}} - \frac{\partial \omega}{\partial \mathbf{X}} \nu_{\Lambda} \\ \omega \frac{\partial \nu_{\Phi}}{\partial \hat{\varphi}} &= \frac{\partial K}{\partial \lambda} - \frac{\partial \bar{K}}{\partial \bar{\lambda}} + \frac{\partial \omega}{\partial \mathbf{X}} \nu_X \\ \omega \frac{\partial \nu_{\Lambda}}{\partial \hat{\varphi}} &= -\frac{\partial K}{\partial \hat{\varphi}} \end{aligned} \quad (9) $$

with the conditions to resolve integration constants

$$ \bar{\nu}_X = 0, \quad \bar{\nu}_P = 0, \quad \bar{\nu}_\Phi = 0, \quad \bar{\nu}_\Lambda = 0 \quad (10) $$

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. Formal solution of Eqs. (9) and (10) is

$$ \begin{aligned} \nu_X &= -\frac{i}{\omega} \sum_{k \in \mathbb{Z}^0} \frac{1}{k} \frac{\partial K^{(k)}}{\partial \mathbf{P}} e^{ik\hat{\varphi}} \\ \nu_P &= \frac{i}{\omega} \sum_{k \in \mathbb{Z}^0} \frac{1}{k} \left( \frac{\partial K^{(k)}}{\partial \mathbf{X}} - \frac{K^{(k)}}{\omega} \frac{\partial \omega}{\partial \mathbf{X}} \right) e^{ik\hat{\varphi}} \\ \nu_\Phi &= -\frac{i}{\omega} \sum_{k \in \mathbb{Z}^0} \frac{1}{k} \left( \frac{\partial K^{(k)}}{\partial \lambda} - \frac{i}{k\omega} \frac{\partial \omega}{\partial \mathbf{X}} \cdot \frac{\partial K^{(k)}}{\partial \mathbf{P}} \right) e^{ik\hat{\varphi}} \\ \nu_\Lambda &= -\frac{K(\mathbf{X}, \mathbf{P}, \hat{\varphi}, 0) - \bar{K}(\mathbf{X}, \mathbf{P}, 0)}{\omega(\mathbf{X})} \end{aligned} \quad (11) $$

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

Equations of motion of $\delta X$, $\delta P$, and $\delta\Phi$ are obtained by expanding Eq. (3) 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{align}
\frac{d\delta\mathbf{X}}{d\tau} &= \overline{\frac{\partial^2 K}{\partial\mathbf{X}\partial\mathbf{P}}(\delta\mathbf{X} + \boldsymbol{\nu}_{\mathbf{X}})} + \overline{\frac{\partial^2 K}{\partial\mathbf{P}^2}(\delta\mathbf{P} + \boldsymbol{\nu}_{\mathbf{P}})} + \overline{\frac{\partial^2 K}{\partial\varphi\partial\mathbf{P}}\nu_{\Phi}} + \overline{\frac{\partial^2 K}{\partial\lambda\partial\mathbf{P}}(\Lambda + \nu_{\Lambda})} \nonumber \\
\frac{d\delta\mathbf{P}}{d\tau} &= -\overline{\frac{\partial^2 K}{\partial\mathbf{X}^2}(\delta\mathbf{X} + \boldsymbol{\nu}_{\mathbf{X}})} - \overline{\frac{\partial^2 K}{\partial\mathbf{P}\partial\mathbf{X}}(\delta\mathbf{P} + \boldsymbol{\nu}_{\mathbf{P}})} - \overline{\frac{\partial^2 K}{\partial\varphi\partial\mathbf{X}}\nu_{\Phi}} - \overline{\frac{\partial^2 K}{\partial\lambda\partial\mathbf{X}}(\Lambda + \nu_{\Lambda})} \tag{12} \\
&\quad - \Lambda \overline{\frac{\partial^2\omega}{\partial\mathbf{x}}\delta\mathbf{X}} - \overline{\frac{\partial\omega}{\partial\mathbf{x}}\delta\Lambda} \nonumber \\
\frac{d\delta\Phi}{d\tau} &= \overline{\frac{\partial\omega}{\partial\mathbf{X}}\cdot\delta\mathbf{X}} + \overline{\frac{\partial\mathbf{K}}{\partial\lambda}} \nonumber
\end{align}
$$

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

Expansion of the hamiltonian up to the second order in ε enables an explicit expression of δΛ as
a function of the other variables, namely

$$
\delta\Lambda = -\frac{1}{\omega} \left[ \left( \frac{\partial K}{\partial \mathbf{X}} + (\Lambda + \nu_{\Lambda}) \frac{\partial \omega}{\partial \mathbf{X}} \right) \cdot (\delta \mathbf{X} + \boldsymbol{\nu}_{\mathbf{X}}) \right. \\
\left. \qquad + \frac{\partial K}{\partial \mathbf{P}} \cdot (\delta \mathbf{P} + \boldsymbol{\nu}_{\mathbf{P}}) + \frac{\partial K}{\partial \Phi} \nu_{\Phi} + \frac{\partial K}{\partial \lambda} (\Lambda + \boldsymbol{\nu}_{\Lambda}) \right]
$$

Denoting

$$
\boldsymbol{\nu}_{\mathbf{X}}^{(k)} := -\frac{i}{\omega k} \frac{\partial K^{(k)}}{\partial \mathbf{P}}, \quad \boldsymbol{\nu}_{\mathbf{P}}^{(k)} := \frac{i}{\omega k} \left( \frac{\partial K^{(k)}}{\partial \mathbf{X}} - \frac{K^{(k)}}{\omega} \frac{\partial \omega}{\partial \mathbf{X}} \right), \quad \boldsymbol{\nu}_{\Phi}^{(k)} := \frac{1}{\omega} \left( \frac{\partial \omega}{\partial \mathbf{X}} \cdot \boldsymbol{\nu}_{\mathbf{X}}^{(k)} + \frac{\partial K^{(k)}}{\partial \lambda} \right)
$$

the k-th coefficients of the series in Eq. (11), the averaged value of δΛ (which appears in the time derivative of δP ) is given by

$$
\begin{align*}
\overline{\delta\Lambda} ={}& -\frac{1}{\omega} \left[ \left( \frac{\partial K^{(0)}}{\partial \mathbf{X}} + \Lambda \frac{\partial \omega}{\partial \mathbf{X}} \right) \cdot \delta \mathbf{X} + \frac{\partial K^{(0)}}{\partial \mathbf{P}} \cdot \delta \mathbf{P} + \frac{\partial K^{(0)}}{\partial \lambda} \Lambda \right] \\
& - \frac{1}{\omega} \sum_{\mathbb{Z}^0} \left[ \left( \frac{\partial K^{(-k)}}{\partial \mathbf{X}} + \nu_{\Lambda}^{(-k)} \frac{\partial \omega}{\partial \mathbf{X}} \right) \cdot \boldsymbol{\nu}_{\mathbf{X}}^{(k)} + \frac{\partial K^{(-k)}}{\partial \mathbf{P}} \cdot \boldsymbol{\nu}_{\mathbf{P}}^{(k)} \right. \\
& \qquad \left. - \left( \nu_{\Phi}^{(-k)} + \frac{1}{\omega} \frac{\partial K^{(-k)}}{\partial \lambda} \right) K^{(k)} \right]
\end{align*}
$$

where we extensively used the integration rule of the average product of two periodic functions
$f(\varphi)$ and $g(\varphi)$

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

<footnote>*The convention $\left[ \frac{\partial^2 K}{\partial \mathbf{X} \partial \mathbf{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 \mathbf{X} \partial \mathbf{P}} = \left( \frac{\partial^2 K}{\partial \mathbf{P} \partial \mathbf{X}} \right)^T$.</footnote>
Replacing Eq.(11) into Eq. (12), integrating by parts terms that include derivatives of $\varphi$, and again using the rule of Eq. (13), yields

$$ 
\begin{aligned}
\frac{d}{d\tau} \delta \mathbf{X} &= \frac{\partial^2 K^{(0)}}{\partial \mathbf{X} \partial \mathbf{P}} \delta \mathbf{X} + \frac{\partial^2 K^{(0)}}{\partial \mathbf{P}^2} \delta \mathbf{P} + \frac{\partial^2 K^{(0)}}{\partial \lambda \partial \mathbf{P}} \Lambda \\
&\quad + \sum_{k \in \mathbb{Z}^0} \left( \frac{\partial^2 K^{(-k)}}{\partial \mathbf{X} \partial \mathbf{P}} \nu_{\mathbf{X}}^{(k)} + \frac{\partial^2 K^{(-k)}}{\partial \mathbf{P}^2} \nu_{\mathbf{P}}^{(k)} - \frac{\partial K^{(-k)}}{\partial \mathbf{P}} \nu'_{\Phi}^{(k)} - \frac{\partial^2 K^{(-k)}}{\partial \lambda \partial \mathbf{P}} \frac{K^{(k)}}{\omega} \right) \\
\frac{d}{d\tau} \delta \mathbf{P} &= - \left( \frac{\partial^2 K^{(0)}}{\partial \mathbf{X}^2} + \Lambda \frac{\partial^2 \omega}{\partial \mathbf{X}^2} \right) \delta \mathbf{X} - \frac{\partial^2 K^{(0)}}{\partial \mathbf{P} \partial \mathbf{X}} \delta \mathbf{P} - \frac{\partial^2 K^{(0)}}{\partial \lambda \partial \mathbf{X}} \Lambda - \frac{\partial \omega}{\partial \mathbf{X}} \overline{\delta \Lambda} \\
&\quad - \sum_{k \in \mathbb{Z}^0} \left( \frac{\partial^2 K^{(-k)}}{\partial \mathbf{X}^2} \nu_{\mathbf{X}}^{(k)} + \frac{\partial^2 K^{(-k)}}{\partial \mathbf{P} \partial \mathbf{X}} \nu_{\mathbf{P}}^{(k)} - \frac{\partial K^{(-k)}}{\partial \mathbf{X}} \nu'_{\Phi}^{(k)} - \frac{\partial^2 K^{(-k)}}{\partial \lambda \partial \mathbf{X}} \frac{K^{(k)}}{\omega} \right) \\
\frac{d}{d\tau} \delta \Phi &= \frac{\partial \omega}{\partial \mathbf{X}} \cdot \delta \mathbf{X} + \frac{\partial K^{(0)}}{\partial \lambda}
\end{aligned}
\quad (14) $$

Equation (14) is a linear system in $\delta \mathbf{X}$, $\delta \mathbf{P}$, and $\delta \Phi$. 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⁶ and FFT.

## ALGORITHM FOR THE GENERATION OF BACON PLOTS

We propose an algorithm for the generation of approximated solutions of Problem (1).

1. Find a one-parameter family of solutions of the average problem

$$ 
\begin{aligned}
\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{x}_0, \mathbf{x}_f, \Lambda) = 0
\end{aligned}
$$

We note that the same initial conditions of the osculating problem, $\mathbf{x}_0$ and $\mathbf{x}_f$, are used here. After finding the solution for $\Lambda = 0$, a numerical continuation scheme can be used ($\Lambda$ serves as 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 parameter to streamline the notation.

2. Compute sensitivities

$$ \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 \mathbf{x}_0}, \frac{\partial \bar{S}}{\partial \mathbf{x}_f}, \frac{\partial \Phi_f}{\partial \tau_f}, \frac{\partial \Phi_f}{\partial \mathbf{P}_0}, \frac{\partial \Phi_f}{\partial \mathbf{x}_0} $$

where $\Phi_f = \Phi(\tau_f^\Lambda | \mathbf{x}_0, \mathbf{P}_0^\Lambda, \Lambda)$ is the final re-scaled fast variable of the average system.

3. Evaluate perturbations of the boundary conditions

find $\delta \tau_f^\Lambda$, $\delta \mathbf{P}_0^\Lambda$, $\delta \Phi_f^\Lambda$ s.t.:

$$ 
\begin{aligned}
\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{x}_0} \delta \mathbf{X}_0 - \frac{\partial \bar{S}}{\partial \mathbf{x}_f} \delta \mathbf{X}_f \\
\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{x}_0} \delta \mathbf{X}_0 + \delta \Phi (\tau_f^\Lambda)
\end{aligned}
$$
where $\delta X_0$ and $\delta X_f$ include both short-periodic variations and second-order terms, namely

$$ \delta X_0 = -\nu_X (x_0, P_0^\Lambda, \varphi_0(t_0)) $$

$$ \delta X_f = -\nu_X (x_f, P(\tau_f^\Lambda), \varphi_f(t_f)) - \delta X(\tau_f^\Lambda) $$

Here, $\delta X(\tau_f^\Lambda)$ and $\delta \Phi(\tau_f^\Lambda)$ are evaluated by integrating Eq. (14) with initial conditions $\delta X = \delta P = 0$ and $\delta \Phi = 0$.

### 4. Evaluate minimum $\varepsilon$ and identify the good members of the family of average solutions

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

where

$$ \varepsilon^\Lambda = \frac{\tau_f^\Lambda}{t_f - t_0 - \delta \tau_f^\Lambda} $$

is such that the duration of the trajectory equals $t_f - t_0$.

Concerning the complexity of the methodology: Step 1 implies first the solution of a single average problem for $\Lambda = 0$ (this is the only shooting problem that has to be solved without a priori knowledge). Then, differential continuation can be used to obtain the solution in the range $[\Lambda^{min}, \Lambda^{max}]$. Step 2 requires the numerical integration of the average system's linearization. AD can be used to evaluate derivatives. Specifically, we used the software Hampath* described in⁶ to compute derivatives of both $\bar{H}$ and $\bar{S}$ which are necessary in Steps 1 and 2, respectively. Hampath only requires in 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 $(P_0^\Lambda, \tau_f^\Lambda)$. Step 3 requires the integration of second-order terms, the evaluation of short-periodic corrections (a single FFT is sufficient to evaluate $\nu(\cdot)$ for any $\varphi$ and, as such, for any $t_0$ or $t_f$), and the solution of n-dimensional linear systems. Step 4 only requires elementary algebraic operations.

We stress that all steps only use the average 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 $t_0$ and $t_f$.

## CASE STUDY

The proposed case study consists of an Earth-Venus low-thrust transfer. Table 1 lists simulation parameters necessary for repeatability.

The relevance of second-order terms is disclosed in Figure 1. Contributions of average, second-order, and short-periodic variations are pointed out in Figure 1(a), where a sample trajectory of the adjoint $p_4$ is drawn. Drift between $p_4$ and $P_4$ is of order $\varepsilon$ at the end of the maneuver, and (for this particular trajectory) including $\delta P_4$ is more important than short-periodic variations $\nu_{P_4}$ to reduce

<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. Earth-Venus transfer. Elements of departure and arrival orbits.</caption><thead><tr><th></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 · 10<sup>-2</sup></td><td>6.78 · 10<sup>-3</sup></td></tr><tr><td>Ecliptic inclination [deg]</td><td>0</td><td>3.3947</td></tr><tr><td>Argument of perifocus [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_2.png" /><figcaption>(a) Adjoint $p_4$. Here, $\Lambda = -1$ and $\varepsilon = 0.05$ mN/kg.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>(b) Error between true fast variable and first-order approximation, $\frac{\Phi}{\varepsilon}$. We note that such error converges to the second-order correction, $\delta\Phi$, as $\varepsilon$ approaches 0. Here, $\Lambda = -1$.</figcaption></figure>

<caption>Figure 1. Importance of second-order terms.</caption>
the error with the original trajectory. In general, contributions of second-order terms and short-periodic variations are of the same order of magnitude. Necessity to account for second-order terms stems from the need of a first-order estimation of the fast variable. As depicted in Figure 1(b), error between $\varphi$ and $\Phi$ converges to $\delta\Phi$ (which is of order one) as $\varepsilon$ approaches zero. Indeed, a drawback of the proposed framework is that quality of the approximation is tied to the small magnitude of $\varepsilon$. Hence, the methodology is not adequate to approximate short transfers, which are characterized by large values of $\varepsilon$.

An example of bacon plot is illustrated in Figure 2(a). Solutions along red-dashed lines have the same maneuvering time, which is larger for upper lines. Variations of $\varepsilon$ along these lines are small when long transfers are envisaged, i.e., importance of the position on departure and arrival orbits dwindles when several revolutions are necessary to complete the maneuver. Bacon plots are not smooth. 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 Figure 2(a), occur when two existing local solutions have the same $\varepsilon$. This can be better appreciated by inspecting Figure 2(b), which depicts the value of $\Lambda$. In fact, each point of the bacon plot is associated to a precise member of the one-parameter family of average 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 Figure 3, which shows the $\varepsilon$ of local solutions along the red-dashed line of Figure 2(b). Two distinct solutions are identified for some values of the departure date, but only the one with lower $\varepsilon$ is of interest. Red dot 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 Figure 2 (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.

## CONCLUSION

This paper offers an efficient numerical solution to the low-thrust Lambert problem which may serve as a tool for preliminary mission design. Overall, the methodology requires the solution of a single simplified TPBVP to draw a full chart of minimum thrust magnitude for any combination of departure and arrival dates. A major advantage of the methodology is the immediate identification of cut points enabling the straightforward discarding of non-minimizing local solutions.

Formal solution of short-periodic variations and an implementable expression of time derivatives of second-order terms are provided. These expressions can be efficiently evaluated by means of the FFT algorithm.

Variable thrust profiles were neglected in this study, but their inclusion in the methodology is straightforward. Perspectives of the work encompass solution of minimum-fuel transfers, and the possibility to account for $\Delta V$ at the departure and arrival dates, yielding a "hybrid" Lambert problem.

## REFERENCES

[1] D. Izzo, "Lambert's problem for exponential sinusoids," *Journal of guidance, control, and dynamics*, Vol. 29, No. 5, 2006, pp. 1242-1245.

[2] G. Avanzini, A. Palmas, and E. Vellutini, "Solution of Low-Thrust Lambert Problem with Perturbative Expansions of Equinoctial Elements," *Journal of Guidance, Control, and Dynamics*, Vol. 38, sep 2015, pp. 1585-1601, 10.2514/1.g001018.
<figure><img src="image_4.png" /><figcaption>(a) Minimum thrust magnitude.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>(b) Continuation parameter, $\Lambda$, of best solutions identified in Step 4.</figcaption></figure>

<caption>**Figure 2.** Earth-Venus transfer for various deparutre and arrival dates.</caption>
<figure><img src="image_6.png" /><figcaption>Figure 3. Value of the small parameter $\epsilon$ for local solutions along the red-dashed line of Figure 2(b). The red dot points out where discontinuities occur in the maps of Figure 2.</figcaption></figure>

[3] R. C. Woolley, F. Laipert, A. K. Nicholas, and Z. Olikara, "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.

[4] L. Dell'Elce, J.-B. Caillau, and J.-B. Pomet, "On the convergence of time-optimal maneuvers of fast-oscillating control systems," *2021 European Control Conference (ECC)*, IEEE, 2021, pp. 2008-2013.

[5] V. Jurdjevic, *Geometric Control Theory*. Cambridge University Press, dec 1996, 10.1017/cbo9780511530036.

[6] J.-B. Caillau, O. Cots, and J. Gergaud, "Differential continuation for regular optimal control problems," *Optimization Methods and Software*, Vol. 27, No. 2, 2012, pp. 177-196.