Machine learning and evolutionary techniques in interplanetary trajectory design

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 [28, 15, 22, 31, 48, 52, 35, 23], machine learning techniques [2, 4, 33, 27], evolutionary neuro-controllers [7, 8], tree search methods [17, 18, 50, 27, 23] 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 [47] 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 [43, 42, 44].

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) [28] 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 [40] and impulsive thrust strategies [15] were deployed. Studies on the use of some form of GA, also leveraging on the ever increasing computational power available, continued (see for example [41, 3, 14, 9, 12]) 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 [34, 10, 22]. 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 [52, 55, 35, 48, 32, 37]. 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 [23] 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) [32] and Ant Colony Optimisation (ACO) [39, 5]. More recently, Covariance Matrix Evolutionary Strategy (CMA-ES) was shown to be potentially outperforming self-adaptive DE on a class of transfers [24].

While the no free lunch theorem [54] 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 [53] (a sort of open trajectory gym) was created and made available to the scientific community and is still the subject of active research [51, 1, 45, 20, 4, 49]. 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 [12] 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 [31] followed by deeper studies on Non Dominated Sorting Genetic Algorithm (NSGA-II) that included also, as integer variables, the planetary fly-by sequence [9]. Multi-objective versions of PSO have also been considered early on [30]. The more modern Multi Objective Evolutionary Algorithm by Decomposition (MOEA/D) was later identified as a most performing technique in applicable cases [24], 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 [46, 5], 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 [26, 25, 13] or ACO [50] have been used in bi-level optimisation schemes, often coupled with smart tree search strategies such as Beam Search (BS) [17, 50], Monte Carlo Tree Search (MCTS) [18], Multi Objective Beam Search (MOBS) [27] or Lazy Race Tree Search (LRTS) [23] 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.

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 [4] on some of the trajectory problems in the GTOP database [53] 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 [2]. 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 [19, 33]. 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 [27], or to define new box bounds and hence focus successive evolutionary runs in promising areas of the search space (cluster pruning [21]).

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 [7] 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 [42, 43, 44]. 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.

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 vehicles landing accuracy [6, 11]. 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 [11] 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 [43, 42]. 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 [43]:

• 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. $\boldsymbol{\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. $\hat{\mathbf{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^{*}$).

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:

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_{\mathrm{sp}}g_{0}$. We have denoted with $I_{\mathrm{sp}}$ the propulsion specific impulse and with $g_{0}$ the Earth gravity constant at sea level. The control variables $u(t)$ and $\hat{\mathbf{i}}_{u}(t)$ (or $\hat{\mathbf{u}}(t)$ for brevity) describe the thrust level (here also called throttle) and its direction and are constrained as follows $\lvert u(t)\rvert\le1$ and $\lvert\hat{\mathbf{i}}_{u}(t)\rvert=1,\forall t\in[t_{1},t_{2}]$.

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

and able to steer the spacecraft from some initial point $\mathbf{x}_{1}\in\mathcal{S}_{1}$ to some final point $\mathbf{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 [38]).

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 $\boldsymbol{\lambda}_{r},\boldsymbol{\lambda}_{v}$ and $\lambda_{m}$, or, for brevity, $\boldsymbol{\lambda}$. We then introduce the Hamiltonian function:

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 $\hat{\mathbf{i}}_u$ we have that $\mathcal{H}=...+c_1\frac{\lambda_v}{m}\cdot\hat{\mathbf{i}}_u 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 $\boldsymbol{\lambda}_v$ for $\mathcal{H}$ to be minimised, more formally:

Substituting this expression back into the Hamiltonian we get:

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:

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^*$:

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,

as:

The costate equations

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

while for the first three equations it is trivial to get back the expressions in (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

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

We may thus write the corresponding Hamilton equation as:

The remaining two equations are then also easily obtained as:

and:

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

The two-point boundary value problem

Following the maximum principle, we know that there exist initial values for the co-states $\boldsymbol{\lambda}(t_1)=\boldsymbol{\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,\boldsymbol{\lambda}_1$, for a time $t_2-t_1$ leading to a final condition $\mathbf{x}_2,\boldsymbol{\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:

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.

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 $\boldsymbol{\theta}_1$ and $\boldsymbol{\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:

and equivalently,

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, $\boldsymbol{\theta}_2=[0,0,\theta_2]$ and the transversality condition for a free final mass is thus:

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 $\boldsymbol{\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:

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

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:

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

Similar expressions are derived in [36] using true anomaly and not the eccentric anomaly and using a different method.

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 $\boldsymbol{\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).

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:

Vector (PDF)

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.

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 1 which is run starting from 10 different points equally spaced along the nominal trajectory and using $\alpha= 0$. Let the final data-set

contain $I$ trajectories

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 $\boldsymbol{\lambda}$, and the times $t$. We refer to this data set as QOC (quadratic optimal control).

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}_{\mathrm{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 2 and results in a new data set $\mathcal{T}_{\mathrm{MOC}}$ containing mass optimal trajectories and visualised in Figure 2.

Algorithm 1 Random Walk

1. procedure $\operatorname{RANDOM\ WALK}(x_{\mathrm{nom},j},\lambda_{\mathrm{nom},j},\Delta t_{\mathrm{nom},j},E_{f_{\mathrm{nom}}},\alpha,\gamma,n)$

2. $\mathcal{T} \gets\{\}$ ▷ Set of optimal trajectories

3. $\mathbf{x}_0 \gets\mathbf{x}_{\mathrm{nom},j}$ ▷ Set initial state as nominal

4. $\boldsymbol{\lambda}_0 \gets\boldsymbol{\lambda}_{\mathrm{nom},j}$ ▷ Set initial co-state as nominal

5. $\Delta t \gets\Delta t_{\mathrm{nom},j}$ ▷ Set time of flight as nominal

6. $E_f \gets E_{f_{\mathrm{nom}}}$ ▷ Set arrival eccentric anomaly as nominal

7. $\gamma\gets\bar{\gamma}$ ▷ Nominal perturbation percentage

8. for $i \gets1,\ldots,n$ do ▷ $n$ random perturbations

9. $\boldsymbol{\beta} \gets U(-1,1) \in\mathbb{R}^{7}$ ▷ Vector of random uniformly distributed numbers

10. $\mathbf{x}_1 \gets\mathbf{x}_0+\mathbf{x}_0\boldsymbol{\beta}\gamma$ ▷ Perturb state in random direction by $\gamma$

11. $(\boldsymbol{\lambda}_1,\Delta t_1,E_{f_1}) \gets\operatorname{TPBVP}(\mathbf{x}_1,\boldsymbol{\lambda}_0,\Delta t,E_f,\alpha)$

12. if successful then ▷ If the TPBVP is successfully solved

13. Update $\mathcal{T}$ ▷ Save successful trajectory

14. $\mathbf{x}_0 \gets\mathbf{x}_1$ ▷ Accept perturbed state

15. $\boldsymbol{\lambda}_0 \gets\boldsymbol{\lambda}_1$ ▷ Accept new computed co-states

16. $\Delta t_0 \gets\Delta t_1$ ▷ Accept new computed time of flight

17. $E_{f_0} \gets E_{f_1}$ ▷ Accept new computed eccentric anomaly

18. $\gamma\gets(\gamma+\bar{\gamma})/2$ ▷ Increase perturbation size

19. else ▷ If solution is unfeasible

20. $\gamma\gets\gamma/2$ ▷ Decrease perturbation size and solve again

21. end if

22. end for

23. return $\mathcal{T}$ ▷ Return set of optimal trajectories

24. end procedure

Algorithm 2 Homotopy

1. procedure $\operatorname{HOMOTOPY}(\mathbf{x}_{\mathrm{QOC},0},\boldsymbol{\lambda}_{\mathrm{QOC},0},\Delta t_{\mathrm{QOC}},E_{f_{\mathrm{QOC}}},\alpha_{\mathrm{tol}})$

2. $\alpha\gets1$ ▷ Try to solve for mass-optimal at first

3. loop

4. $(\boldsymbol{\lambda},\Delta t,E_f) \gets\operatorname{TPBVP}(\mathbf{x}_{\mathrm{QOC},0},\boldsymbol{\lambda}_{\mathrm{QOC},0},\Delta t_{\mathrm{QOC}},E_{f_{\mathrm{QOC}}},\alpha)$

5. if successful then ▷ If the TPBVP is successfully solved

6. $\alpha^{\star} \gets\alpha$ ▷ Current best $\alpha$

7. $(\boldsymbol{\lambda}_{\mathrm{QOC},0},\Delta t_{\mathrm{QOC}},E_{f_{\mathrm{QOC}}}) \gets(\boldsymbol{\lambda},\Delta t,E_f)$ ▷ $\mathbf{x}_0^{\star}$ doesn’t change

8. if $\alpha<\alpha_{\mathrm{tol}}$ then

9. $\alpha\gets(1+\alpha)/2$ ▷ Increase $\alpha$

10. else if $\alpha\ge\alpha_{\mathrm{tol}}$ then ▷ If $\alpha$ is close to 1

11. $\alpha\gets1$

12. else if $\alpha=1$ then ▷ If trajectory is feasible and mass-optimal

13. return $\mathcal{T}$ ▷ Break the loop and return trajectory

14. end if

15. else

16. $\alpha\gets(\alpha+\alpha^{\star})/2$ ▷ Decrease $\alpha$ and retry

17. end if

18. end loop

19. end procedure

Learning details and network architecture

We consider for both of the data-sets, $\mathcal{J}_{\mathrm{MOC}}$ and $\mathcal{J}_{\mathrm{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:

where $\theta$ is the polar angle, $\phi$ is the azimuth angle and $u$ is the throttle magnitude. In previous work [42], 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 [16] 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 [29], 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.

Results

We train different networks to predict the three outputs both concurrently or separately and on both $\mathcal{J}_{\mathrm{MOC}}$ and $\mathcal{J}_{\mathrm{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{J}_{\mathrm{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.

Network	MSE	$u$	$\phi$	$\theta$	Epochs	Network	MSE	$u$	$\phi$	$\theta$	Epochs
$\mathcal{N}^{QOC}_{u,\phi,\theta}$		0.0178	0.0690	0.0870	453	$\mathcal{N}^{MOC}_{u,\phi,\theta}$		0.7250	0.0255	0.0337	794
$\mathcal{N}^{QOC}_{\phi,\theta}$		-	0.0755	0.0907	461	$\mathcal{N}^{MOC}_{\phi,\theta}$		0.0212	0.0289		316
$\mathcal{N}^{QOC}_{u}$		0.0075	-		247	$\mathcal{N}^{MOC}_{u}$		0.7740	-	-	721
$\mathcal{N}^{QOC}_{\phi}$		-	0.0698	-	394	$\mathcal{N}^{MOC}_{\phi}$		-	0.0204	-	248
$\mathcal{N}^{QOC}_{\theta}$		-	-	0.0976	503	$\mathcal{N}^{MOC}_{\theta}$		-	-	0.0305	324

Vector (PDF)

Vector (PDF)

Vector (PDF)

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 [42], by integrating forward the whole spacecraft dynamics considering the deep network predictions as a control and will be part of a separate future publication.