# Particle Swarm Optimization-Based Co-State Initialization for Low-Thrust
Minimum-Fuel Trajectory Optimization

Grant R. Hecht<EMAIL_ADDRESS>Eleonora M. Botta<EMAIL_ADDRESS>Department of Mechanical and Aerospace Engineering, University at Buffalo,
Buffalo, NY, United States

###### Abstract

In this paper, Particle Swarm Optimization with energy-to-fuel continuation is
proposed for initializing the co-state variables for low-thrust minimum-fuel
trajectory optimization problems in the circular restricted three-body
problem. Particle Swarm Optimization performs a search of the solution space
by minimizing the weighted sum of squares of the two-point boundary-value
problem final boundary condition residuals for the minimum-energy problem.
Next, an energy-to-fuel homotopy is employed to transition the minimum-energy
trajectory to a minimum-fuel trajectory, starting from the generated guess.
The proposed methodology is applied to two low-thrust transfer problems in the
Earth-Moon system: a transfer from a geostationary transfer orbit to an L1
halo orbit, and a transfer from an L2 halo orbit to an L1 halo orbit. The
resulting minimum-fuel trajectories are validated with the literature. It is
demonstrated that the methodology can successfully generate guesses for the
initial co-state variables which converge to a solution for both scenarios. A
strategically chosen particle swarm size is shown to improve the convergence
rate of the methodology. The proposed approach is of simple implementation,
can be easily extended to other trajectory optimization problems, facilitates
the discovery of multiple candidate trajectories, and does not require a user-
provided guess, all of which are advantageous features for the preliminary
phase of mission design.

###### keywords:

Low-Thrust Trajectory Optimization , Minimum-Fuel Trajectory Optimization,
Homotopic Continuation , Particle Swarm Optimization , Optimal Control

††journal: Acta Astronautica

CR3BP

    Circular Restricted Three Body Problem
TPBVP

    Two-Point Boundary-Value Problem
PSO

    Particle Swarm Optimization
STM

    State Transition Matrix
PMP

    Pontryagin’s Maximum Principle
COV

    Calculus of Variations
ACT

    Adjoint Control Transformation

##  1 Introduction

Low-thrust propulsion has gained much attention in recent years and is now
often considered the best option for a wide range of space missions. The
associated reduction in required propellant mass compared to traditional
chemical propulsion results in cost savings and an increase in the feasible
payload ratio. Low-thrust propulsion does, however, require long, continuous
thrusting arcs and trajectories often spanning many revolutions about a large
body, which can lead to high-dimensional optimization problems with many local
minima. Low-thrust spacecraft trajectory optimization problems most commonly
focus on minimizing the time of flight, fuel required, or combinations of
these conflicting objectives. As the name suggests, minimum-time trajectory
optimization is necessary when mission requirements specify that the
spacecraft must reach its destination in a timely manner, regardless of the
fuel required, as is often the case when transporting human life or trying to
acquire time-sensitive scientific data. On the other hand, minimum-fuel
trajectory optimization is useful when time is not such an important
constraint and it is desired to reduce the required fuel, which allows for
increasing the feasible payload ratio and decreasing mission costs. Many
approaches for optimizing low-thrust spacecraft trajectories have been
proposed and applied successfully. Typically, approaches can be separated into
two categories known as direct and indirect methods. Hybrid approaches also
exist which combine the characteristics of direct and indirect methods [1]. A
broad review of the state-of-the-art techniques employed for spacecraft
trajectory optimization can be found in work by Chai et al. [2].

Direct trajectory optimization methods formulate the trajectory optimization
problem by parameterizing the state and control variables, forming a nonlinear
programming (NLP) problem [1]. These methods can be further categorized into
shooting, collocation, dynamic programming, or differential inclusion-based
methods [2]. When considering trajectories spanning tens to hundreds of days
with many separate, continuous thrusting arcs – as is common for low-thrust
propulsion –, this parameterization can lead to very high-dimensional problems
which may be intractable if not strategically formulated, due to the
approaching curse of dimensionality. Additionally, this parameterization
generally enforces the control to take some predetermined form, which can
result in suboptimal solutions. Nevertheless, direct methods continue to be an
attractive approach for solving low-thrust trajectory optimization problems
for several reasons. When compared to indirect methods, direct methods can
more easily incorporate complicated constraints and tend to be much less
sensitive to the quality of the initial guess supplied to the algorithm [1].
Furthermore, physical intuition can be exploited when an initial guess for the
solution is determined and supplied to an NLP solver.

Indirect methods apply Calculus of Variations (COV) and Pontryagin’s Maximum
Principle (PMP) to analytically derive the necessary conditions for
optimality, introducing additional co-state variables and formulating the
problem as a Two-Point Boundary-Value Problem (TPBVP) [3]. Collocation or
finite element-based methods can then be used to approximate the solution of
the TPBVP, or a shooting-based method may be used to solve the TPBVP directly
[2]. Collocation or finite element-based approaches can lead to useful
physical insight about the problem while providing approximate solutions [1].
Alternatively, shooting-based methods result in a fully continuous optimal
control law for the low-thrust spacecraft trajectory. Unfortunately, the
convergence of the indirect shooting-based approach is typically not
guaranteed and is highly dependent on the selection of the initial co-state
variables used to initialize the algorithm [1, 3].

Typically, the co-state variables are abstract quantities with no clear
physical significance and the convergence radius of indirect shooting methods
is often very small. Therefore, determining an initial choice of co-state
variables to initialize a shooting routine can be a daunting task. To mitigate
some of this difficulty, several methods have been proposed to assist in
initializing the co-state variables. One such approach, first proposed by
Dixon and Biggs [4, 5] and termed Adjoint Control Transformation (ACT),
involves relating the co-state variables to alternative variables of physical
significance, thereby decreasing the sensitivity of the problem and allowing
for the use of physical insight when selecting a guess. Application of ACT to
problems in astrodynamics is discussed, for example, in References [6, 7, 8,
9]. An alternative approach involves formulating an approximation of the
initial co-state variables to avoid the need for a user-provided guess. For
example, within the field of astrodynamics, Thorne and Hall approximated the
initial co-state variables by solving a similar problem with an analytical
solution by assuming zero gravity and a constant spacecraft mass [10], whereas
Lee et al. derived an approximate relationship between the co-state variables
and the initial constrained physical state based on the behavior of the
initial co-states for a range of optimal trajectories [11, 12]. Unfortunately,
these techniques are restricted to planar motion or transfers of less than one
revolution. Another possibility to estimate the initial co-state variables
involves employing a solution acquired via direct trajectory optimization.
Examples of such methods include a technique that exploits approximate
solutions acquired via shape-based methods and nonlinear least squares [3],
and an approach that allows one to relate the Lagrange multipliers of the
Karush-Kuhn-Tucker equations associated with various collocation methods to
the co-states of the optimal control problem [13, 14, 15, 16, 17]. Of course,
a primary difficulty of such approaches is that a solution must be found
through the appropriate direct trajectory optimization method first, a task
that can itself be difficult. Furthermore, once a solution is discovered via
direct trajectory optimization, successful estimation of the co-state
variables is not guaranteed.

A particularly interesting approach, proposed by Pontani and Conway, employs
Particle Swarm Optimization (PSO) to select co-state variables which minimize
the sum of the absolute value of the final boundary condition residuals [18,
19], solving the TPBVP directly. Such an approach is highly desirable as it
can, in principle, be applied to any TPBVP and does not require a user to
provide any guess (although box constraints must be defined to limit the
search space). Unfortunately, due to the extreme sensitivity of the final
boundary condition residuals to perturbations in the initial co-state
variables for many indirect trajectory optimization problems, this initial
application of PSO was limited to planar minimum-time transfers in the two-
body problem, which reduces the number of unknown variables and constraints.

In this paper, a method is proposed which builds upon the work by Pontani and
Conway [18, 19] and allows for solving non-planar low-thrust minimum-fuel
trajectory optimization problems with the indirect shooting approach by
leveraging PSO to initialize the co-state variables. PSO co-state
initialization and an energy-to-fuel homotopy technique are coupled [20] to
widen the convergence radius of the TPBVP and reduce the difficulty of the
parameter optimization problem to be solved by PSO. PSO generates co-states
for the easier-to-solve minimum-energy problem, before transitioning the
solution to minimum-fuel through continuation, all without the need of a user-
provide guess. Herein, the proposed methodology is presented in the framework
of the Circular Restricted Three Body Problem (CR3BP), and important
considerations in regard to solution optimality and computational performance
are discussed. Case studies, based on previous works [21, 22], are presented
which shed light on the promising performance of the proposed methodology
while allowing for its validation. Guidance on the selection of the particle
swarm size, a heuristic parameter of the optimization algorithm which has a
significant effect on co-state initialization performance, is also provided.

The remainder of this paper is organized as follows. In Section 2, the problem
is formulated, beginning with the CR3BP dynamics, and then deriving the
indirect minimum-fuel TPBVP with an energy-to-fuel homotopy through the
application of COV and PMP. The proposed method of co-state initialization is
discussed in Section 3, along with the single shooting algorithm used to
robustly transition PSO-initialized co-state values, corresponding to the
minimum-energy problem, to a solution of the minimum-fuel problem. In Section
4, two minimum-fuel trajectory optimization scenarios (i.e., a transfer from a
geostationary transfer orbit to L1 halo orbit, and a transfer from an L2 halo
orbit to an L1 halo orbit) are employed to validate the methodology,
demonstrate the performance of the proposed method, and analyze the effect of
the particle swarm size on convergence. Finally, conclusions are made which
summarize the achievements of this work and future steps.

##  2 Indirect Optimal Control Formulation

In this Section, the indirect optimal control TPBVP is formulated for the
minimum-fuel problem with terminal state constraints and a fixed time of
flight. This is done using the CR3BP dynamics model, which is provided and
discussed. Throughout the TPBVP formulation, an energy-to-fuel homotopy
technique is incorporated by introducing a perturbed energy term to the
minimum-fuel cost function, before the Hamiltonian is formed and the
analytical necessary conditions of optimality are derived.

###  2.1 Circular Restricted Three Body Problem Dynamics

The CR3BP dynamics describe the motion of a body in a synodic (rotating)
reference frame under the influence of two primary bodies which travel along
co-planar circular paths centered about the barycenter of the system [23]. The
mass of the third body is indicated with $m$ and is assumed to be much lower
than those of the primary bodies $m_{1}$ and $m_{2}$; therefore it has no
effect on their motion. In the following, all variables are non-
dimensionalized such that the angular velocity of both primary bodies, the sum
of the masses of the primary bodies, and the distance between the primary
bodies are set to unity to improve numerical stability. If the third body is a
spacecraft with the capability to apply propulsive thrust in any direction,
the equations of motion describing the evolution of the spacecraft’s position
and velocity in a synodic reference frame, along with its mass are given by
[9]

$\displaystyle\dot{\mathbf{x}}=\mathbf{f}(\mathbf{x},\hat{\boldsymbol{\alpha}},u)=\begin{bmatrix}\dot{\mathbf{r}}\\\
\dot{\mathbf{v}}\\\ \dot{m}\end{bmatrix}=\begin{bmatrix}\mathbf{v}\\\
\mathbf{g}(\mathbf{r})+\mathbf{h}(\mathbf{v})+uT_{max}\hat{\boldsymbol{\alpha}}/m\\\
-uT_{max}/c\end{bmatrix}$ (1)

where $\mathbf{r}$ and $\mathbf{v}$ are the non-dimensionalized position and
velocity vectors of the spacecraft, $T_{max}$ is its maximum capable thrust,
and $c$ is the propellant exhaust velocity, which is assumed to be constant
and given by $c=I_{sp}g_{0}$, with $I_{sp}$ the specific impulse of the
propellant and $g_{0}=9.81$ m/s2. The control variables are $u$ and
$\hat{\boldsymbol{\alpha}}$, where $u\in[0,1]$ is the thrust throttling factor
and $\hat{\boldsymbol{\alpha}}$ is a unit vector representing the direction at
which thrust is applied. Functions $\mathbf{g}(\mathbf{r})$ and
$\mathbf{h}(\mathbf{v})$ are given by [21]

$\displaystyle\mathbf{g}(\mathbf{r})$
$\displaystyle=\begin{bmatrix}r_{x}-\frac{(1-\mu)(r_{x}+\mu)}{r_{1}^{3}}-\frac{\mu(r_{x}+\mu-1)}{r_{2}^{3}}\\\
r_{y}-\frac{(1-\mu)r_{y}}{r_{1}^{3}}-\frac{\mu r_{y}}{r_{2}^{3}}\\\
-\frac{(1-\mu)r_{z}}{r_{1}^{3}}-\frac{\mu r_{z}}{r_{2}^{3}}\end{bmatrix}$ (2)
$\displaystyle\mathbf{h}(\mathbf{v})$
$\displaystyle=\begin{bmatrix}2v_{y}&-2v_{x}&0\end{bmatrix}^{T}$ (3)

where subscripts $x$, $y$, and $z$ denote the components of the position and
velocity vector along the corresponding Cartesian direction of the synodic
reference frame and $\mu=m_{2}/(m_{1}+m_{2})$ is the non-dimensionalized mass
of the second primary body. Additionally, $r_{1}$ and $r_{2}$ denote the non-
dimensionalized distance of the spacecraft from the first and second primary
bodies respectively and are given by

$\displaystyle r_{1}$
$\displaystyle=\sqrt{(r_{x}+\mu)^{2}+r_{y}^{2}+r_{z}^{2}}$ (4) $\displaystyle
r_{2}$ $\displaystyle=\sqrt{(r_{x}+\mu-1)^{2}+r_{y}^{2}+r_{z}^{2}}$ (5)

###  2.2 Indirect Minimum-Fuel Trajectory Optimization using Energy-to-Fuel
Homotopy

For the minimum-fuel problem, it is desired to find a solution of Eq. (1)
which minimizes the cost function

$\displaystyle J_{mf}=\frac{T_{max}}{c}\int_{t_{i}}^{t_{f}}udt$ (6)

where $t_{i}$ and $t_{f}$ correspond to the initial and final time of the
trajectory. A known characteristic of minimum-fuel trajectories is the bang-
bang nature of their control profile [9, 20, 21], which results in
discontinuities in the throttling factor $u$ along the trajectory and can
severely restrict the convergence radius of numerical methods. To mitigate
this issue, Bertrand and Epenoy proposed the use of an energy-to-fuel homotopy
continuation for transforming solutions from the discontinuity-free minimum-
energy problem to the minimum-fuel problem [20]. A perturbed energy term is
introduced into the cost function by means of the perturbation parameter
$\epsilon$, which is gradually reduced until its effect is nulled. Employing
this approach, the new cost function is given by

$\displaystyle J=\frac{T_{max}}{c}\int_{t_{i}}^{t_{f}}\left[u-\epsilon
u(1-u)\right]dt\hskip 14.22636pt\epsilon\in[0,1]$ (7)

To begin the derivation of the TPBVP, we proceed by forming the Hamiltonian as
[21]

$\displaystyle H$
$\displaystyle=\boldsymbol{\lambda}^{T}\mathbf{f}(\mathbf{x},\hat{\boldsymbol{\alpha}},u)+\frac{T_{max}}{c}\left[u-\epsilon
u(1-u)\right]$ (8)
$\displaystyle=\boldsymbol{\lambda}_{r}^{T}\mathbf{v}+\boldsymbol{\lambda}_{v}^{T}\left[\mathbf{g}(\mathbf{r})+\mathbf{h}(\mathbf{v})+\frac{uT_{max}}{m}\hat{\boldsymbol{\alpha}}\right]-\lambda_{m}\frac{uT_{max}}{c}+\frac{T_{max}}{c}\left[u-\epsilon
u(1-u)\right]$ (9)

where
$\boldsymbol{\lambda}^{T}=\left[\boldsymbol{\lambda}_{r}^{T},\boldsymbol{\lambda}_{v}^{T},\lambda_{m}\right]$
is the introduced co-state vector, which evolves in time according to

$\displaystyle\dot{\boldsymbol{\lambda}}=\begin{bmatrix}\dot{\boldsymbol{\lambda}}_{r}\\\
\dot{\boldsymbol{\lambda}}_{v}\\\
\dot{\lambda}_{m}\end{bmatrix}=-\left(\frac{\partial
H}{\partial\mathbf{x}}\right)^{T}=\begin{bmatrix}-\mathbf{G}^{T}\boldsymbol{\lambda}_{v}\\\
-\boldsymbol{\lambda}_{r}-\mathbf{H}^{T}\boldsymbol{\lambda}_{v}\\\
\frac{uT_{max}}{m^{2}}\boldsymbol{\lambda}_{v}^{T}\hat{\boldsymbol{\alpha}}\end{bmatrix}$
(10) $\displaystyle\mathbf{G}$
$\displaystyle=\frac{\partial\mathbf{g}(\mathbf{r})}{\partial\mathbf{r}}$ (11)
$\displaystyle\mathbf{H}$
$\displaystyle=\frac{\partial\mathbf{h}(\mathbf{v})}{\partial\mathbf{v}}$ (12)

where the non-zero terms of $\mathbf{G}$ and $\mathbf{H}$ are

$\displaystyle G_{1,1}$
$\displaystyle=1-\frac{1-\mu}{r_{1}^{3}}+\frac{3(1-\mu)(x+\mu)^{2}}{r_{1}^{5}}-\frac{\mu}{r_{2}^{3}}+\frac{3\mu(x+\mu-1)^{2}}{r_{2}^{5}}$
(13) $\displaystyle G_{2,2}$
$\displaystyle=1-\frac{1-\mu}{r_{1}^{3}}+\frac{3(1-\mu)y^{2}}{r_{1}^{5}}-\frac{\mu}{r_{2}^{3}}+\frac{3\mu
y^{2}}{r_{2}^{5}}$ (14) $\displaystyle G_{3,3}$
$\displaystyle=-\frac{1-\mu}{r_{1}^{3}}+\frac{3(1-\mu)z^{2}}{r_{1}^{5}}-\frac{\mu}{r_{2}^{3}}+\frac{3\mu
z^{2}}{r_{2}^{5}}$ (15) $\displaystyle G_{1,2}$
$\displaystyle=G_{2,1}=\frac{3(1-\mu)(x+\mu)y}{r_{1}^{5}}+\frac{3\mu(x+\mu-1)y}{r_{2}^{5}}$
(16) $\displaystyle G_{1,3}$
$\displaystyle=G_{3,1}=\frac{3(1-\mu)(x+\mu)z}{r_{1}^{5}}+\frac{3\mu(x+\mu-1)z}{r_{2}^{5}}$
(17) $\displaystyle G_{2,3}$
$\displaystyle=G_{3,2}=\frac{3(1-\mu)yz}{r_{1}^{5}}+\frac{3\mu yz}{r_{2}^{5}}$
(18) $\displaystyle H_{1,2}$ $\displaystyle=-H_{2,1}=2$ (19)

Applying the weak form of PMP, the optimal thrust direction unit vector is
characterized by [3]

$\displaystyle\hat{\boldsymbol{\alpha}}^{*}\in\text{arg}\hskip
2.84526pt\underset{\hat{\boldsymbol{\alpha}}}{\text{min}}\hskip 2.84526ptH$
(20)

which, noting that $uT_{max}/m\geq 0$ and $\lambda_{v}\neq 0$ in general,
gives

$\displaystyle\hat{\boldsymbol{\alpha}}^{*}=-\frac{\boldsymbol{\lambda}_{v}}{\lambda_{v}}$
(21)

where $\lambda_{v}$ is the Euclidean norm of $\boldsymbol{\lambda}_{v}$, and
$-\boldsymbol{\lambda}_{v}$ is defined as the primer vector [24]. Substituting
this result into Eq. (9) and defining a switching function $S$ as [21]

$\displaystyle S=-\frac{c}{m}\lambda_{v}-\lambda_{m}+1,$ (22)

the Hamiltonian can now be written as

$\displaystyle
H=\boldsymbol{\lambda}_{r}^{T}\mathbf{v}+\boldsymbol{\lambda}_{v}^{T}\left[\mathbf{g}(\mathbf{r})+\mathbf{h}(\mathbf{v})\right]+\frac{uT_{max}}{c}\left(S-\epsilon+\epsilon
u\right)$ (23)

Again applying the weak form of PMP, the optimal thrust throttling factor is
given by[21]

$u^{*}=\left\\{\begin{array}[]{ccc}0&S>\epsilon\\\
\frac{\epsilon-S}{2\epsilon}&-\epsilon\leq S\leq\epsilon\\\
1&S<-\epsilon\end{array}\right.$ (24)

Note that the relations between the co-states and optimal control variables
$\hat{\boldsymbol{\alpha}}^{*}$ and $u^{*}$, given in Eqs. (21) and (24), are
equivalent to Lawden’s primer vector control law when $\epsilon=0$.111Lawden’s
primer vector control law is a well-known relationship between the co-states,
optimal throttling factor and thrust direction when the indirect optimal
control approach is employed to solve minimum-fuel problems [24].

To complete the formulation of the trajectory optimization problem, terminal
constraints placed on the physical state variables are given by

$\begin{array}[]{ccc}\mathbf{r}(t_{i})-\mathbf{r}_{i}=0&\mathbf{v}(t_{i})-\mathbf{v}_{i}=0&m(t_{i})-1=0\\\
\mathbf{r}(t_{f})-\mathbf{r}_{f}=0&\mathbf{v}(t_{f})-\mathbf{v}_{f}=0\end{array}$
(25)

where $\mathbf{r}_{i}$, $\mathbf{v}_{i}$, $\mathbf{r}_{f}$, and
$\mathbf{v}_{f}$ are the desired initial and final position and velocity
vectors, and the initial mass has been scaled to one. Noting that the final
mass of the spacecraft is unconstrained, a transversality condition provides
the last constraint, which is given by [25]

$\displaystyle\lambda_{m}(t_{f})=0$ (26)

Substituting Eq. (21) into Eq. (1) and combining with Eq. (10), a
14-dimensional system of first order ordinary differential equations is
obtained

$\displaystyle\dot{\mathbf{y}}=\begin{bmatrix}\dot{\mathbf{x}}\\\
\dot{\boldsymbol{\lambda}}\end{bmatrix}=\begin{bmatrix}\mathbf{v}\\\
\mathbf{g}(\mathbf{r})+\mathbf{h}(\mathbf{v})-\boldsymbol{\lambda}_{v}uT_{max}/(\lambda_{v}m)\\\
-uT_{max}/c\\\ -\mathbf{G}^{T}\boldsymbol{\lambda}_{v}\\\
-\boldsymbol{\lambda}_{r}-\mathbf{H}^{T}\boldsymbol{\lambda}_{v}\\\
-\lambda_{v}uT_{max}/m^{2}\end{bmatrix}$ (27)

which, when paired with the optimal throttling factor given in Eq. (24) and
the constraints given in Eqs. (25) and (26), fully define the TPBVP. The
solution of the TPBVP is the set of initial co-state variables which satisfy
the terminal constraints. Here, it is important to note that solutions of the
derived TPBVP only imply satisfaction of the first-order necessary conditions
of optimality [9]. Therefore, without further investigation of the second-
order sufficient conditions of optimality, solutions are only guaranteed to be
stationary or extremal solutions of the optimal control cost function given in
Eq. (7), also known as candidates of optimality [26]. Furthermore, one should
note that multiple local extrema may exist, resulting in multiple solutions of
the TPBVP.

##  3 TPBVP Solution Methodology

Solving the derived TPBVP requires a good guess of the initial co-state
variables $\boldsymbol{\lambda}(t_{i})$. The methodology proposed here
involves the use of the PSO algorithm to first initialize the co-state
variables for the minimum-energy problem, corresponding to the perturbation
parameter $\epsilon=1$ in Eq. (24). The PSO-initialized co-states are then
used to seed a single shooting procedure with homotopy continuation, where a
trust-region nonlinear solver [27] iteratively discovers the solution to the
TPBVP for decreasing values of the perturbation parameter until a solution
corresponding to $\epsilon=0$ is found.

In this Section, the methodology employed to integrate the state, co-state,
and State Transition Matrix (STM) differential equations while performing
switching detection is discussed. The proposed method for initializing the co-
state variables with PSO is also presented, along with an approach to perform
single shooting with a homotopy continuation. Along the way, important
considerations in regard to programming language selection and computational
efficiency are made.

###  3.1 Numerical Integration with Switching Detection

When integrating the state and co-state differential equations, in addition to
the STM differential equations when single shooting is performed, it is
important to ensure that a sufficiently high degree of accuracy is maintained,
so that convergence can be achieved. Due to the piecewise continuous function
which defines the optimal throttling factor (see Eq. (24)), integration error
can grow about the points at which switching of the throttling factor occurs,
and can severely restrict the convergence of a shooting algorithm. However,
this can be avoided if the points at which the switching occurs are explicitly
determined and numerically integrated up to. Additionally, both the PSO co-
state initialization and single shooting phases require numerous repeated
numerical integrations of the differential equations. Performing numerical
integration in a programming language which exhibits high computational
performance with a differential equation solving suite that provides high-
order Runge-Kutta algorithms and robust event detection is thus desirable. For
these reasons, the Julia [28] programming language and the Julia-implemented
differential equation solving package, DifferentialEquations.jl [29], were
chosen for all numerical computation and integration in this work.

Julia is a dynamically typed language with syntax similar to that of MATLAB or
Python, thereby lending itself to quick development. Historically, a known
characteristic of this category of languages is the slow execution time due to
the required use of an interpreter, an issue which Julia circumvents through
the use of a just-in-time (JIT) LLVM-based compiler. Therefore, although Julia
is a dynamic language with easy-to-read syntax, it allows producing software
with performance that rivals statically typed languages like C and FORTRAN.
Furthermore, DifferentialEquations.jl is highly optimized and provides a
plethora of modern differential equation solvers and features, including
robust event detection schemes, termed callbacks, which can employ high-order
interpolants to check whether an event has occurred at multiple points
backwards in time. Through the application of DifferentialEquations.jl and its
continuous callback feature (i.e. continuous event detection), numerical
integration can be performed efficiently by exploiting an adaptive time-
stepping method while robustly detecting throttling factor switching to within
64-bit floating point number precision.

###  3.2 Co-State Initialization with PSO

Beginning with the inception of the PSO algorithm detailed in the work of
Kennedy and Eberhart [30], many have proposed modifications with improved
convergence characteristics [31, 32]. The term PSO has thus evolved to
encompass a class of heuristic global optimization algorithms rather than any
one strictly defined algorithm. For the sake of reproducibility, in this work,
an implementation of PSO in the Julia programming language is employed, which
is based on the MATLAB PSO algorithm [33].

To employ PSO to initialize the co-state variables for the minimum-energy
problem ($\epsilon=1$), a nonlinear, unconstrained optimization problem is
formulated as

$\displaystyle\underset{\boldsymbol{\lambda}(t_{i})}{\text{Minimize}}\hskip
8.53581ptJ_{PSO}(\boldsymbol{\lambda}(t_{i}))=\mathbf{e}^{T}(t_{f})\mathbf{W}\mathbf{e}(t_{f})$
(28)

where $\mathbf{W}\in\mathbb{R}^{7\times 7}$ is a diagonal weighting matrix and
$\mathbf{e}(t_{f})$ is the final time boundary condition residual given by

$\displaystyle\mathbf{e}(t_{f})=\begin{bmatrix}\mathbf{r}(t_{f})-\mathbf{r}_{f}\\\
\mathbf{v}(t_{f})-\mathbf{v}_{f}\\\ \lambda_{m}(t_{f})\end{bmatrix}$ (29)

which is computed by integrating the state and co-state differential equations
from $t_{i}$ to $t_{f}$ as described before, starting from an initial state
vector constructed as a concatenation of the constrained initial state and the
current guess for the co-state values at the initial time

$\displaystyle\mathbf{y}(t_{i})=\begin{bmatrix}\mathbf{r}_{i}^{T}&\mathbf{v}_{i}^{T}&1.0&\boldsymbol{\lambda}^{T}(t_{i})\end{bmatrix}^{T}$
(30)

Through computational investigation, the implemented single shouting routine,
discussed in the following Section, was found to be most sensitive to large
residuals in the final position. Therefore, the weighting matrix was defined
as

$\displaystyle\mathbf{W}=\text{Diag}\left(\begin{bmatrix}10&10&10&1&1&1&1\end{bmatrix}\right)$
(31)

such that the squared position residuals are weighted by a factor of ten,
while the remaining weights are left at unity.

Due to the stochastic nature of the PSO algorithm, some particles may travel
to positions that result in trajectories that pass below the surface of a
primary body, which is physically non-realizable. Therefore, the integration
routine was set up to halt the integration of a given trajectory if the
distance between the spacecraft and the surface of a primary body was smaller
than zero at any time.

When employing PSO for co-state initialization, especially when a large swarm
size is required and/or the cost function is expensive to evaluate, it is
important to take advantage of the parallelizable nature of the algorithm to
reduce compute time. In fact, for each iteration of PSO, the cost function
must be re-evaluated for each particle. In general, the evaluation of the cost
function for any given particle is independent of all other particles in the
swarm, which makes the process well suited for shared or distributed memory
parallelism, which was taken advantage of in this work.

###  3.3 Single Shooting with Continuation

After co-state initialization with PSO, a single shooting routine based on a
trust-region nonlinear solver provided by NLsolve.jl [27] is employed for
solving the TPBVP. To do so, the nonlinear 7-dimensional vector function which
is desired to be solved is defined as

$\displaystyle\mathbf{Z}(\boldsymbol{\lambda}(t_{i}),\epsilon)=\mathbf{e}(t_{f})$
(32)

where a solution corresponds to $\mathbf{e}(t_{f})=\mathbf{0}_{7\times 1}$. In
other words, it is desired to find the initial co-state vector
$\boldsymbol{\lambda}(t_{i})$ which, for a given value of $\epsilon\in[0,1]$ –
integrating Eq. (27) from $t_{i}$ to $t_{f}$ as described before – results in
a final time boundary condition residual of
$\mathbf{e}(t_{f})=\mathbf{0}_{7\times 1}$. Henceforth,
$\mathbf{Z}(\boldsymbol{\lambda}(t_{i}),\epsilon)$ in Eq. (32) is referred to
as the shooting function. Seeded by the PSO-initialized co-states, the trust-
region algorithm is used to iteratively solve Eq. (32) for values of the
perturbation parameter decreasing from one to zero according to the
continuation law [21]

$\displaystyle\epsilon_{j}=\frac{j^{2}-1}{N^{2}-1}\hskip
14.22636ptj=N,N-1,\dots,2,1$ (33)

where $N$ is the number of iterations that are desired to be performed for
transitioning the perturbation parameter from one to zero. This value can be
chosen somewhat arbitrarily, but it is important to consider that higher
values of $N$ result in smaller changes in the perturbation parameter at each
iteration, providing a smoother transition from the minimum-energy to the
minimum-fuel cost function, albeit at a higher computational cost. Throughout
this work, a value of $N=25$ was chosen to balance numerical stability and
computational efficiency.

It is important to note that the trust-region algorithm requires that the
Jacobian of the shooting function, i.e.,
$\partial\mathbf{Z}(\boldsymbol{\lambda}(t_{i}),\epsilon)/\partial\boldsymbol{\lambda}(t_{i})$,
also be computed. Multiple methods exist to compute this Jacobian (e.g.,
automatic differentiation and finite difference), among which the STM was
selected for this work. The STM is defined as

$\displaystyle\boldsymbol{\Phi}(t_{f},t_{i})=\frac{\partial\mathbf{y}(t_{f})}{\partial\mathbf{y}(t_{i})}$
(34)

To compute the STM, the $14^{2}$ STM differential equations must also be
integrated along with the state and co-state differential equations in Eq.
(27) according to

$\displaystyle\dot{\boldsymbol{\Phi}}(t,t_{i})=\mathbf{F}\boldsymbol{\Phi}(t,t_{i})=\frac{\partial\dot{\mathbf{y}}}{\partial\mathbf{y}}\Big{|}_{t}\boldsymbol{\Phi}(t,t_{i}),\hskip
14.22636pt\boldsymbol{\Phi}(t_{i},t_{i})=\mathbf{I}_{14\times 14}$ (35)

where $\mathbf{F}$ is the Jacobian of Eq. (27) evaluated at time $t$, the
analytical expression of which is provided by Zhang et al. [21]. Additionally,
when discontinuities are present in the shooting function, as is the case for
$\epsilon=0$ when switching occurs (see Eq. (24)), the STM across the
discontinuity must also be computed as [9]

$\displaystyle\boldsymbol{\Psi}(t_{n})=\frac{\partial\mathbf{y}(t_{n}^{+})}{\partial\mathbf{y}(t_{n}^{-})}=\mathbf{I}_{14\times
14}+\left(\dot{\mathbf{y}}\big{|}_{t_{n}^{+}}-\dot{\mathbf{y}}\big{|}_{t_{n}^{-}}\right)\left(\frac{\partial
S}{\partial\mathbf{y}}\frac{1}{\dot{S}}\right)$ (36)

where $t_{n}^{-}$ and $t_{n}^{+}$ represent the time immediately before and
after the discontinuity occurred, respectively, and

$\displaystyle\frac{\partial S}{\partial\mathbf{y}}$
$\displaystyle=\begin{bmatrix}\mathbf{0}_{1\times
6}&\frac{c}{m^{2}}\lambda_{v}&\mathbf{0}_{1\times
3}&-\frac{c}{m}\frac{\boldsymbol{\lambda}_{v}^{T}}{\lambda_{v}}&-1\end{bmatrix}$
(37) $\displaystyle\dot{S}$
$\displaystyle=\frac{c}{m}\left(\boldsymbol{\lambda}_{r}+\mathbf{H}^{T}\boldsymbol{\lambda}_{v}\right)^{T}\frac{\boldsymbol{\lambda}_{v}}{\lambda_{v}}$
(38)

For a trajectory with $M$ discontinuities between $t_{i}$ and $t_{f}$, the
state transition matrix is then given by [21]

$\displaystyle\boldsymbol{\Phi}(t_{f},t_{i})=\boldsymbol{\Phi}(t_{f},t_{M}^{+})\boldsymbol{\Psi}(t_{M})\boldsymbol{\Phi}(t_{M}^{-},t_{M-1}^{+})\boldsymbol{\Psi}(t_{M-1})\dots\boldsymbol{\Phi}(t_{2}^{-},t_{1}^{+})\boldsymbol{\Psi}(t_{1})\boldsymbol{\Phi}(t_{1}^{-},t_{i})$
(39)

Once the STM has been propagated from $t_{i}$ to $t_{f}$ along with the state
and co-state variables, the Jacobian of the shooting function is constructed
from components of the STM, such that

$\displaystyle\frac{\partial\mathbf{Z}(\boldsymbol{\lambda}(t_{i}),\epsilon)}{\partial\boldsymbol{\lambda}(t_{i})}=\begin{bmatrix}\frac{\partial\mathbf{r}(t_{f})}{\partial\boldsymbol{\lambda}_{r}(t_{i})}&\frac{\partial\mathbf{r}(t_{f})}{\partial\boldsymbol{\lambda}_{v}(t_{i})}&\frac{\partial\mathbf{r}(t_{f})}{\partial\lambda_{m}(t_{i})}\\\
\frac{\partial\mathbf{v}(t_{f})}{\partial\boldsymbol{\lambda}_{r}(t_{i})}&\frac{\partial\mathbf{v}(t_{f})}{\partial\boldsymbol{\lambda}_{v}(t_{i})}&\frac{\partial\mathbf{v}(t_{f})}{\partial\lambda_{m}(t_{i})}\\\
\frac{\partial\lambda_{m}(t_{f})}{\partial\boldsymbol{\lambda}_{r}(t_{i})}&\frac{\partial\lambda_{m}(t_{f})}{\partial\boldsymbol{\lambda}_{v}(t_{i})}&\frac{\partial\lambda_{m}(t_{f})}{\partial\lambda_{m}(t_{i})}\\\
\end{bmatrix}$ (40)

A diagram describing the algorithm from co-state initialization to minimum-
fuel solution is provided in the Appendix.

##  4 Results

Table 1: Constant Parameters Constant | Value | Units
---|---|---
Mass Parameter $\mu$ | $1.21506038\times 10^{-2}$ | -
Grav. Constant $g_{0}$ | 9.81 | m/$\text{s}^{2}$
Time Unit (TU) | $3.75162997\times 10^{5}$ | s
Length Unit (LU) | $3.84400000\times 10^{5}$ | km
Velocity Unit (VU) | $1.02462131$ | km/s
Mass Unit (MU) | $m_{i}$ | kg

In this Section, case studies are presented to validate the methodology,
analyze the performance of the proposed method for co-state initialization, as
well as to assess the convergence of the PSO-initialized co-states to
solutions of the minimum-fuel TPBVP for varying particle swarm sizes. The
proposed methodology is applied to determine low-thrust minimum-fuel transfers
in the following two scenarios in the Earth-Moon system: 1) from a
geostationary transfer orbit (GTO) to an $L_{1}$ halo orbit and 2) from an
$L_{2}$ halo orbit to a $L_{1}$ halo orbit. Constant parameters used for both
scenarios, including the time, length, speed, and mass units employed in the
non-dimensionalization of the CR3BP dynamics, are provided in Table 1. Note
that the mass unit (which is a unique value for a given scenario) is defined
as the initial mass of the spacecraft $m_{i}$. The chosen scenarios provide
problems of moderate difficulty to benchmark our approach while allowing for
validation and comparison to previous works [21, 22, 34, 35, 36, 37].

Due to the stochastic nature of PSO, 100 co-state initialization trials are
performed for swarm sizes ranging from 250 to 5000 particles in order for
statistical conclusions to be made about the performance of the algorithm. The
position of the $k$-th particle in the swarm is defined in the 7-dimensional
co-state space as

$\displaystyle\boldsymbol{\lambda}_{k}(t_{i})=\begin{bmatrix}\boldsymbol{\lambda}_{r,k}^{T}(t_{i})&\boldsymbol{\lambda}_{v,k}^{T}(t_{i})&\lambda_{m,k}(t_{i})\end{bmatrix}^{T}$
(41)

The swarm is initialized with particle positions uniformly distributed between
upper and lower bounds such that $\boldsymbol{\lambda}_{r}(t_{i})\in[-40,40]$,
$\boldsymbol{\lambda}_{v}(t_{i})\in[-2,2]$, and $\lambda_{m}(t_{i})\in[0,2]$.
It is important to note that these are not constraints placed on the positions
of particles, but are rather bounds used to define the distribution which
starting particle positions are sampled from. Once the swarm is initialized,
particles are then free to travel anywhere in the search space, which is
restricted to $\boldsymbol{\lambda}_{r}(t_{i})\in[-100,100]$,
$\boldsymbol{\lambda}_{v}(t_{i})\in[-10,10]$, and
$\lambda_{m}(t_{i})\in[0,10]$. Additionally, convergence criteria are chosen
such that a PSO trial is halted if the objective function value does not
improve by more than $1\times 10^{-6}$ over the course of 50 iterations, or if
a maximum run time of 30 minutes is met. Throughout the following case
studies, the effect of different PSO swarm sizes only is analyzed for the
discussion to remain tractable. The remaining PSO parameters are held at their
default values (as defined by MathWorks [33]), with the exception of the
minimum adaptive neighborhood size fraction, which is set to 0.05.

For all numerical integration performed throughout the case studies, Verner’s
“most efficient” Runge-Kutta 9(8) method was employed, with an absolute and
relative tolerance of $10^{-14}$, along with its corresponding $9$-th order
interpolant for switching detection [38]. All computations were performed on a
workstation PC with specifications relevant to this work provided in Table 2.
At each iteration of PSO, the evaluations of the cost function for each
particle were computed in parallel using shared memory parallelism on 32 of
the available 64 physical cores within the CPU.

Table 2: PC Specifications Specification | Value
---|---
CPU | AMD Ryzen Threadripper Pro 3995WX
RAM | 32 GB DDR4 1600 MHz
Operating System | Windows 10 Enterprise
Julia | Version 1.7.0

###  4.1 Scenario 1: GTO to L1 Halo Orbit

For the GTO to $L_{1}$ halo orbit transfer problem, the spacecraft was chosen
to have an initial mass $m_{i}=1500$ kg, a maximum thrust $T_{max}=10$ N, and
a specific impulse $I_{sp}=3000$ s, as in reference [21]. In the minimum-fuel
problem, the time of flight (TOF) was fixed at 8.6404 days. The terminal
physical boundary conditions for the TPBVP, which identify the initial
position of the spacecraft and the final point targeted on the halo orbit,
were defined as [21]

$\displaystyle\mathbf{r}_{i}$
$\displaystyle=\begin{bmatrix}-0.0194885115&-0.0160334798&0.0\end{bmatrix}^{T}\text{
LU}$ $\displaystyle\mathbf{v}_{i}$
$\displaystyle=\begin{bmatrix}8.9188819237&-4.0817936888&0.0\end{bmatrix}^{T}\text{
VU}$ $\displaystyle\mathbf{r}_{f}$
$\displaystyle=\begin{bmatrix}0.8233851820&0.0&-0.0222775563\end{bmatrix}^{T}\text{
LU}$ $\displaystyle\mathbf{v}_{f}$
$\displaystyle=\begin{bmatrix}0.0&0.1341841703&0.0\end{bmatrix}^{T}\text{ VU}$

such that the GTO had periapsis and apoapsis altitudes of $h_{p}=400$ km and
$h_{a}=35864$ km, respectively, and the $L_{1}$ halo orbit had an out-of-plane
amplitude of 8000 km.

(a) Trajectory A

(b) Trajectory B

(c) Trajectory C

(d) Trajectory D

(e) Trajectory E

(f) Trajectory F

Figure 1: Minimum-Fuel Trajectories for Scenario 1 Table 3: Minimum-Fuel TPBVP Solutions for Scenario 1 Trajectory | $\boldsymbol{\lambda}_{r}^{T}(t_{i})$ | $\boldsymbol{\lambda}_{v}^{T}(t_{i})$ | $\lambda_{m}(t_{i})$ |  $\Delta m$ (kg) | # Revs.
---|---|---|---|---|---
A | $\left[23.2524,\hskip 0.28453pt50.6272,-0.08489\right]$ | $\left[-0.1546,\hskip 0.28453pt0.0706,\hskip 0.28453pt-0.0002\right]$ | $0.1385$ | 140.0 | 4
B | $\left[15.6850,\hskip 0.28453pt33.0313,\hskip 0.28453pt-0.0938\right]$ | $\left[-0.1020,\hskip 0.28453pt0.0450,\hskip 0.28453pt-0.0002\right]$ | 0.1334 | 134.4 | 5
C | $\left[7.7397,\hskip 0.28453pt14.5669,\hskip 0.28453pt-0.1239\right]$ | $\left[-0.0466,\hskip 0.28453pt0.0180,\hskip 0.28453pt-0.0001\right]$ | 0.1535 | 139.1 | 6
D | $\left[1.1404,\hskip 0.28453pt-0.6176,\hskip 0.28453pt-0.1535\right]$ | $\left[-0.0006,\hskip 0.28453pt-0.0042,\hskip 0.28453pt-0.0002\right]$ | 0.1935 | 153.5 | 7
E | $\left[-9.4937,\hskip 0.28453pt-25.2095,\hskip 0.28453pt-0.2092\right]$ | $\left[0.0738,\hskip 0.28453pt-0.0401,\hskip 0.28453pt-0.0002\right]$ | 0.2916 | 178.5 | 8
F | $\left[-24.3124,\hskip 0.28453pt-60.4677,\hskip 0.28453pt-0.2883\right]$ | $\left[0.1801,\hskip 0.28453pt-0.0928,\hskip 0.28453pt-0.0002\right]$ | 0.5613 | 219.7 | 9

####  4.1.1 Minimum-Fuel TPBVP Solutions

As discussed before, multiple local extrema of the optimal control cost
function, and therefore multiple solutions of the TPBVP may exist. A total of
six solutions of the minimum-fuel TPBVP (corresponding to $\epsilon=0$) were
found. These are labeled Trajectories A through F, and are shown in Figure 1,
with thrusting arcs depicted in red and coasting arcs depicted in blue. It is
clear that each of the trajectories follows different paths and thrusting and
coasting arcs are placed at different locations along the transfers. It is
observed that Trajectories A, E, and F (in Figures 1(a), 1(e), and 1(f))
exhibit much shorter coasting arcs during either the initial spiraling about
the Earth or the final leg of the trajectory when approaching the $L_{1}$ halo
orbit, compared to the remaining trajectories. It is also observed that each
trajectory exhibits a unique number of revolutions about the Earth before the
final transfer to the halo orbit.

Table 3 displays the initial co-state variables corresponding to each of the
solutions, along with the total change in mass, i.e. fuel required, which is
computed as $\Delta m=[1-m(t_{f})]\times 1500$ kg. As can be seen by observing
the $\Delta m$ column of Table 3, Trajectory B corresponds to the most optimal
minimum-fuel trajectory found, while the remaining solutions correspond to
other extremal solutions of the minimum-fuel cost function. The total number
of revolutions about the Earth is shown in the final column of Table 3. The
unique number of revolutions for each solution suggests that each trajectory
may be locally optimal for a given number of revolutions about the Earth. It
is noted that the initial co-states corresponding to Trajectory B, along with
$\Delta m$, align with a previously found solution by Zhang et al. [21].

####  4.1.2 Co-State Initialization

In this Section, the performance of the proposed method for converging to
solutions of the derived TPBVP which satisfy all boundary conditions is
analyzed. For the purpose of this discussion, a trial is considered to have
converged if the $\infty$-norm of the TPBVP final boundary condition residual
is reduced to less than $1\times 10^{-10}$, i.e.,
$||\mathbf{e}(t_{f})||_{\infty}=\max\\{|e_{i}(t_{f})|:i=1,2,\dots,7\\}<1\times
10^{-10}$.

Table 4: Scenario 1 Convergence Characteristics for Different Swarm Sizes Convergence Characteristic | Swarm Size
---|---
| 250 | 500 | 750 | 1000 | 2000 | 3000 | 4000 | 5000
Min. Fuel % Converged | 22 | 42 | 41 | 50 | 63 | 58 | 64 | 54
Avg. Final Obj. Func. | 7.93 | 5.74 | 5.10 | 4.28 | 3.84 | 3.50 | 2.45 | 2.97
Avg. Time to Converge (min) | 2.33 | 7.75 | 13.48 | 15.66 | 22.88 | 25.91 | 28.04 | 28.84
Avg. Function Evaluations | 3.02E5 | 9.24E5 | 1.56E6 | 2.00E6 | 2.82E6 | 3.06E6 | 3.71E6 | 3.58E6

For each of the investigated swarm sizes, Table 4 displays the percentage of
trials for which the initialized co-states converged to a solution of the
minimum-fuel ($\epsilon=0$) TPBVP, the average final value of the PSO
objective function (see Eq. (28)), and the average time and objective function
evaluations before PSO convergence. The percentage of converged trials was
found to be lowest for the smallest number of particles and to increase with
swarm size, reaching a maximum of 64% converged trials with a swarm size of
4000 particles, before reducing slightly for the largest swarm size. A similar
trend of improvement was observed in the average value of the final PSO
objective function and the number of function evaluations. The average time to
converge increases with the swarm size, approaching the maximum time allocated
to PSO (i.e., 30 minutes per trial) for larger swarm sizes.

These results suggest that for swarm sizes smaller than 1000 particles, early
stagnation of objective function improvement (due to the limited number of
particles exploring the solutions space) results in premature convergence of
the PSO algorithm. For larger swarm sizes, ranging from 2000 to 4000
particles, PSO does not suffer from premature convergence and the ability of
the algorithm to produce guesses for the initial co-states which converge to
solutions of the TPBVP is improved greatly. These larger swarm sizes are more
effective at exploring the solution space and are able to reduce the PSO
objective function further, before improvement stalls and the algorithm
converges. It is also important to note that these larger swarm sizes do come
with an additional computational cost: the objective function must be
evaluated many more times per PSO iteration resulting in longer periods spent
waiting for a guess to be generated or prematurely halting optimization
depending on the maximum time limit imposed. For a swarm size of 5000
particles, we begin to see the effect of this computational cost: the average
time to converge is nearly at the maximum time allocated to PSO per trial,
indicating that in a majority of cases PSO did not converge and was instead
halted after 30 minutes of run time (i.e., the PSO objective function was
still improving when the algorithm was halted). It is expected that this trend
would continue for even larger swarm sizes unless more time was allocated for
PSO co-state initialization.

####  4.1.3 Illustration and Analysis of the Solution Generation Process

Within this Section, the process of generating a solution of the minimum-fuel
TPBVP, from co-state initialization to single shooting with homotopy
continuation, is illustrated and analyzed. The discussion is focused on a
single trial that converged to the most fuel-optimal solution found throughout
the Scenario 1 case study, i.e., Trajectory B.

(a) After 5 iterations

(b) At convergence

Figure 2: Distribution of PSO particles Figure 3: PSO objective function
versus iteration number

By way of example, Figure 2 displays the distribution of 2000 particles (each
represented by a blue dot) in the $x$-$y$ velocity co-state space after 5 PSO
iterations and at PSO convergence, along with the location of the PSO-
generated guess (represented by a yellow dot) and the minimum-energy solution
which the generated guess converged to after performing single shooting
(represented by a red dot). Additionally, Figure 3 displays the PSO objective
function value as the number of PSO iterations increases. Observing Figure
2(a), it can be seen that the particles quickly distributed themselves across
the space from their initial placement of
$\lambda_{v_{x}},\lambda_{v_{y}}\in[-2,2]$. The particles then spent the next
few hundred iterations traveling across the full search space, rapidly
decreasing the PSO objective function as is shown in Figure 3. After
approximately 200 iterations, once significantly better solutions (i.e., co-
states which resulted in a significantly smaller objective function value)
were discovered, the rate at which the objective function improved is seen to
decrease significantly and the swarm began to contract, spending more time
searching the area near the better-known solutions. This slow improvement of
the objective function continued until one of the stopping criteria of PSO was
met. From Figure 2(b), it can be observed that, upon halting of PSO, the
particles had traveled to lie within a small region, which happens to be
centered about the generated guess. It can also be seen that the generated
guess lies near the solution to the minimum-energy problem. This confirms the
capability of the algorithm of generating a guess for the co-state variables
which lies near a solution of the minimum-energy TPBVP.

Figure 4: Trajectory B evolution through homotopy continuation Figure 5:
Switching function and throttling factor through Trajectory B continuation

After PSO-based co-state initialization, the next step in the solution
generation process requires employing single shooting to first converge to a
solution of the minimum-energy TPBVP (i.e., with $\epsilon=1$) starting from
the PSO-generated guess, before repeatedly resolving the TPBVP while
decreasing $\epsilon$ according the Eq. (33). Figure 4 displays the evolution
of Trajectory B through the continuation process as it is transitioned from
minimum-energy to minimum-fuel with $\epsilon=1.0,0.72,0.49,0.25$ and 0.0. One
should note that 25 steps in total are taken during the continuation process,
and only 5 are selected here for the sake of illustration. Also, Figure 5
displays the switching function and throttling factor in time for each of the
trajectories displayed in Figure 4. As can be seen in Figure 4, all
trajectories start from and arrive at the same positions as expected. The
minimum-energy and minimum-fuel trajectories follow nearly identical paths
until the final spiral and transfer to the halo orbit. However, it is observed
that coasting arcs appear to increase as the continuation approaches a
minimum-fuel trajectory. Figure 5 illustrates the gradual introduction of the
bang-bang throttling discontinuity throughout the continuation, which is the
core strength of the energy-to-fuel homotopy, which widens the convergence
radius and improves the probability of PSO-generated guesses producing a
solution of the minimum-fuel cost function.

###  4.2 Scenario 2: L2 Halo Orbit to L1 Halo Orbit

For the $L_{2}$ halo orbit to $L_{1}$ halo orbit transfer problem, the
spacecraft was chosen to have an initial mass $m_{i}=2000$ kg, a maximum
thrust $T_{max}=1.5$ N, and a specific impulse $I_{sp}=2000$ s, as in
reference [22]. The transfer duration was fixed to 12.7 days and the terminal
physical boundary conditions for the TPBVP, identifying the initial state of
the spacecraft on an $L_{2}$ halo orbit and the target point on an $L_{1}$
halo orbit, were defined as [22]

$\displaystyle\mathbf{r}_{i}$
$\displaystyle=\begin{bmatrix}1.1599795702248494&0.009720428035815552&-0.12401864915284157\end{bmatrix}^{T}\text{
LU}$ $\displaystyle\mathbf{v}_{i}$
$\displaystyle=\begin{bmatrix}0.008477705130550553&-0.20786307954141953&-0.010841912833115475\end{bmatrix}^{T}\text{
VU}$ $\displaystyle\mathbf{r}_{f}$
$\displaystyle=\begin{bmatrix}0.8484736688482315&0.00506488863463682&0.17343680487577373\end{bmatrix}^{T}\text{
LU}$ $\displaystyle\mathbf{v}_{f}$
$\displaystyle=\begin{bmatrix}0.005241131023638693&0.26343491250951045&-0.008541420325316247\end{bmatrix}^{T}\text{
VU}$

These correspond to initial and target halo orbits with periods of 14.2 and
11.2 days respectively.

(a) Trajectory $\alpha$ three-dimensional view

(b) Trajectory $\alpha$ planar projection

(c) Trajectory $\beta$ three-dimensional view

(d) Trajectory $\beta$ planar projection

(e) Trajectory $\gamma$ three-dimensional view

(f) Trajectory $\gamma$ planar projection

Figure 6: Minimum-Fuel Trajectories for Scenario 2 Table 5: Minimum-Fuel TPBVP Solutions Trajectory | $\boldsymbol{\lambda}_{r}^{T}(t_{i})$ | $\boldsymbol{\lambda}_{v}^{T}(t_{i})$ | $\lambda_{m}(t_{i})$ |  $\Delta m$ (kg)
---|---|---|---|---
$\alpha$ | $\left[0.12603,-0.07665,-0.05635\right]$ | $\left[0.03999,-0.00518,-0.06410\right]$ | $0.02236$ | 35.34
$\beta$ | $\left[-0.01486,0.01215,-0.07936\right]$ | $\left[0.01015,0.04457,0.01256\right]$ | $0.07632$ | 81.28
$\gamma$ | $\left[-0.02195,0.00659,0.07490\right]$ | $\left[-0.04314,0.03615,0.03842\right]$ | $0.03489$ | 61.27

####  4.2.1 Minimum-Fuel TPBVP Solutions

Similarly to Scenario 1, multiple solutions of the TPBVP were discovered by
applying the proposed method to solve the $L_{2}$ halo orbit to $L_{1}$ halo
orbit transfer problem. Figure 6 displays the resultant trajectories, both in
a three-dimensional view and projected onto the $x$-$y$ plane of the synodic
reference frame with the origin shifted to lie at the center of the Moon. It
is evident that the three discovered solutions differ drastically. Both
Trajectory $\alpha$ and $\beta$ (shown in Figures 6(a), 6(b), 6(c), and 6(d))
appear to take advantage of the lunar gravitational attraction with a close
approach of Moon while transferring between halo orbits. Whereas Trajectory
$\alpha$ efficiently performs the transfer, visibly requiring the least amount
of thrusting among the found solutions, traversing Trajectory $\beta$ requires
thrust to be applied for nearly the full duration of the transfer. Employing
an entirely different strategy, Trajectory $\gamma$ (shown in Figures 6(e) and
6(f)) does not use a close approach to the lunar surface, and instead thrusts
away from the Moon, then follows a coasting arc placed at the greatest
distance from the Moon, before thrusting for the final leg of the transfer to
insert into the $L_{1}$ halo orbit.

Table 5 reports the initial co-state variables corresponding to each of the
solutions, along with the fuel required to traverse each trajectory. As
expected from the previous discussion, Trajectory $\alpha$ is the most fuel-
efficient, requiring just 35.34 kg of fuel, while Trajectory $\gamma$ requires
nearly double the fuel mass (i.e., 61.27 kg), and Trajectory $\beta$ is the
most fuel expensive (at 81.28 kg). It should be noted that the fuel-optimal
trajectory discovered through the application of the proposed method (i.e.,
Trajectory $\alpha$) agrees with results published by Aziz et al [22]. A
slightly lower fuel requirement is found here, which is explained by the fact
that the Hybrid Differential Dynamic Programming (HDDP) algorithm applied by
Aziz et al. enforced a constant direction of thrust over each integration step
and can lag in switching thrust on or off [22], whereas indirect methods allow
thrust to be applied in the optimal direction of the primer vector at every
point along the trajectory and enforce throttling on or off at the optimal
times.

####  4.2.2 Co-State Initialization

Table 6: Convergence Characteristics for Different Swarm Sizes in Scenario 2 Convergence Characteristic | Swarm Size
---|---
| 250 | 500 | 750 | 1000 | 2000 | 3000 | 4000 | 5000
Min. Fuel % Converged | 80 | 86 | 79 | 81 | 67 | 79 | 73 | 75
Avg. Final Obj. Func. ($10^{-3}$) | 14.6 | 8.89 | 9.57 | 7.37 | 9.33 | 7.56 | 7.61 | 8.30
Avg. Time to Converge (min) | 2.30 | 3.76 | 4.77 | 6.20 | 10.46 | 11.28 | 12.69 | 14.61
Avg. Function Evaluations | 7.78E5 | 1.31E6 | 1.58E6 | 2.02E6 | 3.21E6 | 3.25E6 | 3.24E6 | 3.43E6

The performance of the proposed method at converging to solutions of the
minimum-fuel $L2$ halo orbit to $L1$ halo orbit transfer TPBVP is analyzed in
the following. For each of the investigated swarm sizes, Table 6 displays the
percentage of trials that converge to a solution of the minimum-fuel TPBVP,
along with the final PSO objective function value, time to converge, and the
number of function evaluations averaged over the 100 trials.

Overall, the proposed methodology performs more favorably at generating a
guess of the initial co-state values for the $L_{2}$ to $L_{1}$ transfer
problem, compared to Scenario 1. In this case, the best performance was
observed when a swarm size of 500 particles was employed, with 86% of trials
converging to a solution of the TPBVP. A swarm size of 250 particles performed
nearly as well, with 80% of trials converging to a TPBVP solution while only
requiring 2.3 minutes on average to generate the guess. Swarm sizes increasing
up to 1000 particles continued to produce quality guesses frequently,
converging in 79% and 81% of trials for 750 or 1000 particles respectively,
albeit at a higher computational cost. Interestingly, for a swarm size greater
than 1000 particles (aside from 3000 particles), co-state initialization
performance appeared to taper off, whereas the average time to converge in all
cases was well below the time limit of 30 minutes, indicating that allocation
of time for co-state initialization was not the issue. This contradicts
results from Scenario 1, which suggested that a larger swarm size should be
preferred as long as enough time is allocated for co-state initialization.
Furthermore, on average, larger swarm sizes successfully reduced the PSO
objective function further compared to the smaller, more often successful
swarm sizes. This is counter-intuitive, as a smaller average final objective
function value indicates that the generated co-state guesses are closer to
satisfying the final time boundary condition constraints and are therefore
expected to converge to a solution of the TPBVP with a higher frequency, as
was observed in Scenario 1.

It was discovered that multiple local minima of the PSO cost function exist
for the problem investigated in Scenario 2. The necessary condition for a
minimum of the PSO cost function is given by

$\left(\frac{\partial
J_{PSO}}{\partial\boldsymbol{\lambda}(t_{i})}\right)^{T}=2\left(\frac{\partial\mathbf{e}(t_{f})}{\partial\boldsymbol{\lambda}(t_{i})}\right)^{T}\mathbf{W}\mathbf{e}(t_{f})=\mathbf{0}_{7\times
1}$ (42)

which is satisfied either 1) when the final boundary condition residual
$\mathbf{e}(t_{f})$ is reduced to the zero vector or 2) when
$\left(\partial\mathbf{e}(t_{f})/\partial\boldsymbol{\lambda}(t_{i})\right)^{T}$
becomes rank deficient and $\mathbf{W}\mathbf{e}(t_{f})$ lies in its null
space. Clearly, global minima of the PSO cost function correspond to a final
boundary condition residual of zero, and local minima of the PSO cost function
exist in regions where the Jacobian of the final boundary condition residual
with respect to the initial co-state vector is rank deficient.

It is expected that the counter-intuitive phenomenon of larger PSO swarm sizes
resulting in lower convergence rates observed in Scenario 2 is due to PSO
converging to regions near local minima of the PSO cost function more
frequently. In trials when PSO produced a guess which did not converge to a
solution through single shooting, the PSO-generated guess was confirmed to lie
in a region where
$\partial\mathbf{e}(t_{f})/\partial\boldsymbol{\lambda}(t_{i})$, and therefore
the Jacobian of the shooting function with respect to the initial co-states
(see Eq. (40)), is rank deficient. This rank deficiency then results in
updates of Newton’s method in the trust-region algorithm which produce little
to no improvement during single shooting and therefore fail to converge to a
solution of the TPBVP.

(a) Three-dimensional view

(b) Planar projection

Figure 7: Initial guess and stalled solution for minimum-energy trajectory

An example of a trial that encountered this difficulty is shown in Figure 7.
Represented with a dashed line is the trajectory corresponding to the PSO-
generated guess of

$\boldsymbol{\lambda}_{PSO}(t_{i})=\begin{bmatrix}-0.75848&0.28231&-9.76567&-1.27696&2.08707&-3.41462&4.57662\end{bmatrix}^{T}$

After performing single shooting and failing to converge, the trajectory shown
as a solid red line corresponds to the values of the initial co-states at
which single shooting stalled, which are given by

$\boldsymbol{\lambda}(t_{i})=\begin{bmatrix}-0.15992&-0.23058&-5.42275&-0.57888,0.97406&-1.98589&2.11896\end{bmatrix}^{T}$

It is important to note that both trajectories are entirely shown in red
because no coasting arcs were exploited by either; however, it should be
recalled that the minimum-energy problem does not require optimal trajectories
to maintain a bang-bang thrust profile (see Eq. (24) with $\epsilon=1.0$). The
$\infty-$norm of the TPBVP final boundary condition residual was
$||\mathbf{e}(t_{f})||_{\infty}=3.95\times 10^{-2}$ and
$||\mathbf{e}(t_{f})||_{\infty}=9.72\times 10^{-3}$ for the generated guess
and the co-states returned after stalled single shooting respectively, which
does confirm slight improvement during single shooting. Furthermore, it was
found that both the generated guess and single shooting-returned co-states
result in a Jacobian of the shooting function with respect to the initial co-
states of rank six, which corresponds to rank deficiency. This example
illustrates a trial where the final boundary condition residuals were nearly
satisfied for both cases but rank deficiency of the Jacobian of the shooting
function prevented convergence to a solution of the TPBVP.

##  5 Conclusion

This paper proposed a method for initializing the co-state variables for
solving low-thrust minimum-fuel trajectory optimization problems with the
indirect shooting approach by employing PSO paired with an energy-to-fuel
homotopy technique. The method was applied successfully to solve two low-
thrust transfer problems in the Earth-Moon system. It was demonstrated that
the methodology is able to successfully generate guesses for the initial co-
state variables which converge to a solution for both investigated scenarios.
The resulting minimum-fuel trajectories were validated by comparison with
published solutions of the optimal low-thrust transfer for both scenarios.
Analysis was performed to determine the effect that varying the number of
particles in the swarm has on the performance of the proposed method. It was
found that for problems involving several low-thrust spirals about a primary
body, a larger swarm size of 2000 to 4000 particles produced the highest
quality guesses for the initial co-state variables, whereas, for problems
involving short transfers between orbits in cislunar space, a smaller swarm
size of 500 particles is preferred and is capable of rapidly generating
solutions in less than 4 minutes. It should be noted that, due to the
heuristic nature of the PSO algorithm, multiple solutions to the derived
minimum-fuel TPBVP were found, which correspond to extremal solutions of the
minimum-fuel objective function, providing interesting physical insight into
the investigated trajectory optimization scenarios. Considerations on the
existence of local minima and Jacobian rank deficiency were also discussed.

The primary benefits of applying PSO for minimum-fuel co-state initialization
are the simplicity and versatility of the approach. Once the TPBVP is derived
through the application of COV and PMP, little additional work is required to
employ PSO to initialize the co-state variables for a range of trajectory
optimization problems in similar dynamical environments. Additionally, PSO
removes the need for a guess to be provided by the end-user of the algorithm,
an advantageous feature for preliminary mission design efforts, when a large
search for many candidate trajectories is required. A large search for
candidate trajectories could also be accelerated drastically through the
exploitation of a high-performance computing environment, as the methodology
developed herein lends itself well to parallel computing and could be
distributed across many CPUs to generate solutions for many initial and target
states, flight times, and spacecraft parameters. Furthering this cause, the
proposed methodology was shown to discover multiple solutions of the derived
TPBVP, most of which could be considered candidate trajectories during
preliminary mission design efforts once local optimality is verified through
an investigation of the second-order sufficient conditions of optimality.
Finally, while not investigated in this work, the multiple solutions
discovered through the application of the proposed methodology could be
improved through further homotopy continuation on the fixed flight time or
thrust magnitude to produce even more candidate trajectories.

## Acknowledgments

This material is based upon work supported by the National Science Foundation
Graduate Research Fellowship under Grant No. 88983. Grant Hecht also
acknowledges reception of the NASA Space Grant Fellowship.

## References

  * [1] R. Woollands, E. Taheri, J. Junkins, Efficient Computation of Optimal Low Thrust Gravity Perturbed Orbit Transfers, in: 29th AAS/AIAA Spaceflight Mechanics Meeting, Univelt Inc., 2019, pp. 899–918.
  * [2] R. Chai, A. Savvaris, A. Tsourdos, S. Chai, Y. Xia, A Review of Optimization Techniques in Spacecraft Flight Trajectory Design, Progress in Aerospace Sciences 109 (2019) 100543.  doi:https://doi.org/10.1016/j.paerosci.2019.05.003.
  * [3] E. Taheri, N. I. Li, I. Kolmanovsky, Co-State Initialization for the Minimum-Time Low-Thrust Trajectory Optimization, Advances in Space Research 59 (9) (2017) 2360–2373.  doi:https://doi.org/10.1016/j.asr.2017.02.010.
  * [4] L. C. W. Dixon, M. C. Biggs, The Advantages of Adjoint–Control Transformations When Determining Optimal Trajectories by Pontryagin’s Maximum Principle, The Aeronautical Journal (1968) 76 (735) (1972) 169–174.  doi:10.1017/S0001924000042731.
  * [5] L. C. W. Dixon, M. C. Biggs, Adjoint–Control Transformations for Solving Practical Optimal Control Problems, Optimal Control Applications and Methods 2 (4) (1981) 365–381.  doi:https://doi.org/10.1002/oca.4660020405.
  * [6] J. Senent, C. Ocampo, A. Capella, Low-Thrust Variable-Specific-Impulse Transfers and Guidance to Unstable Periodic Orbits, Journal of Guidance, Control, and Dynamics 28 (2) (2005) 280–290.  doi:10.2514/1.6398.
  * [7] C. L. Ranieri, C. A. Ocampo, Optimization of Roundtrip, Time-Constrained, Finite Burn Trajectories via an Indirect Method, Journal of Guidance, Control, and Dynamics 28 (2) (2005) 306–314.  doi:10.2514/1.5540.
  * [8] C. L. Ranieri, C. A. Ocampo, Indirect Optimization of Spiral Trajectories, Journal of Guidance, Control, and Dynamics 29 (2006) 1360–1366.  doi:10.2514/1.19539.
  * [9] R. P. Russell, Primer Vector Theory Applied to Global Low-Thrust Trade Studies, Journal of Guidance, Control, and Dynamics 30 (2) (2007) 460–472.  doi:10.2514/1.22984.
  * [10] J. D. Thorne, C. D. Hall, Approximate Initial Lagrange Costates for Continuous–Thrust Spacecraft, Journal of Guidance, Control, and Dynamics 19 (1996) 283–288.  doi:10.2514/3.21616.
  * [11] D. Lee, H. Bang, Efficient Initial Costates Estimation for Optimal Spiral Orbit Transfer Trajectories Design, Journal of Guidance, Control, and Dynamics 32 (2009) 1943–1947.  doi:10.2514/1.44550.
  * [12] D. Lee, H. Bang, H.-D. Kim, Optimal Earth-Moon Trajectory Design Using New Initial Costate Estimation Method, Journal of Guidance, Control, and Dynamics 35 (2012) 1671–1676.  doi:10.2514/1.55863.
  * [13] F. Fahroo, I. M. Ross, Costate Estimation by a Legendre Pseudospectral Method, Journal of Guidance, Control, and Dynamics 24 (2) (2001) 270–277.  doi:10.2514/2.4709.
  * [14] D. A. Benson, G. T. Huntington, T. P. Thorvaldsen, A. V. Rao, Direct Trajectory Optimization and Costate Estimation via an Orthogonal Collocation Method, Journal of Guidance, Control, and Dynamics 29 (2006) 1435–1440.  doi:10.2514/1.20478.
  * [15] D. Garg, M. A. Patterson, C. Francolin, C. L. Darby, G. T. Huntington, W. W. Hager, A. V. Rao, Direct Trajectory Optimization and Costate Estimation of Finite-Horizon and Infinite-Horizon Optimal Control Problems Using a Radau Pseudospectral Method, Computational Optimization and Applications 49 (2011) 335–358.  doi:10.1007/s10589-009-9291-0.
  * [16] C. C. Françolin, D. A. Benson, W. W. Hager, A. V. Rao, Costate Approximation in Optimal Control Using Integral Gaussian Quadrature Orthogonal Collocation Methods, Optimal Control Applications and Methods 36 (2015) 381–397.  doi:https://doi.org/10.1002/oca.2112.
  * [17] J. D. Eide, W. W. Hager, A. V. Rao, Modified Radau Collocation Method for Solving Optimal Control Problems with Nonsmooth Solutions Part II: Costate Estimation and the Transformed Adjoint System, in: 2018 IEEE Conference on Decision and Control (CDC), 2018, pp. 1651–1656.  doi:10.1109/CDC.2018.8619426.
  * [18] M. Pontani, B. A. Conway, Particle Swarm Optimization Applied to Space Trajectories, Journal of Guidance, Control, and Dynamics 33 (2010) 1429–1441.  doi:10.2514/1.48475.
  * [19] M. Pontani, B. Conway, Optimal Low-Thrust Orbital Maneuvers via Indirect Swarming Method., Journal of Optimization Theory and Applications 162 (2014) 272–292.  doi:https://doi.org/10.1007/s10957-013-0471-9.
  * [20] R. Bertrand, R. Epenoy, New Smoothing Techniques for Solving Bang–Bang Optimal Control Problems – Numerical Results and Statistical Interpretation, Optimal Control Applications and Methods 23 (4) (2002) 171–197.  doi:https://doi.org/10.1002/oca.709.
  * [21] C. Zhang, F. Topputo, F. Bernelli-Zazzera, Y.-S. Zhao, Low-Thrust Minimum-Fuel Optimization in the Circular Restricted Three-Body Problem, Journal of Guidance, Control, and Dynamics 38 (8) (2015) 1501–1510.  doi:10.2514/1.G001080.
  * [22] J. D. Aziz, D. J. Scheeres, G. Lantoine, Hybrid Differential Dynamic Programming in the Circular Restricted Three-Body Problem, Journal of Guidance, Control, and Dynamics 42 (5) (2019) 963–975.  doi:10.2514/1.G003617.
  * [23] V. Szebehely, Chapter 1 – Description of the Restricted Problem, in: V. Szebehely (Ed.), Theory of Orbit, Academic Press, 1967, pp. 7–41.  doi:https://doi.org/10.1016/B978-0-12-395732-0.50007-6.
  * [24] D. F. Lawden, Optimal Trajectories for Space Navigation, The Mathematical Gazette 48 (366) (1964) 478–479.
  * [25] A. E. Bryson, Applied Optimal Control : Optimization, Estimation, and Control, rev. print. Edition, Hemisphere Pub. Corp., Washington, 1975.
  * [26] J. E. Prussing, S. L. Sandrik, Second-Order Necessary Conditions and Sufficient Conditions Applied to Continuous-Thrust Trajectories, Journal of Guidance, Control, and Dynamics 28 (4) (2005) 812–816.  doi:10.2514/1.4040.
  * [27] P. K. Mogensen, K. Carlsson, S. Villemot, S. Lyon, M. Gomez, C. Rackauckas, T. Holy, D. Widmann, T. Kelman, D. Karrasch, A. Levitt, A. N. Riseth, C. Lucibello, C. Kwon, D. Barton, J. TagBot, M. Baran, M. Lubin, S. Choudhury, S. Byrne, S. Christ, T. Arakaki, T. A. Bojesen, Benneti, M. R. G. Macedo, JuliaNLSolvers/NLsolve.jl: v4.5.1, accessed: 2022-10-04 (Dec. 2020).
URL https://github.com/JuliaNLSolvers/NLsolve.jl

  * [28] J. Bezanson, A. Edelman, S. Karpinski, V. B. Shah, Julia: A Fresh Approach to Numerical Computing, SIAM Review 59 (1) (2017) 65–98.  doi:10.1137/141000671.
  * [29] C. Rackauckas, Q. Nie, DifferentialEquations.jl – A Performant and Feature-Rich Ecosystem for Solving Differential Equations in Julia, Journal of Open Research Software 5 (1) (2017) 15.  doi:http://doi.org/10.5334/jors.151.
  * [30] J. Kennedy, R. Eberhart, Particle Swarm Optimization, in: Proceedings of ICNN’95-International Conference on Neural Networks, Vol. 4, IEEE, 1995, pp. 1942–1948.
  * [31] E. Mezura-Montes, C. A. C. Coello, Constraint-Handling in Nature-Inspired Numerical Optimization: Past, Present and Future, Swarm and Evolutionary Computation 1 (4) (2011) 173–194.  doi:https://doi.org/10.1016/j.swevo.2011.10.001.
  * [32] M. E. H. Pedersen, Good Parameters for Particle Swarm Optimization, Hvass Lab., Copenhagen, Denmark, Tech. Rep. HL1001 (2010) 1551–3203.
  * [33] Particle swarm optimization algorithm, accessed: 2022-05-11 (2022).
URL https://www.mathworks.com/help/gads/particle-swarm-optimization-
algorithm.html

  * [34] G. Mingotti, F. Topputo, F. Bernelli‐Zazzera, Combined Optimal Low‐Thrust and Stable‐Manifold Trajectories to the Earth‐Moon Halo Orbits, AIP Conference Proceedings 886 (1) (2007) 100–112.  doi:10.1063/1.2710047.
  * [35] R. Armellin, F. Topputo, A Sixth-Order Accurate Scheme for Solving Two-Point Boundary Value Problems in Astrodynamics, Celestial Mechanics and Dynamical Astronomy 96 (3-4) (2006) 289–309.  doi:10.1007/s10569-006-9047-4.
  * [36] C. Martin, B. A. Conway, Optimal Low-Thrust Trajectories Using Stable Manifolds, Cambridge Aerospace Series, Cambridge University Press, 2010, p. 238–262.  doi:10.1017/CBO9780511778025.010.
  * [37] M. Rasotto, R. Armellin, P. D. Lizia, F. Bernelli-Zazzera, Optimal Low-Thrust Transfers in Two-Body and Three-Body Dynamics, in: 64th International Astronautical Congress, Vol. 7, 2014, pp. 5230 – 5244.
URL https://epubs.surrey.ac.uk/840753/

  * [38] J. H. Verner, Numerically Optimal Runge–Kutta Pairs with Interpolants, Numerical Algorithms 53 (2-3) (2010) 383–396.  doi:https://doi.org/10.1007/s11075-009-9290-3.

## Appendix: Algorithm for PSO Co-State Initialization

## Vitae

Grant R. Hecht is a Ph.D. candidate in the Department of Mechanical and
Aerospace Engineering at the University at Buffalo. He is also an NSF GRFP
fellow and a NASA Pathways Student. He holds a B.Sc. in Aerospace Engineering
from the Missouri University of Science and Technology. His main research
interests are in spacecraft trajectory optimization, formation guidance, and
inertial navigation.

Dr. Eleonora M. Botta is an assistant professor in the Department of
Mechanical and Aerospace Engineering at the University at Buffalo, NY. She
holds a Ph.D. in Mechanical Engineering from McGill University, a B.Sc. in
Aerospace Engineering and a M.Sc. in Space Engineering from Politecnico di
Milano, and a M.Sc. in Aerospace Engineering from Politecnico di Torino. Her
main research interests are the dynamics and control of tethered space
systems, spacecraft trajectory optimization, and space situational awareness.
She is a member of the AAS Space Flight Mechanics Committee and the vice-chair
of the AIAA Space Tethers Technical Committee.

  *[COV]: Calculus of Variations
  *[PMP]: Pontryagin’s Maximum Principle
  *[TPBVP]: Two-Point Boundary-Value Problem
  *[ACT]: Adjoint Control Transformation
  *[PSO]: Particle Swarm Optimization
  *[CR3BP]: Circular Restricted Three Body Problem
  *[STM]: State Transition Matrix