# Multi-objective low-thrust spacecraft trajectory design using reachability
analysis

Nikolaus Vertovec<EMAIL_ADDRESS>Sina Ober-Blöbaum
<EMAIL_ADDRESS>Kostas Margellos<EMAIL_ADDRESS>Department of Engineering Science, University of Oxford, Oxford, OX1 3PJ,
United Kingdom  Department of Mathematics, Paderborn University, 33100,
Paderborn, Germany

###### Abstract

One of the fundamental problems in spacecraft trajectory design is finding the
optimal transfer trajectory that minimizes the propellant consumption and
transfer time simultaneously. We formulate this as a multi-objective optimal
control (MOC) problem that involves optimizing over the initial or final
state, subject to state constraints. Drawing on recent developments in
reachability analysis subject to state constraints, we show that the proposed
MOC problem can be stated as an optimization problem subject to a constraint
that involves the sub-level set of the viscosity solution of a quasi-
variational inequality. We then generalize this approach to account for more
general optimal control problems in Bolza form. We relate these problems to
the Pareto front of the developed multi-objective programs. The proposed
approach is demonstrated on two low-thrust orbital transfer problems around a
rotating asteroid.

###### keywords:

Optimal control; Reachability analysis; Multi-objective optimization; Pareto
optimality; Hamilton-Jacobi equations.

††journal: European Journal of Control

##  1 Introduction

Since the Galileo mission in 1991 we have seen a steady increase in proposed
missions to asteroids and comets, as they might hold the key to many
scientific questions including the origins of life on earth [3]. The Dawn
mission to Vesta and Ceres proved the viability of low-thrust electric
propulsion for asteroid exploration [1, 32], and it is expected that many
upcoming missions will rely on similar low-thrust propulsion. While there has
been a significant study of interplanetary transfer trajectories using low-
thrust propulsion, comparatively little research has been conducted on the
trajectory design in the vicinity of asteroids. We investigate a spacecraft
trajectory design problem around an asteroid, where the objective is to use
minimal amounts of propellant to raise an orbit while keeping flight times as
short as possible. This is a multi-objective optimal control (MOC) problem,
whereby one seeks to find the optimal way a dynamical system can perform a
certain task, while minimizing or maximizing a set of, usually contradictory
and incommensurable, objective functions [30]. Conventional optimization
techniques for spacecraft trajectory design often fall into two categories,
indirect methods, based on the calculus of variations, and direct methods,
whereby the optimal control problem is reformulated as a nonlinear program.
Direct approaches rely on parametrization and while a candidate solution is
found, there are no guarantees on the optimality of the solution. Indirect
methods, meanwhile rely on necessary analytic conditions for optimality using
Lagrange multipliers. Yet while optimality of the obtained solutions may be
guaranteed, indirect approaches, such as multiple shooting methods, rely on a
good initial approximation of the optimal trajectory [39, 41]. A third
approach is dynamic programming, whereby the optimality conditions are
formulated in continuous time based on the Hamilton-Jacobi-Bellman (HJB)
equation [39], however, it is hampered by the so-called curse of
dimensionality. Despite this, unlike direct approaches, optimality is
guaranteed, and unlike indirect approaches, the solution does not rely on an
initial approximation of the optimal trajectory. We expand on this third
approach by taking advantage of recent developments in reachability analysis.

Reachability analysis aims to find the set of points from which a target can
be reached within a given time, subject to constraints. It forms a fundamental
part of the dynamics and control literature and has been used extensively for
controller synthesis of complex systems [4, 24, 29]. In recent years we have
seen considerable research being conducted into computing reachable sets using
Hamilton-Jacobi (HJ) reachability analysis, whereby the reachable set is
derived from the viscosity solution of a HJB equation accounting also for the
presence of state constraints. Such a HJB framework is presented in [11, 25]
with more general value problems bypassing previous regularity issues
presented in [2]. In [17] an extension of the HJB framework to time-varying
targets and constraints is considered and in [15] the approach is extended to
multi-objective control problems. HJ reachability has also been successfully
applied to various aerospace applications including air traffic control [26],
the climbing problem of multi-stage launchers [8], payload optimization [9],
as well as most recently to the complete model of the ascent problem of multi-
stage launchers [7]. One of the advantages of using HJ reachability is that
the optimal trajectory can easily be constructed once the reachable set has
been computed. This makes HJ reachability attractive for problems that require
computing trajectories for various different initial states.

For the spacecraft trajectory design problem considered in this paper, there
are two possible formulations for minimizing the burnt propellant. The first
assumes that the initial mass is a free optimization variable. This approach
is common during mission design where the total required fuel budget is being
calculated. To this end, we formulate the spacecraft trajectory design problem
as a MOC problem and show that it can be equivalently stated as an
optimization problem subject to a constraint that involves the sub-level set
of a certain value function. The latter is shown to be the unique continuous
viscosity solution of a quasi-variational inequality that involves a HJB
equation. Such value functions have been defined in [11, 25] to account for
the presence of state constraints. This formulation allows characterizing the
Pareto front of the formulated MOC problem and also facilitates its
computation by means of available numerical tools.

The second formulation of the spacecraft trajectory design problem assumes a
fixed initial mass and the objective is to maximize the remaining mass after
completing a given orbital maneuver. This approach is more common when the
maneuver needs to be added to a given mission and the available fuel is non-
negotiable. This formulation had been previously investigated by the authors
in [40]. To solve this second formulation of the spacecraft trajectory design
problem, we draw on research from [15] and extend our formulation of the MOC
problem to introduce an auxiliary state, allowing us to account for arbitrary
problems in Bolza form, and, together with appropriate normalization and
approximations allowing for a reduction of the state space, greatly improving
on the method presented in [40].

Thus our contributions can be summarized as

  1. 1.

the formulation of an efficient constrained MOC problem for low-thrust
spacecraft trajectory design that optimizes only over the set of admissible
initial states and transition times,

  2. 2.

the reduction of the state space through the use of appropriate
approximations,

  3. 3.

the expansion of the MOC problem to allow for a generalization of the proposed
methodology for arbitrary multi-objective problems in Bolza form.

This paper is organized into six sections. Section 2 contains details
regarding the derivations of the spacecraft dynamics as well as the
definitions of the constraints pertaining to its behavior. In Section 3 the
optimal control problem is formulated while Section 4 describes how the set of
admissible initial states is derived from the viscosity solution of a quasi-
variational inequality. Section 5 is dedicated to the numerical computation
and case study of an orbital transfer around a rotating asteroid. Finally,
Section 6 provides concluding remarks and directions for future work.

##  2 Mathematical description and physical modeling

###  2.1 Spacecraft equations of motion

We begin by modeling the dynamics of the spacecraft. The spacecraft thrust is
defined in spherical coordinates as

$\mathbf{u}(t)\coloneqq[\boldsymbol{\alpha}(t),\boldsymbol{\delta}(t),\textbf{T}(t)]\in\mathcal{U},$
(1)

where $\boldsymbol{\alpha}(t)\in[-\pi,\pi]$ is the incidence angle,
$\boldsymbol{\delta}(t)\in[-\frac{\pi}{2},\frac{\pi}{2}]$ is the sideslip
angle and $\textbf{T}(t)\in[0,T_{\max}]$ is the variable thrust, with
$T_{\max}$ denoting the maximal allowable thrust. The Cartesian transformation
of the thrust vector is denoted by $\mathbf{u}_{x}(t)$, $\mathbf{u}_{y}(t)$
and $\mathbf{u}_{z}(t)$, respectively. The compact set
$\mathcal{U}=[-\pi,\pi]\times[-\frac{\pi}{2},\frac{\pi}{2}]\times[0,T_{\max}]$
is the set of possible control input values while
$\mathbf{u}\in\mathcal{U}_{ad}$ denotes the control policy and
$\mathcal{U}_{ad}$ denotes the set of admissible policies which is the set of
Lebesgue time measurable functions from $[-\infty,0]$ to $\mathcal{U}$. Note
that time here is considered to be non-positive to facilitate the reachability
problem exposition in Section 3. Boldface notation is used to denote time
varying functions such as trajectories and policies, while non-boldface
notation is used to denote scalars and vectors. The equations of motion of the
spacecraft around a rotating body can be expressed in 3-dimensional Euclidean
space as a second-order ordinary differential equation (see eg., [21])

$2\mathbf{\Omega}(t)\times\frac{d\textbf{R}(t)}{dt}+\mathbf{\Omega}(t)\times(\mathbf{\Omega}(t)\times\textbf{R}(t))+\frac{d\mathrm{U}(\textbf{R}(t))}{d\textbf{R}}\\\
+\frac{d\mathbf{\Omega}(t)}{dt}\times\textbf{R}(t)-\frac{\textbf{u}(t)}{m(t)}=-\frac{d^{2}\textbf{R}(t)}{dt^{2}},$
(2)

where R(t) is the radius vector from the asteroid’s center of mass to the
particle, the first and second time derivatives of R(t) are with respect to
the body-fixed coordinate system, $\mathrm{U}(\textbf{R}(t))$ is the
gravitational potential of the asteroid and $\Omega$ is the rotational angular
velocity vector of the asteroid relative to inertial space. The term
$2\mathbf{\Omega}(t)\times\frac{d\textbf{R}(t)}{dt}$ describes the Coriolis
forces, $\mathbf{\Omega}(t)\times(\mathbf{\Omega}(t)\times\textbf{R}(t))$, the
centrifugal forces and $\frac{d\mathbf{\Omega}(t)}{dt}\times\textbf{R}(t)$ the
Euler forces. We consider an asteroid rotating uniformly with constant
magnitude $\omega$ around the z-axis. Therefore, the Euler forces can be
neglected and we can express the rotation vector as $\Omega\coloneqq\omega
e_{z}$, where $e_{z}$ is the unit vector along the z-axis. Following [18], the
radius vector and its derivatives are given by

$\mathbf{R}(t)\coloneqq\begin{bmatrix}\mathbf{x}(t)\\\ \mathbf{y}(t)\\\
\mathbf{z}(t)\end{bmatrix},\quad\frac{d\textbf{R}(t)}{dt}=\begin{bmatrix}\mathbf{v}_{x}(t)\\\
\mathbf{v}_{y}(t)\\\ \mathbf{v}_{z}(t)\end{bmatrix}.$ (3)

The Coriolis and centrifugal forces (the first two terms in (2)) acting on the
spacecraft are thus

$\displaystyle
2\Omega\times\frac{d\textbf{R}(t)}{dt}=\begin{bmatrix}-2\>\mathbf{v}_{y}(t)\\\
2\omega\textbf{v}_{x}(t)\\\ 0\end{bmatrix},$ (4)
$\displaystyle\Omega\times(\Omega\times\textbf{R}(t))=\begin{bmatrix}-\omega^{2}\mathbf{x}(t)\\\
-\omega^{2}\mathbf{y}(t)\\\ 0\end{bmatrix}.$ (5)

To model the current position, velocity, and available propellant, we define
the state vector

$r\coloneqq\begin{bmatrix}x,y,z,v_{x},v_{y},v_{z},\Delta
m\end{bmatrix}^{T}\in\mathbb{R}^{7},$ (6)

where $\Delta m\in\mathbb{R}_{\geq 0}$ denotes the available propellant. The
total spacecraft mass can be expressed as $m(t)=m_{0}+\Delta m(t)$, where
$m_{0}$ denotes the dry mass of the spacecraft. Following our derivations from
(2), we can formulate the dynamics of the spacecraft, $\dot{r}=f(r,u)$, as

$f(r,u)=\begin{bmatrix}v_{x}\\\ v_{y}\\\ v_{z}\\\
U_{x}(x,y,z)+\omega^{2}x+2\omega v_{y}+\frac{u_{x}}{m_{0}+\Delta m}\\\
U_{y}(x,y,z)+\omega^{2}y-2\omega v_{x}+\frac{u_{y}}{m_{0}+\Delta m}\\\
U_{z}(x,y,z)+\frac{u_{z}}{m_{0}+\Delta m}\\\
-\frac{\sqrt{u_{x}^{2}+u_{y}^{2}+u_{z}^{2}}}{v_{\mathrm{exhaust}}}\end{bmatrix},$
(7)

where $v_{\mathrm{exhaust}}\in\mathbb{R}_{\geq 0}$ is the exhaust velocity
used to express the depletion of mass as propellant is burned, $U_{x}$,
$U_{y}$ and $U_{z}$ are the derivatives of the gravitational potential in the
direction of the unit vectors $e_{x}$, $e_{y}$ and $e_{z}$, respectively, and
where for brevity we neglect the time dependence by denoting
$r=\textbf{r}(t)$, $v_{x}=\textbf{v}_{x}(t)$, and similarly for the other
states.

###  2.2 State constraints

Since the dynamics of the spacecraft were derived for orbits in the vicinity
of the asteroid, we need to enforce state constraints on $x,y,z$. We naturally
also need to ensure that we bound the amount of propellant available. Assuming
that the burnout mass of the spacecraft is the same as the dry mass, we set
$m_{\min}\coloneqq 0$ and $m_{\max}\coloneqq m_{\mathrm{propellant}}$ and
impose $m_{\min}\leq\Delta m\leq m_{\max}$.

Due to particles ejected from the asteroid, we do not want to fall below a
circular orbit with radius $\rho\coloneqq\sqrt{x^{2}+y^{2}+z^{2}}$ of
approximately $\rho_{\min}=1$ km. Furthermore, in order for the two-body
problem under discussion to be valid and the influence of other bodies in the
solar system to be negligible, we need to stay within the sphere of influence
(SOI) of the asteroid. The SOI can be approximated as in [37] by
$\rho_{SOI}\approx a\left(\frac{M_{1}}{M_{2}}\right)^{\frac{2}{5}}$, where $a$
is the semi-major axis of the asteroid’s orbit around the sun ($1.5907\cdot
10^{8}$ km), $M_{1}$ is the Mass of the asteroid ($1.4091\cdot 10^{12}$ kg)
and $M_{2}$ is the mass of the sun ($1.9890\cdot 10^{30}$ kg). Therefore, the
sphere of influence of the asteroid is approximately
$\rho_{\max}=\rho_{SOI}\approx 8.74$ km. The set of states that satisfy the
aforementioned restrictions is given by

$\mathcal{K}\coloneqq\left\\{r\in\mathbb{R}^{7}:\rho\in[\rho_{\min},\rho_{\max}],m\in[m_{\min},m_{\max}]\right\\}.$

The target orbit that we would like to transfer to is denoted by the closed
target set $\mathcal{C}\subset\mathcal{K}$.

The initial orbit that we start at is denoted by the closed initial set
$\mathcal{I}\subset\mathcal{K}$. Note that the initial and the target orbit
restrict only the position and the velocity, but allow the mass to take any
admissible value within $[m_{\min},m_{\max}]$.

While Cartesian coordinates are useful for modeling the behavior of an object
around a rotating body, since we restrict all admissible states to lie within
the set $\mathcal{K}$, which constrains the radius $\rho$, it is more
efficient to recast our problem in spherical coordinates. To this end, define
$a_{\rho},a_{\theta},a_{\psi}$ as the transformations of

$\displaystyle a_{x}\coloneqq$
$\displaystyle{}U_{x}(x,y,z)+\omega^{2}x+2\omega v_{y},$ $\displaystyle
a_{y}\coloneqq$ $\displaystyle U_{y}(x,y,z)+\omega^{2}y-2\omega v_{x},$
$\displaystyle a_{z}\coloneqq$ $\displaystyle U_{z}(x,y,z).$

The tangential velocity in the x-y plane, $v_{t}$, and its perpendicular
counterpart, $v_{\perp}$, can then be defined as follows:

$\displaystyle\begin{bmatrix}v_{t}\\\ v_{\perp}\end{bmatrix}=$
$\displaystyle{}\begin{bmatrix}\rho\dot{\theta}\sin{\psi}\\\
\rho\dot{\psi}\end{bmatrix},$ (8) $\displaystyle\begin{bmatrix}a_{t}\\\
a_{\perp}\end{bmatrix}=$
$\displaystyle\begin{bmatrix}\sin{\psi}[a_{\theta}\rho+\dot{\theta}v_{\rho}]+\dot{\theta}v_{\perp}\cos{\psi}\\\
\rho\dot{\psi}+a_{\psi}\rho\end{bmatrix}.$ (9)

Then we can restate the system dynamics in spherical coordinates as

$\displaystyle r$
$\displaystyle=\begin{bmatrix}\rho,\theta,\psi,v_{\rho},v_{t},v_{\perp},\Delta
m\end{bmatrix}^{T}\in\mathbb{R}^{7},$ (10) $\displaystyle f(r,u)$
$\displaystyle=\begin{bmatrix}v_{\rho}\\\ \frac{v_{t}}{\rho\sin{\psi}}\\\
\frac{v_{\perp}}{\rho}\\\ a_{\rho}+\frac{T}{m_{0}+\Delta m}\cos{\alpha}\\\
a_{t}+\frac{T}{m_{0}+\Delta m}\sin{\alpha}\sin{\delta}\\\
a_{\perp}+\frac{T}{m_{0}+\Delta m}\sin{\alpha}\cos{\delta}\\\
-\frac{T}{v_{\mathrm{exhaust}}}\end{bmatrix}.$ (11)

With a slight abuse of notation, we also redefine $\mathcal{I}$,
$\mathcal{C}$, and $\mathcal{K}$ in spherical coordinates. Finally, we impose
the following assumptions on the spacecraft dynamics.

###### Assumption 2.1

For every $r\in\mathcal{K}$ the set $\Big{\\{}f(r,u):u\in\mathcal{U}\Big{\\}}$
is a compact convex subset of $\mathbb{R}^{7}$.

###### Assumption 2.2

$f:\mathbb{R}^{7}\times\mathcal{U}\rightarrow\mathbb{R}^{7}$ is bounded and
there exists an $L_{f}>0$ such that for every $u_{1},u_{2}\in\mathcal{U}$,

$||f(r_{1},u_{1})-f(r_{2},u_{2})||\leq L_{f}||r_{1}-r_{2}||.$

Due to norm equivalence, the choice of norm is irrelevant and not further
discussed. Using Assumptions 2.1 and 2.2, for any control policy
$\textbf{u}\in\mathcal{U}_{ad}$, any initial state $r_{0}\in\mathcal{K}$ and
transfer time $t_{f}>0$, the system admits a unique, absolutely continuous
solution on $[-t_{f},0]$ (see [34]).

##  3 Problem Statement

###  3.1 Multi-objective optimal control problem

Having defined the system dynamics, we are now in a position to discuss how to
find trajectories that start on an initial orbit, $\mathcal{I}$, and take the
spacecraft to some final orbit, $\mathcal{C}$. Additionally, the objective is
to keep the flight time and required propellant as small as possible. Thus,
the multi-objective optimal control problem can be formulated as a
minimization problem whereby the first goal is to minimize the required
propellant, $\Delta m$, and the second is to minimize the required time for
the orbit change, i.e., the transfer time, denoted by $t_{f}$. The trajectory,
$\mathbf{r}$, which is the solution of (11), belongs to the Sobolev space
$\mathbb{W}^{1,1}(\mathbb{R}^{7})$. The set of trajectory-control pairs on
$[-t_{f},0]$ starting at $r_{0}$ with transfer time $t_{f}$ is denoted as:

$\Pi_{r_{0},t_{f}}\coloneqq\big{\\{}(\mathbf{r},\mathbf{u}):\dot{\mathbf{r}}(t)=f(\mathbf{r}(t),\mathbf{u}(t)),\quad\forall
t\in[-t_{f},0];\\\
\mathbf{r}(-t_{f})=r_{0}\big{\\}}\subset\mathbb{W}^{1,1}(\mathbb{R}^{7})\times\mathcal{U}_{ad}.$

Note that as in [14] we adopt the convention that $0$ denotes the terminal
time hence the transfer time $t_{f}$ denotes the time duration. Under
Assumption 2.1 and by Filippov’s Theorem [23, pg.121], we can conclude, that
$\Pi_{r_{0},t_{f}}$ is compact.

###### Remark 3.3

For similar applications as the one discussed in this paper, where Assumption
2.1 might not hold, we refer to [7] where convexification of the dynamics is
considered in order to ensure that the set of absolutely continuous solutions
of the problem is closed.

The set of admissible (in the sense of satisfying the state constraints)
trajectory-control pairs on $[-t_{f},0]$ starting at $r_{0}$ with transfer
time $t_{f}$ is denoted as:

$\Pi_{r_{0},t_{f}}^{\mathcal{K},\mathcal{C}}\coloneqq\big{\\{}(\mathbf{r},\mathbf{u})\in\Pi_{r_{0},t_{f}}:\mathbf{r}(t)\in\mathcal{K},\quad\forall
t\in[-t_{f},0];\\\
\mathbf{r}(0)\in\mathcal{C}\big{\\}}\subset\mathbb{W}^{1,1}(\mathbb{R}^{7})\times\mathcal{U}_{ad}.$

Finally, the set of admissible initial state and transfer time pairs is
denoted as

$\Gamma\coloneqq\big{\\{}(r_{0},t_{f})\in\mathbb{R}^{7}\times[0,+\infty)\quad\text{such
that}\quad\Pi_{r_{0},t_{f}}^{\mathcal{K},\mathcal{C}}\neq\emptyset\big{\\}}.$

For a given initial state $r_{0}\in\mathbb{R}^{7}$ and transfer time
$t_{f}\in[0,+\infty)$, we can define the cost functions as
$J_{1}(r_{0},t_{f})\coloneqq\Delta m$ and $J_{2}(r_{0},t_{f})\coloneqq t_{f}$,
where $\Delta m$ is the $7$-th element of the state vector $r_{0}$. The
2-dimensional objective function
$J:\mathbb{R}^{7}\times[0,+\infty)\rightarrow\mathbb{R}^{2}$ can then be
written as

$J(r_{0},t_{f})\coloneqq\left[J_{1}(r_{0},t_{f}),J_{2}(r_{0},t_{f})\right]^{T}.$
(12)

We are now in a position to formulate the multi-objective optimal control
problem under study as

$\displaystyle\operatorname*{minimize}_{(r_{0},t_{f})\in\mathcal{I}\times[0,\infty)}$
$\displaystyle J(r_{0},t_{f})$ (13) subject to
$\displaystyle(r_{0},t_{f})\in\Gamma$

###  3.2 Pareto optimality

The solution of (13) in general does not consist of a single isolated point,
but rather a set of optimal compromises between the objectives $J_{1}$ and
$J_{2}$ [27].

###### Definition 3.4

A solution $(r_{0},t_{f})$ is considered Pareto optimal if
$\nexists(\hat{r_{0}},\hat{t}_{f})\in\Gamma$ such that
$J(\hat{r_{0}},\hat{t}_{f})<J(r_{0},t_{f})$,

where a vector $a$ is considered less than $b$ (denoted $a<b$) if for every
element $a_{i}$ and $b_{i}$ the relation $a_{i}<b_{i}$ holds. The relations
$\leq,\geq,>$ are defined in an analogous way. Following Definition 3.4, a
solution $(r_{0},t_{f})$ is considered Pareto optimal if it is not possible to
improve all its performance metrics $J_{1}(r_{0},t_{f}),J_{2}(r_{0},t_{f})$
simultaneously. The set of Pareto optimal solutions is called the Pareto set
$\mathcal{P}_{S}$, while its image is the Pareto front $\mathcal{P}_{F}$.
Therefore, the solution of (13), i.e., the set of minimizing $(r_{0},t_{f})$
pairs, is the desired Pareto set, while the cost function corresponding to the
minimizing $(r_{0},t_{f})$ pairs is the Pareto front.

To allow for mission designers to determine a compromise between minimizing
required propellant and transfer times, we wish to compute the Pareto front.
However, the unconventional constraint in (13) ensuring a solution
$(r_{0},t_{f})$ is feasible, prevents us from solving (13) with standard MOC
problem solvers. Therefore, we will next discuss how we can recast the
constraint $(r_{0},t_{f})\in\Gamma$ to a standard nonlinear inequality
constraint, which will allow us to compute the Pareto front by means of
conventional MOC problem solvers.

##  4 Solution to multi-objective optimal control problems

To find an equivalent formulation for the constraint $(r_{0},t_{f})\in\Gamma$,
let $g(r)$ and $\nu(r)$ be two Lipschitz functions (with Lipschitz constants
$L_{g}$ and $L_{\nu}$, respectively) chosen such that

$\displaystyle g(r)\leq$ $\displaystyle{}0\iff r\in\mathcal{K},$
$\displaystyle\nu(r)\leq$ $\displaystyle 0\iff r\in\mathcal{C}.$

This can be achieved by choosing $g(r)$ and $\nu(r)$ as the signed distance to
the set $\mathcal{K}$ and $\mathcal{C}$, respectively.

Next, we consider the value function $\omega$:

$\omega(r_{0},t_{f})\coloneqq\inf_{(\mathbf{r},\textbf{u})\in\Pi_{r_{0},t_{f}}}\Big{\\{}\nu(\mathbf{r}(0))\bigvee\max_{\tau\in[-t_{f},0]}g(\mathbf{r}(\tau))\Big{\\}},$
(14)

where $a\bigvee b$ denotes $\max(a,b)$. We are now in a position to use the
value function to decide if, for a given initial state and transfer time,
there exists a corresponding admissible trajectory. Thus we can introduce an
equivalent formulation of (13).

###### Theorem 4.5

The constrained MOC problem, (13), is equivalent to

$\displaystyle\operatorname*{minimize}_{(r_{0},t_{f})\in\mathcal{I}\times[0,\infty)}$
$\displaystyle J(r_{0},t_{f})$ (15) subject to
$\displaystyle\omega(r_{0},t_{f})\leq 0.$

###### Proof 1

We show that $(r_{0},t_{f})\in\Gamma\iff\omega(r_{0},t_{f})\leq 0$.

Case A: Consider $(r_{0},t_{f})\in\Gamma$. For the sake of contradiction
assume that $\omega(r_{0},t_{f})>0$. This then implies that for all
$(\mathbf{r},\mathbf{u})\in\Pi_{r_{0},t_{f}}^{\mathcal{K},\mathcal{C}}$ either
$\nu(\mathbf{r}(0))>0\iff\mathbf{r}(0)\notin\mathcal{C}$ or there exists
$\tau\in[-t_{f},0]$ such that
$g(\mathbf{r}(t))>0\iff\mathbf{r}(t)\notin\mathcal{K}$. This contradicts the
fact that $(r_{0},t_{f})\in\Gamma$ establishing that $(r_{0},t_{f})\in\Gamma$
implies $\omega(r_{0},t_{f})\leq 0$

Case B: Consider $(r_{0},t_{f})\in\mathcal{I}\times[0,\infty)$, such that
$\omega(r_{0},t_{f})\leq 0$. Under Assumption 2.1, applying Weierstrass’
Theorem on the existence of minima for compact sets [2, 6], we can conclude
that the infimum over $\Pi_{r_{0},t_{f}}$ exists, and thus,
$\omega(r_{0},t_{f})\leq 0$ implies the existence of a trajectory-control pair
$(\mathbf{r},\mathbf{u})\in\Pi_{r_{0},t_{f}}$, such that for all
$t\in[-t_{f},0],\>g(\mathbf{r}(t))\leq 0$ and $\nu(\mathbf{r}(0))\leq 0$. By
definition of the function $g$ and $\nu$, we thus have
$\mathbf{r}(t)\in\mathcal{K}$ for all $t\in[-t_{f},0]$ and
$\mathbf{r}(0)\in\mathcal{C}$, which in turns implies
$(\mathbf{r},\mathbf{u})\in\Pi_{r_{0},t_{f}}^{\mathcal{K},\mathcal{C}}$.
Therefore,

$\omega(r_{0},t_{f})\leq 0\Rightarrow(r_{0},t_{f})\in\Gamma,$

thus concluding the proof.

Theorem 4.5 implies that the Pareto front can be computed from the solution of
(15). To achieve this we discuss how to compute $\omega$.

###  4.1 Value function computation

To begin to discuss how $\omega$ can be obtained, we introduce the Hamiltonian
$H:\mathbb{R}^{7}\times\mathbb{R}^{7}\rightarrow\mathbb{R}$,

$H(r,q)\coloneqq-\min_{u\in\mathcal{U}}\left(q^{T}f(r,u)\right),$ (16)

where $q\in\mathbb{R}^{7}$ is the costate vector.

###### Theorem 4.6

The value function $\omega$ is the unique continuous viscosity solution of the
following quasi-variational inequality

$\begin{cases}\begin{aligned}
0=\max\big{\\{}g(r)-\omega(r,t),\partial_{t}\omega+H(r,\nabla_{r}\omega)\big{\\}}\\\
\text{for all}\quad t\in[0,\infty),\>r\in\mathbb{R}^{7},\end{aligned}\\\
\begin{aligned} \omega(r,0)={}&\left(\nu(r)\bigvee g(r)\right)\quad\text{for
all}\quad r\in\mathbb{R}^{7},\end{aligned}\end{cases}$

Since the Dynamic Programming Principle [6] holds for $h\in[0,t_{f}]$,
$(r_{0},t_{f})\in\mathcal{K}\times\mathbb{R}$ with $t_{f}\geq 0$:

$\omega(r_{0},t_{f})=\inf_{(\mathbf{r},\textbf{u})\in\Pi_{r_{0},t_{f}}}\Big{\\{}\omega(\mathbf{r}(h-t_{f}),t_{f}-h)\bigvee\max_{s\in[-t_{f},h-t_{f}]}g(\mathbf{r}(s))\Big{\\}},$

the proof of Theorem 4.6 follows standard arguments for viscosity solutions,
as shown in [2, 25]. Note that the infimum should be understood to be over the
restriction of $\Pi_{r_{0},t_{f}}$ over $[-t_{f},h-t_{f}]$.

In order to solve the quasi-variational inequality in Theorem 4.6, we employ a
finite differences scheme. As in [11, 25], a consequence following from
Theorem 4.6 is the Lipschitz continuity of the value function.

###### Proposition 4.7

The value function $\omega$ is Lipschitz continuous.

The proof is similar to the more general proof of Proposition 4.13 that is
introduced in the sequel. Proposition 4.7 allows us to make statements about
the discrete-continuous error estimate, for which we refer to Theorem 5 in
[11], as well as the convergence of the value function, for which we refer to
Proposition 6 in [16].

The Hamiltonian admits an explicit form. To this end, consider the term

$C(r,q)\coloneqq
q_{1}v_{\rho}+q_{2}v_{t}+q_{3}v_{\perp}+q_{4}a_{\rho}+q_{5}a_{t}+q_{6}a_{\perp}.$

Then we can write the Hamiltonian as

$H(r,q)\coloneqq-\min_{u\in\mathcal{U}}\Big{(}\frac{T}{m_{0}+\Delta
m}\big{(}q_{4}\cos{\alpha}\\\
+\sin{\alpha}\left(q_{5}\sin{\delta}+q_{6}\cos{\delta}\right)\big{)}-q_{7}\frac{T}{v_{\mathrm{exhaust}}}\Big{)}-C(r,q).$
(17)

As $T$ is always positive, the thrust angles can be optimized separately; see
Appendix for more details. To this end, the derivation of $\alpha^{*}$ and
$\delta^{*}$ follows a similar procedure as in [7] and is listed in the
Appendix. After applying the optimal thrust angles, the Hamiltonian becomes
affine in $T$ and we are able to find the optimal thrust magnitude. In
particular, $H(r,q)$ becomes

$H(r,q)=-\min_{T\in[0,T_{\max}]}\Big{(}-\frac{T}{m_{0}+\Delta
m}\sqrt{q_{4}^{2}+q_{5}^{2}+q_{6}^{2}}-q_{7}\frac{T}{v_{\mathrm{exhaust}}}\Big{)}-C(r,q)$
(18) $\Rightarrow
T^{*}\coloneqq\begin{cases}T_{\max}\quad\text{if}\>\frac{q_{7}}{v_{\mathrm{exhaust}}}+\frac{\sqrt{q_{4}^{2}+q_{5}^{2}+q_{6}^{2}}}{m_{0}+\Delta
m}\geq 0\\\ 0\quad\text{otherwise}\end{cases}.$ (19)

###### Remark 4.8

The singular control case,
$\frac{q_{7}}{v_{\mathrm{exhaust}}}+\frac{\sqrt{q_{4}^{2}+q_{5}^{2}+q_{6}^{2}}}{m_{0}+\Delta
m}=0$, is negligible for the consideration of the optimal Hamiltonian. For the
optimal control computation during the trajectory calculation, we have
numerically investigated the occurrence of singular arcs and found no
instances in which $\frac{\sqrt{q_{4}^{2}+q_{5}^{2}+q_{6}^{2}}}{m_{0}+\Delta
m}=-\frac{q_{7}}{v_{\mathrm{exhaust}}}$ over an extended interval, thus we do
not further consider the singular control case.

Finally, applying $T^{*}$ and rewriting the minimum as the maximum of the
negation of the associated function, the Hamiltonian takes the following
analytic form

$H(r,q)=-C(r,q)+\max{\Big{(}q_{7}\frac{T_{\max}}{v_{\mathrm{exhaust}}}+\frac{T_{\max}}{m_{0}+\Delta
m}\sqrt{q_{4}^{2}+q_{5}^{2}+q_{6}^{2}},0\Big{)}}.$

###  4.2 Extension to problems in Bolza form

We will now generalize our approach to problems where the objective functions
do not rely only on the initial state and are written in Bolza form. To
achieve this we need to introduce auxiliary states. As in [2, 15], we show how
problems in Bolza form are reformulated into Mayer form, and then show how
problems in Mayer form are solved in a similar fashion as in Section 3.
Consider the $p$-dimensional objective function defined as:

$J_{\mathrm{Bolza}}(\mathbf{r},\mathbf{u},t_{f})\coloneqq
J_{\mathrm{t}}(\mathbf{r}(0))+\\\
\int_{-t_{f}}^{0}J_{\mathrm{r}}(\mathbf{r}(s),\mathbf{u}(s))ds,$ (20)

where $J_{\mathrm{t}}$ denotes the terminal cost and $J_{\mathrm{r}}$ denotes
the running cost. We impose the following assumptions, as in [15].

###### Assumption 4.9

$J_{\mathrm{t}}$ is locally Lipschitz continuous on $\mathbb{R}^{7}$ with
Lipschitz constant $L_{\mathrm{t}}(R)$ for every neighborhood
$R\subset\mathbb{R}^{7}$.

###### Assumption 4.10

$J_{\mathrm{r}}$ is continuous on $\mathbb{R}^{7}\times\mathcal{U}$. Moreover,
$J_{\mathrm{r}}$ is locally Lipschitz continuous on the first variable with
Lipschitz constant $L_{\mathrm{r}}(R)$ for every neighborhood
$R\subset\mathbb{R}^{7}$.

###### Remark 4.11

To ease notation, we omit the dependents on the neighborhood for the Lipschitz
constants and instead assume the existence of a global Lipschitz constant
$L_{\mathrm{t}}$ and $L_{\mathrm{r}}$, respectively.

Next, we define an auxiliary state $z\in\mathbb{R}^{p}$ as:

$\begin{cases}\dot{\mathbf{z}}(s)=-J_{\mathrm{r}}(\mathbf{r}(s),\mathbf{u}(s)),&\forall
s\in[-t_{f},0]\\\ \mathbf{z}(0)=z_{0},\end{cases}$ (21)

where $z_{0}$ becomes an optimization parameter and
$\mathbf{z}\in\mathbb{W}^{1,1}(\mathbb{R}^{p})$. The auxiliary state captures
the cumulative running cost and thus is treated as an additional state. In the
same manner we previously ensured a trajectory, $\mathbf{r}$, stayed within
the set $\mathcal{K}$, we bound $\mathbf{z}$ and ensure that the integrated
running cost, added with the terminal cost, stays below some value $z_{0}$. To
capture all possible trajectories, we introduce the set

$\mathcal{Z}_{r_{0},t_{f},z_{0}}\coloneqq\big{\\{}(\mathbf{r},\mathbf{u},\mathbf{z}):(\mathbf{r},\mathbf{u})\in\Pi_{r_{0},t_{f}};\dot{\mathbf{z}}(s)=\\\
-J_{\mathrm{r}}(\mathbf{r}(s),\mathbf{u}(s)),\forall
s\in[-t_{f},0];\mathbf{z}(0)=z_{0}\big{\\}},$ (22)

and make the following assumption.

###### Assumption 4.12

For every $r\in\mathbb{R}^{7}$,the set

$\left\\{\begin{bmatrix}f(r,u)\\\
-J_{\mathrm{r}}(r,u)\end{bmatrix}:u\in\mathcal{U},\right\\}$

is a compact convex subset of $\mathbb{R}^{7}\times\mathbb{R}^{p}$.

We now introduce the auxiliary value function $\vartheta$:

$\vartheta(r_{0},t_{f},z_{0})\coloneqq\inf_{(\mathbf{r},\textbf{u},\mathbf{z})\in\mathcal{Z}_{r_{0},t_{f},z_{0}}}\Big{\\{}\\\
\bigvee_{i}\left[J^{i}_{\mathrm{t}}(\mathbf{r}(0))-\mathbf{z}^{i}(-t_{f})\right]\bigvee\nu(\mathbf{r}(0))\bigvee\max_{s\in[-t_{f},0]}g(\mathbf{r}(s))\Big{\\}},$
(23)

where $\bigvee_{i}x^{i}$ denotes the maximum element of the vector $x$. As
with $\omega$, the term $\max_{s\in[-t_{f},0]}g(\mathbf{r}(s))$ and
$\nu(\mathbf{r}(0))$ ensures that any trajectory $\mathbf{r}$ remains in
$\mathcal{K}$ and terminates in $\mathcal{C}$. The additional term
$J^{i}_{\mathrm{t}}(\mathbf{r}(0))-\mathbf{z}^{i}(-t_{f})$ ensures that the
integrated running cost, $J_{\mathrm{r}}$, combined with the terminal cost,
$J_{\mathrm{t}}$, never grows larger than $z_{0}$. Thus, in addition to
ensuring that $(\mathbf{r},\mathbf{u})$ are admissible trajectory control
pairs, the sub-zero level set of $\vartheta$ bounds the terminal and
integrated running cost. Therefore,

$\vartheta(r_{0},t_{f},z_{0})\leq
0\iff\Big{[}\exists(\mathbf{r},\mathbf{u})\in\Pi_{r_{0},t_{f}}^{\mathcal{K},\mathcal{C}},J_{\mathrm{Bolza}}(\mathbf{r},\mathbf{u},t_{f})\leq
z_{0}\Big{]}.$ (24)

We are now in a position to introduce the generalized multi-objective optimal
control problem for objective functions in Bolza form:

$\displaystyle\operatorname*{minimize}_{(r_{0},t_{f})\in\mathcal{I}\times[0,\infty)}$
$\displaystyle z_{0}$ (25) subject to
$\displaystyle\vartheta(r_{0},t_{f},z_{0})\leq 0,$

where $z_{0}$ represents an upper bound for the term
$J_{\mathrm{Bolza}}(\mathbf{r},\mathbf{u},t_{f})$, without explicit knowledge
of $\mathbf{r}$ or $\mathbf{u}$.

The generalized value function can again be obtained as the unique continuous
viscosity solution of a quasi-variational inequality

$\begin{cases}\begin{aligned}
\max\big{\\{}g(r)-&{}\vartheta(r,t,z),\partial_{t}\vartheta+H(r,\nabla_{r}\vartheta,\nabla_{z}\vartheta)\big{\\}}=0\\\
&\text{for
all}\>t\in[0,\infty),\>r\in\mathbb{R}^{7},z\in\mathbb{R}^{p},\end{aligned}\\\
\begin{aligned}
\vartheta(r,0,z)=&{}\bigvee_{i}\left[J^{i}_{\mathrm{t}}(r)-z^{i}\right]\bigvee\nu(r)\bigvee
g(r)\\\ &\text{for all}\>r\in\mathbb{R}^{7},\end{aligned}\end{cases}$ (26)

where the Hamiltonian is defined as

$H(r,q_{r},q_{z})\coloneqq\min_{u\in\mathcal{U}}\left(q_{r}^{T}f(r,u)-q_{z}^{T}J_{\mathrm{r}}(r,u)\right).$

###### Proposition 4.13

The value function $\vartheta$ is Lipschitz continuous.

The proof can be found in the Appendix. Proposition 4.13 can be used to show
that a numerical solution of (26) (in the viscosity sense) can always be
determined.

Under Assumptions 2.2, 4.9, 4.10 and 4.12, by Filippov’s Theorem [23], the
problem (23) admits an optimal solution, which implies the existence of an
admissible $\mathbf{r}$ and $\mathbf{z}$ [15]. This yields the following
relationship due to (21)

$\bigvee_{i}\left[J^{i}_{\mathrm{t}}(\mathbf{r}(0))-\mathbf{z}^{i}(-t_{f})\right]=\\\
\bigvee_{i}\left[J^{i}_{\mathrm{t}}(\mathbf{r}(0))+\int_{-t_{f}}^{0}J^{i}_{\mathrm{r}}(\mathbf{r}(s),\mathbf{u}(s))ds-z^{i}_{0}\right].$
(27)

##  5 Numerical Approximation and Results

We will now discuss how the value functions can be obtained numerically, prior
to discussing how the spacecraft trajectory design problem is solved.
Following Proposition 4.7, a numerical solution to (26) can be found. To this
end, we employ the Level Set Methods toolbox of [28]. For the computation of
$\omega$, we use a Lax-Friedrich Hamiltonian

$\mathcal{H}(r,p^{-},p^{+})=H(r,\frac{p^{-}+p^{+}}{2})-\sum_{k=1}^{7}\frac{\alpha_{k}}{2}(p^{+}-p^{-}),$
(28)

where $p^{+}$ and $p^{-}$ are the right and left derivatives computed using an
appropriate fifth-order weighted essentially non-oscillatory (WENO) scheme.
The Lax-Friedrich Hamiltonian consists of an analytic expression of the
Hamiltonian (derived previously), as well as a dissipation term, which is
scaled by the dissipation coefficients $\alpha_{k}$. The dissipation
coefficients $\alpha_{k}$ needs to satisfy

$\alpha_{k}\geq\left|\frac{\partial H}{\partial q_{k}}\right|.$ (29)

Since we need $\alpha_{k}$ to serve as an upper bound, we consider the control
input that maximizes the Hamiltonian:

$\left|\frac{\partial H}{\partial q_{k}}\right|=\Big{|}-\frac{\partial
C(r,q)}{\partial q_{k}}+\frac{\partial}{\partial
q_{k}}\max\Big{(}q_{7}\frac{T_{\max}}{v_{\mathrm{exhaust}}}\\\
+\frac{T_{\max}}{m_{0}+\Delta
m}\sqrt{q_{4}^{2}+q_{5}^{2}+q_{6}^{2}},0\Big{)}\Big{|}$ (30)
$\alpha_{k}=\begin{cases}|v_{\rho}|&k=1\\\
|\frac{v_{t}}{\rho\sin{\psi}}|&k=2\\\ |\frac{v_{\perp}}{\rho}|&k=3\\\
|a_{\rho}-\max(0,\frac{T_{\max}q_{3}}{(m_{0}+\Delta
m)\sqrt{q_{3}^{2}+q_{4}^{2}+q_{5}^{2}}})|&k=4\\\
|a_{t}-\max(0,\frac{T_{\max}q_{4}}{(m_{0}+\Delta
m)\sqrt{q_{3}^{2}+q_{4}^{2}+q_{5}^{2}}})|&k=5\\\
|a_{\perp}-\max(0,\frac{T_{\max}q_{5}}{(m_{0}+\Delta
m)\sqrt{q_{3}^{2}+q_{4}^{2}+q_{5}^{2}}})|&k=6\\\
\frac{T_{\max}}{v_{\mathrm{exhaust}}}&k=7,\end{cases}$ (31)

For a further discussion of the Lax-Friedrich Hamiltonian and WENO scheme, we
refer to [31], while for a discussion of the convergence of $\omega$ and the
derivation of a necessary Courant-Friedrichs-Lewy condition, we refer to [11,
19, 28].

###  5.1 Implementation

To illustrate the theoretical results of the previous sections, we consider a
spacecraft on an initial near circular orbit around asteroid Castalia 4769.
The goal is to compute an efficient transfer trajectory that raises the
spacecraft to a stable orbit at an altitude of $6117.5$ m above the asteroid.
For the derivation of a stable orbit around Castalia 4769 we refer to [22] and
references therein. The gravity of Castalia 4769 was modeled by means of a
spherical harmonic expansion as discussed in [20, 35, 36]. Even though the
proposed theoretical framework allows us to tackle problems of any state
dimension, the available numerical tools and computational power limit us to
only study the lower-dimensional planar case of the application. We,
therefore, omit the states $\psi$ and $v_{\perp}$ to consider only

$r=\begin{bmatrix}\rho,\theta,v_{\rho},v_{t},\Delta
m\end{bmatrix}^{T}\in\mathbb{R}^{5}.$ (32)

To avoid ill-conditioning when solving the HJB equation, the state vector is
normalized using the constants introduced in Table 1. This results in the
following dynamics:

$f(r,u)=\begin{bmatrix}v_{\rho}\\\ \frac{v_{t}}{\rho}\\\
a_{\rho}+\frac{cT}{m_{0}+\Delta m}\cos{\alpha}\\\ a_{t}+\frac{cT}{m_{0}+\Delta
m}\sin{\alpha}\\\ -\frac{cT}{v_{\mathrm{exhaust}}}\end{bmatrix},$ (33)

where $c=T_{\max}\rho_{0}/(m_{0}V_{0}^{2})$ is a normalization constant.
$\rho_{0}$, $m_{0}$, and $v_{0}$ denote the initial radius, mass, and velocity
of the initial orbit, respectively.

Table 1: Normalization Scale | Values
---|---
Distance |  Initial radius $\rho_{0}$
Velocity |  $2$ m/s
Time | $\rho_{0}/v_{0}$
Mass |  $1$ kg
Force |  Maximum thrust $T_{\max}$

The optimal control policy and trajectory
$(\mathbf{r},\mathbf{u})\in\Pi_{r_{0},t_{f}}^{\mathcal{K},\mathcal{C}}$ can be
constructed efficiently using the numerical approximation of $\omega$. For a
given $N\in\mathbb{N}$ we consider the time step $h=\frac{1}{N}$ and a uniform
grid of $[-t_{f},0]$ with spacing $s^{k}=\frac{k}{N}$. Let us define the state
$\\{r^{k}\\}_{k=0}^{N}$ and control $\\{u^{k}\\}_{k=0}^{N-1}$ for the
numerical approximation of the optimal trajectory and control policy. Setting
$r_{0}$ as the initial orbit, we proceed by iteratively computing the control
value

$u^{k}(r^{k})\in\operatorname*{arg\,min}_{u\in\mathcal{U}}\omega(r^{k}+hf(r^{k},u),s^{k})\bigvee
g(r^{k}).$

For a given $\omega(r^{k},s^{k})$, this is done by numerically taking the
partial derivatives along each grid direction to estimate the costate vector,
$q$, and then determining the optimal control values as the minimizer of the
Hamiltonian, $H$. After $u^{k}$ is determined we compute $r^{k+1}$ using the
Matlab ode113 function, a variable-order Adams-Bashforth-Moulton method of
order 1 to 13 [38], and increment $k$. For the implementation, we discretized
the interval $[-t_{f},0]$ using $N=4000$ grid points.

Figure 1: Gravitational acceleration comparing $U_{\rho}$ and $U_{\theta}$
around Castalia 4769.

As shown in Figure 1, when considering orbits further than $4$ km away from
the surface of the asteroid, the variation of the gravitational acceleration
along $\theta$ becomes negligible. It is, therefore, possible to approximate
the gravitational terms in spherical coordinates as

$\displaystyle U_{\rho}(\rho,\theta)$ $\displaystyle\approx U_{\rho}(\rho)$
$\displaystyle U_{\theta}(\rho,\theta)$ $\displaystyle\approx 0.$

Using this approximation makes $a_{\rho}$ and $a_{t}$ independent of $\theta$.
This allows us to omit a grid dimension while numerically solving the quasi-
variational inequality in Theorem 4.6, greatly reducing the computational
cost. Thus the final set of states used for the computation of the value
function is

$r=\begin{bmatrix}\rho,v_{\rho},v_{t},\Delta
m\end{bmatrix}^{T}\in\mathbb{R}^{4}.$ (34)

During the calculation of the optimal trajectory, $\theta$ can easily be
reconstructed by forward integrating the dynamics at each time step $s^{k}$,
i.e.,

$\theta^{k+1}=\theta^{k}+\int_{s^{k}}^{s^{k+1}}\frac{v_{t}^{k}}{\rho^{k}}ds.$
(35)

###  5.2 Simulation results

The spacecraft is modeled with $750$ kg of dry mass, $600$ mN of maximum
thrust, and an exhaust velocity of $40$ km/s. Using an initial orbit with
radius $5.1$ km and tangential velocity of $-2.4$ m/s, we are able to compute
the numerical approximation of $\omega$ using the spatial grid described in
Table 2 in combination with a temporal grid using $N=1000$ grid points. The
propagation of the zero level set over time is shown in Figure 2 and 3. Since
the temporal grid for the calculation of the value function is more coarse
than that used for the trajectory calculation, we need to interpolate the
value function while computing the final trajectory.

Table 2: Spatial grid configuration | $\rho$ | $v_{\rho}$ | $v_{t}$ | $\Delta m$
---|---|---|---|---
Points | 50 | 40 | 40 | 32
Spacing | 0.0045 | 0.0416 | 0.0088 | 0.0067
Minimum | 0.8067 | -0.2495 | -1.4154 | -0.0533
Maximum | 1.0270 | 1.3722 | -1.0704 | 0.1533

The final computed trajectory, for an initial propellant mass of $24.89$ g and
transfer time $2757$ s is shown in Figure 4. The asteroid rendering for Figure
4 was computed as in [22].

Figure 2: The evolution of the value function, $\omega$, projected in the
$\rho$-$v_{t}$ plane for $v_{\rho}=0$ and $\Delta m=100$ g. Figure 3: The
evolution of the zero level set of the value function, $\omega$, projected in
the $\rho$-$v_{t}$ plane for $v_{\rho}=0$ and $\Delta m=100$ g.

The accuracy of the final orbit is within $41$ meters of the target orbit.
Calculating $\omega$ took 9 hours using a 3 GHz 8-Core Intel Core i7-9700
processor running Matlab with an extension of the Level Set Methods toolbox
[13]. The run time and accuracy can be significantly improved upon when using
optimized code such as [10, 12]. To show how the accuracy of the solution, the
memory usage as well as the CPU time varies with the grid size, we recompute
the value function and the trajectory with coarser spatial grids, as presented
in Table 3. We use only single-precision arrays to store the value function,
yet utilized double-precision arrays for all numerical calculations on the
value function. The subscript final denotes the values of the final point of
the computed trajectory, i.e. $\mathbf{r(0)}$, while target refers to the
target orbit used to define $\mathcal{C}$ and subsequently initialize the
computation of the value function.

Table 3: Comparison of the accuracy, CPU time, and memory usage as the size of the grid changes | Low | Medium | High
---|---|---|---
Grid Size | [24 18 18 12] | [32 24 24 16] | [50 40 40 32]
CPU Time | 20 min | 44 min | 9.28 hours
Memory to store $\omega$ | 374 MB | 1181 MB | 10.25 GB
$|\rho_{\mathrm{final}}-\rho_{\mathrm{target}}|$ | 82.72 m | 64.72 m | 40.83 m
$|v_{\rho,\mathrm{final}}-v_{\rho,\mathrm{target}}|$ |  6.80e-5 $\frac{m}{s}$ |  4.26e-5 $\frac{m}{s}$ |  8.31e-7 $\frac{m}{s}$
$|v_{t,\mathrm{final}}-v_{t,\mathrm{target}}|$ |  2.40e-4 $\frac{m}{s}$ |  4.01e-4 $\frac{m}{s}$ |  2.41e-4 $\frac{m}{s}$

Once $\omega$ is computed, it is incorporated into (15), which is solved using
Matlab’s paretosearch function. Solving the MOC problem took 120 seconds and
the resulting Pareto front is shown in Figure 5.

Figure 4: Initial orbit and transfer trajectory to a circular orbit $\approx
1$ km further away from the asteroid. Figure 5: Pareto front of the objective
functions $J_{1}=\Delta m$ in g and $J_{2}=t_{f}$ in seconds. The trajectory
in Figure 4 is derived using the point marked by the black diamond, while the
time-optimal solution used in Figure 6 uses the point denoted by the magenta
pentagram.

A comparison of the thrust magnitude of the smoothed optimal control policy is
shown in Figure 6. As can be seen, the time-optimal solution uses near
continuous thrust to reach the target, at the cost of using a large amount of
fuel. The control policy of the trajectory presented in Figure 4 meanwhile,
has noticeable cruising phases where no fuel is consumed. As expected, the
thrust magnitude of both policies follows a bang-bang structure.

Figure 6: A comparison of the thrust magnitude of the smoothed optimal control
policy of the time-optimal solution and the trajectory presented in Figure 4.

To illustrate the results of problems in Bolza form, we consider the case of
optimizing the remaining propellant in oppose to the initial propellant.
Therefore, let us fix the initial propellant to $100$ g. We use a similar
setup as in 5, with the addition of the auxiliary state, $z$, defined as a
terminal cost. We consider the uniform spaced grid over $r$ and $z$, defined
in Table 4. Storing the value function using single precision requires 15.2 GB
of memory.

Table 4: Spatial grid configuration for the Bolza problem | $\rho$ | $v_{\rho}$ | $v_{t}$ | $\Delta m$ | $z$
---|---|---|---|---|---
Points | 32 | 24 | 24 | 32 | 32
Spacing | 0.007 | 0.069 | 0.016 | 0.0067 | 0.0067
Minimum | 0.805 | -0.277 | -1.446 | -0.0533 | -0.1533
Maximum | 1.029 | 1.317 | -1.068 | 0.1533 | 0.0533

Since the objective is to maximize the remaining propellant, the optimization
problem needs to minimize $-\Delta m$. Therefore, for a given final state
$r_{f}\in\mathbb{R}^{7}$ and transfer time $t_{f}\in[0,+\infty)$, we define
the cost functions as $J_{1}(r_{f},t_{f})\coloneqq-\Delta m$ and
$J_{2}(r_{f},t_{f})\coloneqq t_{f}$, where $\Delta m$ denotes the $7$-th
element of the state vector $r_{f}$ (the mass in our case). The 2-dimensional
objective function
$J:\mathbb{R}^{7}\times[0,+\infty)\rightarrow\mathbb{R}^{2}$ can then be
written as

$J(r_{f},t_{f})\coloneqq\left[J_{1}(r_{f},t_{f}),J_{2}(r_{f},t_{f})\right]^{T}.$
(36)

The resulting Pareto front is shown in Figure 7. Calculating the reachable set
took approximately $32$ hours of CPU time. Using the reachable set,
calculating the Pareto front took approximately $258$ seconds. As expected,
the Pareto front looks similar to that of Figure 5, yet not identical due to
numerical inaccuracy and change in the initial mass, resulting in modified
dynamics.

Using an initial propellant mass of $100$ g, z of $-75$ g, and transfer time
of $2698$ s, the transfer orbit is computed in the same way as before. As
expected, the transfer trajectory is similar to the trajectory shown in Figure
4, and the accuracy of the transfer orbit is within $66$ meters of the target
orbit. For comparison to the first approach, the accuracy of the second
approach is shown in Table 5. From comparing the Pareto front, it can be seen
that both methods have a minimum transfer time of around 2500 seconds as well
as a minimum propellant requirement of just over 20 g.

Figure 7: Pareto front of the objective functions $J_{1}=z$ in g and $J_{2}=t_{f}$ in seconds, for the Bolza problem. The point used for trajectory calculations is marked by the black diamond. Table 5: Comparison of the accuracy, cpu time, and memory usage for the second formulation in Bolza form. | Bolza-formulation
---|---
Grid Size | [32 24 24 32 32]
CPU Time | 32 hours
Memory to store $\omega$ | 15.2 GB
$|\rho_{\mathrm{final}}-\rho_{\mathrm{target}}|$ | 66.14 m
$|v_{\rho,\mathrm{final}}-v_{\rho,\mathrm{target}}|$ |  1.8e-6 $\frac{m}{s}$
$|v_{t,\mathrm{final}}-v_{t,\mathrm{target}}|$ |  0.0018 $\frac{m}{s}$

##  6 Conclusion

We have presented a novel method of using the value function of a quasi-
variational inequality to compute the decision space of multi-objective
optimization problems. The feasibility and effectiveness of the proposed
approach was demonstrated by applying it to the problem of low-thrust
trajectory design. The approach is applicable to arbitrary multi-objective
optimization problems where the control variable is required to lie within a
reachable set.

Since the Hamiltonian becomes affine with respect to the thrust magnitude,
once the thrust angles have been fixed, future research aims to exploit this
fact by integrating classification based approaches [33]. Furthermore,
utilizing approximations of the reachable set as in [5], as well as
decomposing the reachable sets as in [14], seems promising.

## References

  * [1] Matthew J. Abrahamson, Alessandro Ardito, Dongsuk Han, Robert Haw, Brian Kennedy, Nick Mastrodemos, Sumita Nandi, Ryan Park, Brian Rush, and Andrew Vaughan.  Dawn orbit determination team: Trajectory modeling and reconstruction processes at Vesta.  Advances in the Astronautical Sciences, 148(February):2271–2290, 2013.
  * [2] Albert Altarovici, Olivier Bokanowski, and Hasnaa Zidani.  A general Hamilton-Jacobi framework for non-linear state-constrained control problems.  ESAIM - Control, Optimisation and Calculus of Variations, 19(2):337–357, 2013.
  * [3] P. G. Antreasian, F. T. Nicholson, P. H. Kallemeyn, S. Bhaskaran, R. J. Haw, and P. Halamek.  Galileo orbit determination for the Ida encounter.  Advances in the Astronautical Sciences, 87(2):1027–1048, 1994.
  * [4] Jean Pierre Aubin, John Lygeros, Marc Quincampoix, Shankar Sastry, and Nicolas Seube.  Impulse differential inclusions: A viability approach to hybrid systems.  IEEE Transactions on Automatic Control, 47(1):2–20, 2002.
  * [5] Somil Bansal and Claire J. Tomlin.  Deepreach: A deep learning approach to high-dimensional reachability.  CoRR, abs/2011.02082, 2020.
  * [6] Bardi, Martino and Capuzzo-Dolcetta, Italo.  Optimal Control and Viscosity Solutions of Hamilton-Jacobi-Bellman Equations.  Birkhäuser Basel, 1997.
  * [7] O. Bokanowski, E. Bourgeois, A. Désilles, and H. Zidani.  Global optimization approach for the ascent problem of multi-stage launchers.  In Hans Georg Bock, Willi Jäger, Ekaterina Kostina, and Hoang Xuan Phu, editors, Modeling, Simulation and Optimization of Complex Processes HPSC 201, pages 1–42, Cham, 2021. Springer International Publishing.
  * [8] Olivier Bokanowski, Eric Bourgeois, A Désilles, and Hasnaa Zidani.  Global optimization approach for the climbing problem of multi-stage launchers .  working paper or preprint, February 2015.
  * [9] Olivier Bokanowski, Eric Bourgeois, Anna Désilles, and Hasnaa Zidani.  Payload optimization for multi-stage launchers using HJB approach and application to a SSO mission.  IFAC-PapersOnLine, 50(1):2904–2910, 2017.
  * [10] Olivier Bokanowski, Anna Désilles, Hasnaa Zidani, and J Zhao.  ROC-HJ solver.
  * [11] Olivier Bokanowski, Nicolas Forcadel, and Hasnaa Zidani.  Reachability and minimal times for state constrained nonlinear problems without any controllability assumption.  SIAM Journal on Control and Optimization, 48(7):pp. 4292–4316, 2010.
  * [12] Mo Chen, Jaime Fisac, and Claire Tomlin.  Berkeley efficient api in c++ for level set methods.
  * [13] Mo Chen, Sylvia Herbert, Somil Bansal, and Claire Tomlin.  Optimal control helper toolbox.
  * [14] Mo Chen, Sylvia L. Herbert, Mahesh S. Vashishtha, Somil Bansal, and Claire J. Tomlin.  Decomposition of Reachable Sets and Tubes for a Class of Nonlinear Systems.  IEEE Transactions on Automatic Control, 63(11), 2018.
  * [15] Anna Désilles and Hasnaa Zidani.  Pareto front characterization for multiobjective optimal control problems using Hamilton-Jacobi approach.  SIAM Journal on Control and Optimization, 57(6):3884–3910, 2019\.
  * [16] I.J. Fialho and T.T. Georgiou.  Worst case analysis of nonlinear systems.  IEEE Transactions on Automatic Control, 44(6):1180–1196, 1999.
  * [17] Jaime F. Fisac, Mo Chen, Claire J. Tomlin, and S. Shankar Sastry.  Reach-avoid problems with time-varying dynamics, targets and constraints.  In Proceedings of the 18th International Conference on Hybrid Systems: Computation and Control, HSCC ’15, page 11–20, New York, NY, USA, 2015\. Association for Computing Machinery.
  * [18] Donald T Greenwood.  Principles of dynamics.  Prentice-Hall, Englewood Cliffs, N.J, 2nd ed. edition, 1988.
  * [19] C. Hermosilla, P. R. Wolenski, and H. Zidani.  The Mayer and Minimum Time Problems with Stratified State Constraints.  Set-Valued and Variational Analysis, 26(3):643–662, 2018.
  * [20] R. Scott Hudson and Steven J. Ostro.  Shape of asteroid 4769 Castalia (1989 PB) from inversion of radar images.  Science, 1994.
  * [21] Yu Jiang and Hexi Baoyin.  Orbital Mechanics near a Rotating Asteroid.  Journal of Astrophysics and Astronomy, 2014.
  * [22] Shankar Kulumani and Taeyoung Lee.  Low-thrust trajectory design using reachability sets near asteroid 4769 castalia.  In AIAA/AAS Astrodynamics Specialist Conference, September 2016\.
  * [23] Daniel Liberzon.  Calculus of variations and optimal control theory: A concise introduction.  Princeton University Press, 2011.
  * [24] John Lygeros, Claire Tomlin, and Shankar Sastry.  Controllers for reachability specifications for hybrid systems.  Automatica, 35(3):349–370, 1999.
  * [25] Kostas Margellos and John Lygeros.  Hamilton-jacobi formulation for reach-avoid differential games.  IEEE Transactions on Automatic Control, 56(8):1849–1861, aug 2011\.
  * [26] Kostas Margellos and John Lygeros.  Toward 4-D trajectory management in air traffic control: A study based on monte carlo simulation and reachability analysis.  IEEE Transactions on Control Systems Technology, 21(5):1820–1833, 2013.
  * [27] Kaisa Miettinen.  Nonlinear Multiobjective Optimization, volume 12 of International Series in Operations Research & Management Science.  Springer US, Boston, MA, 1998.
  * [28] Ian M. Mitchell.  The flexible, extensible and efficient toolbox of level set methods.  Journal of Scientific Computing, 35(2-3):300–329, 2008.
  * [29] Ian M. Mitchell, Alexandre M. Bayen, and Claire J. Tomlin.  A time-dependent Hamilton-Jacobi formulation of reachable sets for continuous dynamic games.  IEEE Transactions on Automatic Control, 50(7):947–957, 2005.
  * [30] S. Ober-Blöbaum, M. Ringkamp, and G. zum Felde.  Solving multiobjective optimal control problems in space mission design using discrete mechanics and reference point techniques.  In 2012 IEEE 51st IEEE Conference on Decision and Control (CDC), pages 5711–5716, 2012.
  * [31] S. Osher and R. Fedkiw.  Level Set Methods and Dynamic Implicit Surfaces.  Applied Mathematical Sciences. Springer New York, 2003.
  * [32] Daniel W. Parcher.  Low-thrust orbit transfer design for Dawn operations at Vesta.  Advances in the Astronautical Sciences, 140(818):1213–1232, 2011\.
  * [33] Vicenç Rubies-Royo, David Fridovich-Keil, Sylvia Herbert, and Claire J. Tomlin.  A classification-based approach for approximate reachability.  In 2019 International Conference on Robotics and Automation (ICRA), pages 7697–7704, 2019.
  * [34] Shankar Sastry.  Nonlinear systems : analysis, stability, and control.  Interdisciplinary applied mathematics ; v. 10. Springer, New York ; London, 1999.
  * [35] D. J. Scheeres, J. K. Miller, and D. K. Yeomans.  The Orbital Dynamics Environment of 433 Eros: A Case Study for Future Asteroid Missions.  BBA - Bioenergetics, 1987.
  * [36] D. J. Scheeres, S. J. Ostro, R. S. Hudson, and R. A. Werner.  Orbits close to asteroid 4769 Castalia.  Icarus, 121(1):67–87, 1996.
  * [37] W. Seefelder.  Lunar Transfer Orbits Utilizing Solar Perturbations and Ballistic Capture.  Luft- und Raumfahrt. Utz, Wiss., 2002.
  * [38] Lawrence F. Shampine and Mark W. Reichelt.  The MATLAB ODE Suite.  SIAM Journal on Scientific Computing, 18(1):1–22, jan 1997.
  * [39] Abolfazl Shirazi, Josu Ceberio, and Jose A. Lozano.  Spacecraft trajectory optimization: A review of models, objectives, approaches and solutions.  Progress in Aerospace Sciences, 102:76–98, 2018.
  * [40] Nikolaus Vertovec, Sina Ober-Blöbaum, and Kostas Margellos.  Multi-objective minimum time optimal control for low-thrust trajectory design.  In 2021 European Control Conference (ECC), pages 1975–1980, 2021\.
  * [41] O. von Stryk and R. Bulirsch.  Direct and indirect methods for trajectory optimization.  Annals of Operations Research, 37:357–373, 12 1992.

##  Appendix A Proposition A.14

###### Proposition A.14

Under Assumption 2.2, any two trajectories $\mathbf{r}$ and $\mathbf{\hat{r}}$
reconstructed from $f$, with $\mathbf{r}(-t_{f})=r_{0}$ and
$\mathbf{\hat{r}}(-t_{f})=\hat{r}_{0}$, respectively, are such that
$\big{|}\big{|}\mathbf{r}(\tau)-\mathbf{\hat{r}}(\tau)\big{|}\big{|}\leq\big{|}\big{|}r_{0}-\hat{r}_{0}\big{|}\big{|}e^{(t_{f}+\tau)L_{f}}$
for all $\tau\in[-t_{f},0]$.

###### Proof 2

Let $r_{0},\hat{r}_{0}\in\mathbb{R}^{7}$ be two initial states, and
$t_{f}\in[0,\infty)$. For the same $t_{f}$, we choose two trajectory control
pairs $(\mathbf{r},\mathbf{u})\in\Pi_{r_{0},t_{f}}$ and
$(\mathbf{\hat{r}},\mathbf{\hat{u}})\in\Pi_{\hat{r}_{0},t_{f}}$. Then by
Carathéodory’s existence of solutions [34], the following relation holds:

$\displaystyle\big{|}\big{|}\mathbf{r}(-t)-\mathbf{\hat{r}}(-t)\big{|}\big{|}$
$\displaystyle\leq\big{|}\big{|}r_{0}-\hat{r}_{0}\big{|}\big{|}+\int_{-t_{f}}^{-t}\big{|}\big{|}f(\mathbf{r}(s),\mathbf{u}(s))-f(\mathbf{\hat{r}}(s),\mathbf{\hat{u}}(s))\big{|}\big{|}ds$
$\displaystyle\leq\big{|}\big{|}r_{0}-\hat{r}_{0}\big{|}\big{|}+L_{f}\int_{-t_{f}}^{-t}\big{|}\big{|}\mathbf{r}(s)-\mathbf{\hat{r}}(s)\big{|}\big{|}ds$
$\displaystyle\leq\big{|}\big{|}r_{0}-\hat{r}_{0}\big{|}\big{|}e^{(t_{f}-t)L_{f}},$

where the second inequality is due to Assumption 2.2, while the last
inequality is due to the Bellman-Gronwall Lemma [34].

##  Appendix B Proof of Proposition 4.13

###### Proof 3

Fix
$(r_{0},z_{0}),(\hat{r}_{0},\hat{z}_{0})\in\mathbb{R}^{7}\times\mathbb{R}^{p}$,
$t_{f}\in[0,\infty)$ and let $\epsilon>0$. We choose
$(\hat{\mathbf{r}},\hat{\mathbf{u}},\hat{\mathbf{z}})\in\mathcal{Z}_{\hat{r}_{0},t_{f},\hat{z}_{0},}$
such that

$\vartheta(\hat{r}_{0},t_{f},\hat{z}_{0})\geq\bigvee_{i}\left[J^{i}_{\mathrm{t}}(\mathbf{\hat{r}}(0))-\mathbf{\hat{z}}^{i}(-t_{f})\right]\bigvee\nu(\mathbf{\hat{r}}(0))\bigvee\max_{s\in[-t_{f},0]}g(\mathbf{\hat{r}}(s))-\epsilon.$

By definition of $\vartheta$, for any
$(\mathbf{r},\mathbf{u})\in\Pi_{{r}_{0},{t}}$, this yields the following
relation

$\vartheta(r_{0},t_{f},z_{0})-\vartheta(\hat{r}_{0},t_{f},\hat{z}_{0})\leq\bigvee_{i}\left[J^{i}_{\mathrm{t}}(\mathbf{r}(0))-\mathbf{z}^{i}(-t_{f})\right]\bigvee\nu(\mathbf{r}(0))\bigvee\max_{s\in[-t_{f},0]}g(\mathbf{r}(s))\\\
-\bigvee_{i}\left[J^{i}_{\mathrm{t}}(\mathbf{\hat{r}}(0))-\mathbf{\hat{z}}^{i}(-t_{f})\right]\bigvee\nu(\mathbf{\hat{r}}(0))\bigvee\max_{s\in[-t_{f},0]}g(\mathbf{\hat{r}}(s))+\epsilon.$

Let $\kappa\in[-t_{f},0]$ be such that

$g(\mathbf{r}(\kappa))=\max_{s\in[-t_{f},0]}g(\mathbf{r}(s)).$

We then have

$\vartheta(r_{0},t_{f},z_{0})-\vartheta(\hat{r}_{0},t_{f},\hat{z}_{0})\leq\bigvee_{i}\left[J^{i}_{\mathrm{t}}(\mathbf{r}(0))-\mathbf{z}^{i}(-t_{f})\right]\bigvee\nu(\mathbf{r}(0))\bigvee
g(\mathbf{r}(\kappa))\\\
-\bigvee_{i}\left[J^{i}_{\mathrm{t}}(\mathbf{\hat{r}}(0))-\mathbf{\hat{z}}^{i}(-t_{f})\right]\bigvee\nu(\mathbf{\hat{r}}(0))\bigvee
g(\mathbf{\hat{r}}(\kappa))+\epsilon.$

Using Proposition A.14, we distinguish the following cases.

Case A:
$g(\mathbf{r}(\kappa))\geq\bigvee_{i}\left[J^{i}_{\mathrm{t}}(\mathbf{r}(0))-\mathbf{z}^{i}(-t_{f})\right]\bigvee\nu(\mathbf{r}(0))$

$\vartheta(r_{0},t_{f},z_{0})-\vartheta(\hat{r}_{0},t_{f},\hat{z}_{0})\leq\\\
g(\mathbf{r}(\kappa))-\bigvee_{i}\left[J^{i}_{\mathrm{t}}(\mathbf{\hat{r}}(0))-\mathbf{\hat{z}}^{i}(-t_{f})\right]\bigvee\nu(\mathbf{\hat{r}}(0))\bigvee
g(\mathbf{\hat{r}}(\kappa))+\epsilon\\\ \leq
g(\mathbf{r}(\kappa))-g(\mathbf{\hat{r}}(\kappa))+\epsilon\leq
L_{g}e^{(t_{f}+\kappa)L_{f}}\big{|}\big{|}r_{0}-\hat{r}_{0}\big{|}\big{|}+\epsilon,$

where the last inequality is due to the fact that $g$ is Lipschitz continuous.

Case B:
$\nu(\mathbf{r}(0))\geq\bigvee_{i}\left[J^{i}_{\mathrm{t}}(\mathbf{r}(0))-\mathbf{z}^{i}(-t_{f})\right]\bigvee
g(\mathbf{r}(\kappa))$

$\vartheta(r_{0},t_{f},z_{0})-\vartheta(\hat{r}_{0},t_{f},\hat{z}_{0})\leq\\\
\nu(\mathbf{r}(0))-\bigvee_{i}\left[J^{i}_{\mathrm{t}}(\mathbf{\hat{r}}(0))-\mathbf{\hat{z}}^{i}(-t_{f})\right]\bigvee\nu(\mathbf{\hat{r}}(0))\bigvee
g(\mathbf{\hat{r}}(\kappa))+\epsilon\\\
\leq\nu(\mathbf{r}(0))-\nu(\mathbf{\hat{r}}(0))+\epsilon\leq
L_{\nu}e^{t_{f}L_{f}}\big{|}\big{|}r_{0}-\hat{r}_{0}\big{|}\big{|}+\epsilon$

Case C:
$\bigvee_{i}\left[J^{i}_{\mathrm{t}}(\mathbf{{r}}(0))-\mathbf{{z}}^{i}(-t_{f})\right]\geq
g(\mathbf{r}(\kappa))\bigvee\nu(\mathbf{r}(0))$

Recall (27), then under Assumption 2.2, any two trajectories $\mathbf{z}$ and
$\mathbf{\hat{z}}$ reconstructed from $J_{\mathrm{r}}$ with $z_{0}$ and
$\hat{z}_{0}$, respectively, are bounded within a given time interval
$[-t,0]$. To see this,

$\displaystyle\big{|}\big{|}\mathbf{z}(-t)-\mathbf{\hat{z}}(-t)\big{|}\big{|}$
$\displaystyle\leq\big{|}\big{|}z_{0}-\hat{z}_{0}\big{|}\big{|}+\int_{-t}^{0}\big{|}\big{|}J_{\mathrm{r}}(\mathbf{r}(s),\mathbf{u}(s))-J_{\mathrm{r}}(\mathbf{\hat{r}}(s),\mathbf{\hat{u}}(s))\big{|}\big{|}ds$
$\displaystyle\leq\big{|}\big{|}z_{0}-\hat{z}_{0}\big{|}\big{|}+L_{\mathrm{r}}\int_{-t}^{0}\big{|}\big{|}\mathbf{r}(s)-\mathbf{\hat{r}}(s)\big{|}\big{|}ds$
$\displaystyle\leq\big{|}\big{|}z_{0}-\hat{z}_{0}\big{|}\big{|}+L_{\mathrm{r}}\int_{-t}^{0}\big{|}\big{|}r_{0}-\hat{r}_{0}\big{|}\big{|}e^{(t_{f}+s)L_{f}}ds$
$\displaystyle\leq\big{|}\big{|}z_{0}-\hat{z}_{0}\big{|}\big{|}+\big{|}\big{|}r_{0}-\hat{r}_{0}\big{|}\big{|}L_{\mathrm{r}}e^{t_{f}L_{f}}\frac{1-e^{-tL_{f}}}{L_{f}},$

where the third inequity is due to Proposition A.14, and the last one follows
by performing the integration.

Next, let $j\in[1,\ldots,p]$ be such that

$J^{j}_{\mathrm{t}}(\mathbf{{r}}(0))-\mathbf{{z}}^{j}(-t_{f})=\bigvee_{i}\left[J^{i}_{\mathrm{t}}(\mathbf{{r}}(0))-\mathbf{{z}}^{i}(-t_{f})\right].$

Then it follows, that

$\displaystyle\vartheta(r_{0},t_{f},z_{0})-\vartheta(\hat{r}_{0},t_{f},\hat{z}_{0})$
$\displaystyle\leq\bigvee_{i}\left[J^{i}_{\mathrm{t}}(\mathbf{{r}}(0))-\mathbf{{z}}^{i}(-t_{f})\right]$
$\displaystyle-\bigvee_{i}\left[J^{i}_{\mathrm{t}}(\mathbf{\hat{r}}(0))-\mathbf{\hat{z}}^{i}(-t_{f})\right]\bigvee\nu(\mathbf{\hat{r}}(0))\bigvee
g(\mathbf{\hat{r}}(\kappa))+\epsilon$
$\displaystyle\leq\bigvee_{i}\left[J^{i}_{\mathrm{t}}(\mathbf{{r}}(0))-\mathbf{{z}}^{i}(-t_{f})\right]$
$\displaystyle-\bigvee_{i}\left[J^{i}_{\mathrm{t}}(\mathbf{\hat{r}}(0))-\mathbf{\hat{z}}^{i}(-t_{f})\right]+\epsilon$
$\displaystyle\leq\left[J^{j}_{\mathrm{t}}(\mathbf{{r}}(0))-\mathbf{z}^{j}(-t_{f})\right]-\left[J^{j}_{\mathrm{t}}(\mathbf{\hat{r}}(0))-\mathbf{\hat{z}}^{j}(-t_{f})\right]+\epsilon$
$\displaystyle\leq\left[J^{j}_{\mathrm{t}}(\mathbf{{r}}(0))-J^{j}_{\mathrm{t}}(\mathbf{\hat{r}}(0))\right]-\left[\mathbf{z}^{j}(-t_{f})-\mathbf{\hat{z}}^{j}(-t_{f})\right]+\epsilon$

By Proposition A.14 and under Assumption 4.9

$\left[J^{j}_{\mathrm{t}}(\mathbf{{r}}(0))-J^{j}_{\mathrm{t}}(\mathbf{\hat{r}}(0))\right]\leq
L_{\mathrm{t}}\big{|}\big{|}r_{0}-\hat{r}_{0}\big{|}\big{|}e^{t_{f}L_{f}}.$

Finally, this yields the relationship

$\left[J^{j}_{\mathrm{t}}(\mathbf{{r}}(0))-J^{j}_{\mathrm{t}}(\mathbf{\hat{r}}(0))\right]-\left[\mathbf{z}^{j}(-t_{f})-\mathbf{\hat{z}}^{j}(-t_{f})\right]+\epsilon\\\
\leq\big{|}\big{|}z_{0}-\hat{z}_{0}\big{|}\big{|}+\big{|}\big{|}r_{0}-\hat{r}_{0}\big{|}\big{|}\Big{[}L_{\mathrm{t}}e^{t_{f}L_{f}}+L_{\mathrm{r}}\frac{e^{t_{f}L_{f}}-1}{L_{f}}\Big{]}+\epsilon.$

Thus in every case, on an interval $[0,t_{f}]$ there exists a set of constants
$C_{r}$ and $C_{z}$, such that

$\vartheta(r_{0},t_{f},z_{0})-\vartheta(\hat{r}_{0},t_{f},\hat{z}_{0})\\\ \leq
C_{r}\big{|}\big{|}r_{0}-\hat{r}_{0}\big{|}\big{|}+C_{z}\big{|}\big{|}z_{0}-\hat{z}_{0}\big{|}\big{|}+\epsilon$

The same argument conducted with $(r_{0},t_{f},z)$ and
$(\hat{r}_{0},t_{f},\hat{z})$ reversed establishes that

$\vartheta(\hat{r}_{0},t_{f},\hat{z}_{0})-\vartheta(r_{0},t_{f},z_{0})\\\ \leq
C_{r}\big{|}\big{|}r_{0}-\hat{r}_{0}\big{|}\big{|}+C_{z}\big{|}\big{|}z_{0}-\hat{z}_{0}\big{|}\big{|}+\epsilon.$

Since $\epsilon$ is arbitrary, we conclude that

$\big{|}\big{|}\vartheta(\hat{r}_{0},t_{f},\hat{z}_{0})-\vartheta(r_{0},t_{f},z_{0})\big{|}\big{|}\\\
\leq
C_{r}\big{|}\big{|}r_{0}-\hat{r}_{0}\big{|}\big{|}+C_{z}\big{|}\big{|}z_{0}-\hat{z}_{0}\big{|}\big{|},$

thus concluding the proof.

The proof for $\omega$ is similar to that of $\vartheta$ and will therefore be
omitted in the interest of space.

##  Appendix C Derivation of the optimal thrust angles

Notice that since the applied thrust, $T$, is always positive, the term

$\left(q_{4}\cos{\alpha}+\sin{\alpha}\left(q_{5}\sin{\delta}+q_{6}\cos{\delta}\right)\right)$

in the Hamiltonian can be minimized separately from $T$.

To this end, we introduce the auxiliary variables

$\displaystyle\chi(\delta)\coloneqq$
$\displaystyle{}\sqrt{q_{5}^{2}+q_{6}^{2}}\cos{(\delta-\arctan{\frac{q_{5}}{q_{6}}})},$
$\displaystyle A(\delta)\coloneqq$
$\displaystyle\sqrt{q_{4}^{2}+\chi(\delta)^{2}},$

using the trigonometric identity

$a\cos{x}+b\sin{x}=R\cos{(x-\arctan{\frac{b}{a}})},$

with $R=\sqrt{a^{2}+b^{2}}$. First optimizing over $\alpha$, and subsequently
over $\delta$ (notice that this sequential minimization is exact since
$A(\delta)\geq 0$) results in

$\min_{\alpha,\delta\in[-\pi,\pi]\times[-\frac{\pi}{2},\frac{\pi}{2}]}\left(q_{4}\cos{\alpha}+\sin{\alpha}\left(q_{5}\sin{\delta}+q_{6}\cos{\delta}\right)\right)\\\
=\min_{\delta\in[-\frac{\pi}{2},\frac{\pi}{2}]}A(\delta)\min_{\alpha\in[-\pi,\pi]}\cos{(\alpha-\arctan{\frac{\chi(\delta)}{q_{4}}})}.$

Thus, the optimal thrust angles are given by

$\alpha^{*}(\delta)\coloneqq\pm\pi+\arctan{\frac{\chi(\delta)}{q_{4}}}.$

Since $\cos{(\alpha^{*}(\delta)-\arctan{\frac{\chi(\delta)}{q_{4}}})}=-1$,
after applying $\alpha^{*}(\delta)$, it follows that

$\delta^{*}\in\operatorname*{arg\,min}_{\delta\in[-\frac{\pi}{2},\frac{\pi}{2}]}-A(\delta)=\arctan{\frac{q_{5}}{q_{6}}}\pm\pi.$

Subsequently,

$q_{4}\cos{\alpha^{*}}+\sin{\alpha^{*}}\left(q_{5}\sin{\delta^{*}}+q_{6}\cos{\delta^{*}}\right)\\\
=-\sqrt{q_{4}^{2}+q_{5}^{2}+q_{6}^{2}}.$ (37)

Substituting (37) into (17) results in (18) which depends in an affine fashion
on $T$.