11institutetext: Dario Izzo 22institutetext: European Space Agency, Noordwijk,
2201 AZ, The Netherlands 22email<EMAIL_ADDRESS>33institutetext:
Christopher Iliffe Sprague 44institutetext: KTH Royal Institute of Technology,
Teknikringen 14, SE-100 44 Stockholm, Sweden 44email<EMAIL_ADDRESS>55institutetext: Dharmesh Vijay Tailor 66institutetext: European Space Agency,
Noordwijk, 2201 AZ, The Netherlands 66email<EMAIL_ADDRESS>
# Machine learning and evolutionary techniques in interplanetary trajectory
design

Dario Izzo  Christopher Sprague  Dharmesh Tailor

###### Abstract

After providing a brief historical overview on the synergies between
artificial intelligence research, in the areas of evolutionary computations
and machine learning, and the optimal design of interplanetary trajectories,
we propose and study the use of deep artificial neural networks to represent,
on-board, the optimal guidance profile of an interplanetary mission. The
results, limited to the chosen test case of an Earth-Mars orbital transfer,
extend the findings made previously for landing scenarios and quadcopter
dynamics, opening a new research area in interplanetary trajectory planning.

##  1 Introduction

The use of Artificial Intelligence (AI) techniques for the design of
interplanetary trajectories, in particular in their early design phases, has
been proposed and studied by numerous scientists over the past few decades.
Among the AI techniques deployed to help the optimisation of spacecraft
trajectories are evolutionary algorithms Janin and Gomez-Tierno (1985); Gage
et al (1995); Izzo et al (2007); Lee et al (2005); Sentinella and Casalino
(2009); Vasile et al (2010); Olds et al (2007); Izzo et al (2013), machine
learning techniques Ampatzis and Izzo (2009); Cassioli et al (2012); Mereta et
al (2017); Izzo et al (2016), evolutionary neuro-controllers Dachwald (2004);
Dachwald and Ohndorf (2019), tree search methods Grigoriev and Zapletin
(2013); Hennes and Izzo (2015); Simões et al (2017); Izzo et al (2016, 2013)
to only name a few widely studied methods. While the synergies between
research in AI and in trajectory design were becoming increasingly apparent,
the field of AI as a whole experienced a renaissance in the second decade of
the millennium, delivering exciting developments that are of significance to
many scientific fields. Certainly the area of machine learning is one that is
powering such an AI renaissance, in particular with Deep Learning techniques
Schmidhuber (2015) becoming fundamental to define new benchmarks in the most
diverse applications. It thus seems likely that also the optimisation of
interplanetary trajectories, already receptive of AI methods in general, is
going to benefit from the vast amount of new knowledge being produced in the
context of AI research.

In this chapter we first briefly provide an overview on the state-of-the art
of the use of evolutionary techniques (Section §2) and machine learning
techniques (Section §3) when it comes to optimising interplanetary
trajectories, and we then introduce a novel idea concerning the application of
Deep Learning for the on-board representation of the optimal guidance profile
for an interplanetary probe, extending previous work made on planetary landing
scenarios Sánchez-Sánchez et al (2016); Sánchez-Sánchez and Izzo (2016);
Schiavone et al (2012).

##  2 Evolutionary algorithms for trajectory planning, an overview

Evolutionary algorithms are a class of global optimisation techniques that
make use of heuristic rules, often inspired but not limited to natural
paradigms such as Darwinian evolution, to search for optimal solutions in
discontinuous, rugged landscapes. As such they are a good match to solve
interplanetary trajectory optimisation problems where the planetary geometry
defines a quite complex solution landscape already in simple cases. Already in
1985 a Genetic Algorithm (GA) Janin and Gomez-Tierno (1985) was proposed and
studied in the context of interplanetary trajectory optimisation concluding
that

> … a considerable effort is still needed for developing efficient schemes
> using genetic algorithms. However, they appear to offer an entirely original
> way for solving a large class of [interplanetary trajectory] global
> optimisation problems …

Given the CPU power available at the time, it is only natural that this
pioneering work complained about efficiency. In the following decade the
“entirely original way” was consolidated and more successful applications of
GAs for both low-thrust propelled spacecraft Rauwolf and Coverstone-Carroll
(1996) and impulsive thrust strategies Gage et al (1995) were deployed.
Studies on the use of some form of GA, also leveraging on the ever increasing
computational power available, continued (see for example Rogata et al (2003);
Biesbroek and Ancarola (2002); Gad and Abdelkhalik (2011); Deb et al (2007);
Elsayed et al (2011)) showing how the original intuition, back in 1985, was
one rich of consequences. It was indeed relatively consequential to substitute
the GA with any other form of stochastic optimiser, so that in 2004 the
European Space Agency performed a series of studies, in collaboration with
academia, to test and compare several global optimisers on interplanetary
trajectory problems Myatt et al (2004); Di Lizia and Radice (2004); Izzo et al
(2007). It then became clear how other evolutionary algorithms were offering
alternative, and in many ways superior, choices to help the preliminary phases
of trajectory design. Among them, Differential Evolution (DE) and Particle
Swarm Optimisation (PSO) were identified and benchmarked on several chemical
propulsion problems opening the way to further independent confirmations of
their performances Vasile et al (2010); Yao et al (2017); Olds et al (2007);
Sentinella and Casalino (2009); Luo et al (2005); Pontani and Conway (2010).
Interestingly, in 2013, a self-adaptive Differential Evolution algorithm was
used to plan the grand tour of the Jupiter moons that was awarded the golden
Humies medal Izzo et al (2013) for “human-competitive results produced by
genetic and evolutionary computation”. Many other evolutionary approaches have
been proved to be of use by various researchers during the early 2000s,
including Simulated Annealing (SA) Luo et al (2005) and Ant Colony
Optimisation (ACO) Radice and Olmo (2006); Ceriotti and Vasile (2010). More
recently, Covariance Matrix Evolutionary Strategy (CMA-ES) was shown to be
potentially outperforming self-adaptive DE on a class of transfers Izzo et al
(2014a).

While the no free lunch theorem Wolpert and Macready (1997) guarantees that no
evolutionary approach is better than any other on average, restricting the
class of optimisation problems to those representing interplanetary transfers
allows to identify the most useful approaches. To this end, a set of problems
called GTOP database Vinkó and Izzo (2008) (a sort of open trajectory gym) was
created and made available to the scientific community and is still the
subject of active research Stracquadanio et al (2011); Addis et al (2011);
Schlueter (2014); Islam et al (2012); Cassioli et al (2012); Simões et al
(2014). Some of the the interplanetary trajectory problems (i.e. Messenger and
Cassini2) in the GTOP database were also used during the CEC2011 competition
attracting the attention of the larger scientific community of evolutionary
computations (see Elsayed et al (2011) for the competition winner).

Most of the research mentioned so far, considers continuous optimisation
problems with a single objective and box constraints. In its most general
case, though, the problem encountered in interplanetary trajectory design are
multi-objective, with nonlinear constraints and, possibly, integer decision
variables. Genetic approaches to multi-objective trajectory design were
benchmarked already in 2005 Lee et al (2005) followed by deeper studies on Non
Dominated Sorting Genetic Algorithm (NSGA-II) that included also, as integer
variables, the planetary fly-by sequence Deb et al (2007). Multi-objective
versions of PSO have also been considered early on Lavagna (2007). The more
modern Multi Objective Evolutionary Algorithm by Decomposition (MOEA/D) was
later identified as a most performing technique in applicable cases Izzo et al
(2014a), in the same work, several constraint handling techniques including
co-evolution and immune systems were tested.

Whenever the the representation of an interplanetary trajectory requires
integer and continuous variables, the resulting optimisation problem (a MINLP)
typically becomes intractable also for evolutionary approaches and while some
results have been obtained, for example using the ACO paradigm Schlueter et al
(2013); Ceriotti and Vasile (2010), a more convincing approach is to consider
the continuous part and the integer part of the problem separately and
architect some optimisation scheme tackling the two problems with different
nested techniques (bi-level optimisation). Evolutionary approaches based on
GAs Izzo et al (2015, 2014b); Englander (2013) or ACO Simões et al (2017) have
been used in bi-level optimisation schemes, often coupled with smart tree
search strategies such as Beam Search (BS) Grigoriev and Zapletin (2013);
Simões et al (2017), Monte Carlo Tree Search (MCTS) Hennes and Izzo (2015),
Multi Objective Beam Search (MOBS) Izzo et al (2016) or Lazy Race Tree Search
(LRTS) Izzo et al (2013) to take care of the integer part.

Evolutionary approaches to interplanetary trajectory planning have proven
their worth beyond any criticism and while they are still not as widely used
by the aerospace industry as they could, it is likely that we will see a
larger penetration in the industrial sector in the upcoming years, as the
interest in Artificial Intelligence methods powered by Deep Learning and
Machine Learning will put also evolutionary techniques in the spotlight.

##  3 Machine Learning and interplanetary trajectories

The use of machine learning (ML) algorithms to aid the design of
interplanetary trajectory is not as widely researched as that of evolutionary
techniques and is limited to fewer works. The reasons are to be found in the
less obvious applicability of these methods to the problems encountered in the
design of interplanetary trajectories and in the lack of data sets produced
and made available by the aerospace community. In this section we try to
summarise the ideas that have so far been proposed.

During the optimisation of an interplanetary transfer, as in any optimisation
task, a large number of solutions are computed and assessed to inform the
search for better candidates. Typically all these design points are discarded
and lost after they have been used to produce new promising search directions.
The idea of applying supervised learning on such a data set in order to build
a model improving the further selection of initial guesses to guide successive
evolutionary runs was proposed and tested Cassioli et al (2012) on some of the
trajectory problems in the GTOP database Vinkó and Izzo (2008) using Support
Vector Machines (SVM). Following a similar reasoning, a ML model can be
trained, using the points sampled during the evolution, in order to construct
a surrogate model of the trajectory worth, which then avoids expensive
evaluations of the objective function Ampatzis and Izzo (2009). Building
surrogate models is particularly relevant when the interplanetary mission
fitness requires a high degree of computational resources such as in the case
of optimal low-thrust transfers. In that case a surrogate model approximating
the final optimal transfer mass enables to quickly search for good launch and
arrival epochs, as well as planetary body sequences, e.g. in the case of
multiple asteroid or debris rendezvous missions Hennes et al (2016); Mereta et
al (2017). Unsupervised learning techniques such as clustering or nearest
neighbours have also been used to select the target of transfers in multiple
asteroid rendezvous missions, upon proper definition of a metric coping with
the orbital non linearities Izzo et al (2016), or to define new box bounds and
hence focus successive evolutionary runs in promising areas of the search
space (cluster pruning Izzo (2010)).

A second field where ML algorithms have been used in the context of research
in interplanetary trajectory design methods, is that of the on-board
representation of the optimal guidance profile. Already in 2004 Dachwald
(2004) the idea was put forward of using machine learning to learn a
representation of the optimal spacecraft guidance profile. In later years the
technique was studied exclusively in the context of evolutionary neuro-
controllers and optimal control. More recently, neuro evolution has been
substituted with supervised learning (back-propagation), the artificial neural
network structure has become deeper, and focus has been shifted to the on-
board real time computation of guidance profiles Sánchez-Sánchez and Izzo
(2016); Sánchez-Sánchez et al (2016); Schiavone et al (2012). These last works
have all been limited, to landing problems and the applicability of these
ideas to interplanetary trajectories has not been studied. In the next section
we will try to contribute closing this gap showing how the optimal guidance
profile of a phase-less Earth-Mars transfer can also be satisfactorily
represented by a deep artificial neural network.

##  4 On-board optimal guidance via a deep network

In planetary landing problems, the possibility to generate optimal guidance
profiles on-board has been studied and suggested to significantly enhance a
vehicle’s landing accuracy de Croon and Izzo (2012); Dueri et al (2016). The
resulting algorithms need to run on a radiation-hardened flight processor,
that is on a significantly slower processor than those available on modern
desktops and one having significant architectural differences. As a rule of
thumb, a space qualified processor is an order of magnitude slower than the
average desktop processor. The resulting CPU load on the spacecraft has been
estimated in a practical scenario Dueri et al (2016) to be 0.7s for the
computation of one optimal action. To the same end, an alternative is offered
by deep feed-forward neural networks trained on the ground to approximate the
optimal control and used on-board in real time. The on-board computational
effort associated to this architecture is that of one forward pass of the
network which, though deep, is not necessarily large Sánchez-Sánchez et al
(2016); Sánchez-Sánchez and Izzo (2016). Extending on previous studies on
landing and quadcopter dynamics, the feasibility of such a scheme for
interplanetary low-thrust transfers is studied here. The same overall scheme
is used Sánchez-Sánchez et al (2016):

  * •

Step 1: We solve thousands of optimal control problems using Pontryagin’s
maximum principle. We store, with some time sampling, the obtained solutions
in a data set.

  * •

Step 2: We train deep feed-forward neural networks on the data set to learn
the optimal control structure.

  * •

Step 3: We use, on-board, the trained network to compute the optimal feedback.

With respect to the previously studied case of planetary landings the
methodology is, essentially, unchanged and it is of great interest to see how
it can be applied on a radically different problem having larger
dimensionality, different non linearities and, arguably, a more complex
structure.

#### A short note on notation

We use boldface symbols (e.g. $\bm{\lambda},\mathbf{r}$, etc.) to denote
vector quantities, while scalars or the vector norm will be indicated by
normal symbols (e.g. $\lambda,r$, etc.). Unit vectors will be indicated with a
small hat (e.g. $\mathbf{\hat{i}}$) and time derivatives with a dot (e.g.
$\dot{x}$). We will sometime use the asterisk as a superscript to indicate
optimal quantities (e.g. $u^{*}$).

###  4.1 Spacecraft Dynamics

The motion of a spacecraft equipped with a constant specific impulse nuclear
electric low-thrust propulsion system can be described in some heliocentric,
inertial, frame using Cartesian coordinates by the equations:

$\begin{array}[]{l}\dot{\mathbf{r}}=\mathbf{v}\\\
\dot{\mathbf{v}}=-\frac{\mu}{r^{3}}\mathbf{r}+c_{1}\frac{u(t)}{m}\mathbf{\hat{i}_{u}}(t)\\\
\dot{m}=-c_{2}u(t)\end{array}$ (1)

where $\mathbf{r}$, $\mathbf{v}$ and $m$ denote the spacecraft position,
velocity and mass (also denoted as spacecraft state $\mathbf{x}$), $\mu$ is
the gravitational parameter of the Sun, $c_{1}$ the maximum thrust achievable
by the onboard propulsion system and $c_{2}=c_{1}/I_{sp}g_{0}$. We have
denoted with $I_{sp}$ the propulsion specific impulse and with $g_{0}$ the
Earth gravity constant at sea level. The control variables $u(t)$ and
$\mathbf{\hat{i}_{u}}(t)$ (or $\mathbf{u}(t)$ for brevity) describe the thrust
level (here also called throttle) and its direction and are constrained as
follows $|u(t)|\leq 1$ and $|\mathbf{\hat{i}_{u}}(t)|=1,\forall
t\in[t_{1},t_{2}]$.

We are interested in the problem of finding the controls $u(t)$ and
$\mathbf{\hat{i}_{u}}(t)$ to minimise the functional:

$J(t_{1},t_{2},u(t))=\alpha\int_{t_{1}}^{t_{2}}u(t)dt+(1-\alpha)\int_{t_{1}}^{t_{2}}u^{2}(t)dt$
(2)

and able to steer the spacecraft from some initial point
$x_{1}\in\mathcal{S}_{1}$ to some final point $x_{2}\in\mathcal{S}_{2}$ where
the sets $\mathcal{S}_{1}$ and $\mathcal{S}_{2}$ are some closed subsets
(hypersurfaces) of the state space. Note that the functional $J$ is
parameterised by the continuation parameter $\alpha\ \in[0,1]$ which weights
two contributions corresponding to mass optimality ($\alpha=1$) and quadratic
control optimality ($\alpha=0$).

We now derive the, known, necessary optimality conditions for the above stated
problem applying the maximum principle from Pontryagin (note that we have
stated a minimisation problem, hence the conditions are actually slightly
different from the ones originally derived in Pontryagin work Pontryagin et al
(1962)).

###  4.2 Minimisation of the Hamiltonian $\mathcal{H}$

We start by introducing seven auxiliary functions defined in
$t\in[t_{1},t_{2}]$, the co-states, which we will denote with
$\bm{\lambda}_{r}$,$\bm{\lambda}_{v}$ and $\lambda_{m}$, or, for brevity,
$\bm{\lambda}$. We then introduce the Hamiltonian function:

$\mathcal{H}(\mathbf{x},\bm{\lambda},\mathbf{u})=\bm{\lambda}_{r}\cdot\mathbf{v}+\bm{\lambda}_{v}\cdot\left(-\frac{\mu}{r^{3}}\mathbf{r}+c_{1}\frac{u}{m}\mathbf{\hat{i}_{u}}\right)-\lambda_{m}c_{2}u+\alpha
u+(1-\alpha)u^{2}$ (3)

which, following Pontryagin theory, needs to be minimised by our controls
along an optimal trajectory. Isolating the relevant part of the Hamiltonian
that depends on the control $\mathbf{\hat{i}_{u}}$ we have that
$\mathcal{H}=...+c_{1}\frac{u}{m}\bm{\lambda}_{v}\cdot\mathbf{\hat{i}_{u}}+...$,
which allows to conclude that, since $m,u$ and $c_{1}$ are all positive
numbers, the thrust direction must be in the opposite direction of
$\bm{\lambda}_{v}$ for $\mathcal{H}$ to be minimised, more formally:

$\mathbf{\hat{i}_{u}}^{*}=-\frac{\bm{\lambda}_{v}}{\lambda_{v}}$ (4)

Substituting this expression back into the Hamiltonian we get:

$\mathcal{H}(\mathbf{x},\bm{\lambda},\mathbf{u})=\bm{\lambda}_{r}\cdot\mathbf{v}-\frac{\mu}{r^{3}}\bm{\lambda}_{v}\cdot\mathbf{r}-c_{1}\frac{u}{m}\lambda_{v}-\lambda_{m}c_{2}u+\alpha
u+(1-\alpha)u^{2}$

which we now need to minimise with respect to $u$. Isolating the relevant part
of $\mathcal{H}$, i.e. the terms that depend on $u$, we have:

$\mathcal{H}=...+(1-\alpha)u^{2}+u\left(\alpha-\frac{c_{1}}{m}\lambda_{v}-\lambda_{m}c_{2}\right)+...$

which is a convex parabola that will take its minimal value in
$u=\frac{\frac{c_{1}}{m}\lambda_{v}+\lambda_{m}c_{2}-\alpha}{2(1-\alpha)}$.
Since $u\in[0,1]$, we get the final expression for an optimal $u^{*}$:

$u^{*}=\min\left[\max\left(\frac{\frac{c_{1}}{m}\lambda_{v}+\lambda_{m}c_{2}-\alpha}{2(1-\alpha)},0\right),1\right]$
(5)

Note that in the corner case $\alpha=1$, which correspond to a mass optimal
control, the above expression results to be singular and the minimiser of the
Hamiltonian can be conveniently rewritten introducing the switching function,

$S_{sw}=c_{1}\lambda_{v}+mc_{2}\lambda_{m}-m\alpha$

as:

$u^{*}=\left\\{\begin{array}[]{ll}1&\mbox{if}\quad S_{sw}>0\\\
0&\mbox{if}\quad S_{sw}<0\\\ \end{array}\right.$ (6)

###  4.3 The costate equations

The states $\mathbf{x}(t)$ and co-states $\bm{\lambda}(t)$ must be a solution
to the set of differential equations that are elegantly written using the
Hamiltonian formalism as:

$\begin{array}[]{ll}\dot{\mathbf{r}}=\frac{\partial
H}{\partial\bm{\lambda}_{r}},&\dot{\bm{\lambda}}_{r}=-\frac{\partial
H}{\partial\mathbf{r}}\\\ \dot{\mathbf{v}}=\frac{\partial
H}{\partial\bm{\lambda}_{v}},&\dot{\bm{\lambda}}_{v}=-\frac{\partial
H}{\partial\mathbf{v}}\\\ \dot{m}=\frac{\partial
H}{\partial\lambda_{m}},&\dot{\lambda}_{m}=-\frac{\partial H}{\partial m}\\\
\end{array}$

while for the first three equations it is trivial to get back the expressions
in Eq.(1), the co-states equations do require some extra work, in particular
when deriving the gravity gradient. Let us then start deriving the co-state
equations by computing $\frac{\partial H}{\partial\mathbf{r}}$ starting from
the expression in Eq.(4.2). We have

$\frac{\partial
H}{\partial\mathbf{r}}=-\bm{\lambda}_{v}\cdot\nabla\left(\frac{\mu}{r^{3}}\mathbf{r}\right)$

where the symbol $\nabla$ denotes the gradient operator. Its easier to compute
the expression regrouping as follows:

$\displaystyle-\frac{1}{\mu}\frac{\partial H}{\partial\mathbf{r}}$
$\displaystyle=$
$\displaystyle\nabla\left(\bm{\lambda}_{v}\cdot\frac{1}{r^{3}}\mathbf{r}\right)=$
$\displaystyle=$
$\displaystyle\frac{1}{r^{3}}\nabla\left(\bm{\lambda}_{v}\cdot\mathbf{r}\right)+\left(\bm{\lambda}_{v}\cdot\mathbf{r}\right)\nabla\left(\frac{1}{r^{3}}\right)=$
$\displaystyle=$
$\displaystyle\frac{\bm{\lambda}_{v}}{r^{3}}-3\left(\bm{\lambda}_{v}\cdot\mathbf{r}\right)\frac{\nabla
r}{r^{4}}=$ $\displaystyle=$
$\displaystyle\frac{\bm{\lambda}_{v}}{r^{3}}-3\left(\bm{\lambda}_{v}\cdot\mathbf{r}\right)\frac{\mathbf{r}}{r^{5}}$

We may thus write the corresponding Hamilton equation as:

$\dot{\bm{\lambda}}_{r}=-\frac{\partial
H}{\partial\mathbf{r}}=\mu\frac{\bm{\lambda}_{v}}{r^{3}}-3\mu\left(\bm{\lambda}_{v}\cdot\mathbf{r}\right)\frac{\mathbf{r}}{r^{5}}$

The remaining two equations are then also easily obtained as:

$\dot{\bm{\lambda}}_{v}=-\frac{\partial
H}{\partial\mathbf{v}}=-\bm{\lambda}_{r}$

and:

$\dot{\lambda}_{m}=-\frac{\partial H}{\partial
m}=c_{1}u\frac{\bm{\lambda}_{v}\cdot\mathbf{\hat{i}_{u}}}{m^{2}}$

The final system of differential equations describing the state and co-states
evolution along an optimal trajectory may be now summarised:

$\left\\{\begin{array}[]{l}\dot{\mathbf{r}}=\mathbf{v}\\\
\dot{\mathbf{v}}=-\frac{\mu}{r^{3}}\mathbf{r}-c_{1}\frac{u^{*}(t)}{m}\frac{\bm{\lambda}_{v}}{\lambda_{v}}\\\
\dot{m}=-c_{2}u^{*}(t)\\\
\dot{\bm{\lambda}}_{r}=\mu\frac{\bm{\lambda}_{v}}{r^{3}}-3\mu\left(\bm{\lambda}_{v}\cdot\mathbf{r}\right)\frac{\mathbf{r}}{r^{5}}\\\
\dot{\bm{\lambda}}_{v}=-\bm{\lambda}_{r}\\\
\dot{\lambda}_{m}=-c_{1}u^{*}(t)\frac{\lambda_{v}}{m^{2}}\end{array}\right.$
(7)

###  4.4 The two-point boundary value problem

Following the maximum principle, we know that there exist initial values for
the co-states $\bm{\lambda}(t_{1})=\bm{\lambda}_{1}$ such that an optimal
trajectory steering the system from a starting state
$\mathbf{x}_{1}\in\mathcal{S}_{1}$ to a final state
$\mathbf{x}_{2}\in\mathcal{S}_{2}$ will be found integrating the system of
equations stated in Eq.(7) from the initial condition
$\mathbf{x}_{1},\bm{\lambda}_{1}$, for a time $t_{2}-t_{1}$ leading to a final
condition $\mathbf{x}_{2},\bm{\lambda}_{2}$. Since we also want to find
optimal values for $\mathbf{x}_{1}\in\mathcal{S}_{1}$,
$\mathbf{x}_{2}\in\mathcal{S}_{2}$ and $t_{2}$ we need some additional
conditions. The problem of finding the optimal $t_{2}$ ($t_{1}$ can be assumed
to be 0 as our system is autonomous) is solved by adding a condition on the
Hamiltonian:

$\mathcal{H}(\mathbf{x}(t_{2}),\bm{\lambda}(t_{2}),\mathbf{u}^{*}(t_{2}))=0$
(8)

To optimally select $\mathbf{x}_{1}\in\mathcal{S}_{1}$ and
$\mathbf{x}_{2}\in\mathcal{S}_{2}$, instead, we will need to add some
conditions, called transversality conditions and derived in the next section.

###  4.5 Transversality conditions

Transversality conditions ensure that the initial and final states are
selected optimally in the allowed sets $\mathcal{S}_{1}$ and
$\mathcal{S}_{2}$. They can be elegantly written introducing the vectors
$\bm{\theta}_{1}$ and $\bm{\theta}_{2}$ belonging to the hyperplane tangent to
the hypersurfaces $\mathcal{S}_{1}$ and $\mathcal{S}_{2}$ in $\mathbf{x}_{1}$
and $\mathbf{x}_{2}$. For all such vectors it must be that:

$\bm{\lambda}(t_{1})\cdot\bm{\theta}_{1}=0$

and equivalently,

$\bm{\lambda}(t_{2})\cdot\bm{\theta}_{2}=0$

Let us now derive explicitly these relations in the case of the interplanetary
transfer dynamics. We start focusing on the terminal condition
$\mathbf{x}_{2}\in\mathcal{S}_{2}$. Typically the final mass of the spacecraft
is not specified and we may, in this case, express $\mathcal{S}_{2}$ as
$\mathbf{x}_{2}=[\mathbf{r}_{2},\mathbf{v}_{2},c]$ with $c$ being a free
parameter and $\mathbf{r}_{2},\mathbf{v}_{2}$, being at this stage considered
as fixed. Trivially, in this case,
$\bm{\theta}_{2}=[\mathbf{0},\mathbf{0},\theta_{2}]$ and the transversality
condition for a free final mass is thus:

$\bm{\lambda}_{m}(t_{2})=0$ (9)

More work is needed to derive transversality conditions in the case where also
$\mathbf{r}_{2},\mathbf{v}_{2}$ are not fixed but constrained to belong to
some manifold. A typical situation would be, for example, to leave some of the
final Keplerian orbital elements as free. In this case, in order to find
$\bm{\theta}_{2}$, i.e. the vector tangent to the defined manifold, we need to
write a parameterisation of such manifold, that is an expression of
$\mathbf{r}_{2}$ and $\mathbf{v}_{2}$ as a function of the desired orbital
parameter. Using the Keplerian elements $a,e,i,\omega,\Omega,E$ we may use
well know relations to establish such a parameterisation:

$\begin{array}[]{l}\mathbf{r}_{2}=\mathbf{R}\left[\begin{array}[]{c}a(\cos
E-e)\\\ a\sqrt{1-e^{2}}\sin E\\\ 0\end{array}\right]\\\
\mathbf{v}_{2}=\sqrt{\frac{\mu}{a^{3}}}\frac{1}{1-e\cos
E}\mathbf{R}\left[\begin{array}[]{c}-a\sin E\\\ a\sqrt{1-e^{2}}\cos E\\\
0\end{array}\right]\end{array}$ (10)

where $\mathbf{R}$ is the rotation matrix from the orbital reference frame to
the inertial defined as:

$\mathbf{R}=\left[\begin{array}[]{ccc}\cos\Omega\cos\omega-\sin\Omega\sin\omega\cos
i&-\cos\Omega\sin\omega-\sin\Omega\cos\omega\cos i&\sin\Omega\sin i\\\
\sin\Omega\cos\omega+\cos\omega\sin\omega\cos
i&-\sin\Omega\sin\omega+\cos\Omega\cos\omega\cos i&-\cos\Omega\sin i\\\
\sin\omega\sin i&\cos\omega\cos i&\cos i\end{array}\right]$

the tangent vector to the resulting manifold will then be simply obtained
deriving the above expression with respect to the chosen parameter.

For example, let us take the case of a free final anomaly. This corresponds to
a transfer to some final target orbit regardless of the phase. Deriving the
above expressions with respect to $E$, or in this case equivalently $t$, we
get the transversality condition:

$\bm{\lambda}_{r}\cdot\mathbf{v}_{2}-\bm{\lambda}_{v}\cdot\frac{\mu}{r_{2}^{3}}\mathbf{r}_{2}=0$
(11)

Similarly if we leave the semi-major axis $a$ as free, the conditions is
easily obtained as:

$2\bm{\lambda}_{r}\cdot\mathbf{r}_{2}-\bm{\lambda}_{v}\cdot\mathbf{v}_{2}=0$
(12)

Similar expressions are derived in Pan et al (2013) using true anomaly and not
the eccentric anomaly and using a different method.

##  5 Test case description (nominal trajectories)

As a test case, we consider a low-thrust Earth-Mars orbital transfer. A
spacecraft having mass $m_{0}=1000$ (kg) and equipped with a propulsion system
able to deliver a constant $T_{max}=0.3$ (N) with a specific impulse of
$I_{sp}=2500$ (s) starts its transfer from Earth’s orbit to reach, after a
time of flight $\Delta t$, Mars’ orbit. Both starting and target orbits are
considered Keplerian. The optimal transfer strategy is found by solving the
resulting TPBVP where we seek the time of flight $\Delta t$, departure
eccentric anomaly $E_{0}$, arrival eccentric anomaly $E_{f}$ and departure co-
state variables $\lambda_{0}$ that satisfy the boundary conditions, the
dynamics Eq.(7), the transversality condition on free time Eq.(8), free mass
Eq.(9) and free anomalies Eq.(12).

[width=0.8]nominal_trajs.pdf

Figure 1: The two nominal trajectories considered.

The TPBVP is solved first for quadratic control optimality (QOC - $\alpha=0$).
Once the optimal trajectory is found, the problem is then solved with
gradually increasing $\alpha$ until mass-optimal control (MOC) is achieved
$\alpha=1$, at which point both solutions (nominal trajectories visualised in
Figure 1) are stored as:

$\mathbf{T}_{nom}=\left\\{(\mathbf{x},\bm{\lambda},t)_{nom,j}\right\\},j=0,\dots,J_{nom}$

where the grid points are determined by the adaptive integrator used. This
will result in more points allocated in areas where the dynamics gradient is
higher, which are also areas where we would want to have more training samples
and is thus considered as an appropriate mechanism.

##  6 Generating the QOC data set

As we want to train a deep neural architecture to represent the optimal
control around our nominal trajectories, we need to generate a large data set
describing the functional relationship that is to be learned. Such a
functional relationship is the optimal state feedback
($\mathbf{x}\to\mathbf{u}^{*}$) which can be computed along single
trajectories neighbouring the nominal ones at the cost of solving the
resulting TPBVPs. In order to make use of previously computed solutions
(starting from the nominal trajectories) to help the convergence of the TPBVP
solver, a continuation technique is then employed where new TPBVPs are created
using as initial states those of the nominal trajectories perturbed along a
random walk and as initial co-states those computed from previous steps. This
allows to generate efficiently the database continuing one starting solution,
i.e. the nominal trajectory. The TPBVP is solved by means of a simple single
shooting method. The pseudo algorithm used to fill the data-set is given in
Algorithm LABEL:alg:randomwalk which is run starting from 10 different points
equally spaced along the nominal trajectory and using $\alpha=0$. Let the
final data-set

$\mathcal{T}_{\mbox{QOC}}=\left\\{\mathbf{T}_{i}\right\\},i=0,\dots,I$ (13)

contain $I$ trajectories

$\mathbf{T}_{i}=\left\\{(\mathbf{x},\bm{\lambda},t)_{i,j}\right\\},j=0,\dots,J_{i}$
(14)

sampled in $J_{i}$ nodes (as determined by the outcome of an adaptive step
numerical integration) each recording the states $\mathbf{x}$, the co-states
$\bm{\lambda}$, and the times $t$. We refer to this data set as QOC (quadratic
optimal control).

[width=0.8]moc_data2.png

Figure 2: Visualisation of the MOC data set containing $308,579$ optimal
state-control pairs around the nominal trajectory.

###  6.1 Generating the MOC data set

A second data set containing mass optimal trajectories is obtained via
continuation over the parameter $\alpha$ (homotopy). Iterating through the
initial states and co-states in $\mathcal{T}_{\mbox{QOC}}$, an attempt is made
to solve, starting from that guess, directly a mass-optimal control
($\alpha=1$). If the resulting trajectory is feasible, it is stored and the
next data set entry is considered. If the trajectory is not feasible, the
homotopy parameter is decreased, and the optimisation is reattempted. If
successful, the homotopy parameter is saved as the current best and its
current value is increased. If, instead, the TPBVP solver is not converging,
the homotopy parameter is decreased, and the algorithm continues. The pseudo
code describing the homotopy method is outlined in Algorithm LABEL:alg:homo
and results in a new data set $\mathcal{T}_{\mbox{MOC}}$ containing mass
optimal trajectories and visualised in Figure 2.

##  7 Learning details and network architecture

We consider for both of the data-sets, $\mathcal{T}_{\mbox{MOC}}$ and
$\mathcal{T}_{\mbox{QOC}}$, the optimal state-action pairs
($\mathbf{x},\mathbf{u}$) and we build a deep model of the relation
$\mathbf{u}(\mathbf{x})$, i.e. the optimal state feedback. The dimension of
the values we have to model is $D=3$ as they represent a three dimensional
vector (i.e. the thrust vector). We choose to represent such a vector in its
polar coordinates so that:

$\mathbf{u}=u\left[\sin{\theta}\cos{\phi},\sin{\theta}\sin{\pi},\cos{\theta}\right]$

where $\theta$ is the polar angle, $\phi$ is the azimuth angle and $u$ is the
throttle magnitude.

In previous work Sánchez-Sánchez and Izzo (2016), it was found that a
feedforward, fully-connected neural network with 3 hidden layers and 32
units/layer could satisfactorily represent the optimal guidance profile for a
problem with simpler dynamics. Starting from that knowledge, and from the fact
that we are here considering a higher dimensional case with a higher degree of
complexity, we experimented with deeper and wider networks. The result of our
investigation indicated that peak performance could be obtained by a neural
network with four hidden layers and 200 units per layer. Larger networks
either matched this performance or performed worse. We settled on the
Rectified Linear Unit (ReLU) as the activation function for the hidden layers.
In addition, we used a hyperbolic tangent activation function for the output
layer. Thus we scaled the targets in the dataset to the range $[-1,1]$. We
initialised the network’s weights using the heuristic in Glorot and Bengio
(2010) in which random values are sampled from a uniform distribution close to
zero. Biases were initialised to zero. We also normalised our input data such
that features had a mean of zero and standard deviation of one; this helped to
speed up the optimisation process. We used the Adam training algorithm Kingma
and Ba (2014), a variant of stochastic gradient descent. Weights were updated
based on a minibatch of 64 training examples. For the performance metric we
use the commonly used mean squared error loss function (MSE). A 10% split is
used to define training and validation data. We used an adaptive learning rate
starting from an initial value of $10^{-3}$. The learning rate was reduced by
a factor of 10 when we observed a plateau in the training loss for a period
greater than 10 epochs. Furthermore, we stopped training when we observed a
plateau in the training loss for greater than 50 epochs (‘early stopping’).
The plateau was defined to be a decrease in loss less than $10^{-4}$. We
experimented with different values for the hyperparameters described above
before settling on the final configuration.

Network | MSE | Epochs
---|---|---
u | $\phi$ | $\theta$
$\mathcal{N}^{QOC}_{u,\phi,\theta}$ | 0.0178 | 0.0690 | 0.0870 | 453
$\mathcal{N}^{QOC}_{\phi,\theta}$ | _ | 0.0755 | 0.0907 | 461
$\mathcal{N}^{QOC}_{u}$ | 0.0075 | _ | _ | 247
$\mathcal{N}^{QOC}_{\phi}$ | _ | 0.0698 | _ | 394
$\mathcal{N}^{QOC}_{\theta}$ | _ | _ | 0.0976 | 503

Network | MSE | Epochs
---|---|---
u | $\phi$ | $\theta$
$\mathcal{N}^{MOC}_{u,\phi,\theta}$ | 0.7250 | 0.0255 | 0.0337 | 794
$\mathcal{N}^{MOC}_{\phi,\theta}$ | _ | 0.0212 | 0.0289 | 316
$\mathcal{N}^{MOC}_{u}$ | 0.7740 | _ | _ | 721
$\mathcal{N}^{MOC}_{\phi}$ | _ | 0.0204 | _ | 248
$\mathcal{N}^{MOC}_{\theta}$ | _ | _ | 0.0305 | 324

Table 1: MSE for neural networks trained on the MOC and QOC data sets. The
networks differ on the output dimension.

[width=]dt_learning_curves.pdf

Figure 3: Loss history for the networks $\mathcal{N}^{QOC}_{u,\phi,\theta}$
and $\mathcal{N}^{MOC}_{u,\phi,\theta}$. Note the jumps correspond to the
learning rate being reduced.

[width=]dt_qoc.pdf

Figure 4: Network $\mathcal{N}^{QOC}_{u,\phi,\theta}$ predictions on four
different transfers from the QOC test set.

[width=]dt_moc.pdf

Figure 5: Network $\mathcal{N}^{MOC}_{u,\phi,\theta}$ predictions on four
different transfers from the MOC test set.

##  8 Results

We train different networks to predict the three outputs both concurrently or
separately and on both $\mathcal{T}_{\mbox{MOC}}$ and
$\mathcal{T}_{\mbox{QOC}}$ data sets. We indicate each network using
superscripts and subscripts. For example, $\mathcal{N}^{MOC}_{\phi,\theta}$
indicates the network trained to model the variables $\phi$ and $\theta$ in
the $\mathcal{T}_{\mbox{MOC}}$ data set. We report in Table 1 the obtained
validation loss in the cases tested. We note that the best results for the
output variables is obtained from different models indicating the absence of a
clear trend to be exploited, while one could expect that predicting the entire
thrust vector at once would bring advantages as the network would be able to
share some of the weights as part of the necessary computations. We also
report in Figure 3 the loss during the training for the representative case of
the $\mathcal{N}^{MOC}_{u,\phi,\theta}$ and
$\mathcal{N}^{QOC}_{u,\phi,\theta}$ networks showing a relatively standard
trend. We finally show four different optimal transfers and plot alongside
them the optimal values against the predicted network values. The results are
visualised in Figure 4 and Figure 5.

It appears that a deep network is able to represent the optimal guidance
structure of the Earth-Mars transfer quite satisfactorily introducing errors
that are, on average, rather small. One should also keep in mind that the
actual quantitative value of the losses can be improved by further fine-tuning
of the hyperparameters used in training the models. In general, we observe the
prediction of the throttle $u$ is better attained in the QOC case. This is not
surprising as the structure of the optimal control in the MOC case is highly
discontinuous being a bang-bang control. We further note that data points
close to the switching points are very difficult to predict correctly and when
erroneously predicted result in big contributions to the overall loss. These
points form a majority in the datasets; a consequence of the time grid defined
by the numerical integrator for the TPBVP shooting method solver. This is the
reason for the MSE on the $u$ prediction to be of a greater magnitude in the
case of the MOC networks, while the visualised plot of the predictions for the
same cases are actually returning a much better scenario.

A definitive assessment on the capability of the trained deep models to steer
correctly the spacecraft has to be conducted, in a similar fashion to what
done in Sánchez-Sánchez and Izzo (2016), by integrating forward the whole
spacecraft dynamics considering the deep network predictions as a control and
will be part of a separate future publication.

## References

  * Addis et al (2011) Addis B, Cassioli A, Locatelli M, Schoen F (2011) A global optimization method for the design of space trajectories. Computational Optimization and Applications 48(3):635–652
  * Ampatzis and Izzo (2009) Ampatzis C, Izzo D (2009) Machine learning techniques for approximation of objective functions in trajectory optimisation. In: Proceedings of the IJCAI-09 workshop on artificial intelligence in space, pp 1–6
  * Biesbroek and Ancarola (2002) Biesbroek RG, Ancarola BP (2002) Optimization of Launcher Performance and Interplanetary Trajectories for Pre-Asessment Studies. In: IAF Abstracts, 34th COSPAR Scientific Assembly
  * Cassioli et al (2012) Cassioli A, Di Lorenzo D, Locatelli M, Schoen F, Sciandrone M (2012) Machine learning for global optimization. Computational Optimization and Applications 51(1):279–303
  * Ceriotti and Vasile (2010) Ceriotti M, Vasile M (2010) MGA trajectory planning with an ACO-inspired algorithm. Acta Astronautica 67(9):1202–1217
  * de Croon and Izzo (2012) de Croon G, Izzo D (2012) Real-time landing based on optimality principles and vision. In: 23rd International Symposium on Space Flight Dynamics (ISSFD)
  * Dachwald (2004) Dachwald B (2004) Low-thrust trajectory optimization and interplanetary mission analysis using evolutionary neurocontrol. PhD thesis, Doctoral thesis, Universität der Bundeswehr München Fakultät für Luft-und Raumfahrttechnik
  * Dachwald and Ohndorf (2019) Dachwald B, Ohndorf A (2019) Global Optimization of Continuous-Thrust Trajectories Using Evolutionary Neurocontrol. In: Fasano G, Pinter J (eds) Modeling and Optimization in Space Engineering - 2018, Springer International Publishing, Switzerland (2019)
  * Deb et al (2007) Deb K, Padhye N, Neema G (2007) Interplanetary trajectory optimization with swing-bys using evolutionary multi-objective optimization. In: International Symposium on Intelligence Computation and Applications, Springer, pp 26–35
  * Di Lizia and Radice (2004) Di Lizia P, Radice G (2004) Advanced global optimisation for mission analysis and design. Final Report Ariadna id 04/4101
  * Dueri et al (2016) Dueri D, Açıkmeşe B, Scharf DP, Harris MW (2016) Customized real-time interior-point methods for onboard powered-descent guidance. Journal of Guidance, Control, and Dynamics
  * Elsayed et al (2011) Elsayed SM, Sarker RA, Essam DL (2011) GA with a new multi-parent crossover for solving IEEE-CEC2011 competition problems. In: Evolutionary Computation (CEC), 2011 IEEE Congress on, IEEE, pp 1034–1040
  * Englander (2013) Englander J (2013) Automated trajectory planning for multiple-flyby interplanetary missions. University of Illinois at Urbana-Champaign
  * Gad and Abdelkhalik (2011) Gad A, Abdelkhalik O (2011) Hidden genes genetic algorithm for multi-gravity-assist trajectories optimization. Journal of Spacecraft and Rockets 48(4):629–641
  * Gage et al (1995) Gage P, Braun R, Kroo I (1995) Interplanetary trajectory optimization using a genetic algorithm. Journal of the Astronautical Sciences 43(1):59–76
  * Glorot and Bengio (2010) Glorot X, Bengio Y (2010) Understanding the difficulty of training deep feedforward neural networks. In: Proceedings of the thirteenth international conference on artificial intelligence and statistics, pp 249–256
  * Grigoriev and Zapletin (2013) Grigoriev I, Zapletin M (2013) Choosing promising sequences of asteroids. Automation & Remote Control 74(8)
  * Hennes and Izzo (2015) Hennes D, Izzo D (2015) Interplanetary Trajectory Planning with Monte Carlo Tree Search. In: IJCAI, pp 769–775
  * Hennes et al (2016) Hennes D, Izzo D, Landau D (2016) Fast approximators for optimal low-thrust hops between main belt asteroids. In: Computational Intelligence (SSCI), 2016 IEEE Symposium Series on, IEEE, pp 1–7
  * Islam et al (2012) Islam SM, Das S, Ghosh S, Roy S, Suganthan PN (2012) An adaptive differential evolution algorithm with novel mutation and crossover strategies for global numerical optimization. IEEE Transactions on Systems, Man, and Cybernetics, Part B (Cybernetics) 42(2):482–500
  * Izzo (2010) Izzo D (2010) Global optimization and space pruning for spacecraft trajectory design. Spacecraft Trajectory Optimization 1:178–200
  * Izzo et al (2007) Izzo D, Becerra VM, Myatt DR, Nasuto SJ, Bishop JM (2007) Search space pruning and global optimisation of multiple gravity assist spacecraft trajectories. Journal of Global Optimization 38(2):283–296
  * Izzo et al (2013) Izzo D, Simões LF, Märtens M, De Croon GC, Heritier A, Yam CH (2013) Search for a grand tour of the Jupiter galilean moons. In: Proceedings of the 15th annual conference on Genetic and evolutionary computation, ACM, pp 1301–1308
  * Izzo et al (2014a) Izzo D, Hennes D, Riccardi A (2014a) Constraint handling and multi-objective methods for the evolution of interplanetary trajectories. Journal of Guidance, Control, and Dynamics
  * Izzo et al (2014b) Izzo D, Simoes LF, Yam CH, Biscani F, Di Lorenzo D, Addis B, Cassioli A (2014b) GTOC5: results from the European Space Agency and University of Florence. Acta Futura 8:45–55
  * Izzo et al (2015) Izzo D, Getzner I, Hennes D, Simões LF (2015) Evolving solutions to TSP variants for active space debris removal. In: Proceedings of the 2015 Annual Conference on Genetic and Evolutionary Computation, ACM, pp 1207–1214
  * Izzo et al (2016) Izzo D, Hennes D, Simões LF, Märtens M (2016) Designing complex interplanetary trajectories for the global trajectory optimization competitions. In: Space Engineering, Springer, pp 151–176
  * Janin and Gomez-Tierno (1985) Janin G, Gomez-Tierno M (1985) The genetic algorithms for trajectory optimization. In: Stockholm International Astronautical Federation Congress
  * Kingma and Ba (2014) Kingma DP, Ba J (2014) Adam: A method for stochastic optimization. CoRR abs/1412.6980, URL http://arxiv.org/abs/1412.6980, 1412.6980
  * Lavagna (2007) Lavagna MR (2007) Multi-objective pso for interplanetary trajectory design. In: Proceedings of the 9th annual conference on Genetic and evolutionary computation, ACM, pp 175–175
  * Lee et al (2005) Lee S, von Allmen P, Fink W, Petropoulos A, Terrile R (2005) Multi-objective evolutionary algorithms for low-thrust orbit transfer optimization. In: Genetic and Evolutionary Computation Conference (GECCO 2005)
  * Luo et al (2005) Luo YZ, Tang GJ, Zhou LN (2005) Simulated annealing for solving near-optimal low-thrust orbit transfer. Engineering Optimization 37(2):201–216
  * Mereta et al (2017) Mereta A, Izzo D, Wittig A (2017) Machine Learning of Optimal Low-Thrust Transfers Between Near-Earth Objects. In: International Conference on Hybrid Artificial Intelligence Systems, Springer, pp 543–553
  * Myatt et al (2004) Myatt D, Becerra VM, Nasuto SJ, Bishop J (2004) Advanced global optimisation for mission analysis and design. Final Report Ariadna id 04/4101
  * Olds et al (2007) Olds AD, Kluever CA, Cupples ML (2007) Interplanetary mission design using differential evolution. Journal of Spacecraft and Rockets 44(5):1060–1070
  * Pan et al (2013) Pan B, Chen Z, Lu P, Gao B (2013) Reduced transversality conditions in optimal space trajectories. Journal of Guidance, Control, and Dynamics
  * Pontani and Conway (2010) Pontani M, Conway BA (2010) Particle swarm optimization applied to space trajectories. Journal of Guidance, Control and Dynamics 33(5):1429–1441
  * Pontryagin et al (1962) Pontryagin LS, Boltyanskii V, Gamkrelidze R (1962) Ef mishchenko the mathematical theory of optimal processes. New York: Interscience
  * Radice and Olmo (2006) Radice G, Olmo G (2006) Ant colony algorithms for two-impulse interplanetary trajectory optimization. Journal of Guidance Control and Dynamics 29(6):1440
  * Rauwolf and Coverstone-Carroll (1996) Rauwolf GA, Coverstone-Carroll VL (1996) Near-optimal low-thrust orbit transfers generated by a genetic algorithm. Journal of Spacecraft and Rockets 33(6):859–862
  * Rogata et al (2003) Rogata P, Di Sotto E, Graziano M, Graziani F (2003) Guess value for interplanetary transfer design through genetic algorithms. Advances in the Astronautical Sciences 114:613–627
  * Sánchez-Sánchez and Izzo (2016) Sánchez-Sánchez C, Izzo D (2016) Real-time optimal control via Deep Neural Networks: study on landing problems. arXiv preprint arXiv:161008668
  * Sánchez-Sánchez et al (2016) Sánchez-Sánchez C, Izzo D, Hennes D (2016) Learning the optimal state-feedback using deep networks. In: Computational Intelligence (SSCI), 2016 IEEE Symposium Series on, IEEE, pp 1–8
  * Schiavone et al (2012) Schiavone G, Izzo D, Simões LF, De Croon GC (2012) Autonomous spacecraft landing through human pre-attentive vision. Bioinspiration & biomimetics 7(2):025,007
  * Schlueter (2014) Schlueter M (2014) MIDACO software performance on interplanetary trajectory benchmarks. Advances in Space Research 54(4):744–754
  * Schlueter et al (2013) Schlueter M, Erb SO, Gerdts M, Kemble S, Rückmann JJ (2013) MIDACO on MINLP space applications. Advances in Space Research 51(7):1116–1131
  * Schmidhuber (2015) Schmidhuber J (2015) Deep learning in neural networks: An overview. Neural networks 61:85–117
  * Sentinella and Casalino (2009) Sentinella MR, Casalino L (2009) Hybrid evolutionary algorithm for the optimization of interplanetary trajectories. Journal of Spacecraft and Rockets 46(2):365
  * Simões et al (2014) Simões LF, Izzo D, Haasdijk E, Eiben AE (2014) Self-adaptive genotype-phenotype maps: neural networks as a meta-representation. In: International Conference on Parallel Problem Solving from Nature, Springer, pp 110–119
  * Simões et al (2017) Simões LF, Izzo D, Haasdijk E, Eiben A (2017) Multi-rendezvous Spacecraft Trajectory Optimization with Beam P-ACO. In: European Conference on Evolutionary Computation in Combinatorial Optimization, Springer, pp 141–156
  * Stracquadanio et al (2011) Stracquadanio G, La Ferla A, De Felice M, Nicosia G (2011) Design of Robust Space Trajectories. In: SGAI Conf., Springer, pp 341–354
  * Vasile et al (2010) Vasile M, Minisci E, Locatelli M (2010) Analysis of some global optimization algorithms for space trajectory design. Journal of Spacecraft and Rockets 47(2):334
  * Vinkó and Izzo (2008) Vinkó T, Izzo D (2008) Global optimisation heuristics and test problems for preliminary spacecraft trajectory design. Eur Space Agency, Adv Concepts Team, ACT Tech Rep, id: GOHTPPSTD
  * Wolpert and Macready (1997) Wolpert DH, Macready WG (1997) No free lunch theorems for optimization. IEEE transactions on evolutionary computation 1(1):67–82
  * Yao et al (2017) Yao W, Luo J, Macdonald M, Wang M, Ma W (2017) Improved Differential Evolution Algorithm and Its Applications to Orbit Design. Journal of Guidance, Control, and Dynamics


# Analysis of Independent Learning in Network Agents: A Packet Forwarding Use
Case

Abu Saleh Md Tayeen  Dept. of Computer Science
New Mexico State University
Las Cruces, NM, U.S.A.
<EMAIL_ADDRESS>Milan Biswal  Dept. of Computer Science
New Mexico State University
Las Cruces, NM, U.S.A.
<EMAIL_ADDRESS>Abderrahmen Mtibaa  Dept. of Computer Science
University of Missouri–St Louis
St. Louis, Missouri, U.S.A.
<EMAIL_ADDRESS>Satyajayant Misra  Dept. of Computer Science
New Mexico State University
Las Cruces, NM, U.S.A.
<EMAIL_ADDRESS>
###### Abstract

Multi-Agent Reinforcement Learning (MARL) is nowadays widely used to solve
real-world and complex decisions in various domains. While MARL can be
categorized into independent and cooperative approaches, we consider the
independent approach as a simple, more scalable, and less costly method for
large-scale distributed systems, such as network packet forwarding. In this
paper, we quantitatively and qualitatively assess the benefits of leveraging
such independent agents learning approach, in particular IQL-based algorithm,
for packet forwarding in computer networking, using the Named Data Networking
(NDN) architecture as a driving example. We put multiple IQL-based forwarding
strategies (IDQF) to the test and compare their performances against very
basic forwarding schemes and simple topologies/traffic models to highlight
major challenges and issues. We discuss the main issues related to the poor
performance of IDQF, and quantify the impact of these issues on isolation when
training and testing the IDQF models under different model tuning parameters
and network topologies/characteristics.

###### Index Terms:

NDN Forwarding, Reinforcement Learning, Multi-agents, Independent Learning,
Forwarding Strategy

##  I Introduction

Reinforcement learning (RL) [1] helps explore and interact with local
environment to acquire knowledge to solve real-world problems, especially
complex multi-agent systems such as traffic signal control [2], controlling
drones [3], robot soccer [4], machines in factories [5], and automated trading
[6]. RL-based solutions often deploy (i) a single centralized agent [7], or
(ii) multiple cooperative agents [8] to acquire the knowledge of such
distributed systems. While these solutions may not scale for large distributed
systems due to the amount of data exchange (i.e., overhead), and additional
command and control delays, we investigate in this paper the efficacy of non-
cooperative/autonomous agents approach, Independent Q-Learning (IQL) [9, 2]
and quantitatively and qualitatively assess its performance using packet
forwarding in an ISP-scale network as a use-case.

ISP networks are often large-scale and complex, which renders using (i) a
single centralized agent (e.g., SDN), or (ii) multiple cooperative agents
(e.g., message exchange among routers) very costly–large overhead,
single/multiple point of failure, and hard to maintain/upgrade. We assess the
potential benefits and efficacy of using independent agents that learn to make
forwarding decisions using only localized information (e.g., incoming,
outgoing traffic state). We choose to implement all our IQL-based models using
Named Data Networking (NDN), which provides many “useful” features such as an
easy to use forwarding plane, and context awareness decisions at the routers,
to test our models, configurations, and tuning parameters. We provide a short
NDN primer and a motivation behind our NDN choice in Section II. However, we
believe that the paper’s discussions, outcomes, and conclusions remain valid
for IP architecture networks.

We compare the performance of a neural network equipped IQL-based forwarding
strategy, IDQF (Independent Deep Q-Network Forwarding) against the most basic
(i.e., the simplest) forwarding scheme in NDN, namely best_route (BR) [10]. We
show that BR outperforms IDQF in most of the scenarios by simply choosing one
path, while IDQF often fails to select the most appropriate path towards the
producer.

This comparison mainly aims at highlighting the challenges of IQL-based
approach in such environment. We ran multiple experiments with various
parameters, tuning, and setups to outline the challenges that led to
unsuccessful decisions made by IDQF. The identified challenges (not a complete
list) include: (i) the non-stationarity of the system due to the lack of
global synchronization leading to triggered changes as soon as one agent
changes its policy; (ii) the generalization of the training setup that is
difficult due to network delays and frequent changes in the environment; (iii)
the inability to assess the reward of an action immediately due to the
variations of network delays; and (iv) the partial observability that prevents
the routers to see overlapping paths, conflicting flows, etc. and make optimal
decisions. We discuss these challenges in detail and present experiments
and/or examples to support the impact of these challenges on the performance
of IDQF in particular, and distributed autonomous multi-agent RL schemes in
general.

The rest of this paper is organized as follows. Section II includes background
and related works for this work. Section III summarizes the design of the IDQF
forwarding mechanism and describe the simulation results. In Section IV, we
discuss the challenges regarding the IDQF forwarding scheme and autonomous
multi-agent RL schemes. Finally, Section V concludes the paper and suggest
future work.

##  II Background And Related Work

###  II-A ICN/NDN Background

Different from IP networks that use IP addresses (hosts identities) to inform
intermediate routers’ forwarding decisions, the fundamental idea of Named Data
Networking (NDN) architecture is to retrieve the named-pieces of information
(named network-layer packets), from any node that can provide it. In NDN,
routers are equipped with a content store (CS), a pending interest table
(PIT), and a forwarding information base (FIB). The FIB (similar to the
forwarding table in IP routers) gets populated using a routing algorithm. Any
node/router that receives a special packet (known as Interest packet)
requesting for a content (Data packet) performs a CS-lookup on the content
name. If the content is not available in the CS, the router performs a lookup
in its PIT to check whether there is an existing request for the content that
the router has forwarded into the network. If the PIT lookup is successful,
the router adds the incoming Interest packet’s interface to the PIT entry
(interest aggregation) and drops the Interest. If no PIT match is found, the
router creates a new PIT entry for the Interest. It then utilizes a forwarding
strategy to select the best interface to upstream router(s) from multiple
interfaces obtained from the FIB and forwards the Interest packet in the
direction of the data source(s).

###  II-B Related Work

There exists a good number of literature which applied RL to routing in
networks [11]. We mainly outline a few research works that applied RL for
designing forwarding strategies in NDN. In [12, 13], authors proposed
Q-routing based strategies that alternates between exploration and
exploitation phases to learn the best interface in terms of a cost (e.g.
delay) during dynamic network conditions. However, the switching of the phases
are dependent on some parameters and hampers the ability of the strategies to
react quickly to network changes such as congestion. Moreover, it requires
altering the default NDN packet formats and data structures to calculate
experience information for Q-routing, which violates the principle of data
immutability in NDN. The works [14] and [15] exploited neural networks to
learn the mapping from the state of the network condition to the action of
choosing the optimal interface for forwarding interests. However, there is no
balancing of load across interfaces. Moreover, the strategies were verified
using a single agent topology, which makes it difficult to assess the efficacy
in larger/real networks.

##  III IQL-based Stateful Forwarding

In this section, we illustrate our IQL-based NDN forwarding scheme, IDQF where
each node makes autonomous forwarding decisions based on their own
observations and rewards. We also assess the performance of this scheme by
comparing it to the best_route (BR) [10].

###  III-A Design of IDQF Forwarding Scheme

The goal of the agent in each router is to autonomously learn a decision
policy to forward Interest packets through the optimal next-hop interface. The
schematics of the IDQF forwarding scheme is shown in Figure 1. The agent
learns/discovers the optimal forwarding policy indirectly using the action-
value function. In the vanilla Q-learning approach, the state-action pairs are
mapped to corresponding Q-values (maximum expected rewards) by maintaining a
Q-table. For a communication network with numerous routers and communication
flows, the size and complexity of state-action pairs or the corresponding
Q-table will be very large, which poses a significant challenge for IQL to
learn the optimal decision policy. To resolve this, we approximate the action-
value function using neural networks and use the Deep-Q-Network (DQN) training
algorithm [16].

Figure 1: High level design of the IDQF forwarding strategy implemented at any
router in the network. Periodically, routers update their states when
receiving incoming Interests, and measure rewards when Data packets are
received from the producer.

As shown in Figure 1, the router $RA$ is equipped with neural network DQN. At
decision time $t$, $RA$, calculates the state of its available next hop
interfaces, Face1 and Face2 by observing some local parameters such as round
trip times (RTTs), interest satisfaction ratios, number of re-transmitted
packets collected between time $t$ and $t-\delta t$. Here, $\delta t$ is the
time period for sensing and updating the weights of the DQN. The state is fed
into the DQN as input. DQN then predicts either Face1 or Face2 as an action.
The $RA$ agent then forwards the current Interest packet and all subsequent
Interest packets through the predicted interface until next decision time
occurs at $t+\delta t$. The chosen action results in a change of state for the
corresponding interface. $RA$ also gets a feedback for its action, which is
the reward. We define the reward, $RW(t)$ at time $t$ by the following
equation:

$RW(t)=-\Big{[}\Big{(}\frac{1}{M}\sum_{k=(t-\delta
t)}^{t}{RTT(k)}\Big{)}+C\times R(t)\Big{]}$ (1)

Here, $M$ is the number of Data packets received between $t$ and $t-\delta t$,
$RTT(k)$ is the Round-Trip Time (RTT) of the received Data packet at time $k$,
$R(t)$ is the number of re-transmitted packets which were actually transmitted
through the interface chosen for the $\delta t$ period, and $C$ is a penalty
constant multiplier that aims at discouraging choosing interface with high
drop ratios (we choose $C=4$ seconds in our case). Note that in case of drops
of Interests/Data packets, there will be no RTT measurements, thus the penalty
constant helps in assigning a maximum delay for every packet dropped. The $RA$
agent uses the reward to improve its decision policy.

###  III-B Evaluation of IDQF

We use _ndnSIM_ [17], the NDN network simulator, to evaluate the forwarding
scheme performance. To implement the RL agents, we used a concurrently running
daemon process that interfaces with the _ndnSIM_ and exchanges training data
and prediction of the next hop interface for each node. In this paper, we will
share the results using the Sprint topology shown in Figure 2. We have tested
the performance using multiple network topologies (grid, tree), and performed
tuning of different parameters of DQN.

Figure 2: Sprint topology used in our simulation; propagation delays (in
milliseconds) are shown as link weights; link capacities vary from 1Mbps at
the core (black links), 2Mbps (grey) and 5Mbps at the edge (blue)

We consider a simple scenario, where only one consumer ($C_{0}$) requests data
from only one producer ($P_{0}$) at a fixed rate $r$. Note that we used such a
simple scenario to highlight challenges that may arise even in a very
predictable network traffic with no loss, errors, contention, etc. The Sprint
topology in Figure 2 consists of $11$ router nodes and $18$ links among them.
Each router node has a queue size of 100 packets. In our simple scenario, we
deploy the agents in nodes $N_{3}$, $N_{4}$, $N_{6}$, and $N_{9}$ only. Each
of these agents have their own IDQF-strategy. In our results, we have
considered the size of the state vector to be two and chose the best two
interfaces at each agent (those with the least delays to the producer). Note
that we have used and tested different sizes which result in performance
degradation when we scale the exploration beyond two interfaces.

In our simulations, the DQN of each agent has one hidden layer, one input
layer and one output layer. Each hidden layer consists of 32 nodes, and the
layers are fully connected with each other. While training the node agents we
tuned several parameters of DQN. The first parameter is the decay rate, which
controls the number of times the agents will perform exploration by randomly
choosing actions. Second, we adjusted the learning rate to control the
convergence speed of DQN. Third, we varied the capacity of the replay buffer
which governs the diversity of samples available to the agent at any moment.

We trained the DQN of agents with decay rate=$0.001$ and learning rate=$1.0$
for 50 episodes where an episode corresponds to one instance of simulation of
the network environment with duration 60 seconds. The agents were trained
separately for five different request rate starting from low Interest rate
(41% of network capacity) to high Interest rate (123% of network capacity) in
packets/second (pkts/s). The interest rates are 100 pkts/s (0.82 Mbps), 150
pkts/s (1.23 Mbps), 200 pkts/s (1.64 Mbps), 250 pkts/s (2.05 Mbps), and 300
pkts/s (2.46 Mbps). We deployed the agents for each Interest rate and ran
simulations five times to get average results. We compare the throughput and
average delay measured in $C_{0}$ for both BR and IDQF scheme and present the
results in Figure 3. The height of the bar for each point in Figure 3
represents the standard deviation of throughput and delay measurements at
every second of simulation.

(a) Throughput measurements at the consumer node $C_{0}$

(b) End-to-end delay measurements at the consumer node $C_{0}$

Figure 3: Throughput and delay comparison of _best_route_ (BR) and IDQF
performance; Error bars indicate the variation of the measurements during the
simulation

Figure 3 shows that the IDQF strategy does not outperform even the most basic
strategy, best_route (BR). BR always uses the interface with the shortest
delay regardless of congestion. It will react and switch an interface only if
the best route fails (after at least one timeout event). Note that we share
the IDQF configuration that provides the best results for this scenario. As
shown in Figure 3(a), BR outperforms IDQF by 10%-36% in terms of throughput
for all Interest rates. We notice that IDQF often fails to choose the best
interface, provided the cost (e.g., latency from the router to the producer)
of using each interface. Similarly, in Figure 3(b), we find that IDQF strategy
fails to deliver data faster than BR in most cases, especially at lower
Interest rates. Interests are often forwarded to interfaces that lead to
congestion and drops, resulting in an average of 2 to 3 re-transmissions per
Interest.

To summarize, IDQF strategy, while tested even with a simple scenario with a
single flow and limited network variability, fails to converge into a stable
state that learns the non-changing network characteristics and makes decisions
that can often be solved by simple utility functions. In this study, we have
identified a list of challenges (discussed in the next section) which advocate
that reinforcement learning based approaches driven by autonomous decision
making paradigm will result in similar non-optimal or sub-optimal solutions.
These challenges are generalized beyond NDN, and may apply to IP and other
network architectures/systems.

##  IV Challenges of using autonomous agents in distributed systems

We will build on the results shown in the previous section that highlight the
sub-optimal performance of the IDQF strategy, and autonomous multi-agent RL
algorithms in general. In this section, we will discuss and analyze
qualitatively and quantitatively few challenges that contribute to this “low”
performance.

###  IV-A Training and Generalization of Agents

The exhaustiveness of the sample training set and the method of training play
a significant role in proper generalization of the neural network of an agent.
In a multi-agent setting, the goal of every agent is to maximize its
cumulative reward or the expected return, which is simply the sum of all
rewards the agent received starting from the initial state to a terminal state
of the environment. An agent can reach a terminal state if certain score is
obtained or the task is complete. To learn an optimized policy the agents are
iteratively trained through multiple episodes where each episode is a sequence
of states, actions, and rewards starting from an initial state to a terminal
state. Often, it is difficult to determine a terminal state or a target score
for the network environment. Therefore, we argue that setting the episode as a
period of time, which is long enough to incorporate actions, and their
rewards, will depend on the network topology and its characteristics (traffic,
congestion). Based on the choice of episode length (simulation duration),
there is a possibility of over-fitting or under-fitting of the neural networks
of the agents, that can often lead to bad generalization and sub-optimal
policy learning for the agents.

To better illustrate this challenge, let us consider the scenario depicted in
Figure 2. Agents closer to the consumer node would receive much more Interest
packets, resulting in more data points, when compared to others closer to the
producer node. Indeed, agents explore their interfaces by spreading received
Interests across their next hop nodes, resulting in fewer number of Interests
received by core routers which render learning more challenging at these
nodes. For instance, nodes such as $N_{4}$ will not get enough number of
Interests to build its training data within the short simulation duration. In
our controlled experiments with a fixed network topology, we tested and tuned
different simulation periods, yet failed to find optimal forwarding strategies
for the nodes.

The generalization of the agents also depends on the size of a memory buffer
called _experience replay_ [16] that is used to store experience data (i.e.,
state transition record as tuple $<$state, action, reward, next state$>$).
This replay buffer helps the neural networks to break correlation in
observation sequences and learn from more independently and identically
distributed past experiences to stabilise the training. Without replay buffer,
agents will update the weights of their neural networks based on only current
state, action, and reward. However, sometimes state transition tuples
generated at earlier stage of exploration constitute a good policy for the
agent. These records need to be stored and can be utilized for training the
neural network to achieve better generalization. During training with replay
buffer, the agents update the weights of their neural networks or DQN by
randomly sampling mini-batches of experience data from it.

To demonstrate the effect of replay buffer, we ran two experiments: one
involving agents trained without any replay buffer and another by training
agents with replay buffer. We show the cumulative rewards obtained by node
agent $N_{3}$ per episode of training in Figure 4(a).

(a) Cumulative reward of node agent $N_{3}$ with different sizes of replay
buffer

(b) Cumulative reward of node agent $N_{3}$ when trained with different number
of agents

(c) CDF of delays of node agent $N_{3}$

Figure 4: Demonstration of non-stationarity, effect of replay buffer size and
delay measurements during training

We show that the absence of replay buffer leads to higher oscillation which
makes the system more susceptible to variations in the input states resulting
in issues in convergence across all agents in the network. Deploying a replay
buffer (we have tested different buffer sizes) often exhibits less variations,
however, at the same time the cumulative reward does not move toward higher
values as training progresses through every episode. This results in
additional delays in reacting to any sudden change in the network
characteristic (e.g., input state, reward). As a result, experience replay is
difficult to tune in such autonomous multi-agent RL scenario.

###  IV-B Non-stationarity of the Environment

One of the major challenges for distributed autonomous multi-agent RL is the
_non-stationarity_. In a multi-agent setting of Markov Decision Process, the
state and reward of each agent depend on the actions of all agents and every
agent adapts its behaviors in the context of other co-adapting learning
agents. Since all agents are learning concurrently, the environment is not
stationary and constantly changing from the perspective of any given agent.
This phenomenon actually rules out any convergence guarantees.

To demonstrate the effect of non-stationarity we performed two experiments on
topology in Figure 2. In one experiment, we put an agent only in node $N_{3}$
and train the single agent which will learn to choose between interfaces to
node $N_{6}$ and $N_{9}$. In another experiment, we put agents in both $N_{3}$
and $N_{9}$ and train the two agents concurrently. Here, $N_{9}$ will learn to
select between its next hops $N_{8}$ and $N_{4}$. We assume that in both cases
all the other nodes without any agents know the forwarding next hop node to
get optimal performance. In Figure 4(b), we show the cumulative reward of node
agent $N_{3}$ for the two cases. From Figure 4(b), we observe that training
with two agents cause the cumulative rewards of $N_{3}$ to be less and to
fluctuate more compared to that of $N_{3}$ when training was done with one
agent. This is because with two agents the input state and reward of agent
$N_{3}$ depends on the actions (face selections) of agent $N_{9}$. Similarly,
the input state of agent $N_{9}$ is conditioned upon $N_{3}$’s distribution of
interests to $N_{9}$. As a result, during learning both agents keep changing
their policy of choosing an interface based on the actions taken by each
other. These actions can be conflicting (as shown by the large oscillations in
the two agent cumulative reward) and result in degradation of performance as
shown in our previous results. On the other hand, with only one agent $N_{3}$,
its state and rewards solely depend on its actions, resulting in faster
convergence and better performance. This phenomenon is amplified when: (i) all
agents are running simultaneously, (ii) the number of agents in the network
scales, (iii) the network experience multiple concurrent flows (e.g., many
consumers and producers). One of the ways to deal with non-stationarity is to
adopt a centralized architecture where agents will know the actions of one
another. However, it is infeasible for forwarding packets in large distributed
computer networks.

###  IV-C Partial Observability of the Environment

Another key factor in autonomous multi-agent RL is the _partial
observability_. In many real-world environments, agents cannot perceive the
complete state of the environment. For example, in our context of packet
forwarding, the agents at each router observe the state per interface per
producer in terms of some parameters such as RTT or drop ratio which indicate
the congestion status of their respective interfaces. However, agents cannot
distinguish certain events, which may occur beyond their local knowledge,
including if two or more of their available interfaces lead to overlapping
paths towards a given producer, or if a given interface will compete with
other distant flows (i.e., using portions of the path towards the producer).
All these situations may result in similar measured rewards, but require
different mitigation (e.g., drop packets, choose different interface, etc.).
Besides, the node agents neither know the actions (next hop choices) other
agents have taken nor can compute the global reward (overall delay in
consumer) by themselves. In order to find a globally optimal policy, the
agents need to take into account the states and behavior of other agents and
adapt to the joint behavior accordingly. Global observability, such as in SDN,
or message-exchange-based coordination is needed to overcome this problem. We
argue that even message-exchange-based coordination is challenging due to
delays in message exchanged across agents unless a mechanism such as separate
control channel is in place to transmit the messages out of band.

###  IV-D Inability to Assess Actual Reward

Most of the RL algorithms have been designed for environments where the
feedback or reward is instantly available upon execution of an action. These
algorithms cannot adequately solve the problems where actual rewards can not
be assessed or are not available immediately [18]. The delayed feedback
phenomenon is also visible in our packet forwarding problem. In any network,
delays occur due to many factors including transmission, propagation, queuing,
etc. Thus, measuring the state parameters and true reward value for the last
action can not be collected before the next decision step occurs. Since RTTs
could take from few milliseconds to seconds (especially due to drop and re-
transmissions), estimating the current state and reward, e.g., RTT or drops,
for any forwarded Interest packet is pushed by additional and often variable
delays. Throughout these waiting times, the agents will have to make
forwarding decisions for incoming Interest packets and assign some
approximated rewards against its state-action observations. These approximated
rewards make it difficult to determine which of the state-action pairs
contributed to the achievement of the optimization goal. This can lead to
agents obtaining sub-optimal policies. For example, the topology in Figure 2
shows that the theoretical RTT of the optimal path between router $N_{3}$ and
the producer is 49.51ms. However, as shown in Figure 4(c), based on our
experiments we measured delays ranging from 61ms (i.e., counting queuing and
other processing delays) to 1s (maximum delay before dropping PIT entries),
where more than 50% of the delay measurements exceeded 300ms. These large
variations in the delay measurements in a controlled experiment, where number
of flows, traffic, and routes are static, pushes autonomous multi-agent RL
approaches to make approximate forwarding decisions which are often far from
the optimal policy.

##  V Conclusion

In this paper, we applied IQL, the autonomous and non-cooperative multi-agent
RL approach to perform packet forwarding in NDN. We approximated the
Q-function of IQL using DQN. We performed many experiments in different
topologies with tuning of several parameters to achieve a better
generalization of the DQN. We show that our implemented forwarding scheme,
IDQF does not significantly improve the network performance compared to the
default forwarding strategy available in NDN. We have also presented several
challenges (e.g., non-stationarity, partial observability, etc.) of this
approach using examples and experiments to support its effect on the
performance of IDQF. A possible direction for our future work is to try fully
cooperative multi-agent reinforcement learning approaches to design forwarding
strategies.

## References

  * [1] R. S. Sutton, A. G. Barto _et al._ , _Introduction to reinforcement learning_. MIT press Cambridge, 1998, vol. 135.
  * [2] H. Wei, C. Chen, G. Zheng, K. Wu, V. Gayah, K. Xu, and Z. Li, “Presslight: Learning max pressure control to coordinate traffic signals in arterial network,” in _Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining_, 2019, pp. 1290–1298.
  * [3] S.-M. Hung and S. N. Givigi, “A q-learning approach to flocking with uavs in a stochastic environment,” _IEEE transactions on cybernetics_ , vol. 47, no. 1, pp. 186–197, 2016.
  * [4] D. Schwab, Y. Zhu, and M. Veloso, “Zero shot transfer learning for robot soccer,” in _Proceedings of the 17th International Conference on Autonomous Agents and MultiAgent Systems_ , 2018, pp. 2070–2072.
  * [5] H. Kazmi, J. Suykens, A. Balint, and J. Driesen, “Multi-agent reinforcement learning for modeling and control of thermostatically controlled loads,” _Applied energy_ , vol. 238, pp. 1022–1035, 2019.
  * [6] J. Jin, C. Song, H. Li, K. Gai, J. Wang, and W. Zhang, “Real-time bidding with multi-agent reinforcement learning in display advertising,” in _Proceedings of the 27th ACM International Conference on Information and Knowledge Management_ , 2018, pp. 2193–2201.
  * [7] G. Stampa, M. Arias, D. Sánchez-Charles, V. Muntés-Mulero, and A. Cabellos, “A deep-reinforcement learning approach for software-defined networking routing optimization,” _arXiv preprint arXiv:1709.07080_ , 2017\.
  * [8] J. K. Gupta, M. Egorov, and M. Kochenderfer, “Cooperative multi-agent control using deep reinforcement learning,” in _International Conference on Autonomous Agents and Multiagent Systems_. Springer, 2017, pp. 66–83.
  * [9] M. Tan, “Multi-agent reinforcement learning: Independent vs. cooperative agents,” in _Proceedings of the tenth international conference on machine learning_ , 1993, pp. 330–337.
  * [10] C. Yi, A. Afanasyev, I. Moiseenko, L. Wang, B. Zhang, and L. Zhang, “A case for stateful forwarding plane,” _Computer Communications_ , vol. 36, no. 7, pp. 779–791, 2013.
  * [11] Z. Mammeri, “Reinforcement learning based routing in networks: Review and classification of approaches,” _IEEE Access_ , vol. 7, pp. 55 916–55 950, 2019.
  * [12] R. Chiocchetti, D. Perino, G. Carofiglio, D. Rossi, and G. Rossini, “Inform: a dynamic interest forwarding mechanism for information centric networking,” in _Proceedings of the 3rd ACM SIGCOMM workshop on Information-centric networking_ , 2013, pp. 9–14.
  * [13] M. Zhang, X. Wang, T. Liu, J. Zhu, and Q. Wu, “Afsndn: A novel adaptive forwarding strategy in named data networking based on q-learning,” _Peer-to-Peer Networking and Applications_ , pp. 1–9, 2020.
  * [14] Y. Zhang, B. Bai, K. Xu, and K. Lei, “Ifs-rl: An intelligent forwarding strategy based on reinforcement learning in named-data networking,” in _Proceedings of the 2018 Workshop on Network Meets AI & ML_, 2018, pp. 54–59.
  * [15] Y. A. B. de Sena, K. L. Dias, and C. Zanchettin, “Dqn-af: Deep q-network based adaptive forwarding strategy for named data networking,” in _2020 IEEE Latin-American Conference on Communications (LATINCOM)_. IEEE, 2020, pp. 1–6.
  * [16] V. Mnih, K. Kavukcuoglu, D. Silver, A. A. Rusu, J. Veness, M. G. Bellemare, A. Graves, M. Riedmiller, A. K. Fidjeland, G. Ostrovski _et al._ , “Human-level control through deep reinforcement learning,” _Nature_ , vol. 518, no. 7540, pp. 529–533, 2015.
  * [17] S. Mastorakis, A. Afanasyev, and L. Zhang, “On the evolution of ndnSIM: an open-source simulator for NDN experimentation,” _ACM Computer Communication Review_ , Jul. 2017.
  * [18] S. Shen and M. Chi, “Reinforcement learning: the sooner the better, or the later the better?” in _Proceedings of the 2016 conference on user modeling adaptation and personalization_ , 2016, pp. 37–44.