# A Koopman-Operator Control Optimization for Relative Motion in Space

Simone Servadio111Postdoctoral Associate, Department of Aeronautics and
Astronautics, Massachusetts Institute of Technology, MA 02139, USA. email:
<EMAIL_ADDRESS>Massachusetts Institute of Technology, Cambridge, MA, 02139,
USA  Roberto Armellin222Professor, Mechanical Engineering,
<EMAIL_ADDRESS>Auckland Space Institute, The University of
Auckland, Auckland, NZL  Richard Linares333Boeing Assistant Professor,
Department of Aeronautics and Astronautics, Massachusetts Institute of
Technology, MA 02139, USA. email<EMAIL_ADDRESS>Massachusetts Institute of
Technology, Cambridge, MA, 02139, USA

###### Abstract

A high order optimal control strategy implemented in the Koopman operator
framework is proposed in this work. The new technique exploits the Koopman
representation of the solution of the equations of motion to develop an energy
optimal inverse control methodology. The operator theory can reformulate a
nonlinear dynamical system of finite dimension into a linear system with an
infinite number of dimensions. As a results, the state of any nonlinear
dynamics is represented as a linear combination of high-order orthogonal
polynomials, which creates the state transition polynomial map of the
solution. Since the optimal control technique can be reduced to a two-points
boundary value problem, the Koopman map is used to connect the state and
control variables in time, such that optimal values are obtained through map
inversion and polynomial evaluation. The new technique is applied to
rendezvous applications in space, where the relative motion between two
satellites is modelled with a high-order polynomial series expansion of the
Lagrangian of the system, such that the Clohessy-Wiltshire equations represent
the reduction of the high-order model to a linear truncation.

## Introduction

Trajectory optimization is a key element in every space mission [1].
Trajectory optimization enables complex interplanetary transfers [2],
operative orbit acquisition (e.g., from a geostationary transfer orbit to a
geostationary one [3]), rendezvous and docking [4], formation flying control
and reconfiguration [5] and station-keeping [6], just to mention a few. The
minimization of propellant or transfer time is the typical objective of
trajectory optimization.

Low-thrust propulsion is more frequently used as it allows to reduce the
propellant mass thanks to high specific impulse. Trajectory optimization with
low-thrust propulsion is a challenging task as it involves solving complex
optimal control problems (OCP) [7]. Two main methods are used to solve low-
thrust optimization problems: direct and indirect. In the direct methods, the
control is parametrized and the dynamics are reduced into a set of
constraints. In this way, the OPC becomes a parametric optimization problem
solved with nonlinear programming tools. The large dimension of the problem
and possible convergence to suboptimal solutions are two main issue of direct
methods. Indirect methods can find optimum solutions more accurately and
result in a lower dimension numerical problem. These methods apply calculus of
variation and Pontryagin’s maximum principle to transform the OCP into a two-
point boundary value problem [8]. The main drawback of indirect methods is
that the convergence domain is small, guesses for initial costate are not
easily found, and discontinuity in the control laws are difficult to handle
[9].

In the astrodynamics community there is a strong interest in developing
techniques to solve OCP in an efficient way for set of initial and final
conditions. This interest is fueled by the need for performing trade-off
analyses efficiently, quantifying the effect of uncertainties and better
estimate mission budgets, or use OCP solutions for spacecraft guidance. The
formulation direct methods as convex optimization problems has received a lot
of attention in recent years [10]. Thanks to the convergence properties of
convex optimization and the efficiency of modern solvers, OCP can be
potentially solved onboard and the solutions used for spacecraft guidance [11,
12].

An alternative approach for autonomous guidance is to computing the expansion
of the solution of OCPs with respect to uncertainties in initial and final
conditions, using state transition tensors (STTs) [13, 14] and differential
algebra (DA) [15, 16]. These expansions can allow for the update of the
controls due to the presence of errors in the spacecraft or target states.
Indeed the required update of the reference trajectory is reduced to a simple
evaluation of polynomials.

In the same vein, this work introduces an alternative way to efficiently
compute the solution of OCP formulated in a indirect approach. Our approach is
based on the Koopman operator (KO). The KO theory was initially developed by
Koopman [17] and further developed by Von Neumann in the 1930s [18]. They
transitioned the operator theory initially developed in quantum mechanics [19,
20] to classical mechanics problems. The key result was the development of a
linear operator that can reformulate a nonlinear dynamical system of finite
dimension into a linear system with an infinite number of dimensions. The
linearity of the KO is very appealing, but this benefit is contrasted with the
fact that it is infinite-dimensional. However, this issue can be overcome by
capturing the evolution on a finite subspace spanned by a finite set of basis
functions. In effect, this is a truncation of the KO to a finite subspace.
More recently, the KO has been employed in astrodynamics to model the motion
of a satellite about an oblate planet [21, 22] and to analyze the orbits of a
satellite in the close proximity of libration points [23].

The KO solution of the system represents the time propagation of the state as
a linear combination of well selected polynomial basis functions. This works
exploit this representation to create a high-order state polynomial transition
map to solve the TPBVP, following the inverse control approach. Indeed, the
map can be inverted [15] to obtain the optimal control that connects, in the
desired amount of time, the state of the system between the boundary
conditions. The final result is a high-order inverse control technique that
does not require the evaluation of a reference solution, contrary to previous
approaches [15, 16].

The problem we aim to solve with the KO-based approach is the one of a
rendezvous with a chief spacecraft moving in a circular Keplerian orbit. The
transfer is implemented with a low-thrust propulsion system and the objective
is to minimize the energy of the transfer. Differently from the classical
Clohessy-Wiltshire equations [24], we retain terms up to an arbitrary high
order in the Legendre expansion of the gravitational potential. The equations
of motion are derived from the polynomial series expansion of the Lagrangian
of the system [25], where the accuracy of the representation of relative
motion increases due to the presence of high order terms [24]. This results in
a nonlinear control problem that can only be solved numerically due to the
nonlinearities in the Euler-Lagrange equations.

The paper is structured as follows. The Koopman operator theory is presented
in the following section. Afterwards, the energy optimal control technique is
derived in the KO framework, highlighting the advantages of the newly
developed methodology. Then, the numerical application are presented. First,
the Duffing oscillator is controlled to give a representation of the new
technique for a easy problem. Secondly, the linear and high-order controlled
Clohessy-Wiltshire equations are derived, providing comparison between the two
dynamics. The optimal inverse KO control is applied to achieve rendezvous
between a chaser spacecraft and a target for different times of flight and
boundary conditions. Lastly, conclusions are drawn.

## The Koopman Operator

A classical definition of nonlinear dynamical systems is given by the initial
value problem, which can be represented by a set of coupled autonomous
ordinary differential equations in the form:

$\left\\{\begin{tabular}[]{l}$\displaystyle\frac{d}{dt}{\bf x}(t)={\bf f}({\bf
x})$\\\ ${\bf x}(t_{0})={\bf x_{0}}$\end{tabular}\right.$ (1)

where ${\mathbf{x}}\in\mathbb{R}^{m}$ is the state which depends on the time
evolution $t$, ${\bf f}:\mathbb{R}^{m}\rightarrow\mathbb{R}^{m}$ is the
nonlinear dynamics model, $m$ is the number of dimensions in which the problem
is defined, and ${\bf x}_{0}$ is the initial condition of the system at time
$t_{0}$. The Koopman operator (KO), $\mathcal{K}(\cdot)$ is an infinite-
dimensional linear operator that evolves all observable functions
$\mathcal{G}({\bf x})$ of the state and it allows to redefine any problem
written in classical mechanics into an equivalent Koopman formulation.

Let $\mathcal{F}$ be a vector space of observable functions, where
$\mathcal{G}({\bf x})\in\mathcal{F}$. Since the KO is an infinite-dimensional
linear operator, this space of functions $\mathcal{F}$, which the observables
are defined on, is also infinite-dimensional. Therefore, if
$g\subseteq\mathcal{G}({\bf x})$ is a given observable in this space, the
evolution of $g$ in the dynamical system is represented by:

$\mathcal{K}\left(g({\bf x})\right)=\frac{d}{dt}g({\bf x})=\left(\nabla_{{\bf
x}}g({\bf x})\right)\frac{d}{dt}{\bf x}(t)=\left(\nabla_{{\bf x}}g({\bf
x})\right){\bf f}({\bf x}),$ (2)

where $\nabla_{{\bf x}}g=(\partial g/\partial x_{1},\partial g/\partial
x_{2},\dots,\partial g/\partial x_{m})$. That way, the evolution of any
observable subjected to the dynamical system is provided by the Koopman
operator:

$\mathcal{K}\left(\cdot\right)=\left(\nabla_{{\bf x}}\cdot\right){\bf f}({\bf
x}).$ (3)

Note that the evolution of the observables is provided by the application of
the chain rule for the time derivative of $g({\bf x})$. Consequently, the
defined operator is linear [17], such that:

$\mathcal{K}\left(\beta_{1}g_{1}({\bf x})+\beta_{2}g_{2}({\bf
x})\right)=\beta_{1}\mathcal{K}\left(g_{1}({\bf
x})\right)+\beta_{2}\mathcal{K}\left(g_{2}({\bf x})\right),$ (4)

for any pair of observables $g_{1}\subseteq\mathcal{G}({\bf x})$ and
$g_{2}\subseteq\mathcal{G}({\bf x})$ and any arbitrary constants $\beta_{1}$
and $\beta_{2}$. The linearity of the Koopman Operator is very appealing, but
this benefit is contrasted with the fact that it is infinite-dimensional.
However, this issue can be overcome by capturing the evolution of the system
on a finite subspace spanned by a finite set of basis functions instead of all
measurement functions in a Hilbert space. In effect, this is a truncation of
the Koopman operator to a finite subspace $\mathcal{F}_{D}$ of dimension $n$,
where $\mathcal{F}_{D}\in\mathcal{F}$. This subspace $\mathcal{F}_{D}$ can be
spanned by any set of eigenfunctions $\phi_{i}\in\mathcal{F}_{D}$, with
$i\in\\{1,2,\dots,n\\}$, defined as:

$\mathcal{K}\left(\phi_{i}({\bf x})\right)=\frac{d}{dt}\phi_{i}({\bf
x})=\lambda_{i}\phi_{i}({\bf x}),$ (5)

where $\lambda_{i}$ are the eigenvalues associated with the eigenfunctions
$\phi_{i}$, and $n$ is the number of eigenfunctions chosen to represent the
space. Therefore, the Koopman eigenfunctions can be used to form a
transformation of variables that linearizes the system. Particularly, let
${\bf\Phi}({\bf x})=\left(\phi_{1}({\bf x}),\dots,\phi_{n}({\bf
x})\right)^{T}$ be the set of eigenfunctions of the KO in $\mathcal{F}_{D}$.
Then, using the relation in Eq. (5), it is possible to write the evolution of
${\bf\Phi}$ as:

$\mathcal{K}\left(\bf\Phi\right)=\frac{d}{dt}{\bf\Phi}=\Lambda{\bf\Phi},$ (6)

where $\boldsymbol{\Lambda}=\text{diag}([\lambda_{1},\dots,\lambda_{n}])$ is
the diagonal matrix containing the eigenvalues of the system in
$\mathcal{F}_{D}$. This transformation is called the Koopman Canonical
Transform [26]. The solution of Eq. (6) is:

${\bf\Phi}(t)=\exp(\boldsymbol{\Lambda}t){\bf\Phi}(t_{0}),$ (7)

where ${\bf\Phi}(t_{0})$ is the value of the eigenfunctions at the initial
time $t_{0}$. This result will be used later to solve the complete system once
the eigenfunctions of the operator are obtained.

In general, we are interested in the identity observable, that is, ${\bf
g}({\bf x})={\bf x}$. Therefore, it is required to be able to represent these
observables in terms of the KO eigenfunctions. This is achieved using the
Koopman modes, i. e., the projection of the full-state observable onto the KO
eigenfunctions. If this projection can be found, the evolution of the state is
represented by means of the evolution of the KO eigenfunctions, and thus, an
approximate solution to the system can be provided.

### The Galerkin Method

The eigenfunction of the Koopman operator are obtained through the Galerkin
method. The eigenfunctions give rise to a set of linear first-order PDEs for
the eigenfunctions in the form:

$\mathcal{K}(\phi_{i})=\left(\nabla_{x}\phi_{i}({\bf x})\right){\bf f}({\bf
x})=\lambda_{i}\phi_{i}({\bf x}),$ (8)

or in a more expanded notation:

$\frac{d\phi_{i}({\bf x})}{dt}={f}_{1}\left({\bf
x}\right)\frac{\partial}{\partial x_{1}}\phi_{i}({\bf
x})+\cdots+{f}_{d}\left({\bf x}\right)\frac{\partial}{\partial
x_{d}}\phi_{i}({\bf x})=\lambda_{i}\phi_{i}({\bf x}),$ (9)

where ${\bf f}({\bf x})=(f_{1}({\bf x}),f_{2}({\bf x}),\dots,f_{m}({\bf
x}))^{T}$. This equation is a linear first-order PDE and in general has no
closed-form solution. However, it is possible to approximate the solution
using the Galerkin method.

This work makes use of Legendre polynomials, due to the advantages they
provide in the computation of the Koopman matrix [21]. Legendre polynomials
are a set of orthogonal polynomials defined in a Hilbert space that generate a
complete basis. The idea of this methodology is to represent any function of
the space by using this set of basis functions. This is done by the use of
inner products and the correct normalization of the Legendre polynomials. Let
$f$ and $g$ be two arbitrary functions from the Hilbert space considered.
Then, the inner product between these two functions is defined as:

$\langle f,g\rangle=\int_{\Omega}f({\bf x})g({\bf x})w({\bf x})d{\bf x},$ (10)

where $w({\bf x})$ is a positive weighting function defined on the space
domain $\Omega$. This work makes use of the normalized Legendre polynomials as
proposed by Arnas and Linares [27], due to the advantages they provide in the
computation of the inner product and in the accuracy of obtained in the
Koopman matrix. First, normalized Legendre polynomials are a set of
orthonormal polynomials, which simplifies the projection process. Second,
Legendre polynomials are defined in a bounded domain $\Omega=[-1,1]$, as
opposed to other basis functions, like Hermite polynomials, that are defined
in all $\mathbb{R}$. Third, their weighted function is a constant, $w({\bf
x})=1$. This aspect not only simplifies the expression of the inner product,
such that it becomes easier to implement a convenient automatic computation of
the integrals, but also approximates functions uniformly in the domain of
definition [27, 28]. In addition, the normalized Legendre polynomials are
defined such that:

$\langle\mathcal{L}_{i},\mathcal{L}_{j}\rangle=\int_{\Omega}\mathcal{L}_{i}({\bf
x})\mathcal{L}_{j}({\bf x})w({\bf x})d{\bf x}=\delta_{ij},$ (11)

where $\mathcal{L}_{i}$ and $\mathcal{L}_{j}$ with
$\\{i,j\\}\in\\{1,\dots,n\\}$ are two given normalized Legendre polynomials
from the set of basis functions selected, and $\delta_{ij}$ is Kronecker’s
delta.

### Time Evolution of the Basis Functions and the Koopman Matrix

The KO is based on the selection of a set of $n$ orthogonal basis functions
that describe the space of solutions. Legendre polynomials have been selected
as a base for the theory developed in this paper. Therefore, following the KO
approach, it is desired to describe the rate of change of the basis functions,
in time, as a linear combination of the functions themselves:

$\dfrac{d\mathcal{L}}{dt}={\mathbf{K}}\mathcal{L}$ (12)

where $\mathcal{L}$ are the Legendre polynomials and ${\mathbf{K}}$ is the
Koopman matrix. Equation (12) is linear but high dimensional. The original ODE
is expressed as a combination of linear differential equations, which solution
would be exact for an infinite dimensional solution space. However, due to
practicality, the Koopman solution of the ODE is an approximation that works
on a subset of the infinite Hilbert space that correctly describes the
solution. As such, the higher the order of the KO, the more accurate the
Koopman solution, since the observables are described using more
eigenfunctions that better deal with high nonlinearities. Each row of Eq. (12)
describes how each single total derivative of the Legendre polynomial behaves
as a linear combination of the Legendre polynomials themselves. Therefore,
each entry of the Koopman matrix is evaluated through the inner product, which
calculates the projection of the derivative into the basis functions. Thus,
the $i$th row can be written as

$\displaystyle\dfrac{d\mathcal{L}_{i}}{dt}$
$\displaystyle=\langle\dfrac{d\mathcal{L}_{i}}{dt},\mathcal{L}_{0}\rangle\mathcal{L}_{0}+\langle\dfrac{d\mathcal{L}_{i}}{dt},\mathcal{L}_{1}\rangle\mathcal{L}_{1}+\langle\dfrac{d\mathcal{L}_{i}}{dt},\mathcal{L}_{2}\rangle\mathcal{L}_{2}+\dots=\sum_{j=0}^{n}\langle\dfrac{d\mathcal{L}_{i}}{dt},\mathcal{L}_{j}\rangle\mathcal{L}_{j}$
(13)

From this expression, the single entry $\mathbf{K}_{ij}$ of the Koopman matrix
can be evaluated as

$\mathbf{K}_{ij}=\langle\dfrac{d\mathcal{L}_{i}}{dt},\mathcal{L}_{j}\rangle$
(14)

that indicates the projection of the rate of change of the $i$th basis
function onto the $j$th basis function. Therefore, the KO matrix is square
with dimensions $n\times n$. This formulation completely describes the
evolution of the system in time, for each basis function. The inner product
has already been defined, and it evaluates each component of the Koopman
matrix via the Galerkin method. The remaining term that requires handling is
the total derivative of the basis functions. However, using chain rule, we can
decompose each total derivative in terms of its partial derivatives with
respect to the states. Considering the general case of a system with $m$
states, the $i$th total derivative is expanded as

$\displaystyle\dfrac{d\mathcal{L}_{i}}{dt}$
$\displaystyle=\dfrac{\partial\mathcal{L}_{i}}{\partial
x_{1}}\dfrac{dx_{1}}{dt}+\dfrac{\partial\mathcal{L}_{i}}{\partial
x_{2}}\dfrac{dx_{2}}{dt}+\dfrac{\partial\mathcal{L}_{i}}{\partial
x_{3}}\dfrac{dx_{3}}{dt}+\dots=\sum_{j=1}^{m}\dfrac{\partial\mathcal{L}_{i}}{\partial
x_{j}}\dfrac{dx_{j}}{dt}=\sum_{j=1}^{m}\dfrac{\partial\mathcal{L}_{i}}{\partial
x_{j}}{f}_{j}({\bf x}).$ (15)

By looking at the equation, it can be noted that each partial derivative is
always known and it does not depend directly on the problem. On the contrary,
${f}_{j}({\bf x})$ is one equation of motion, which changes depending on the
problem and that describe the selected application.

### The Koopman Solution of the System

The system is described through the Koopman matrix. However, in order to
obtain the solution of the system, we are interested in the time evolution of
the state. Thanks to the Koopman representation, any function of the state can
be approximated as a linear combination of the basis functions. Let
$\mathbf{g}(\mathbf{x})$ be a set of observables that we are interested to
evaluate. Any $i$th observable can be projected on the set of basis functions
using the inner product:

$\mathbf{g}=\sum_{j=0}^{n}\langle\mathbf{g},\mathcal{L}_{j}\rangle\mathcal{L}_{j}$
(16)

These observables can be any function of the original variables $\mathbf{x}$,
including the states themselves, which constitutes identity observables.
Therefore, the set of coefficients that expresses the projection of the
observables into the basis functions can be represented in matrix form,
likewise for the Koopman matrix:

$\mathbf{H}_{ij}=\langle\mathbf{g}_{i},\mathcal{L}_{j}\rangle$ (17)

where $\mathbf{H}_{i,j}$ indicates the projection of the $i$th observable onto
the $j$th basis function. The number of observables dictates the dimensions of
$\mathbf{H}$: in the particular case of the identity observables, where it is
desired to obtain the state of the system, matrix $\mathbf{H}$ has dimensions
$m\times n$.

The observables are represented as a linear combination of the basis
functions, which time evolution is approximated by the Koopman matrix. Under
the assumption of matrix $\mathbf{K}$ being diagonalizable, the
eigendecomposition of the dynamics can be performed;

$\mathbf{V}\mathbf{K}=\boldsymbol{\Lambda}\mathbf{V}$ (18)

where $\mathbf{V}$ is the eigenvectors matrix and $\boldsymbol{\Lambda}$ is
the diagonal matrix of eigenvalues. Rewriting Eq. (12) by highlighting the
variable dependency of the functions

$\dfrac{d}{dt}(\mathcal{L}(\mathbf{x}(t)))=\mathbf{K}\mathcal{L}(\mathbf{x}(t)),$
(19)

thanks to the eigendecomposition of the Koopman matrix, Eq. (18), the
eigenfunctions of the system are obtained from the Legendre polynomials using
the eigenvectors matrix:

$\boldsymbol{\phi}(\mathbf{x}(t))=\mathbf{V}\mathcal{L}(\mathbf{x}(t))$ (20)

The analysis of the eigenfunctions brings us to a simpler description of the
system where each differential equation is decoupled from the others, and,
therefore, easier to solve. Let us take the time derivative of the
eigenfunctions and substitute what has been derived so far:

$\displaystyle\dfrac{d}{dt}\boldsymbol{\phi}(\mathbf{x}(t))$
$\displaystyle=\dfrac{d}{dt}(\mathbf{V}\mathcal{L}(\mathbf{x}(t)))=\mathbf{V}\dfrac{d}{dt}(\mathcal{L}(\mathbf{x}(t)))$
(21a) $\displaystyle=\mathbf{V}\mathbf{K}\mathcal{L}(\mathbf{x}(t))$ (21b)
$\displaystyle=\boldsymbol{\Lambda}\mathbf{V}\mathcal{L}(\mathbf{x}(t))$ (21c)
$\displaystyle=\boldsymbol{\Lambda}\boldsymbol{\phi}(\mathbf{x}(t))$ (21d)

where the substitutions come, respectively, from Eq. (19), Eq. (18), and Eq.
(20). The last relation, Eq. (21d), shows a diagonal system of ODEs whose
solution is known

$\boldsymbol{\phi}(\mathbf{x}(t))=\exp(\boldsymbol{\Lambda}t)\boldsymbol{\phi}(\mathbf{x}(t_{0}))$
(22)

where $\boldsymbol{\phi}(\mathbf{x}(t_{0}))$ indicates the value of the
eigenfunctions at the initial time $t_{0}$.

The evolution of the eigenfunctions with time is used to find the solution of
any observable, and thus, of the state of the system (when the observable is
the identity). Each observable function $\mathbf{g}(\mathbf{x}(t))$ has
already been represented as a linear combination of the basis functions of the
system through matrix $\mathbf{H}$. Therefore, after some manipulations, the
time evolution of the observables can be expressed solely as a function of
time $t$ and of the state initial condition $\mathbf{x}(t_{0})$ :

$\displaystyle\mathbf{g}(\mathbf{x}(t))$
$\displaystyle=\mathbf{H}\mathcal{L}(\mathbf{x}(t))$ (23a)
$\displaystyle=\mathbf{H}\mathbf{V}^{-1}\boldsymbol{\phi}(\mathbf{x}(t))$
(23b)
$\displaystyle=\mathbf{H}\mathbf{V}^{-1}\exp(\boldsymbol{\Lambda}t)\boldsymbol{\phi}(\mathbf{x}(t_{0}))$
(23c)
$\displaystyle=\mathbf{H}\mathbf{V}^{-1}\exp(\boldsymbol{\Lambda}t)\mathbf{V}\mathcal{L}(\mathbf{x}(t_{0}))$
(23d)

where the substitutions come, respectively, from the inversion of Eq. (20),
Eq. (22), and Eq. (20). The solution of the equations of motion has been
found, having picked the state of the system as observables.

A more complete and detailed explanation of the Koopman framework, with a
curated analysis and derivation of the mathematics of the operator, is offered
in Ref. [21, 28, 22, 23].

## Optimal Control via Koopman Polynomial Map Inversion

The KO solution of the system is expressed, for each time step, as a linear
combination of the system basis functions, which get evaluated at the initial
condition. This representation suits the evaluation of optimal control via
indirect method.

### Energy Optimal Costate Evaluation

Assume that we want to link two states, $\mathbf{x}_{0}$ and $\mathbf{x}_{f}$
in a given amount of time, time-of-flight (ToF) $\bar{t}$,

$\mathbf{x}(t_{0})=\mathbf{x}_{0},\quad\quad\mathbf{x}(t_{f})=\mathbf{x}_{f},\quad\quad
t_{f}-t_{0}=\bar{t},$ (24)

consuming minimum energy. The costs function $\mathcal{J}$ related to the
performance index $\mathcal{P}$ to be optimized is

$\mathcal{J}=\int^{t_{f}}_{t_{0}}\mathcal{P}dt=\int^{t_{f}}_{t_{0}}\dfrac{1}{2}\mathbf{\tau}^{2}dt$
(25)

After separating the state of the system into position and velocity,
$\mathbf{x}=[\mathbf{r}\quad\mathbf{v}]^{T}$, the problem statement can be
written as:

$\mathbf{p}(\mathbf{x},\tau,\hat{\boldsymbol{\alpha}})=\begin{cases}\dot{\mathbf{r}}=\mathbf{v}\\\
\dot{\mathbf{v}}=\mathbf{f}\big{(}\mathbf{r},\mathbf{v}\big{)}+\tau\hat{\boldsymbol{\alpha}}\end{cases}$
(26)

where $\hat{\boldsymbol{\alpha}}$ describes the direction of the control,
while $\mathbf{f}\big{(}\mathbf{r},\mathbf{v}\big{)}$ represents the dynamics,
assumed as polynomial, or as the polynomial approximation of any selected
application.

The evaluation of the optimal control, at this step, follows classic
derivations [8]. The solution method consists into defining an auxiliary
function, the control Hamiltonian:

$\displaystyle\mathcal{H}(\mathbf{r},\mathbf{v},\tau,\hat{\boldsymbol{\alpha}})$
$\displaystyle=\mathcal{P}+\boldsymbol{\lambda}^{T}\mathbf{p}\big{(}\mathbf{x},\tau,\hat{\boldsymbol{\alpha}}\big{)}$
(27)
$\displaystyle=\dfrac{1}{2}\mathbf{\tau}^{2}+\boldsymbol{\lambda}_{\mathbf{r}}^{T}\mathbf{v}+\boldsymbol{\lambda}_{\mathbf{v}}^{T}\left(\mathbf{f}\big{(}\mathbf{r},\mathbf{v}\big{)}+\tau\hat{\boldsymbol{\alpha}}\right)$
(28)

where vector $\boldsymbol{\lambda}$ represents the costates, with dimension
equal to the number of states, divided into costates relative to position,
$\boldsymbol{\lambda}_{\mathbf{r}}$, and costates relative to the velocity,
$\boldsymbol{\lambda}_{\mathbf{v}}$. The equations of motion are obtained by
applying the Hamilton’s equations to the Hamiltonian, where optimal control is
found through minimization.

$\begin{cases}\dot{\mathbf{r}}=\dfrac{\partial\mathcal{H}}{\partial\boldsymbol{\lambda}_{\mathbf{r}}}=\mathbf{v}\\\
\dot{\mathbf{v}}=\dfrac{\partial\mathcal{H}}{\partial\boldsymbol{\lambda}_{\mathbf{v}}}=\mathbf{f}\big{(}\mathbf{r},\mathbf{v}\big{)}+\tau\hat{\boldsymbol{\alpha}}\\\
\dot{\boldsymbol{\lambda}_{\mathbf{r}}}=-\dfrac{\partial\mathcal{H}}{\partial\mathbf{r}}=-\left[\dfrac{\partial\mathbf{f}}{\partial\mathbf{r}}\right]^{T}\boldsymbol{\lambda}_{\mathbf{v}}\\\
\dot{\boldsymbol{\lambda}_{\mathbf{v}}}=-\dfrac{\partial\mathcal{H}}{\partial\mathbf{v}}=-\boldsymbol{\lambda}_{\mathbf{r}}-\left[\dfrac{\partial\mathbf{f}}{\partial\mathbf{v}}\right]^{T}\boldsymbol{\lambda}_{\mathbf{v}}\end{cases}$
(29)

This system represents the augmented dynamics, where the time behaviour of the
costates has been stacked onto the original set of ODEs. In order to find the
optimal control, the terms in the Hamiltonian that depend only on $\tau$ and
$\hat{\boldsymbol{\alpha}}$ are isolated:

$\displaystyle\mathcal{H}(\tau,\hat{\boldsymbol{\alpha}})=\dfrac{1}{2}\mathbf{\tau}^{2}+\boldsymbol{\lambda}_{\mathbf{v}}^{T}\tau\hat{\boldsymbol{\alpha}}$
(30)

By inspection, it is clear that the direction of trust that minimizes
$\mathcal{H}$ is

$\hat{\boldsymbol{\alpha}}^{*}=-\dfrac{\lambda_{\mathbf{v}}}{||\lambda_{\mathbf{v}}||}$
(31)

where starred variables indicate optimal values. By substituting the optimal
direction back into the Hamiltonian

$\mathcal{H}=\dfrac{1}{2}\mathbf{\tau}^{2}-\tau||\boldsymbol{\lambda}_{\mathbf{v}}||$
(32)

the necessary condition for the evaluation of $\tau$ is

$\dfrac{\partial\mathcal{H}}{\partial\tau}=0\quad\Rightarrow\quad\tau-||\boldsymbol{\lambda}_{\mathbf{v}}||=0\quad\Rightarrow\quad\tau^{*}=||\boldsymbol{\lambda}_{\mathbf{v}}||$
(33)

Therefore, the optimal control law is explicit in $\mathbf{x}$ and
$\boldsymbol{\lambda}$

$\tau^{*}\hat{\boldsymbol{\alpha}}^{*}=||\boldsymbol{\lambda}_{\mathbf{v}}||\left(-\dfrac{\lambda_{\mathbf{v}}}{||\lambda_{\mathbf{v}}||}\right)=-\boldsymbol{\lambda}_{\mathbf{v}}$
(34)

and it can be substituted into the equations of motion:

$\begin{cases}\dot{\mathbf{r}}=\dfrac{\partial\mathcal{H}}{\partial\boldsymbol{\lambda}_{\mathbf{r}}}=\mathbf{v}\\\
\dot{\mathbf{v}}=\dfrac{\partial\mathcal{H}}{\partial\boldsymbol{\lambda}_{\mathbf{v}}}=\mathbf{f}\big{(}\mathbf{r},\mathbf{v}\big{)}-\boldsymbol{\lambda}_{\mathbf{v}}\\\
\dot{\boldsymbol{\lambda}_{\mathbf{r}}}=-\dfrac{\partial\mathcal{H}}{\partial\mathbf{r}}=-\left[\dfrac{\partial\mathbf{f}}{\partial\mathbf{r}}\right]^{T}\boldsymbol{\lambda}_{\mathbf{v}}\\\
\dot{\boldsymbol{\lambda}_{\mathbf{v}}}=-\dfrac{\partial\mathcal{H}}{\partial\mathbf{v}}=-\boldsymbol{\lambda}_{\mathbf{r}}-\left[\dfrac{\partial\mathbf{f}}{\partial\mathbf{v}}\right]^{T}\boldsymbol{\lambda}_{\mathbf{v}}\end{cases}$
(35)

The optimal control has been resolved and the problem is now a TPBVP where the
two initial states are known,
$\mathbf{x}(t_{0})=[\mathbf{r}(t_{0})\quad\mathbf{v}(t_{0})]^{T}=\mathbf{x}_{0}$
and
$\mathbf{x}(t_{f})=[\mathbf{r}(t_{f})\quad\mathbf{v}(t_{f})]^{T}=\mathbf{x}_{f}$,
with the unknown of the initial costates
$\boldsymbol{\lambda}(t_{0})=[\boldsymbol{\lambda}_{\mathbf{r}}(t_{0})\quad\boldsymbol{\lambda}_{\mathbf{v}}(t_{0})]^{T}$,
that propagated forward in time give exactly $\mathbf{x}_{f}$.

### The Koopman Solution

The equations of motion expressed in system (35) can be solved using the KO
approach. In the KO methodology, the solution of the set of ODEs is expressed
as a linear combination of well-defined basis functions, where the evaluation
of the Koopman matrix describes their time evolution and rate of change.
Therefore, considering the augmented state composed of the original state
variables with the addition of the costates, the solution of the state can be
obtained as described by Eq. (23d). Let us highlight the dependencies on time,
state, and costate variables;

$\mathbf{x}(\mathbf{x}(t_{0}),\boldsymbol{\lambda}(t_{0}),t)=\mathbf{H}\mathbf{V}^{-1}\exp(\boldsymbol{\Lambda}t)\mathbf{V}\mathcal{L}(\mathbf{x}(t_{0}),\boldsymbol{\lambda}(t_{0}))$
(36)

where the matrix of observables $\mathbf{H}$ has been calculated using the
Galerkin method, after selecting the identity observables for the states.
Equation (36) is a polynomial map, function of $\mathbf{x}_{0}$,
$\boldsymbol{\lambda}_{0}$, and $t$, created by the basis functions, for each
component of the state. In the TPBVP, we are interested in the values of the
state at the given final time $t_{f}$, thus, the time dependency on the
polynomials can be removed:

$\mathbf{x}_{f}(\mathbf{x}(t_{0}),\boldsymbol{\lambda}(t_{0}))=\mathbf{H}\mathbf{V}^{-1}\exp(\boldsymbol{\Lambda}t_{f})\mathbf{V}\mathcal{L}(\mathbf{x}(t_{0}),\boldsymbol{\lambda}(t_{0}))=\mathcal{M}_{t_{0}\rightarrow
t_{f}}(\mathbf{x}(t_{0}),\boldsymbol{\lambda}(t_{0}))$ (37)

where $\mathcal{M}_{t_{0}\rightarrow t_{f}}$ represents the polynomial map
that connects any initial condition,
$[\mathbf{x}(t_{0})\quad\boldsymbol{\lambda}(t_{0})]^{T}$, to the final state
at the given time, $[\mathbf{x}(t_{f})\quad\boldsymbol{\lambda}(t_{f})]^{T}$.
Therefore, the KO solution can be interpreted as a polynomial transition map
up to the selected final time $t_{f}$.

### Control Solution Using Map Inversion

The aim of the procedure is to find the value of costate at time $t_{0}$ that
drives the state form $\mathbf{x}_{0}$ to the desired final values
$\mathbf{x}_{f}$ at time $t_{f}$. However, the current KO solution, Eq. (37),
has $\boldsymbol{\lambda}(t_{0})$ as one of the variables of the polynomials.
Therefore, this map needs to be inverted, such that the initial conditions
become a function of the final variables, and not vice-versa. In order for a
map to be invertible, similarly to matrices, it needs to be square and full
rank. The polynomials in $\mathcal{M}$ map
$\mathbb{R}^{m}\rightarrow\mathbb{R}^{2m}$, where $m$ is the number of states.
This issue is overcome by augmenting the solution and introducing the map
corresponding to the identity function
$\mathbf{x}(t_{0})=\mathcal{I}(\mathbf{x}(t_{0}))$ [29, 15]. This identity map
of the state can be written highlighting the dependence on the costate
variables, where all monomials connected to the costates are zero;
$\mathbf{x}(t_{0})=\mathcal{I}_{\mathbf{x}}(\mathbf{x}(t_{0}),\boldsymbol{\lambda}(t_{0}))$.
Thus, the augmented equation can be written as

$\begin{bmatrix}\mathbf{x}(t_{f})\\\
\mathbf{x}(t_{0})\end{bmatrix}(\mathbf{x}(t_{0}),\boldsymbol{\lambda}(t_{0}))=\begin{bmatrix}\mathcal{M}\\\
\mathcal{I}_{\mathbf{x}}\end{bmatrix}(\mathbf{x}(t_{0}),\boldsymbol{\lambda}(t_{0}))$
(38)

such that the augmented polynomial map meets the requirements to be inverted.
Using inversion techniques for high order polynomials [30], the map is
inverted to obtain

$\begin{bmatrix}\boldsymbol{\lambda}(t_{0})\\\
\mathbf{x}(t_{0})\end{bmatrix}(\mathbf{x}(t_{f}),\mathbf{x}(t_{0}))=\begin{bmatrix}\mathcal{M}\\\
\mathcal{I}_{\mathbf{x}}\end{bmatrix}^{-1}(\mathbf{x}(t_{f}),\mathbf{x}(t_{0}))=\mathcal{W}(\mathbf{x}(t_{f}),\mathbf{x}(t_{0}))$
(39)

The inverted map $\mathcal{W}$ is a function of two variables: the state at
the initial time $\mathbf{x}(t_{0})$, and the values of the states at final
time $\mathbf{x}_{f}(t_{f})$. The initial costates are found by considering
the first half of the inverted $\mathcal{W}$ map, since the first $m$
polynomials are related to the control:

$\boldsymbol{\lambda}(t_{0})=\mathcal{W}_{\boldsymbol{\lambda}}(\mathbf{x}(t_{f}),\mathbf{x}(t_{0}))$
(40)

where $\mathcal{W}_{\boldsymbol{\lambda}}$ indicates the map with only the
costate polynomials selected. The optimal numerical values are obtained
through polynomial evaluation

$\boldsymbol{\lambda}_{0}=\mathcal{W}_{\boldsymbol{\lambda}}(\mathbf{x}_{f},\mathbf{x}_{0})$
(41)

where the substitutions $\mathbf{x}(t_{f})=\mathbf{x}_{f}$ and
$\mathbf{x}(t_{0})=\mathbf{x}_{0}$ guarantee to satisfy the boundary
conditions of the TPBVP. The control optimization is concluded and the time
evolution evolution of both the states, and the costates, can be found by
solving system (35) with initial conditions $\mathbf{x}(t_{0})=\mathbf{x}_{0}$
and $\boldsymbol{\lambda}(t_{0})=\boldsymbol{\lambda}_{0}$.

It is important to underline that contrary to other techniques based on map
inversion [29, 15], the KO solution does not require the evaluation of any
nominal control to perform accurately. This aspect is due to the global nature
of the KO approximation, and to the fact that the KO map of the solution is
able to work directly with state variables over deviations. Looking at Eq.
(40), the inverted map creates a function of the optimal costates on the
boundary conditions. Therefore, if we seek a different solution where either
the initial condition and/or the desired final state are changed, the optimal
control can be rapidly obtained by evaluating the map as in Eq. (41). The
Koopman approach has completely described the influence of the boundary
conditions on the optimal initial costate, such that once the Koopman solution
is obtained, any change in the application can be resolved by looking directly
at the state transition polynomial maps.

## Numerical Applications

The proposed control techniques is applied to three different numerical
applications. First, the state of the Duffing oscillator is driven to rest
from multiple configurations, displaying the control effort and the
relationship with the maneuver time. Afterwards, the Clohessy-Wiltshire (CW)
problem is presented and the equations of motion are derived. This application
is used to prove the efficiency of the new technique for rendezvous problems.
Lastly, the assumptions that lead to the classical linear CW equations are
relaxed and a new, high-order, system that better describes the relative
motion between two satellites is presented, where nonlinearites play an
important role in increasing the accuracy of the relative motion.

### The Duffing Oscillator

The application to the Duffing Oscillator system is now proposed, where the KO
inverse control technique is applied to control the state to a set of desired
values for different time requirements. The Duffing oscillator describes the
oscillations of a mass $M$ attached to a nonlinear spring, with stiffness
constant $k$, and a damper. The scalar system can therefore be described by
the following equations of motion, where $q$ represents the position, while
$p$ is the velocity;

$\displaystyle\ \dot{q}$ $\displaystyle=\dfrac{p}{M}$ (42a)
$\displaystyle\dot{p}$ $\displaystyle=-kq-ka^{2}\epsilon q^{3}$ (42b)

where $a$ is a unit transformation constant and $\epsilon$ is a small
parameter. In this scenario, the ODE parameters are defined as $a=1$, $M=1$,
$k=1$, and $\epsilon=0.001$. The system is at equilibrium at rest, for null
values of position and velocity, $(0,0)$. When perturbed, the mass starts to
oscillate in a periodic motion. After adding control to the oscillator, the
controlled system, in the form of Eqs. (35), has the following form:

$\displaystyle\ \dot{q}$ $\displaystyle=p$ (43a) $\displaystyle\dot{p}$
$\displaystyle=-q-\epsilon q^{3}-\lambda_{p}$ (43b)
$\displaystyle\dot{\lambda}_{q}$ $\displaystyle=3\epsilon
q^{2}\lambda_{p}+\lambda_{p}$ (43c) $\displaystyle\dot{\lambda}_{p}$
$\displaystyle=-\lambda_{q}$ (43d)

Figure 1: Bringing the Duffing oscillator to its equilibrium point for
different times.

At first, let us assume that the mass is moved from its equilibrium position
to $q_{0}=1$, and given no initial velocity, $p_{0}=0$. The goal of this
application is to bring the mass back at rest in its equilibrium point given a
fixed time length. Figure 1 shows the system brought to the origin for
different amounts of times. The uncontrolled systems has the mass oscillating
between -1 and 1 indefinitely, displayed in black. The energy optimal approach
has the state of the system rotating in a spiral around the equilibrium point,
at a rate that gets slower the larger the final time $t_{f}$.

Figure 2: Control magnitude requested for each time length.

Thus, the required magnitude of control effort decreases as the system is
allowed more time to reach stasis. Figure 2 reports the amount of control
applied to the system for different times, both for position and velocity,
evaluated as

$\hat{\lambda_{i}}=\dfrac{1}{t_{f}}\sqrt{\int_{0}^{t_{f}}\lambda_{i}^{2}dt}$
(44)

with $i=q,p$. The control effort monotonically decreases as time increases,
with the costate for position and velocity overlapping for large amount of
times.

Figure 3: Controlled pathways to equilibrium for $t_{f}=2$ seconds.

A second analysis is performed with a fixed required time $t_{f}=2$ seconds,
where the initial condition of the state is picked following a grid, and the
system is lead to stasis. Figure 3 shows the path to equilibrium from the 12
different initial conditions, and the relative uncontrolled behaviour in the
background. Therefore, control has been achieved given different states and
time constraints.

### The Clohessy-Wiltshire problem

The Clohessy-Wiltshire (CW) equations represent a linear approximation of the
motion between two, relatively close, satellites in orbit [24]. The problem
can be formulated using different set of coordinates [25]: the approach of
deriving the Lagrangian for the relative motion in Cartesian coordinates is
here chosen. Then, the Hamiltonian of the relative system can be evaluated
using a Legendre transformation, such that it is possible to divide the linear
part from the nonlinear contribution.

The Cartesian coordinates are the most convenient when allowing control and
simulation techniques. Let us consider the rotating Cartesian Euler-Hill
system, where the origin of the frame is set on a circular orbit of radius $a$
about the primary celestial body. Therefore, this frame is rotating with mean
motion $n=\sqrt{\dfrac{\mu}{a^{3}}}$, where $\mu$ is the gravitational
parameter. The classic Clohessy-Wiltshire equations are obtained with the
assumption of relative motion of two spacecrafts with respect to a circular
orbit. Let us define with $\mathbf{r}_{1}$ the position vector of the target
spacecraft, and with $\mathbf{r}_{2}$ the position vector of the chaser, as
shown in Figure 4, with respect to the
$\mathcal{I}=\\{\hat{I},\hat{J},\hat{K}\\}$ reference frame.

Figure 4: Relative motion rotating frame.

Thus, the relative position vector in the rotating frame
$\mathcal{R}=\\{\hat{x},\hat{y},\hat{z}\\}$ is defined as
$\mathbf{\rho}=[x,y,z]^{T}$, and the angular velocity between the two frames
$\mathcal{I}$ and $\mathcal{R}$ is $\boldsymbol{\omega}=[0,0,n]^{T}$. Under
the circular orbit assumption, $||\mathbf{r}_{1}||=a$, and the component-wise
relative velocity in the rotating frame can be calculated by

$\mathbf{v}=\boldsymbol{\omega}\times\mathbf{r}_{1}+\boldsymbol{\omega}\times\boldsymbol{\rho}+\dfrac{d^{\mathcal{R}}}{dt}\boldsymbol{\rho}=\begin{bmatrix}\dot{x}-ny\\\
\dot{y}+nx+na\\\ \dot{z}\end{bmatrix}$ (45)

The equations of motion are obtained by applying the Euler-Lagrange equation
to the Lagrangian of the system $\mathbf{L}$. The Lagrangian is evaluated by
subtracting the potential energy $\mathbf{U}$ to the kinetic energy
$\mathbf{T}$:

$\mathbf{L}=\mathbf{T}-\mathbf{U}$ (46)

Working with quantities per unit of mass, the kinetic energy is:

$\mathbf{T}=\dfrac{1}{2}||\mathbf{v}||=\dfrac{1}{2}(\dot{x}+\dot{y}+\dot{z})^{2}+n(x\dot{y}-y\dot{x}+a\dot{y})+\dfrac{1}{2}n^{2}y^{2}+\dfrac{1}{2}n^{2}x^{2}+n^{2}ax+\dfrac{1}{2}n^{2}a^{2}$
(47)

The potential energy of the chaser is the one of a spacecraft attracted by a
spherical body. This gravitational potential is expanded using Legendre
polynomials following the developments from [25]:

$\mathbf{U}=-\dfrac{\mu}{||\mathbf{r}_{2}||}=-\dfrac{\mu}{||\mathbf{r}_{1}+\boldsymbol{\rho}||}=-\dfrac{\mu}{a\sqrt{1+2\dfrac{\mathbf{r}_{1}\cdot\boldsymbol{\rho}}{a^{2}}+\dfrac{\rho^{2}}{a^{2}}}}=-\dfrac{\mu}{a}\sum_{i=0}^{\infty}\mathcal{P}_{k}(\cos\alpha)\left(\dfrac{\rho}{a}\right)^{k}$
(48)

where $\rho=||\boldsymbol{\rho}||$, and $\mathcal{P}_{k}(\cos\alpha)$ are the
Legendre polynomials. Looking back at Figure 4, angle $\alpha$ is the angle
between the position of the target and the relative position of the chaser.
Therefore, the cosine can be rewritten as

$\cos\alpha=-\dfrac{\mathbf{r}_{1}\cdot\boldsymbol{\rho}}{a\rho}=-\dfrac{x}{\rho}$
(49)

Substituting in the Lagrangian, the formulation becomes

$\mathbf{L}=\dfrac{1}{2}(\dot{x}+\dot{y}+\dot{z})^{2}+n(x\dot{y}-y\dot{x}+a\dot{y})+\dfrac{1}{2}n^{2}y^{2}+\dfrac{1}{2}n^{2}x^{2}+n^{2}ax+\dfrac{1}{2}n^{2}a^{2}+n^{2}a^{2}\sum_{k=0}^{\infty}\mathcal{P}_{k}\left(-\dfrac{x}{\rho}\right)\left(\dfrac{\rho}{a}\right)^{k}$
(50)

which, in practice, needs to be approximated by selecting the maximum order of
the terms in the summation. Therefore, it is convenient to define the set of
polynomials $\mathcal{Q}_{k}$ as

$\mathcal{Q}_{k}=n^{2}a^{2}\mathcal{P}_{k}\left(-\dfrac{x}{\rho}\right)\left(\dfrac{\rho}{a}\right)^{k}\quad\forall\quad
k=0,\dots,\infty$ (51)

such that it eases the contribution of each singular term added to the
Lagrangian by the potential for each order. Table 1 shows the contribution of
the first 6 polynomials.

Table 1:  Potential polynomials for different orders. $k$ | $\mathcal{Q}_{k}$
---|---
0 | $n^{2}a^{2}$
1 | $-n^{2}ax$
2 | $\dfrac{3}{2}n^{2}x^{2}-\dfrac{1}{2}\rho^{2}a^{2}$
3 | $-\dfrac{5}{2}\dfrac{n^{2}}{a}x^{3}+\dfrac{3}{2}\dfrac{\rho^{2}n^{2}}{a}x$
4 | $\dfrac{35}{8}\dfrac{n^{2}}{a^{2}}x^{4}-\dfrac{30}{8}\dfrac{\rho^{2}n^{2}}{a^{2}}x^{2}+3\dfrac{\rho^{4}n^{2}}{a^{2}}$
5 | $-\dfrac{63}{8}\dfrac{n^{2}}{a^{3}}x^{5}+\dfrac{70}{8}\dfrac{\rho^{2}n^{2}}{a^{3}}x^{3}-\dfrac{15}{8}\dfrac{\rho^{4}n^{2}}{a^{3}}x$
$\dots$ | $\dots$

In the Clohessy-Wiltshire equations for relative motions, only small
deviations from the reference orbit are considered, where only the first three
polynomials of the potential expansion are considered. Thus, substituting the
polynomials with $\rho=\sqrt{x^{2}+y^{2}+z^{2}}$, we obtain the low order
Lagrangian

$\mathbf{L}_{LIN}=\dfrac{1}{2}(\dot{x}+\dot{y}+\dot{z})^{2}+n(x\dot{y}-y\dot{x}+a\dot{y})+\dfrac{3}{2}n^{2}a^{2}+\dfrac{3}{2}n^{2}x^{2}-\dfrac{1}{2}n^{2}z^{2}$
(52)

where the subscript “$LIN$" indicates that this Lagrangian provides the linear
equations of motion. In fact, the classic linear CW system of ODE is obtained
by applying the Euler-Lagrange equation

$\dfrac{\partial\mathbf{L}}{\partial
x}-\dfrac{d}{dt}\left(\dfrac{\partial\mathbf{L}}{\partial\dot{x}}\right)=0$
(53)

repeated for variables $y$ and $z$. The usual equations of motion are then
obtained:

$\displaystyle\ddot{x}$ $\displaystyle=2n\dot{y}+3n^{2}x$ (54a)
$\displaystyle\ddot{y}$ $\displaystyle=-2n\dot{x}$ (54b)
$\displaystyle\ddot{z}$ $\displaystyle=n^{2}z$ (54c)

However, it is common procedure to normalize these equations and work with a
simplified notation. Therefore, all distances get normalized by the semi-major
axis $a$, which corresponds to a normalization of the rates of a factor $n$.
Under this transformation, the new variables are dimensionless and the
Lagrangian can be rewritten as

$\tilde{\mathbf{L}}_{LIN}=\dfrac{1}{2}(\dot{x}+\dot{y}+\dot{z})^{2}+x\dot{y}-y\dot{x}+\dot{y}+\dfrac{3}{2}+\dfrac{3}{2}x^{2}-\dfrac{1}{2}z^{2}$
(55)

which gives the normalized system of ODEs:

$\displaystyle\ddot{x}$ $\displaystyle=2\dot{y}+3x$ (56a)
$\displaystyle\ddot{y}$ $\displaystyle=-2\dot{x}$ (56b)
$\displaystyle\ddot{z}$ $\displaystyle=z$ (56c)

It can be noted that the evolution of the system in the $z$ direction is
completely decoupled from the $(x,y)$ plane motion and it could be solved by
considering the contribution separately. In the KO solution, the system gets
further scaled by a factor of 0.05 to best fit the evaluation of each
component of the Koopman matrix with the Galerkin method.

Therefore, consider a target spacecraft orbiting around earth with a semi-
major axis $a=6678000$ meters. In this application, a chaser satellite aims to
rendezvous with the target in one day. The relative (unscaled) initial state
of the chaser with respect to the target is

$\displaystyle\mathbf{r}_{0}$
$\displaystyle=\begin{bmatrix}-2.0772&4.5157&0\end{bmatrix}\text{km}$ (57a)
$\displaystyle\mathbf{v}_{0}$ $\displaystyle=\begin{bmatrix}-8.6074\times
10^{-5}&4.2376\times 10^{-3}&0\end{bmatrix}\text{km/s}$ (57b)

which means that the chaser is getting farther away from the target, and that
it needs to slow down and invert its relative direction of motion. The
uncontrolled motion for one day is displayed in Figure 5, where the left graph
shows the relative position of the chaser with respect to the target in the
$(x,y)$ plane, while the graph on the right portraits the relative velocity
profile.

Figure 5: Relative uncontrolled motion of the chaser with respect to the
target for one day.

The uncontrolled chaser keeps getting away from the target in a periodic
motion. Thus, control is added to the system in an energy optimal manner, such
that the controlled dynamics are written in the formulation provided by
Equations (35):

$\displaystyle\dot{x}$ $\displaystyle=v_{x}$ (58a) $\displaystyle\dot{y}$
$\displaystyle=v_{y}$ (58b) $\displaystyle\dot{z}$ $\displaystyle=v_{z}$ (58c)
$\displaystyle\dot{v}_{x}$ $\displaystyle=3n^{2}x+2nv_{y}-\lambda_{v_{x}}$
(58d) $\displaystyle\dot{v}_{y}$ $\displaystyle=-2nv_{x}-\lambda_{v_{y}}$
(58e) $\displaystyle\dot{v}_{z}$ $\displaystyle=-n^{2}z-\lambda_{v_{z}}$ (58f)
$\displaystyle\dot{\lambda}_{x}$ $\displaystyle=-3n^{2}\lambda_{v_{x}}$ (58g)
$\displaystyle\dot{\lambda}_{y}$ $\displaystyle=0$ (58h)
$\displaystyle\dot{\lambda}_{z}$ $\displaystyle=n^{2}\lambda_{v_{z}}$ (58i)
$\displaystyle\dot{\lambda}_{v_{x}}$
$\displaystyle=-\lambda_{x}+2n\lambda_{v_{y}}$ (58j)
$\displaystyle\dot{\lambda}_{v_{y}}$
$\displaystyle=-\lambda_{y}-2n\lambda_{v_{x}}$ (58k)
$\displaystyle\dot{\lambda}_{v_{z}}$ $\displaystyle=-\lambda_{z}$ (58l)

After solving the system in the KO framework, and inverting the polynomial map
of the augmented state (order three selected), the values of the initial
costates that bring the system to rendezvous are

$\displaystyle\boldsymbol{\lambda}_{\mathbf{r},0}$
$\displaystyle=\begin{bmatrix}-4.3659\times 10^{-11}&1.6400\times
10^{-13}&0\end{bmatrix}^{T}$ (59)
$\displaystyle\boldsymbol{\lambda}_{\mathbf{v},0}$
$\displaystyle=\begin{bmatrix}-1.0017\times 10^{-9}&-1.5900\times
10^{-8}&0\end{bmatrix}^{T}$ (60)

The controlled pathway of the chaser can be appreciated in Figure 6, where it
can be noted that the spacecraft slows down and changes its direction to
approach the origin with smaller oscillations each time. The velocity profile,
the right plot of the figure, shows how the chaser decreases its velocity in a
spiral, reaching stasis at the desired final time.

Figure 6: Rendezvous of the chaser to the target in one day.

The rendezvous is achieved and the control technique has been proven to work
successfully.

In this specific case that system (58) is linear, it is possible to find the
analytical solution of the inverse control problem. Indeed, let us consider
the state transition matrix (STM) $\boldsymbol{\Phi}$ such that

$\begin{bmatrix}\mathbf{x}_{f}\\\
\boldsymbol{\lambda}_{f}\end{bmatrix}=\begin{bmatrix}\boldsymbol{\Phi}_{11}&\boldsymbol{\Phi}_{12}\\\
\boldsymbol{\Phi}_{21}&\boldsymbol{\Phi}_{22}\end{bmatrix}\begin{bmatrix}\mathbf{x}_{0}\\\
\boldsymbol{\lambda}_{0}\end{bmatrix}$ (61)

where the STM is calculated as $\boldsymbol{\Phi}=\exp{(\mathbf{A}t_{f})}$,
with matrix $\mathbf{A}$ obtained by form the set of Equations (58).
Therefore, for the linear case, it is possible to solve for the analytical
solution of the initial costate as

$\boldsymbol{\lambda}_{0}=\boldsymbol{\Phi}_{12}^{-1}\left(\mathbf{x}_{f}-\boldsymbol{\Phi}_{11}\mathbf{x}_{0}\right)$
(62)

When comparing the KO solution, Equations (59) and (60), with the analytical
solution obtained through Equation (62), we obtain the following levels of
relative errors, for the position
$\epsilon_{REL}(\boldsymbol{\lambda}_{\mathbf{r},0})$, and velocity
$\epsilon_{REL}(\boldsymbol{\lambda}_{\mathbf{v},0})$;

$\displaystyle\epsilon_{REL}(\boldsymbol{\lambda}_{\mathbf{r},0})$
$\displaystyle=\begin{bmatrix}0.021&0.001&0\end{bmatrix}^{T}\%$ (63)
$\displaystyle\epsilon_{REL}(\boldsymbol{\lambda}_{\mathbf{v},0})$
$\displaystyle=\begin{bmatrix}0.840&0.033&0\end{bmatrix}^{T}\%$ (64)

proving a correct evaluation of the optimal costates. This comparison
validates the reliability of the KO methodology, that can be now applied to
the nonlinear CW dynamics.

### The High-Order Clohessy-Wiltshire Problem

The relative equations of motion studied so far consider only small deviations
of the chaser from the reference orbit. However, when the magnitude of the
relative position and velocity increases, the linear approximation of the
Clohessy-Wiltshire equations is not sufficient to fully describe the system.
Therefore, in this new application, the effects of the nonlinearities between
chaser and target are taken into account. Thus, looking back at the definition
of the Lagrangian from Equation (50), the potential component can be truncated
at orders higher than two (which gives the linear system). The resulting high-
order Clohessy-Wiltshire system, obtained via the Euler-Lagrange equations, is
fully coupled.

Looking at Table 1, for any order $k$, the polynomials $\mathcal{Q}_{k}$ have
the magnitude of the relative position vector, $\rho$, appearing only with
even powers. Indeed, when looking at the $i$-th polynomial, the monomials of
$\mathcal{Q}_{i}$ are in the form

$\displaystyle\mathcal{Q}_{i}$
$\displaystyle=q_{1}\rho^{i}+q_{2}\rho^{(i-2)}x^{2}+q_{3}\rho^{(i-4)}x^{4}+\dots+q_{n}\rho^{0}x^{i}\quad\forall
i=2,4,6,\dots$ (65) $\displaystyle\mathcal{Q}_{i}$
$\displaystyle=q_{1}\rho^{(i-1)}+q_{2}\rho^{(i-3)}x^{3}+q_{3}\rho^{(i-5)}x^{5}+\dots+q_{n}\rho^{0}x^{i}\quad\forall
i=3,5,7,\dots$ (66)

which means that any $\mathcal{Q}_{k}$ is polynomial in the variables $x$,
$y$, and $z$, since the square root of $\rho=\sqrt{x^{2}+y^{2}+z^{2}}$ is
always cancelled by the even power of its exponent. Therefore, the resulting
Lagrangian for any truncation order will always give polynomial equations of
motion from the Euler-Lagrangian equations, which are suitable to the KO
control technique. While the truncation order of the potential energy can be
chosen arbitrary large, the effects of the nonlinear terms and the accuracy of
the model reaches a plateau for large numbers. Indeed, looking again back at
Table 1, it is worth noticing the the magnitude of the coefficients of the
monomials of $\mathcal{Q}_{k}$ decreases by a factor $a$ for each additional
order. That is, for a large $k$, the contribution of the new polynomials
become negligible when compared to the low order polynomials. The constant
decrease of the magnitude of the coefficients is directly related by the size
of the target’s orbit. Since each high-order polynomial gets divided by the
semi-major axis, the effect of nonlinearities is stronger for low altitude
orbits. Conceptually, a stronger curvature of the orbit, closer to Earth, has
stronger nonlinear effects when compared to far orbits, where the local
linearization is more accurate. Therefore, the high-order Clohessy-Wiltshire
equations present their advantages and improvements in the LEO region.

Figure 7: Comparison between classic and high-order Clohessy-Wiltshire
equations.

Let us consider the initial condition previously selected by Eq. (57a) and Eq.
(57b). For the simulations proposed in this paper, the target is orbiting in
LEO at 307 km altitude. Figure 7 shows the comparison between the classic
(linear) Clohessy-Wiltshire propagation, in black lines, and the high-order
Clohessy-Wiltshire propagation, in blue. The linear propagation is analogous
at Figure 5, which correspsonds to selecting truncation order two for the
potential polynomials. The nonlinear (blue) lines are instead evaluated by
doubling the truncation order to four. Looking at the position graph on the
left, after an initial overlapping between the linear and nonlinear solution,
the blue line starts to drift apart transitioning on the left, while the black
line keeps a perfect periodic behaviour. Therefore, the high-order solution
better approximates the actual relative position of the chaser, as shown in
[24].

The KO control technique is applied to reach rendezvous between chaser and
target when the high-order dynamics is considered. In order to show the
benefits of the new approximation, an initial condition of the type

$\displaystyle\mathbf{r}_{0}$
$\displaystyle=\begin{bmatrix}0&$y$&0\end{bmatrix}\text{km},$ (67a)
$\displaystyle\mathbf{v}_{0}$
$\displaystyle=\begin{bmatrix}0&0&0\end{bmatrix}\text{km/s}$ (67b)

has been selected. This initial condition underlines the benefits of the
nonlinear systems since, in the linear dynamics, the $y$ axis is a set of
equilibrium points. Indeed, looking at system (54), any point in the $y$ axis
is stationary and the forward propagation of such initial condition leads to a
constant vector. Physically, this means that if the chaser is exactly behind
the target with null relative velocity, its relative position is fixed in
time. Such approximation holds only when the curvature of the orbit is
neglected (full linearization). Therefore, the aim of the presented
application is to achieve rendezvous when the chaser is at $y=10$ km ahead of
the target, given the constraint of performing the maneuver in 12 hours.
Figure 8 shows the resulting pathways of the chaser subject to different
dynamics. As stated above, the uncontrolled linear propagation results to a
constant integration, that is a single point in the figure. Therefore, the
uncontrolled linear solution, in green, is not visible since it is static. On
the contrary, the uncontrolled nonlinear dynamics, in black, has been
evaluated considering the first six potential polynomials. This line shows the
actual behaviour related to the uncontrolled relative motion of the
spacecraft, with the presence of the translational drift portrayed in Figure
7. The controlled pathways are reported in the figure with red and blue lines,
for the linear and the high-order dynamics, respectively. Thus, the correct
solution of the controlled chaser can be appreciated by analyzing the blue
curves. As the chaser gains velocity by increasing its $x$ components, it is
able to rapidly approach the target in the $y$ direction. After passing the
half-point, the control starts to slow the chaser down, such that it reaches
rendezvous with null velocity. The controlled linear pathway, in red, is
reported for comparison purposes. The red curve shows the behaviour of the
chaser that is propagated under linear dynamics, which fails to achieve
rendezvous. Thus, the consideration of nonlinear terms is crucial for the high
accuracy and precision evaluation of the costates of the system and the newly
developed energy optimal control technique is able to evaluate such control.

Figure 8: Rendezvous of the chaser to the target in 12 hours under different
models.

The rendezvous application has been repeated with the chaser initial condition
at different relative positions, in order to asses the capability of the KO
technique to evaluate the optimal control to approach the target from every
direction. As such, Figure 9 reports the optimal rendezvous pathways of the
chaser starting from a distance of two kilometers from the target (dotted
circle), with initial condition separated every $\pi/4$. The figure reports
solutions for rendezvous in exactly four hours. The values of the initial
relative velocity $\dot{y}$ is evaluated considering the difference in the
orbit’s tangential velocity, where the reference orbit radius $a$ is increased
by a factor $x$. This aspect best describes the relative initial condition
between chaser and target orbiting in different orbits. Therefore, Figure 9
shows that the KO energy optimal inverse control technique is able to
calculate the exact costates for the chaser approaching from every direction.
These solutions have been calculated using a high-order Lagrangian that
considers the first six polynomials in the expansion of the potential energy.

Figure 9: Rendezvous of the chaser to the target in 4 hours given different
relative positions in orbit.

The newly developed technique has been proven solid against high
nonlinearities in the dynamics. However, the accuracy of the KO control
solution comes at a cost of few limitations and drawbacks. Due to the
similarities of the KO framework to perturbation theory [23], the robustness
of the solution decreases as the relative distance between chaser and target
increases. This drawback is overcome in two separate ways. First, for
extremely large initial conditions, the linear dynamics solution can be used
as a first approach to the desired state, switching to the high-order solution
once closer. Secondly, distances can be furthermore scaled, such that the KO
solves for the costates as the initial relative position is closer to the
desired final state. The optimal values of the costates for the unscaled
system are determined by multiplying the scaled costates by the inverse of the
scaling factor. These techniques are often required to let the KO work with
the scaled states, which are bounded in the [-1,1] domain or close to the
origin, such that the Galerkin method best evaluates the projection of the
dynamics onto the Koopman eigenfunctions.

Lastly, the Koopman approximation of the solution of the system can only be so
accurate up to certain time, before diverging from the actual true pathway.
Therefore, the proposed control technique is not suitable for extremely large
time constraints, like multiple weeks in the high-order Clohessy-Wiltshire
equations.

## Conclusion

A novel energy optimal inverse control technique has been presented in this
paper. The state of the system is augmented to comprehend the costate
variables and the resulting augmented system of ordinary differential equation
is solved in the Koopman framework. Thus, the Koopman solution approximates
the time evolution of the state as a linear combination well-selected basis
functions: the Legendre polynomials have been selected in this work. The
Koopman propagation of the state constitutes the polynomial transition map of
the system. This map can be inverted to evaluate and evaluated to obtain the
initial values of the optimal costates, which lead the system to the desired
state at the final time. The solution of the TPBVP is therefore achieved
considering the separate influence of the constraints in the Koopman solution:
the time requirement is imposed by fixing the coefficients of the polynomials,
while the state requirement is solved by evaluation the Legendre polynomials
at the initial conditions.

The application to the Duffing oscillator offers the explanation of the new
theory in a simple toy problem, and it gives a nice visualization of the
solution. Afterwards, the Clohessy-Wiltshire problem has been selected to
achieve rendezvous between two satellites orbiting in close proximity. The
classic linear system of equations has been presented as a particular case
from the high-order Lagrangian, where the polynomial series expansion of the
potential energy of the system has been truncated at the second order.
Therefore, a comparison between the controlled high-order system, where the
nonlinearities have been taken onto account, with the controlled linear
equations, shows the improved accuracy of the dynamics. The addition of the
nonlinear terms accounts for the drift of the linear periodic solution form
the $y$ axis and, consequently, the optimal costate values calculated with a
large truncation order are more truthful of the real behaviour of the system.
The results show that the KO control technique can evaluate costates inverting
a polynomial map from a Lagrangian composed by the first six polynomials of
the potential. The robustness of the KO technique is tested to prove that the
chaser can approach the target from any direction.

## Acknowledgments

The authors want to acknowledge the support of this work by the Air Force’s
Office of Scientific Research under Contract Number FA9550-18-1-0115.

## References

  * Conway [2010] Conway, B. A., _Spacecraft trajectory optimization_ , Vol. 29, Cambridge University Press, 2010.
  * Jehn et al. [2012] Jehn, R., Garcia, D., Schoenmaekers, J., and Companys, V., “Trajectory design for BepiColombo based on navigation requirements,” _Journal of Aerospace Engineering_ , Vol. 4, No. 1, 2012, p. 1.
  * Haberkorn et al. [2004] Haberkorn, T., Martinon, P., and Gergaud, J., “Low thrust minimum-fuel orbital transfer: a homotopic approach,” _Journal of Guidance, Control, and Dynamics_ , Vol. 27, No. 6, 2004, pp. 1046–1060.  10.2514/1.4022.
  * Luo et al. [2014] Luo, Y., Zhang, J., and Tang, G., “Survey of orbital dynamics and control of space rendezvous,” _Chinese Journal of Aeronautics_ , Vol. 27, No. 1, 2014, pp. 1–11.  10.1016/j.cja.2013.07.042.
  * Di Mauro et al. [2018] Di Mauro, G., Lawn, M., and Bevilacqua, R., “Survey on guidance navigation and control requirements for spacecraft formation-flying missions,” _Journal of Guidance, Control, and Dynamics_ , Vol. 41, No. 3, 2018, pp. 581–602.  10.2514/1.G002868.
  * Gazzino et al. [2019] Gazzino, C., Arzelier, D., Louembet, C., Cerri, L., Pittet, C., and Losa, D., “Long-term electric-propulsion geostationary station-keeping via integer programming,” _Journal of guidance, control, and dynamics_ , Vol. 42, No. 5, 2019, pp. 976–991.  10.2514/1.G003644.
  * Betts [1998] Betts, J. T., “Survey of numerical methods for trajectory optimization,” _Journal of guidance, control, and dynamics_ , Vol. 21, No. 2, 1998, pp. 193–207.  10.2514/2.4231.
  * Bryson and Ho [2018] Bryson, A. E., and Ho, Y.-C., _Applied optimal control: optimization, estimation, and control_ , Routledge, 2018.
  * Jiang et al. [2012] Jiang, F., Baoyin, H., and Li, J., “Practical techniques for low-thrust trajectory optimization with homotopic approach,” _Journal of guidance, control, and dynamics_ , Vol. 35, No. 1, 2012, pp. 245–258.  10.2514/1.52476.
  * Malyuta et al. [2021] Malyuta, D., Yu, Y., Elango, P., and Açıkmeşe, B., “Advances in trajectory optimization for space vehicle control,” _Annual Reviews in Control_ , Vol. 52, 2021, pp. 282–315.  10.1016/j.arcontrol.2021.04.013.
  * Armellin [2021] Armellin, R., “Collision avoidance maneuver optimization with a multiple-impulse convex formulation,” _Acta Astronautica_ , Vol. 186, 2021, pp. 347–362.  10.1016/j.actaastro.2021.05.046.
  * Hofmann and Topputo [2021] Hofmann, C., and Topputo, F., “Rapid low-thrust trajectory optimization in deep space based on convex programming,” _Journal of Guidance, Control, and Dynamics_ , Vol. 44, No. 7, 2021, pp. 1379–1388.  10.2514/1.G005839.
  * Boone and McMahon [2021a] Boone, S., and McMahon, J., “Orbital Guidance Using Higher-Order State Transition Tensors,” _Journal of Guidance, Control, and Dynamics_ , Vol. 44, No. 3, 2021a, pp. 493–504.  10.2514/1.G005493.
  * Boone and McMahon [2021b] Boone, S., and McMahon, J., “Variable Time-of-Flight Spacecraft Maneuver Targeting Using State Transition Tensors,” _Journal of Guidance, Control, and Dynamics_ , Vol. 44, No. 11, 2021b, pp. 2072–2080.  10.2514/1.G005890.
  * Di Lizia et al. [2014a] Di Lizia, P., Armellin, R., Bernelli-Zazzera, F., and Berz, M., “High order optimal control of space trajectories with uncertain boundary conditions,” _Acta Astronautica_ , Vol. 93, 2014a, pp. 217–229.  10.1016/j.actaastro.2013.07.007.
  * Di Lizia et al. [2014b] Di Lizia, P., Armellin, R., Morselli, A., and Bernelli-Zazzera, F., “High order optimal feedback control of space trajectories with bounded control,” _Acta Astronautica_ , Vol. 94, No. 1, 2014b, pp. 383–394.  10.1016/j.actaastro.2013.02.011.
  * Koopman [1931] Koopman, B. O., “Hamiltonian systems and transformation in Hilbert space,” _Proceedings of the National Academy of Sciences_ , Vol. 17, No. 5, 1931, pp. 315–318.  10.1073/pnas.17.5.315.
  * Neumann [1932] Neumann, J. v., “Zur Operatorenmethode in der klassischen Mechanik,” _Annals of Mathematics_ , 1932, pp. 587–642.  10.2307/1968537.
  * Naylor and Sell [2000] Naylor, A., and Sell, G., _Linear Operator Theory in Engineering and Science_ , Applied Mathematical Sciences, Springer, 2000.
  * Prigogine [1962] Prigogine, I., _Non–equilibrium statistical mechanics_ , Dover Publications, 1962.
  * Arnas and Linares [2021] Arnas, D., and Linares, R., “Approximate analytical solution to the zonal harmonics problem using Koopman operator theory,” _Journal of Guidance, Control, and Dynamics_ , Vol. 44, No. 11, 2021, pp. 1909–1923.  10.2514/1.G005864.
  * Arnas et al. [2021] Arnas, D., Linares, R., and Alfriend, K., “An Analysis of Koopman-based Perturbation Theory Applied to the Motion About an Oblate Planet,” _31st AAS/AIAA Space Flight Mechanics Meeting_ , 2021.
  * Servadio et al. [2022] Servadio, S., Arnas, D., and Linares, R., “Dynamics Near the Three-Body Libration Points via Koopman Operator Theory,” _Journal of Guidance Control and Dynamics_ , 2022.  10.2514/1.G006519.
  * Curtis [2013] Curtis, H., _Orbital mechanics for engineering students_ , Butterworth-Heinemann, 2013.  Third Edition, pp 385.
  * Kasdin et al. [2005] Kasdin, N. J., Gurfil, P., and Kolemen, E., “Canonical modelling of relative spacecraft motion via epicyclic orbital elements,” _Celestial Mechanics and Dynamical Astronomy_ , Vol. 92, No. 4, 2005, pp. 337–370.  10.1007/s10569-004-6441-7.
  * Surana and Banaszuk [2016] Surana, A., and Banaszuk, A., “Linear observer synthesis for nonlinear systems using Koopman operator framework,” _IFAC-PapersOnLine_ , Vol. 49, No. 18, 2016, pp. 716–723.  10.1016/j.ifacol.2016.10.250.
  * Arnas and Linares [2021, accepted for publication] Arnas, D., and Linares, R., “Approximate Analytical Solution to the Zonal Harmonics Problem Using Koopman Operator Theory,” _Journal of Guidance, Control, and Dynamics_ , 2021, accepted for publication.  10.2514/1.G005864.
  * Servadio et al. [2021] Servadio, S., Arnas, D., and Linares, R., “A Koopman Operator Tutorial with Othogonal Polynomials,” _arXiv preprint arXiv:2111.07485_ , 2021.
  * Valli et al. [2013] Valli, M., Armellin, R., Di Lizia, P., and Lavagna, M., “Nonlinear mapping of uncertainties in celestial mechanics,” _Journal of Guidance, Control, and Dynamics_ , Vol. 36, No. 1, 2013, pp. 48–63.  10.2514/1.58068.
  * Berz and Hawkes [1999] Berz, M., and Hawkes, P. W., “Modern map methods in particle beam physics,” _Advances in Imaging and Electron Physics_ , Vol. 108, 1999, pp. 1–318.  10.1016/S1076-5670(08)70227-1.