Deep Networks as Approximators of Optimal Transfers Solutions in Multitarget Missions

TL;DR

This paper demonstrates that deep neural networks can accurately estimate optimal transfer solutions in multitarget interplanetary missions, significantly speeding up preliminary mission design processes.

Contribution

The study introduces a deep learning approach to estimate optimal transfer parameters, outperforming traditional methods in speed and accuracy for complex space mission trajectory optimization.

Findings

Deep networks achieve less than 0.5% error for low-thrust transfer estimates.

Deep networks achieve less than 4% error for multi-impulse transfer estimates.

Deep learning models outperform other machine learning algorithms in this context.

Abstract

In the design of multitarget interplanetary missions, there are always many options available, making it often impractical to optimize in detail each transfer trajectory in a preliminary search phase. Fast and accurate estimation methods for optimal transfers are thus of great value. In this paper, deep feed-forward neural networks are employed to estimate solutions to three types of optimization problems: the transfer time of time-optimal low-thrust transfers, fuel consumption of fuel-optimal low-thrust transfers, and the total dv of minimum-dv J2-perturbed multi-impulse transfers. To generate the training data, low-thrust trajectories are optimized using the indirect method and J2-perturbed multi-impulse trajectories are optimized using J2 homotopy and particle swarm optimization. The hyper-parameters of our deep networks are searched by grid search, random search, and the tree-structured Parzen estimators approach. Results show that deep networks are capable of estimating the final mass or time of optimal transfers with extremely high accuracy; resulting into a mean relative error of less than 0.5% for low-thrust transfers and less than 4% for multi-impulse transfers. Our results are also compared with other off-the-shelf machine-learning algorithms and investigated with respect to their capability of predicting cases well outside of the training data.

Tables

a  =	semi-major axis
e  =	eccentricity
i  =	inclination
$\Omega$  =	right ascension of ascending node
$\omega$  =	perigee
$\mu$  =	gravitational constant of central body
$u$  =	engine thrust ratio
$\alpha$  =	unit vector of thrust direction
$g_0$  =	standard value of gravitational acceleration
Tmax  =	maximum thrust magnitude
Isp  =	thruster specific impulse
M  =	transformation matrix of modified equinoctial elements
D  =	gravity vector
$ϵ$  =	homotopy parameter
$m$  =	mass of spacecraft
$\lambda$  =	co-state vector
$L$  =	phase angle
nlayer  =	number of layers
nneuron  =	number of neurons
f  =	activation function
opt  =	optimizer
$\eta$  =	learning rate
$c$  =	decay rate of learning rate
B  =	batch size
F  =	feature

Table 2. Table 2: Search space of grid search

Hyper-parameters	Search space
nlayer	2, 3, 4, 5
nneuron	32, 128
$\eta$	0.01, 0.001
B	200, 800
F	Frv, Fcoe, Fmee

Table 3. Table 3: Search space of random search and TPE

Hyper-parameters	Search space
nlayer	2, 3, 4, 5
nneuron	32, 128
f	sigmoid, tanh, relu
opt	GD, MGD, AGD
$\eta$	0.1 - 0.0001
dm	ED, NED
ds	500 - 1500
$c$	0.8 - 1
B	100 - 1000
F	Frv, Fcoe, Fmee, Farv, Facoe, Famee

Table 4. Table 4: Basic information of low-thrust dataset

	ms (kg)	|$\Delta$a| (AU)	|$\Delta$e|	$|\Delta$i| (deg)	|$\Delta$L| (deg)	Optimal time (d)	ttransfer	Optimal fuel (kg)
mean	1490.2	0.1041	0.0262	1.31	14.85	533.76	719.60	270.15
max	2000	0.3	0.0971	3	30	1350.16	1791.60	873.76
min	1000	1e-6	1e-6	7e-5	8.4e-5	78.63	97.23	37.01

Table 5. Table 5: Basic information of multi-impulse dataset

	|$\Delta$a| (km)	|$\Delta$e|	$|\Delta$i| (deg)	|$\Delta \Omega$| (deg)	ttransfer	Optimal $\Delta$v (m/s)
mean	66.43	0.0055	0.76	1.48	7.18	119.51
max	277.88	0.0192	3.76	3.80	24.94	508.32
min	0.0079	3e-6	3.4e-4	5.3e-5	0.48	3.11

Table 6. Table 6: Best DNN models.

	nlayer	nneuron	f	opt	$\eta$	dm	ds	$\lambda$	B	F
TOP	4	218	sigmoid	AGD	1.99e-3	NED	700	0.85	250	Farv
FOP	5	170	sigmoid	AGD	7.06e-3	NED	700	0.85	1000	Faee
JMP	4	210	sigmoid	AGD	3.41e-3	NED	1000	0.83	850	Fcoe

Table 7. Table 7: Accuracy of test datasets.

	1	2	3
TOP	99.54%	99.43%	99.15%
FOP	99.64%	99.65%	99.66%
JMP	95.77%	96.06%	95.94%

Table 8. Table 8: Search space of other machine-learning algorithms.

Algorithms	Search space
RandomForest	n_estimators: [10, 100], max_depth: [2,100]&None, min_samples_split: [2, 10], min_samples_leaf: [1, 4], bootstrap: [True, False]
Bagging	n_estimators: [10, 100], max_samples: [0.2, 0.8], bootstrap: [True, False]
AdaBoost	n_estimators: [10, 100], learning_rate: [0.01, 2], loss: [linear, square]
ExtraTrees	n_estimators: [10, 100], criterion: [mse, mas], max_depth: [2,100]&None, min_samples_split: [2, 10], min_samples_leaf: [1, 4], bootstrap: [True, False]
GradientBoosting	n_estimators: [10, 100], loss: [ls, lad, huber], max_depth: [2,100], min_samples_split: [2, 10], min_samples_leaf: [1, 4],
DecisionTree	splitter: [best, random], max_depth: [2,100]&None, min_samples_split: [2, 10], min_samples_leaf: [1, 4],
ExtraTree	splitter: [best, random], max_depth: [2,100]&None, min_samples_split: [2, 10], min_samples_leaf: [1, 4],

Table 9. Table 9: Comparison with other machine-learning algorithms.

	TOP	FOP	JMP
RandomForest	94.75%	98.11%	95.10%
Bagging	94.94%	98.12%	95.18%
AdaBoost	86.37%	96.33%	82.92%
ExtraTrees	95.50%	98.29%	95.22%
GradientBoosting	95.15%	98.16%	95.29%
DecisionTree	92.15%	97.56%	93.75%
ExtraTree	91.98%	97.41%	93.68%
DNN	99.54%	99.66%	96.02%

Table 10. Table 10: Generalization capability of low-thrust problems.

	ms	|$\Delta$a|	|$\Delta$e|	|$\Delta$i|	|$\Delta$L|
Mean	2188.5 kg	0.3279 AU	0.1604	3.28 deg	32.96 deg
TOP	98.63%	99.29%	97.37%	99.49%	99.41%
FOP	97.79%	99.33%	97.42%	99.58%	99.49%

Table 11. Table 11: Generalization capability of multi-impulse problems.

	|$\Delta$a|	|$\Delta$e|	|$\Delta$i|	|$\Delta \Omega$|
Mean	300.68 km	0.0209	3.95 deg	4.17 deg
JMP	94.25%	92.10%	96.11%	96.54%

Table 12. Table 12: Results of time-optimal multitarget mission.

Actual rank	Rendezvous sequence	$t_{estimated}$ (d)	$t_{actual}$ (d)	Accuracy (%)	Estimated rank
1	6553 7400 803 12570 5701 11047 2053 3427 11005 693	1287.49	1321.36	97.44	1
2	6553 7400 803 12570 5701 11047 2053 3427 693 11005	1297.48	1322.16	98.13	9
3	7424 803 7400 12570 5701 11047 2053 3427 11005 693	1288.25	1323.62	97.33	2
4	7424 803 7400 12570 5701 11047 2053 3427 693 11005	1298.26	1325.12	97.97	11
5	7424 803 7400 11047 12570 5701 2053 3427 11005 693	1293.04	1330.66	97.17	6
6	6553 803 7400 12570 5701 11047 2053 3427 11005 693	1292.16	1330.75	97.10	4
7	6553 7400 803 11047 12570 5701 2053 3427 11005 693	1290.68	1331.59	96.93	3
8	6553 803 7400 11047 12570 5701 2053 3427 11005 693	1295.14	1332.15	97.22	8
9	7424 803 7400 11047 12570 5701 2053 3427 693 11005	1303.04	1334.50	97.64	17
10	6553 803 7400 12570 5701 11047 2053 3427 693 11005	1302.17	1334.61	97.57	14
11	7424 803 7400 11047 5701 12570 2053 3427 11005 693	1302.96	1335.33	97.58	16
12	6553 7400 803 11047 12570 5701 2053 3427 693 11005	1300.70	1335.70	97.38	13
13	6553 803 7400 11047 12570 5701 2053 3427 693 11005	1305.12	1336.42	97.66	20
14	6553 7400 803 11047 5701 12570 2053 3427 11005 693	1298.34	1336.58	97.14	12
15	803 7400 6553 12570 5701 11047 2053 3427 11005 693	1293.17	1337.27	96.70	7
16	803 7400 6553 12570 5701 11047 2053 3427 693 11005	1303.17	1342.70	97.06	18
17	7400 803 6553 12570 5701 11047 2053 3427 11005 693	1292.47	1368.52	94.44	5
18	3914 9641 10952 4932 3013 8433 7359 8484 7366 9616	1297.95	1370.37	94.72	10
19	7400 803 6553 12570 5701 11047 2053 3427 693 11005	1302.48	1371.74	94.95	15
20	6553 7400 803 12570 11047 5701 2053 3427 11005 693	1304.24	1378.53	94.61	19

Table 13. Table 13: Results of fuel-optimal multitarget mission.

Actual rank	Rendezvous sequence	$m_{estimated}$ (kg)	$m_{actual}$ (kg)	Accuracy (%)	Estimated rank
1	1150 14205 7410 12512 12461 4984 14195 5033 3018 3911	1412.47	1376.76	97.41	1
2	1150 14205 12512 7410 12461 4984 14195 5033 3018 3911	1401.73	1362.17	97.10	2
3	14205 1150 7410 12512 12461 4984 14195 5033 3018 717	1382.29	1360.13	98.37	8
4	1150 14205 7410 12512 12461 4984 14195 3018 5033 717	1382.90	1359.52	98.28	7
5	1150 14205 7410 12512 4984 12461 14195 5033 3018 717	1393.23	1359.36	97.51	4
6	1150 14205 7410 12512 4984 12461 14195 5033 3018 3911	1393.51	1354.31	97.11	3
7	14205 1150 7410 12512 12461 4984 14195 5033 3018 3911	1382.93	1354.13	97.87	6
8	1150 14205 7410 12512 12461 4984 14195 5033 717 7362	1381.12	1351.52	97.81	9
9	1150 14205 7410 12512 12461 4984 7393 7362 717 5033	1377.29	1350.86	98.04	13
10	1150 14205 7410 12512 12461 4984 7362 7393 3911 5033	1374.63	1350.83	98.24	16
11	14205 1150 12512 7410 12461 4984 14195 5033 3018 717	1374.08	1349.14	98.15	18
12	6534 14205 12512 7410 12461 4984 14195 5033 3018 717	1383.67	1348.34	97.38	5
13	1150 14205 7410 12512 12461 4984 14195 5033 3018 7362	1377.42	1347.33	97.77	12
14	1150 14205 7410 12512 4984 12461 14195 5033 717 3018	1378.52	1343.63	97.40	10
15	1150 14205 7410 12512 12461 4984 14195 3018 3911 5033	1375.53	1342.19	97.52	14
16	1150 14205 7410 12512 12461 4984 7362 7393 717 5033	1373.63	1340.00	97.49	19
17	6534 14205 7410 12512 4984 12461 14195 5033 3018 717	1373.53	1337.48	97.30	20
18	1150 14205 7410 12512 12461 4984 14195 5033 3018 2645	1377.83	1337.13	96.96	11
19	1150 14205 12512 7410 12461 4984 14195 3018 5033 3911	1374.71	1336.06	97.11	15
20	6534 14205 7410 12512 4984 12461 14195 5033 3018 3911	1374.27	1332.29	96.85	17

Table 14. Table 14: FOP.

Spacecraft	Removal sequence	$\Delta {\text{<em>v</em>}}_{DNN}$ (m/s)	$\Delta {\text{<em>v</em>}}_{actual}$ (m/s)	Accuracy (%)
1	31 1 5 78 24 23 38 11 32 16	768.00	798.80	96.14
2	52 8 21 18 62 41 9 108 115 59	965.29	965.39	98.86
3	4 63 37 26 109 79 20 61 15 112	1288.36	1285.29	99.76
4	100 29 120 86 55 14 67 46 53 94	1038.73	1014.49	97.61
5	106 85 93 73 39 36 45 50 95 98	1072.01	1064.22	99.27
6	10 91 99 65 22 70 71 33 7 13	1436.10	1434.62	99.90
7	68 102 77 34 54 118 3 113 47 81	1194.01	1195.71	99.86
8	111 58 114 88 27 76 110 101 80 105	719.16	730.01	98.51
9	119 42 0 60 40 43 12 66 83 6	996.17	1012.42	98.40
10	30 17 57 92 56 69 97 2 72 87	908.91	925.45	98.21

Table 15. Table 15: Computation time.

	timeDNN (ms)	timeactual (ms)
TOP	6.1	2061
FOP	3.5	17337
JMP	3.2	399782

Equations

$$\begin{array}[]{l}p=a\left({1-{e^{2}}}\right),\\ {e_{x}}=e\cos\left({\Omega+\omega}\right),\\ {e_{y}}=e\sin\left({\Omega+\omega}\right),\\ {h_{x}}=\tan\left({{i\mathord{\left/{\vphantom{i2}}\right.\kern-1.2pt}2}}\right),\cos\Omega\\ {h_{y}}=\tan\left({{i\mathord{\left/{\vphantom{i2}}\right.\kern-1.2pt}2}}\right),\sin\Omega\\ L=\Omega+\omega+f,\end{array}$$

$$\begin{array}[]{l}\dot{\textbf{\emph{x}}}={\textbf{\emph{M}}}\dfrac{{{T_{\max}}u}}{m}{\textbf{\emph{$\alpha$}}}+{\textbf{\emph{D}}},\\ \dot{m}=-\dfrac{{{T_{\max}}}}{{{I_{sp}}{g_{0}}}}u,\end{array}$$

$$\begin{array}[]{l}J={\lambda_{0}}\int_{{t_{0}}}^{{t_{f}}}{Ldt},\\ {L_{{\rm{TOP}}}}=1,\\ {L_{{\rm{FOP}}}}=\dfrac{{{T_{max}}}}{{{I_{sp}}{g_{0}}}}\left({u-\varepsilon\ln\left[{u\left({1-u}\right)}\right]}\right),\end{array}$$

$$H={\text{λ}}_{\text{x}}^T\text{M}\frac{T_{max}u}{m}\text{α}+{\text{λ}}_{\text{x}}^T\text{D}-{\lambda }_m\frac{T_{\max }u}{I_{sp}{g}_0}+{\lambda }_{0}L.$$

$$\text{α}=-\frac{{\text{M}}^T{\text{λ}}_{\text{x}}}{\Vert {\text{M}}^T{\text{λ}}_{\text{x}}\Vert }.$$

$$\begin{array}[]{l}{u_{{\rm{TOP}}}}=\left\{{\begin{array}[]{*{20}{l}}{0,{\rm{if}}{\rho_{{\rm{TOP}}}}>0},\\ {1,{\rm{if}}{\rho_{{\rm{TOP}}}}<0},\\ {[0,1],{\rm{if}}{\rho_{{\rm{TOP}}}}=0},\end{array}}\right.\\ {\rho_{{\rm{TOP}}}}=-\dfrac{{{I_{{\rm{sp}}}}{g_{0}}\left\|{{\textbf{\emph{M}}^{T}}\textbf{\emph{$\lambda$}}_{\textbf{\emph{x}}}}\right\|}}{{{\lambda_{0}}m}}-\dfrac{{{\lambda_{m}}}}{{{\lambda_{0}}}},\end{array}$$

$$\begin{array}[]{l}{u_{{\rm{FOP}}}}=\dfrac{{2\varepsilon}}{{{\rho_{{\rm{FOP}}}}+2\varepsilon+\sqrt{\rho_{{\rm{FOP}}}^{2}+4{\varepsilon^{2}}}}},\\ {\rho_{{\rm{FOP}}}}=1-\dfrac{{{I_{{\rm{sp}}}}{g_{0}}\left\|{{\textbf{\emph{M}}^{T}}\textbf{\emph{$\lambda$}}_{\textbf{\emph{x}}}}\right\|}}{{{\lambda_{0}}m}}-\dfrac{{{\lambda_{m}}}}{{{\lambda_{0}}}}.\end{array}$$

$$\begin{array}[]{c}\textbf{\emph{x}}({t_{0}})=\textbf{\emph{x}}{{}_{0}},m({t_{0}})={m_{0}},\\ \textbf{\emph{x}}({t_{f}})=\textbf{\emph{x}}{{}_{f}}.\end{array}$$

$$\text{Φ}(\text{λ}(t_0){)}_{\mathrm{TOP}}=[\text{x}(t_f)-{\text{x}}_f,{\lambda }_m(t_f),||\text{λ}(t_0)||-1,H(t_f)-{\text{λ}}_x(t_f)\cdot {\dot{\text{x}}}_f{]}^T=0,$$

$$\text{Φ}(\text{λ}(t_0){)}_{\mathrm{FOP}}=[\text{x}(t_f)-{\text{x}}_f,{\lambda }_m(t_f),||\text{λ}(t_0)||-1{]}^T=0.$$

$$\begin{array}[]{l}\dfrac{{{\rm{d}}{a_{m}}}}{{{\rm{d}}t}}=0,\dfrac{{{\rm{d}}{e_{m}}}}{{{\rm{d}}t}}=0,\dfrac{{{\rm{d}}{i_{m}}}}{{{\rm{d}}t}}=0,\\ \dfrac{{{\rm{d}}{\Omega_{m}}}}{{{\rm{d}}t}}=-\dfrac{3}{2}{J_{2}}{\left({\dfrac{{{R_{e}}}}{p}}\right)^{2}}n\cos\left({{i_{m}}}\right),\\ \dfrac{{{\rm{d}}{\omega_{m}}}}{{{\rm{d}}t}}=\dfrac{3}{4}{J_{2}}{\left({\dfrac{{{R_{e}}}}{p}}\right)^{2}}n\left({5{{\cos}^{2}}\left({{i_{m}}}\right)-1}\right),\\ \dfrac{{{\rm{d}}{M_{m}}}}{{{\rm{d}}t}}=n+\dfrac{3}{4}{J_{2}}{\left({\dfrac{{{R_{e}}}}{p}}\right)^{2}}n\sqrt{1-e_{m}^{2}}\left({3{{\cos}^{2}}\left({{i_{m}}}\right)-1},\right)\end{array}$$

$$\begin{array}[]{l}\Delta{\textbf{\emph{v}}_{\rm{i}}}=\Delta{v_{\rm{i}}}\left[\begin{array}[]{c}\cos\left({{\alpha_{\rm{i}}}}\right)\cos\left({{\beta_{\rm{i}}}}\right)\\ \cos\left({{\alpha_{\rm{i}}}}\right)\sin\left({{\beta_{\rm{i}}}}\right)\\ \sin\left({{\alpha_{\rm{i}}}}\right)\end{array}\right],\\ \Delta{v_{\rm{i}}}=\Delta{v_{\max}}{x_{3{\rm{i}}+1}},\\ {\alpha_{\rm{i}}}=\pi\left({{x_{3{\rm{i}}+2}}-0.5}\right),\\ {\beta_{\rm{i}}}=2\pi{x_{3{\rm{i}}+3}}\,\ i=0,1,\ldots,{\rm{N}}-2,\end{array}$$

$$\begin{array}[]{l}{t_{0}}={t_{\rm{s}}}+{x_{3\left({{\rm{N}}-1}\right)+1}}\left({{t_{\rm{N}}}-{t_{\rm{s}}}}\right),\\ {t_{\rm{i}}}={t_{{\rm{i}}-1}}+{x_{3\left({{\rm{N}}-1}\right)+1+{\rm{i}}}}\left({{t_{\rm{N}}}-{t_{{\rm{i}}-1}}}\right),\\ i=1,2,\ldots,{\rm{N}}-1.\end{array}$$

$$\Delta v=\left(\sqrt{{\left(\frac{\Delta a}{2a_1}\right)}^2+{\left(\frac{\Delta e}{2}\right)}^2+{\left(\Delta i\right)}^2}+\Delta \Omega \sin i_1\right)\sqrt{\frac{2\mu }{a_1}}.$$

$$\text{li+1}=\text{f}(\text{w}{\text{l}}_i+\text{b}),$$

$$f_{sig}(x)=\frac{1}{1+e^{-x}},$$

$$f_{tanh}(x)=\frac{e^x-e^{-x}}{e^x+e^{-x}},$$

$$f_{relu}(x)=\max (0,x).$$

$$E=\frac{1}{n}\sum _{i=1}^n{\left(\hat{y}-y\right)}^2,$$

$$\begin{array}[]{l}{\rm{error}}=\dfrac{1}{n}\sum\limits_{i=1}^{n}{\dfrac{{\left|{\hat{y}-y}\right|}}{y}},\\ {\rm{accuracy}}=1-{\rm{error}}.\end{array}$$

$${\text{w}}^{\prime }=\text{w}-\eta \frac{\partial E}{\partial \text{w}},$$

$${\eta }^{\prime }=\eta c^t,$$

$${\eta }^{\prime }=\eta e^{-ct},$$

$$t=\frac{global\_step}{decay\_step}.$$

$$\begin{array}[]{l}\textbf{\emph{h}}=\textbf{\emph{r}}\times\textbf{\emph{v}},\\ E=-\dfrac{\mu}{{2a}}.\end{array}$$

Full text

Nomenclature

1 Introduction

Multitarget interplanetary missions are a complex domain for optimization. In space missions to multiple targets, substantial benefits are gained since the average expense of exploring each target is lowered. Multitarget missions that are already being conducted are focused mostly on small celestial bodies in our solar system (except multi-gravity-assist missions). Small celestial bodies attract vast research interest because of their great significance in various aspects [1, 2, 3, 4]. NASA’s Near Earth Asteroid Rendezvous (NEAR) Shoemaker performed a flyby of the asteroid 253 Mathilde on the way to the asteroid 433 Eros [5]. In the previous plan for the NEAR mission, an ambitious plan called Small-Body Grand Tour, which aimed to achieve flybys of two comets and two asteroids over a 10-year period, was proposed [6]. Deep Space 1 [7] and Dawn [8] also achieved multitarget visits of small celestial bodies. The Global Trajectory Optimization Competitions (GTOCs), which are one of the most challenging events in space engineering, had several editions greatly focussed on multitarget small celestial body missions555https://sophia.estec.esa.int/gtoc_portal/. In addition to multitarget missions to small celestial bodies, multitarget active debris removal (ADR) missions have also been avidly studied recently. The population of low-earth-orbit (LEO) debris has grown rapidly in the last two decades [9] and will continue to increase even if no new spacecraft is launched, as a consequence of the collision between existing debris, known as Kessler’s syndrome [10]. Many ADR methods are being studied, including contactless interaction [11] and the debris engine666Lan, L., Li, J., and Baoyin, H., “Debris Engine: A Potential Thruster for Space Debris Removal,” 2015, https://www.technologyreview.com/s/544156/junk-eating-rocket-engine-could-clear-space-debris/.. Optimization of multitarget ADR missions is a complex mathematical problem that can be seen as a complex variant to the more famous Travelling Salesman Problem [12] and was the subject of the 9th GTOC [13] and 8th China Trajectory Optimization Competition, of which the best solutions are presented in [14, 15].

The preliminary phases in the design of multitarget space missions often consider these problems as combinatorial-continuous [16]. The combinatorial part is a sequence optimization problem in which many mission options are available, and the continuous part is a single-leg transfer trajectory optimization problem where two types of transfer models, low-thrust and multi-impulse, are mainly considered. The low-thrust trajectories can be optimized using indirect methods [17, 18, 19, 20, 21]. The multi-impulse trajectories can be optimized using particle swarm optimization (PSO) [22, 23] and genetic algorithms [24]. Although these proposed methods are efficient, they are time-consuming when a large number of transfers are required to be optimized. Therefore, it is impractical to optimize each trajectory in the preliminary design of multitarget missions, and fast estimation methods are thus useful to approximate the optimal solutions of transfers between two targets as highlighted in [25].

The famous Edelbaum’s approximation gives the formulas of characteristic velocity requirements for multi-impulse and low-thrust transfers between quasi-circular orbits [26]. In [15], the formulas are extended to J2-perturbed low-thrust problems. Edelbaum’s method is completely analytical, and therefore very suitable for a large-scale sequence search problem. However, for low-thrust optimization problems, these approximation formulas are derived under the assumption that the transfer is a many-revolution transfer. Many-revolution low-thrust transfers are impractical in deep-space multitarget missions because even one revolution will already take years (for example, 4 or 5 years in the main belt). Thus, transfers between two targets in multitarget deep space missions are usually less than one revolution. When the transfer time is short, Edelbaum’s methods become inaccurate, and estimation methods for short transfers are required. An approximate method is presented to estimate the velocity change of low-thrust transfers between co-planar, low-eccentricity orbits [27], and is extended to transfers with small inclination change [28]. These methods are based on Edelbaum’s approximation and are accurate for the evaluation of multitarget missions. In these methods, equations need to be iteratively solved and an additional correction is necessary for rendezvous transfers, making these methods somewhat complex to implement. In addition to the methods based on Edelbaum’s approximation, machine-learning methods [25, 29] are proposed to estimate the maximum final mass of the fuel-optimal rendezvous problem. In their research, after the features are well selected, different regressors are tested, and their performances are compared. Machine-learning methods are found to be vastly superior to other commonly used approximation methods, such as Lambert estimation. Estimation methods of multi-impulse optimization have rarely been studied in the past, but it is also time-consuming to optimize many-revolution J2-perturbed multi-impulse trajectories [23]. Estimation of multi-impulse optimization will gradually become necessary with the increment of studies in, for instance, multitarget ADR missions in low Earth orbit.

Deep neural networks (DNNs) have attracted the interest of many researchers and have had applications in various fields in recent years. DNNs consist of connected layers that are composed of neurons and are able to exhibit desired behavior after learning [30]. DNNs have their applications in the field of astrodynamics as well, including optimal control problems [31, 32], autonomous lunar landing [33], and orbit prediction [34]. In [31], DNNs are trained to learn the optimal control in four different cases of pinpoint landing. Their results, which show that the landings controlled by DNNs are close to optimal ones, allow for the design of on-board real-time optimal control. In [32], DNNs are trained to design the Earth-Mars orbital transfers of mass optimal control. In [33], DNNs, i.e., convolutional neural networks and recurrent neural networks, are trained to perform an autonomous Moon landing using only raw images taken by on-board cameras, with the optimal control directly determined by the images only. In [34], DNNs are trained by historical orbit determination and prediction data to improve orbit prediction accuracy. It can be seen from the above research that DNNs can establish good connections between input and output, which is very difficult to obtain in an analytical way. The estimation of optimal transfers considered in this paper is also a problem for which a connection between the transfer parameters and the optimal solution must be established. Therefore, DNNs are employed in this paper to estimate optimal solutions of transfers.

Datasets that are necessary for the training of DNNs are generated by optimizing transfer trajectories. In this paper, three different cases, time optimal low-thrust, fuel optimal low-thrust, and minimum $\Delta$v J2-perturbed multi-impulse trajectories, are optimized. For the low-thrust optimization, the orbital states are expressed in modified equinoctial elements (MEEs) because they are non-singular and are thus efficient and robust [15]. Co-state normalization [17] is utilized to improve the efficiency of the initial guess of co-states. The logarithmic homotopic method [17] is applied to overcome the difficulties caused by bang-bang control in fuel optimal problems. For J2-perturbed multi-impulse optimization, a many-revolution J2-Lambert problem is solved using J2 homotopy at first. In J2 homotopy, the constant J2 is multiplied by a homotopic parameter $\epsilon$ that varies from 0 to 1 iteratively to approach the J2-perturbed problem from a two-body-problem perspective. PSO is then applied to optimize the total $\Delta$v. Then, 100,000 trajectories are generated for each case, and are divided into a training dataset, validation dataset, and test dataset. The training dataset is used for DNNs to learn. The validation dataset is used for model selection. There are many hyper-parameters in DNN models, such as learning rate and batch size. Hyper-parameters are set prior to learning and not selected by the learning algorithm itself. A hyper-parameter search is necessary because the performance of a DNN is highly affected by hyper-parameters. In this paper, three techniques, grid search [35], random search [36], and the tree-structured Parzen estimators (TPE) approach [37], are applied to search the optimal DNN model. The well-trained DNNs are capable of estimating optimal transfers with extremely high accuracy. DNNs also show excellent performance when compared with other machine-learning algorithms and investigated based on generalization capability. Finally, DNNs are applied to three mult-itarget missions and prove to be highly fast and accurate.

The rest of this paper is organized as follows. In Section II, trajectory optimization methods are presented. In Section III, the datasets are generated using a main-belt asteroid rendezvous mission and multitarget ADR mission as backgrounds. In Section IV, the architecture of DNNs and hyper-parameter search techniques are detailed. In Section V, DNNs are analyzed and applied to multitarget missions. Conclusions are drawn in Section VI.

2 Trajectory Optimization

2.1 Low-Thrust Trajectory Optimization

Modified equinoctial elements (MEEs) are used here as the states x because they are non-singular and thus are efficient and robust. Knowing the classical orbital elements (COEs), the MEEs can be expressed as:

[TABLE]

where a is the semi-major axis, e the eccentricity, i the inclination, $\Omega$ the right ascension of the ascending node, $\omega$ the perigee, and f the true anomaly.

The dynamic equations of low-thrust propelled spacecraft in a two-body model are:

[TABLE]

where x = [p, ex, ey, hx, hy, L], m is the instantaneous mass of the spacecraft, $\mu$ the gravitational constant of the central body, u (0 $\leq$ u $\leq$ 1) the engine thrust ratio, $\alpha$ the unit vector of the thrust direction, g0 = 9.80665 m/s2 the standard value of gravitational acceleration, Tmax the maximum achievable thrust magnitude, and Isp the thruster specific impulse. M is a transformation matrix and D the gravity vector, the details of which can be found in [38].

For time-optimal problems (TOPs) and fuel-optimal problems (FOPs), the performance indexes take the forms:

[TABLE]

where $\lambda_{0}$ is a positive factor that does not change the origin optimal problem and is necessary for the normalization of the co-state variables [17]. For FOPs, the performance index uses the logarithmic homotopic method [39] in which the fuel optimal index is achieved when $\epsilon$ becomes 0 from 1. Introducing the co-state vector $\lambda$(t) = [\lambda$${}_{\textbf{\emph{x}}}, $\lambda$m, $\lambda$0], the Hamiltonian is built as

[TABLE]

According to Pontryagin’s minimum principle, the direction of the optimal thrust is determined as follows:

[TABLE]

The optimal thrust magnitude of TOPs is:

[TABLE]

and the optimal thrust magnitude of FOPs is:

[TABLE]

Given the boundary condition, the optimal control problem can be transformed into a two-point boundary-value problem (TPBVP) and can be solved by shooting methods solving shooting equations. For rendezvous problems, in which the final time is optimized in TOPs and fixed in FOPs, the state variables of the spacecraft satisfy the following boundary conditions:

[TABLE]

The shooting functions of TOPs and FOPs can be deduced according to optimal control theory [20, 15] and are given as follows:

[TABLE]

[TABLE]

Compared with FOPs, TOPs are relatively easy to solve due to the continuity and differentiability of the right-hand sides of the dynamical equations and Euler-Lagrange equations. A number of ideas have been advanced, also recently [17, 19], to improve the FOP convergence.

2.2 J2-perturbed Multi-Impulse Trajectory Optimization

For a higher computing efficiency of J2-perturbed propagation, the dynamic equations are set up using mean orbit elements (MOEs):

[TABLE]

where J2=1.08262668e-3, $\emph{p}=a_{m}(1-e^{2})$,$n=\sqrt[]{\mu/a_{m}^{3}}$, $\mu=3.9860044e14m^{3}/s^{2}$ is the gravitational constant of Earth, and the Earth radius Re=6378137.0 m. The transformation between MOEs and COEs can be found in [40], and then the position r and velocity v can be acquired. In the J2-perturbed multi-impulse problems (JMPs) considered in this paper, the initial and final states, as well as the initial and final epochs ts and tN, are fixed. The entire trajectory of the spacecraft can be split into N+1 legs by N intermediate impulses. The position and velocity of the spacecraft at the epoch ti are denoted ri and vi±, respectively, where the superscripts – and + represent the velocity before and after the impulse, respectively. Given the impulse vectors $\Delta$vi and corresponding epochs ti, i=0,1,2,…,N-2, the first N legs can be propagated sequentially with Eq. (11). The last leg is obtained by solving a J2-Lambert problem using the shooting method to ensure that the rendezvous conditions are met.

J2 homotopy is employed here to overcome the difficulties in the convergence of the shooting process caused by the long transfer duration. The constant J2 is multiplied by a homotopic parameter $\epsilon\in$ [0,1] that varies from 0 to 1 iteratively to approach the J2-perturbed problem from the perspective of a two-body problem. We denote the target position and velocity at the final epoch as rf and vf, respectively. For the last leg, the shooting uses vN-1+ as the shooting variable and $\Delta$r=rN- rf as the shooting function. A two-body Lambert problem in which the spacecraft transfers from rN-1 to rN for a flight time of tN-tN-1 is solved first to obtain the initial guess of vN-1+, and then the homotopy process is implemented. The pseudo-algorithm is given in Algorithm 1. After the J2-Lambert problem is solved, the last two impulses as well as the total $\Delta$v can be calculated.

The optimization of $\Delta$vtotal, which is determined by $\Delta$vi and ti, is a typical parameter optimization problem, and swarm intelligence algorithms can find optimal solutions efficiently [22, 23]. Thus, a PSO is applied to optimize $\Delta$vtotal of the multi-impulse rendezvous problems in this paper. Considering a rendezvous process in which N+1 (N$\geq$2) impulses are employed, the unknown variables contain $\Delta$vi, i=0,1,2,…,N-2, and ti, i=0,1,2,…,N-1, so the dimension of this case is n=(N-1)×3+N=4N-3. All the n unknown variables are denoted by the normalized position vector of a particle x=[x1, x2,…,xn], of which each component is in [0, 1]:

[TABLE]

[TABLE]

When N=1, which means only two impulses are employed, the dimension of x is n=1 and the only unknown variable is t0. The number of impulses (or N) is also unknown and it may differ for optimal solutions with different transfer conditions. To determine the optimal number of impulses, an enumeration for N is applied in this algorithm. Using PSO to obtain the optimal $\Delta$vtotal for the cases where N=1, 2, 3, and 4 separately, the case with the smallest $\Delta$vtotal is considered the optimal solution. Such an enumeration is implemented by parallel computing to improve the computational efficiency; that is, the four cases are computed concurrently through four threads.

3 Dataset Generation

3.1 Dataset Generation of Low-Thrust Trajectories

A main-belt asteroid rendezvous mission is considered as the test case for the low-thrust optimization problem. Asteroid orbital parameters are taken from the GTOC7 competition [41] description. Altogether, 100,000 time-optimal trajectories and 100,000 fuel-optimal trajectories are generated. To make the generation process fast and highly efficient, asteroids with eccentricities larger than 0.1 are excluded. The process of low-thrust database generation is given in Algorithm 2.

A pre-screening procedure is adopted to select asteroid pairs. The a**MAX is 0.3 AU, i**MAX is 3 deg, and L**MAX is 30 deg. This procedure ensures that the two selected asteroids are close enough to each other at the given epoch to improve the convergence rate of low-thrust optimal control problems. The start time is randomly given from [56800, 59800] MJD, and the initial mass is randomly given from [1000, 2000] kg. If the TOP converges, the FOP is then solved with the transfer time transfer randomly given from [1.2, 1.5]×t**top. If also the FOP converges, then the following information will be saved to the database:

A1, A2, start time, initial mass, transfer time of FOP, optimal time of TOP, optimal fuel of FOP.

There are many local optima in complex optimization problems. In the indirect method of solving low-thrust optimization problems, different initial guesses of co-states can lead to different local optima. Therefore, each TOP and FOP is solved 10 times and the best result is saved as the optimal value.

3.2 Dataset Generation of Multi-Impulse Trajectories

A multitarget ADR mission is considered as the test case for the multi-impulse optimization problem. Debris orbital parameters are from the GTOC9 competition [13]. Altogether, 100,000 J2-perturbed multi-impulse trajectories are generated. The process of multi-impulse database generation is given in Algorithm 3.

A $\Delta$v estimation calculated by Eq. (24) [23] is used as pre-screening. If this $\Delta$v is larger than $\Delta$v**max= 500 m/s, than this debris pair will be abandoned:

[TABLE]

The start time is randomly given from [56,800, 59,800] MJD, and the transfer time is randomly given from [5, 30] days. The following information will be saved to the database:

A1, A2, start time, transfer time, optimal $\Delta$v.

The PSO will run three times for each problem and the best result is saved as the optimal value.

4 Neural Architecture

4.1 Architecture of DNNs

The structure of the network considered in this paper is illustrated in Fig.1 and is essentially a feed-forward neural network with multiple hidden layers. The structure is determined by the number of layers n**layer and the number of neurons n**neuron at each hidden layer. The data will flow through these layers and finally produce the output. At each layer, denoting the input li, the output li+1 is calculated as follows:

[TABLE]

where w is the weight matrix, b the bias vector, and f a nonlinear function named the activation function. This mechanism is also shown in Fig.1. There are three most commonly used activation functions for the hidden layers: the sigmoid function expressed in Eq.(16), the hyperbolic tangent (tanh) function expressed in Eq.(17), and the rectified linear (ReLu) function expressed in Eq.(18). For the output layer, a linear function is selected:

[TABLE]

[TABLE]

[TABLE]

The training process consists of adjusting the value of the parameters of each layer to minimize the loss function:

[TABLE]

where y is the actual value, $\hat{y}$ the estimated value, and n the number of data in this training iteration. To evaluate the performance of DNNs, the error function is given as:

[TABLE]

Gradient descent (GD) algorithms are the state-of-art in training the parameters, e.g.,

[TABLE]

where $\eta$ is the learning rate. Some modified gradient descent algorithms are also very effective, such as momentum gradient descent (MGD) [42] and Adam gradient descent (AGD) [43].

A most important hyper-parameter is the learning rate $\eta$. With a small learning rate, the DNN may take a very long time to converge, while with a large learning rate the DNN may oscillate around the optimal solution. Therefore, it is a good strategy to start from a large learning rate and decrease the learning rate gradually. Two decay schedules are applied in this paper: one is exponential decay (ED),

[TABLE]

and the other is natural exponential decay (NED),

[TABLE]

where c is the decay rate and

[TABLE]

In training, the batch size B is also a hyper-parameter that must be chosen. A batch of data is put into the DNN to train it and then the weights are updated. This hyper-parameter effects training time more [35]. Early stopping is utilized here, and the training will stop when there is no improvement in the last N epochs.

The form of input, which is also known as features, is also an important factor that affects the performance of DNNs. For the problem considered in this paper, the input is a set of transfer parameters. The orbit states can be expressed in many different forms [44]. Three forms are considered in this paper: Cartesian form rv (x, y, z, vx, vy, vz), classical orbit elements COEs (a, e, i, $\Omega$, $\omega$, f), and modified equinoctial elements MEEs (p, ex, ey, hx, hy, L), and the corresponding input forms are donated Frv, Fcoe, and Fmee, respectively.

In addition to the basic transfer parameters, other information can also be added into the input. The differences of two integration constants, angular momentum h and mechanical energy E in Eq. (25), between initial and final orbits, can be used as additional features:

[TABLE]

For the fuel-optimal, low-thrust problem, a two-body Lambert solution is used as an additional feature. For the J2-perturbed multi-impulse problem, Eq.(14) is used as an additional feature. The input forms with additional features are donated Fa**rv, Fa**coe, and Fa**mee. All of the inputs are scaled to be between -1 and 1. The actual transfer time (day) of TOPs is multiplied by a factor of 0.0172 as the actual output, the actual final mass (kg) of FOPs is multiplied by a factor of 0.004 as the actual output, and the actual $\Delta\emph{v}_{total}$ (m/s) of JMPs is multiplied by a factor of 0.01 as the actual output.

4.2 Selection of DNN models

The performance of a DNN is highly affected by hyper-parameters, which are set prior to the learning and not selected by the learning algorithm itself. Many model choices are implemented by manual search, which gives the impression of neural network training as an art [35]. There are two drawbacks of manual search. One is that it is not easily reproducible and automated, and the other is that a human may never find the optimal hyper-parameters. Therefore, three techniques, grid search [35], random search [36], and the tree-structured Parzen estimators (TPE) approach [37], are applied in this paper to select the most proper DNN model.

Grid Search: Grid search is an automated search process executed after defining a search space. However, an exponential explosion may occur when the number of hyper-parameters increases. One must therefore define the search space carefully due to the limitation of computing capability.

Random Search: Random search is highly recommended by [36], as they found random sampling very efficient. Random search allows one to explore more hyper-parameters and more values of one hyper-parameter, compared to grid search. The enlarged search space increase the possibility of finding the optimal hyper-parameters, while random sampling does not increase the computational budget. Random search is also simple to implement.

Tree-structured Parzen Estimators Approach: TPE is a sequential model-based optimization (SMBO) algorithm. In grid search and random search, the next step is completely independent of previous steps. SMBO algorithms will consider previous results and decide what to try next. TPE uses random forest regression to provide the prediction of which region of hyper-parameters is optimal. The details of TPE implementation can be found in [37].

The hyper-parameters that must be searched are number of layers nlayer, number of neurons at each hidden layer nneuron, activation function f, batch size B, optimizer opt, initial learning rate $\eta$, decay model dm, decay step ds, and decay rate c. In addition to the hyper-parameters, six feature forms, Frv, Fcoe, Fmee, Fa**rv, Fa**coe, and Fa**mee, also must be searched. The search space of grid search is listed in Table 2 and that of random search and TPE in Table 3. There are altogether 96 trails in grid search; therefore, 100 trails are set in random search and TPE, so the three search methods are comparable. Even though their computing budgets are the same, random search and TPE can explore a larger parameter space. These parameters are uniformly random, except that the learning rate is uniform in the log-domain [40]. The other default hyper-parameters of grid search are f = relu, opt = GD, dm = ED, ds = 1000, and c = 0.9. The step for early stopping is set to 50.

The dataset is divided into 80,000 training datasets, 10,000 validation datasets, and 10,000 test datasets. The training dataset is used to train the network, which is inside the learning algorithm. The validation dataset is used for DNN model selection, which is outside the learning algorithm. The test dataset is finally used to evaluate the performance of selected model. Many take validation and test datasets as one dataset, but when a model-selection process is considered the validation and test datasets must be separated.

5 Simulations

5.1 Dataset Analysis

The basic dataset information is listed in Tables 4 and 5. In Table 4, the phase angle L is calculated at the start time, and t**transfer is the transfer time of FOPs. The maximum value of orbital differences in the low-thrust dataset is the same as the pre-screening settings. In Table 5, the $\Delta\Omega$ is the orbital difference between start debris and target debris at the rendezvous time. Note that the minimum transfer time in Table 5 is less than 1 d, which is the lower boundary when generating the dataset, because the transfer time is the duration between rendezvous and first impulse time, which is an optimization variable.

The distribution of datasets is shown in Fig. 2, where the upper, middle, and lower rows are TOPs, FOPs, and JMPs, respectively. The color bar indicates the amount. In the last sub-figure of JMPs, $\Delta\Omega^{\prime}$ is the difference of $\Omega$ precession rates between start debris and target debris. For low-thrust transfers, it can be found that the phase angle has an apparent relevance to the optimal value. Transfers are concentrated in the center (small orbital differences) of the $\Delta\emph{a}$ and $\Delta\emph{i}$ sub-figures, and in the left (light in mass) of ms sub-figure because it is easy to converge for these transfers. For J2-perturbed multi-impulse transfers, a strong relevance can be found between $\Delta\emph{i}$ and the optimal value because the outer-plane maneuver costs $\Delta\emph{v}$ significantly. $\Delta\Omega^{\prime}$ shows a similar relevance because $\Delta\Omega^{\prime}$ is strongly determined by i, which can be seen in Eq. (11). The $\Delta\Omega$ shows a nearly opposite tendency, although it is the outer-plane change. This is because the precession rate $\Omega^{\prime}$caused by J2 perturbation can be used to adjust the orbital plane instead of using large $\Delta\emph{v}$. The transfer with large $\Delta\Omega$ and small $\Delta\emph{i}$ can pass the pre-screening described in Algorithm 3, and the actual $\Delta\emph{v}$ is small because of the small $\Delta\emph{i}$.

5.2 DNN Results

Three hyper-parameter search techniques, grid search, random search, and TPE, are implemented to select the optimal DNN model. The results of TOPs, FOPs, and JMPs are shown in Figs. 3, 4, and 5, respectively. In these figures, the x axis is the trails and the y axis the error on the validation dataset. It should be noted that the DNN model is trained on the training dataset and selected according to its performance on the validation dataset. Dots with three different colors and shapes represent the three search techniques. The lower sub-figure is the zoom of the bottom of the upper sub-figure. For TOPs, the best DNN model is found by random search, and the corresponding error is 0.61%. For FOPs, the best DNN model is found by TPE, and the corresponding error is 0.35%. For JMPs, the best DNN model is found by random search, and the corresponding error is 4.01%, while the error of the best model found by TPE is 4.06%. Random search and TPE have similar performance, but both are superior to grid search. Random search is recommended because it is highly effective and easy to implement.

Overall, the errors of JMPs are larger than those of TOPs and FOPs. This is because the optimality of low-trust datasets is better than that of the multi-impulse dataset. Among the methods used to optimize low-thrust trajectories, the indirect method utilized in this paper has the advantage of optimality. Although there are many local optima, the approach of optimizing the trajectory 10 times and choosing the best one can make the solution close to the global optimal solution. PSO as an stochastic optimization algorithm does not have an optimality as good as that of the indirect method. Therefore, the sub-optimal solutions in the multitarget dataset make the network difficult to learn.

Some errors are larger than 10% because their learning parameters are not set properly. However, in Fig. 5, errors are larger than 60% in some trails. By investigating their hyper-parameters, it is found that in most of these trails the features are Frv or Fa**rv. This is because, in long-duration multi-impulse transfers, the phase has very little effect on the optimal $\Delta\emph{v}$, but the phase does determine the position and velocity. Even if the spacecraft starts the transfer from two different position along the same orbit, which means the position and velocity are quite different, the optimal $\Delta\emph{v}$ can be nearly the same. Therefore, Frv and Fa**rv are not suitable features for long-duration transfers.

The best DNN models selected by the search algorithms are listed in Table 6. Some general guidance in choosing DNN models can be gleaned from this table, i.e., they are all quite deep and large networks that have several of the same or similar parameters. The activation functions of the three problems are sigmoid, the optimizers are AGD, the initial learning rates are all 1e-3, and the decay models are NED. It should be noted that some of these hyper-parameters may never be selected by manual and grid search.

The performance of the selected DNN models must be checked on the test dataset. Since the initialization of weights and bias are random, each complete training process can lead to different results. Therefore, the DNNs are trained three times and that with the best performance on the test dataset will be saved for the following analysis. The training histories of the best DNNs of the three problems are shown in Fig. 6. Up until now, well-trained DNNs have been created. In the following, accuracy will be used to present the performance of the DNNs. The accuracy on test datasets are listed in Table 7. For low-thrust problems, the accuracy of the DNNs is higher than 99.5%, and for multi-impulse problems it is higher than 96%. These results demonstrate that a well-trained DNN is capable of estimating optimal transfers with extremely high accuracy.

DNNs are compared with other machine-learning algorithms, e.g., those used in [25, 29]. A random search is utilized to tune these algorithms and the search space of the hyper-parameters is shown in Table 8 [45]. Features are the same as those in the best DNN models. The algorithms are trained on training datasets and their accuracy is tested on test datasets. For each case, 100 trails are implemented and the results are listed in Table 9. It is shown that DNNs are superior to all of the other machine-learning algorithms. However, these machine-learning algorithms also have excellent performance. In [25, 29], estimated optimal transfers were already proved to have been improved significantly over those of the Lambert model, which is most commonly used. Our work confirms the effectiveness of machine learning, especially deep learning, in estimating optimal transfers.

To further investigate the performance of well-trained DNNs, generalization capability is studied by using the transfers in which the transfer parameters are outside the boundary of which values are given in Tables 4 and 5. For low-thrust problems, these transfers are generated by a similar procedure described in Algorithm 2, except that the pre-screening is different. In the pre-screening used to study generalization capability, the transfer will be saved if it is 1-1.2 times the boundary value. For instance, to study the generalization capability in terms of |$\Delta$a|, the transfer will be saved only when |$\Delta$a|max < |$\Delta$a| < 1.2|$\Delta$a|max. To study the generalization capability in terms of m**s, the initial mass is given from [2,000, 2,400] kg randomly. When studying |$\Delta$e|, asteroids with eccentricities larger than 0.1 are not excluded as before; otherwise, it will be very difficult for |$\Delta$e| to be in the range 1-1.2 |$\Delta$e|max. The initial mass m**s, |$\Delta$a|, |$\Delta$e|, |$\Delta$i|, and |$\Delta$L| are studied separately, and for each situation 100 transfers are calculated. The generalization capability results of low-thrust problems are given in Table 10. The accuracy is still higher than 99% in terms of |$\Delta$a|, |$\Delta$i|, and |$\Delta$L|. In |$\Delta$e|, the accuracy is 97%, because including asteroids with eccentricities larger than 0.1 also brings in the change of the original problems.

For multi-impulse problems, |$\Delta$a|, |$\Delta$e|, |$\Delta$i|, and |$\Delta\Omega$| are studied separately. For the study of |$\Delta$i| and |$\Delta\Omega$|, transfers are generated by a similar procedure as described in Algorithm 3. In addition to |$\Delta$|max < |$\Delta$| < 1.2|$\Delta$|max being necessary, $\Delta\emph{v}_{max}$ is changed to 700 m/s. Since |$\Delta$a|max and |$\Delta$e|max almost reach the limitation of all the 123 debris, the approach to increasing |$\Delta$a| or |$\Delta$e| is to select one transfer from a dataset randomly and change the COE of the target orbit to make |$\Delta$| in the range 1-1.2 times |$\Delta$|max. The generalization capability results of multi-impulse problems are given in Table 11, and demonstrate that well-trained DNNs have an excellent generalization capability.

5.3 Applications in Multitarget Mission Design

Here, the well-trained DNNs are applied in the design of three multitarget missions. The first two multitarget missions are multitarget main-belt asteroid, low-thrust rendezvous missions. The asteroids’ orbital parameters are from the 7th GTOC [41]. The third multitarget mission is a multitarget ADR J2-perturbed multi-impulse mission. The debris orbital parameters are from the 9th GTOC [13].

In the first mission, 10 asteroids must be visited in the minimum time. The modified Julian date (MJD) of departure is 59,215. After excluding asteroids with eccentricities larger than 0.1, 4,986 asteroids remain. The amount of all possible solutions is 4,986!/4,976! (approximately 1037). Transfers between the two targets are time-optimal, low-thrust transfers, and the optimal transfer time and mass consumption are estimated by the well-trained DNN. Beam search, which is commonly used for this type of problem [15, 46], is applied to search the optimal sequence. The same pre-screening procedure is adopted as described in Algorithm 2. The beam width is set to be 10,000, and the heuristic cost is set to be transfer time. One can tune the options of beam search, but it is not necessary in the present case.

After the searching process is over, the best 20 sequences are optimized using the indirect method. Results are shown in Table 12, where the actual solutions are solved by the indirect method and estimated solutions are acquired by the DNN. Most of the accuracies are higher than 97%, which is quite a high accuracy considering that the length of the sequence is 10. The forward sequences in the actual rank are still forward in the estimated rank. The rank difference is due to the fact that their transfer days are very close to each other.

In the second mission, 10 asteroids must be visited using minimum fuel consumption. Transfers between the two targets are fuel-optimal, low-thrust transfers. The minimum transfer time $t_{min}$ is first estimated by the DNN and the transfer time is given by $t_{fop}=c*t_{min}$. In this simulation, c is set to be 1.35. The fuel consumption is then estimated by the DNN. Beam search is also applied. Other simulation parameters are the same as in the first mission. After the searching process is over, the best 20 sequences are optimized using the indirect method. Results are shown in Table 13. All the accuracies are higher than 96%, and the actual and estimated ranks have a high consistency.

In the third mission, 100 pieces among the total 123 debris of debris must be removed with 10 spacecraft. Each spacecraft will remove 10 debris using minimum $\Delta$v**total. For each spacecraft, the start and transfer times are randomly given, and a pre-screening procedure, described in Algorithm 3, is also utilized. The single leg $\Delta$v is estimated by the DNN. After 10 sequences are obtained, PSO is applied to optimize the actual $\Delta$v**total. Results are shown in Table 14. All the accuracies are higher than 96%. It is interesting that even if the accuracy of the well-trained DNN of JMPs is not as good as that of TOPs and FOPs, an excellent performance in the sequence accuracy is still evident. Results of the three multitarget missions demonstrate that DNNs can provide highly accurate solutions in preliminary mission design.

Another important factor in multitarget mission design is computation time, which is shown in Table 15. The computation time is single-leg computation time calculated by averaging the computation times when optimizing the entire sequence. The actual computation time includes the procedure that solves the problem 10 (for TOPs and FOPs) or three (for JMPs) times and chooses the best solution. There is no doubt that DNNs have an obvious advantage in computation speed.

6 Conclusions

In this paper, three types of optimal transfers, time-optimal low-thrust transfers, fuel-optimal low-thrust transfers, and minimum-$\Delta$v J2-perturbed multi-impulse transfers, are estimated using DNNs. Random search is recommended to search hyper-parameters of DNNs because it is highly effective and easy to implement. After the DNNs are well-trained, the accuracy can be higher than 99.5% for low-thrust trajectories and higher than 96% for multi-impulse trajectories. DNNs are also compared with other machine-learning algorithms that have already been proved to improve accuracy in the estimation of optimal transfers, and are found to be superior to all of them. The investigations of generalization capability show that well-trained DNNs can still estimate optimal transfers accurately even if transfer parameters are outside the area in which the datasets are generated. Finally, DNNs are applied in three multitarget mission designs and are proved to be highly fast and accurate. This work demonstrates the estimation of optimal transfers via a DNN trained on a regression task, and further applications in multitarget mission design can be sought.

Funding Sources

This research is supported by the Chinese National Natural Science Fund for Distinguished Young Scientists of China (No. 11525208).