{'fraction_non_alphanumeric': 0.02717731285188017, 'fraction_numerical': 0.00755311508074419, 'mean_word_length': 4.464927829830337, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
MGA trajectory planning with an ACO-inspired algorithm

Given a set of celestial bodies, the problem of finding an optimal sequence of swing-bys, deep space manoeuvres (DSM) and transfer arcs connecting the elements of the set is combinatorial in nature. The number of possible paths grows exponentially with the number of celestial bodies. Therefore, the design of an optimal multiple gravity assist (MGA) trajectory is a NP-hard mixed combinatorial-continuous problem. Its automated solution would greatly improve the design of future space missions, allowing the assessment of a large number of alternative mission options in a short time. This work proposes to formulate the complete automated design of a multiple gravity assist trajectory as an autonomous planning and scheduling problem. The resulting scheduled plan will provide the optimal planetary sequence and a good estimation of the set of associated optimal trajectories. The trajectory model consists of a sequence of celestial bodies connected by twodimensional transfer arcs containing one DSM. For each transfer arc, the position of the planet and the spacecraft, at the time of arrival, are matched by varying the pericentre of the preceding swing-by, or the magnitude of the launch excess velocity, for the first arc. For each departure date, this model generates a full tree of possible transfers from the departure to the destination planet. Each leaf of the tree represents a planetary encounter and a possible way to reach that planet. An algorithm inspired by Ant Colony Optimization (ACO) is devised to explore the space of possible plans. The ants explore the tree from departure to destination adding one node at the time: every time an ant is at a node, a probability function is used to select a feasible direction. This approach to automatic trajectory planning is applied to the design of optimal transfers to Saturn and among the Galilean moons of Jupiter.


INTRODUCTION
The complete automatic design of multiple gravity assist trajectories (MGA), that is the definition of an optimal sequence of planetary encounters and the definition of one or more locally optimal trajectories for each sequence, has been approached with several different techniques. All of them can be classified in two main categories: two level approaches and integrated approaches.
Two-level approaches split the problem into two sub-problems which lay at two different levels: one sub-problem is to find the optimal sequence of planetary encounters; the other is to find an optimal trajectory for that sequence. Two-level approaches define the planetary sequence independently of the trajectory itself. Once the sequence (or a set of promising sequences) has been selected, then the optimal trajectory can be searched for within the set of selected sequences only [1].
Simplified, low fidelity, models for representing the trajectory [2] are used at the first level: this allows for a quick assessment of many sequences, if not all. At the second level, a full model is used to optimize the trajectory. Each sequence is represented by a string of integer numbers, while the associated trajectory is represented with a string of real and integer numbers defining the time and the characteristics of the events occurring along the trajectory (e.g. launch, deep space manoeuvre, arrival at a celestial body, number of revolutions around the Sun, etc.). Therefore, for each sequence, there is an infinite variety of possible trajectories. The issue with two-level approaches is that it is difficult to assess the optimality of a given planetary sequence without an exhaustive search for all possible trajectories associated with that sequence. Unfortunately, finding an optimal trajectory is a very difficult global optimisation problem in itself. This, combined with the fact that usually there exist a very high number of sequences for a given transfer problem, requires a considerable computational effort. The computational cost can be reduced by discarding non-promising sequences. However, if the low-fidelity model is not accurate enough, either some good sequences are discarded, or many of the retained ones can result to be actually not interesting.
As opposed to the two-level approaches, integrated approaches define a mixed integer-continuous optimization problem, which tackles both the search of the sequence and the optimization of the trajectory by using a single model, at the same time [3]. This kind of problem is known in literature as a hybrid optimization problem [4]. The main difficulty with integrated approaches is that a variation of even a single celestial body in the sequence corresponds to a substantially different set of trajectories. In addition, a variation of the length of the sequence implies varying the number of legs of the trajectory, and thus the total length of the solution vector.
The automatic design of a trajectory with discrete events was recently formulated as a hybrid optimal control problem [5], and a solution was proposed by Conway et al. [6] with a two level approach based on genetic algorithms.
Here we propose to formulate the complete automated design of a multiple gravity assist trajectory as an autonomous planning and scheduling problem. The resulting scheduled plan will provide the planetary sequence for a multiple gravity assist trajectory and a good estimation of the optimality of the associated trajectories.
Although the proposed method can fall in the category of the integrated approaches, the scheduling and the planning of the events are separated at two different levels. A specific MGA trajectory model was developed to automatically schedule the events, if a plan is available, and to provide a good estimation of the feasibility and quality of a trajectory. A novel algorithm, partially inspired by the Ant Colony Optimization (ACO) paradigm [7], was devised to explore the space of possible plans. ACO was originally created to solve the Travelling Salesman Problem [8], and later successfully applied to a number of other discrete optimisation problems. Here the original idea behind ACO was elaborated to solve the planning problem associated to the design of MGA trajectories. The basic ACO paradigm is hybridized with a taboo strategy based on an external archive on infeasible solutions. In the literature, some ACO-derived meta-heuristics exist for the specific solution of different scheduling problems. In particular, Merkle et al. [9] proposed to apply ACO to the solution of the Resource-Constrained Project Scheduling Problem, while Blum, in his work [10], suggested the hybridization of Ant Colony Optimization with a probabilistic version of Beam Search for the solution of the Open Shop Scheduling problem.
In this paper, at first we will present the trajectory model and the integrated scheduling of the events, then the novel ACObased algorithm and how the plan is constructed. Finally, two case studies will demonstrate the effectiveness of the proposed approach at solving known space trajectory design problems.

TRAJECTORY MODEL
The trajectory model is based on a two-dimensional linked-conic approximation of the trajectory and of the orbits of the planets. The trajectory is composed of a sequence of planar conic arcs linked together through discrete, instantaneous events. In particular, the sequence is continuous in position and piecewise continuous in velocity, i.e. each event introduces a discontinuity in the velocity of the spacecraft but not in its position. The discrete events can be: launch, deep space manoeuvre (DSM), swing-by, and brake.
A final assumption of the present implementation is that all the orbits of both spacecraft and celestial bodies are direct, thus no retrograde orbits are allowed.
In summary, the proposed trajectory model is composed of: a launch from the departure celestial body; a series of deep space flight legs connected through gravity assist manoeuvres (modelled through a linked-conic approximation); an arrival at a target celestial body. Each one of these basic components will be explained in the following sections.

Launch
The launch event is modelled as an instantaneous change of the velocity of the spacecraft with respect to the departure planet. The velocity change is given in terms of modulus 0 v (which depends on the capabilities of the launcher) and inplane direction, specified through the angle 0  , measured counter clockwise with respect to the planet's orbital velocity vector P v at time of launch 0 t . 0 t and 0  are free parameters of the model, while launch velocity modulus 0 v will be used to target the next planetary encounter and solve the phasing problem, as explained later.

Swing-by
Gravity assist manoeuvres, or swing-bys, are modelled as instantaneous changes of the velocity vector of the spacecraft due solely to the gravity field of the planet.
Given the relative velocity vector   v before the swing-by, the physical properties of unperturbed hyperbolic orbital motion   , which means that the modulus of the outgoing velocity v   at infinity is known. Its direction can be computed considering the anomaly of the outgoing asymptote (see Fig. 2): Here, P  is the gravity constant of the planet, and p r is the radius of the pericentre of the hyperbola. Once  is computed, the relative outgoing velocity is calculated by rotating   v in the plane of an angle  . As for the launch velocity magnitude, the radius of pericentre ps r is tuned to meet the terminal conditions of the transfer leg following the swing-by.

Deep space flight leg
Each deep space flight leg starts at a departure planet i P and ends at an arrival planet 1 i P  , and is made of two conic arcs linked at a point i M . If the leg contains a deep space manoeuvre, this is applied in this point, and it produces an instantaneous change in the heliocentric velocity vector of the spacecraft, due to an ignition of the engines. In this model, we assume that the DSM is performed either at the apocentre or pericentre of the conic arc preceding the manoeuvre. In addition, the change in velocity is tangential to that arc.
For clarity, in the remainder of this section, we neglect the subscript index i of the leg in all the variables.

First arc
Let us assume that the spacecraft is at a given planet 1 P at time 1 t . Its position coincides with that of the planet, which is known from the 2D ephemeris. The heliocentric velocity of the spacecraft, instead, depends on the preceding launch or swing-by event.
If the transfer leg contains a DSM, the first step is to find the position M and time DSM t of the deep space manoeuvre. The position can either be the pericentre or the apocentre, according to a binary parameter / p a f . The true anomaly of the DSM is given by The time of the DSM DSM t is found by using the Kepler's time law [11]. Note that the use of a tangential DSM at apsidal points is one of the simplifying assumptions of the trajectory model, which allows reducing the number of free parameters. This type of DSM is particularly effective for perigee/apogee raising/decreasing, or period change.
If the leg does not contain any DSM, i.e. 0 DSM m  , the first arc is propagated up to a fictitious point M, defined by adding an angle   to the initial true anomaly of the spacecraft. The reason for using this forced propagation is twofold: first, we want to force the spacecraft to move away from the planet, before computing the second leg. This is done to prevent that, if no full revolutions are considered, the first intersection could happen after a null time. Second, it prevents any event (swingby or DSM) to happen immediately after the swing-by, which would be unfeasible due to spacecraft operation constraints.
Therefore, the quantity   has to be larger than the machine numerical precision but small enough to allow for the modelling of short transfer legs. It is important to underline the value of   is not a design parameter of the model and does

Second arc
The second arc starts at point M and is propagated until the intersection of the orbit of planet 2 P . Given the orbital parameters of the spacecraft at M, and the orbital parameters of planet 2 P , the task is to find the intersection between the two coplanar orbits. If there are no intersections, the leg is unfeasible, and the initial conditions of the leg, or its parameters, have to be changed. Otherwise, one of the two possible intersections is selected according to the binary parameter 1/ 2 f : let us call , int   the true anomalies of the selected intersection, respectively along the orbit of the spacecraft and of the planet.  Figure 3 illustrates a complete leg, including a DSM. The figure highlights that the orbital intersection does not imply, in general, that the planet is at the intersection point at the correct time. This issue will be addressed in the following paragraph.

Solution of the Phasing Problem
In order to perform a gravity assist manoeuvre or a planetary capture, the terminal position of the spacecraft has to match that of the planet. However, at intersection time int t , planet 2 P is at true anomaly 2 P  , which is generally different from  .
The time of intersection is a function of the states at the beginning of the leg, which ultimately depend on 0 v or ps r depending on the starting event. Therefore, if we introduce the symbol  , such that

Complete trajectory
A complete trajectory is made of a sequence of transfer legs connecting 1 P n  celestial bodies 0 1 , , . The trajectory starts from 0 P at time 0 t with a launch event characterized by a departure angle 0  . The first leg goes from 0 P to 1 P , then a number of swing-bys and interplanetary legs follow, until the final planet P n P . Thus, a complete trajectory with 1 P n  planets has P n legs, and 1 P n  swing-bys.
The solution of Eq.
(2) provides a complete scheduling of the trajectory given the initial time 0 t and the five parameters Since these five parameters fully characterize all possible legs from a planet i P to a planet 1 i P  , they are said to define a type of transfer. Conversely, because of the multiplicity of the zeros of Eq. (2), each type of transfer corresponds to a set of trajectories.
Hence, assigning a value to 0 t , 0 which is a tree structure in which every branch, from root to leaves, is a trajectory. Each trajectory is characterised by a different combination of * 0 v and * p r for each leg.
The entire tree is a complete set of trajectories from 0 P to P n P and represents a solution of the MGA trajectory planning problem. Thus, a plan is fully defined by assigning a value to the parameters in Table 1  The whole model was implemented in ANSI C and compiled as a MEX-file for interfacing with MATLAB.

THE ACO-MGA ALGORITHM
The model described in the previous section yields a set of scheduled trajectories provided that a complete or partial plan is available. In this section, we present an optimization procedure, based on the ant colony optimization paradigm, to explore the space of possible plans.

Solution coding
In ACO-MGA, a solution is coded through a string of discrete values assigned to the parameters. However, the set of parameters discussed before is inhomogeneous, as it is made of real, integer and binary variables. In particular, 0 t , 0  and DSM i m are real continuous variables and need to be properly discretised. In the present implementation of ACO-MGA, the values of the departure date 0 t and the departure angle 0  , are assumed to be pre-assigned and therefore the two parameters are removed from the list of the variables. The rationale behind this choice is that, although the launch date has a great impact on the resulting trajectory, if an algorithm exists that is able to efficiently generate a complete plan for a given launch date, a systematic search can be performed along the launch window, with a given time step (as will be shown in one of the test cases). The angle 0  , on the other hand, is often dictated by the launch capabilities and the first planetary encounter [1]. DSM m is a design parameter in the optimisation: the discretisation is chosen a priori and it can be, for example, uniform within a range. The maximum value (in magnitude) depends on the engine capabilities and on the transfer problem itself. It is a task of the user to choose a significant set of values for this parameter.
Using the additional assumptions on 0 t , 0  , and fixing 0 P , each solution can be coded using a vector s of positive integers.
The vector has 2 legs n components. Each pair of consecutive components encodes all the parameters necessary to characterise one leg of the solution (Fig. 6). The first element of the pair is encoding the identification number of the target planet for the corresponding leg according to the following procedure: an ordered list , P i q containing all the planets available as a target for each leg i is predefined (and fixed); then, let q q q q rows and 5 columns (  is the cardinality of a set). The parameters for the j th type of transfer for the i th leg can be obtained as follows:  ). Each taboo list consists of a matrix, which has an arbitrary number of rows and 2i columns.
Dual to the list of taboo partial solutions, the feasible list feas L stores all the solutions, which are completely feasible, i.e. reach P n P . This is once more a matrix with an arbitrary number of rows and 2 P n columns.
Since each solution contained in the feasible list is complete, then it is possible to associate an objective value (or cost) to each one of them. In the following test cases, we will use, as objective value, the v  and a combination of v  and T. Note that, since, in general, there is more than one trajectory for a given solution (i.e. for a given set of free design variables), the objective value of a solution is given by the best trajectory value.

Search engine
The search space is organised as an acyclic oriented tree. Each branch of the tree represents a leg of the problem, while each node (or leave) represents a different destination planet and type of transfer. A population of virtual ants are dispatched to explore the tree, searching for an optimal solution.
The search runs for a given number of iterations

Solution generation
The tree is simultaneously explored, from root to leaves, by m ants. At each iteration, each one of the m ants explores the tree independently of the others, but taking into account the information collected by all the ants, at previous iterations, and associated to the j th element of , P i q is increased as follows: The amount of pheromone which is deposited depends on the objective value This mechanism is analogous to the pheromone deposition of standard ACO and aims at driving the search of the ants toward planetary sequences, which previously led to good solutions. In fact, those planets which generate (partial) sequences that appear either frequently in the feasible list, or rarely, but with low objective function are selected with higher probability. On the other hand, the probability of selecting other planets remains positive, such that one or more ants can probabilistically choose a planet that generates an undiscovered sequence. Note that, if the feasible list is empty, then all the planets have the same probability to be selected.
The parameter planet w controls the learning rate of the ants. A low value of planet w will make the term

5:
For each solution l in feas L that matches temp s 6:

10:
End For

Type of transfer generation
Once an ant has filled in the odd components of a solution s, it proceeds assigning values to the even components (see Algorithm 4). Similarly to the planet sequence generation, for each transfer all the available types of transfer are assigned, one at the time, to the solution s. A vector s for which a value is assigned to both the odd and even components up to leg i represents a partial solution. Similarly to before, a vector , t i τ contains the pheromone associated to the rows of the matrix i G (i.e. to each type of transfer).
For each new partial solution, the taboo list is first checked. If the partial solution appears in the taboo list, then it means that this solution will be infeasible, regardless of the way it is completed. The pheromone of the type of transfer associated to that sequence is set to zero, to avoid future selection of that type of transfer. If the partial solution does not appear in the taboo list, the feasible list is searched for any matching partial solution. For every match found, the pheromone of that type of transfer is modified as follows: where the weight type w is introduced with analogous meaning to planet w . In fact, the higher the coefficient, the higher the chances that solutions similar to the feasible ones are generated. Conversely, a low value of type w will favour the selection of sequences with a different type of transfer, thus increasing the random exploration of the whole solution space. The selection of the type of transfer happens probabilistically, in the same way as for the planet before.
If, at a given leg i, all possible transfer types correspond to partial solutions in the taboo list, the vector of pheromone distribution , t i τ will be full of zeros. As a consequence, the solution s (which can be partial or complete) is discarded, and the ant can stop its exploration of that branch of the tree. At the end of the solution generation step, the solution s is either discarded or completed. Once all the ants have completed their exploration, the result is a number of solutions (less than or equal to the number of ants m) to be evaluated.

Solution evaluation
Once a set of solutions S has been generated by the ants, each solution has to be evaluated to assess its feasibility and its objective value. This is done by calling the trajectory model.

Differences with standard ACO algorithm
Although the proposed approach takes its inspiration from ACO, and substantially applies the same principle, there is a fundamental difference that prevents the use of the standard ACO for the MGA problem. In the MGA problem, as opposed to the TSP tackled with ACO, the pheromone cannot be assigned to single legs of the trajectory: this is due to the fact that each leg (identified by its couple of integers) has no intrinsic value in the trajectory, if disconnected from the previous legs.
In other words, the only two parameters (planet and type of transfer) are not sufficient to fully characterise one leg: the actual value of the leg depends also on its initial conditions, which are in turn determined by the parameters of the previous legs. As a comparison, in the TSP, the distance between pair of cities, that is summed up to build the tour cost, is constant and fixed a-priori, for a given instance of the problem. For this reason, it is not possible to assign pheromone to single legs, independently from the rest of the solution.

CASE STUDIES
The proposed optimisation method was applied to two case studies inspired by the Laplace [12] and Cassini [13] missions, the former of which currently under preliminary study by ESA, NASA and JAXA.
ACO-MGA was tested against genetic algorithms, which are known to perform well on problems with discrete variables. In particular, it was chosen to use the genetic algorithm implemented in MATLAB ® within the Genetic Algorithm and Direct Search Toolbox™ (GATBX), and the Non-dominated Sorting Genetic Algorithm (NSGA-II) [14]. Settings for all the optimisers will be specified for each test case. While NSGA-II can deal with discrete variables, GATBX only uses real variables: a wrapper of the objective function was coded to round the continuous solution vector to the closest integer.
Due to the stochastic nature of the methods involved in the comparative tests, all the algorithms were run 100 times. The performance index used to compare the ACO-MGA against the other global optimisers is the success rate for a given number of function evaluations: according to the theory developed in [15], 100 repetitions give an error in the determination for the exact success rate of less than 6%, with 92% confidence.
Some preliminary tests showed that the best performances of ACO-MGA are achieved if the algorithm is run in 2 steps, using different sets of parameters. In particular, in the first step, the weights , planet type w w are set to null: remembering Eq.
(3) this choice translates into an initial pure random search. In the second step, weights are set to non-null values, to explore around the feasible solutions found.
where obj f  is the expected value for the objective function. In this way, choosing for example 1 w  , a 1 is added to the pheromone of a given element every time a matching solution with objective obj f  appears in the feasible list. The value of pheromone is higher if the objective value of the matching feasible solution is lower.
This two-step procedure can be explained in the following way. The first step allows a random sampling of the solution space, with the aim of finding a good number of feasible solutions. This is done to prevent the algorithm to stagnate around the first feasible solution found. The second step intensifies the search around the feasible solutions which were found at step one. Because of Eqs. (3) and (4), feasible solutions with low objective value are likely to be investigated further. In addition, the random component in the process does not forbid to keep exploring the rest of the search space.
The test cases were run on an Intel ® Pentium ® 4 3 GHz machine running Microsoft ® Windows ® XP.

Laplace case study
In this mission, the spacecraft reaches the sphere of influence of Jupiter after an interplanetary flight, and exploits a swingby of Ganymede to get captured into the Jovian system. At this point, multiple swing-bys of Ganymede and Callisto are used to reduce the relative velocity to Callisto v  , in order to be captured by the moon and start the scientific phase of the mission.
The problem under consideration relates to the second part of the transfer: it is assumed that the interplanetary trajectory has been already optimised, including the first Ganymede gravity assist. The resulting orbit has a 3:1 resonance Three free legs are added to the trajectory. For the first two, the algorithm can choose to target Ganymede or Callisto, while for the third and last, the target must be Callisto:  The radius of pericentre of the swing-bys is bounded between 1 and 3 radii P R of the body. Here a comment is needed: it is usual practice, in preliminary mission design, to consider a safety margin in the closest approach of a planet during a swingby manoeuvre. Usually this translates in considering a radius of pericentre not smaller than 1.1 P R . The main reason for using a slightly lower minimum altitude is that the number of feasible solutions in the solution space increases, since each swing-by can provide a higher deflection of the velocity vector, and thus the search for the optimal solution results more favourable for all the optimisers, but in particular for the population-based ones. We can also assume that the safety margin can be added when the solution is re-optimised with a more complete model, and assume that a DSM can compensate for the lack of swing-by deflection. Note that this trick will not be adopted in the next, more realistic case study. GATBX and NSGA-II were run for 600 function evaluations with the settings shown in Table 2. The mutation and crossover parameters pcross_bin and pmut_bin were tuned through several preliminary runs although very little difference in the quality of the results was registered.
Since the size of the population, popsize, is very important for genetic-based algorithms, and can significantly affect the results, this case study was also run 100 times with different sizes of the population (maintaining the predefined number of total function evaluations by varying the number of generations accordingly ngen). For NSGA-II, it resulted that there was no noticeable change in the quality of the results over 100 runs. This is related to the fact that NSGA-II is not completely converging with only 600 function evaluations. For GATBX, instead, results were changing significantly, and the settings leading to the best solutions were used.
The parameters for ACO-MGA were tuned by running the same test case for different values of the weights, the number of ants and the number of iterations. The maximum number of function evaluations was set to 600 as for the other algorithms although ACO-MGA requires, on average, 545 function evaluations. The best results were obtained with the following settings. 10 ants were used, with a first optimisation step with 30 iterations and , 0 planet type w w  , followed by a second step of 30 more iterations with , 2 0 3 k m / s planet type w w   . Because of the normalisation shown in Eq. (5), the non-dimensional weight w does not depend on the dimensions and value of the objective function of the specific problem: therefore, we will use the same value also for other transfer problems, as will be shown in the next case study. Note that the two weights , planet type w w are tuning parameters of the optimiser. It was found that the performance of ACO-MGA is not very sensitive to the value of the weight (within reasonable variations). Furthermore, it should be easy, for a given problem, to get a value for obj f  , either by experience or simply by random sampling the solution space. In the case obj f  cannot be determined, then the values of , planet type w w shall be selected by tuning the optimiser on the specific problem.
Results in the form of statistical parameters over the 100 runs are presented in Table 3. All the algorithms found at least one feasible solution in every run. The value of 2 km/s as a target value for the v  has been chosen to compute the success rate according to the procedure proposed in [16].
The results in Table 3 point out that, while all the algorithms find feasible solutions in all the runs, the quality of the solution is much better for ACO-MGA. Moreover, GATBX found a good solution only in 31% of the runs, and NSGA-II in 39%.
ACO-MGA, instead, found a good solution in 62% of the runs.
The time for one ACO-MGA run is about 8 min. The simplicity of the test case, together with the implementation of ACO-MGA in a high-level language like MATLAB, makes the use of an optimisation method slower than the systematic scan of the whole solution space. Note that this will not happen in the more complex Cassini test case.  The reference solution for this problem, as chosen by ESA during a preliminary study [17], was re-optimised using a full 3D model with one free deep space manoeuvre per leg [18], and minimising the v  : the resulting trajectory is represented in Fig. 7 (a), starting from the second swing-by of Ganymede. The sequence for this solution is GGCGC, and the objective value, i.e. the final relative velocity, is 1.96 km/s v   . The solution is practically ballistic.
ACO-MGA identified a solution, with an objective value of 1.91 km/s v   , equal flyby sequence and comparable transfer times per leg (see Fig. 7 (b)). Table 4 compares some parameters of the 2D solution with the re-optimised 3D solution.
Slight differences are due to the different models, and mainly to the changes of inclination that are needed in the 3D solution. When re-optimized with the full 3D model, the 2D solution identified by ACO-MGA yields the reference solution in Fig. 7 (a). This demonstrates that the assumptions underneath the 2D model are consistent and allow for identifying first guess solutions that converge locally to higher fidelity solutions with comparable cost and mission characteristics.
Note that the solution chosen by ESA, and used here as a comparison, is not the best from the point of view of the arrival velocity. In fact, this solution was chosen by ESA following a trade off, taking into account not only the v  , but also the presence of DSMs, the total time of flight, the radiation dose, and the arrival velocity vector at Callisto.
ACO-MGA found also solutions with lower v  . The best one has 1.709 km/s v   , and corresponds to the trajectory plotted in Fig. 8. The swing-by sequence is the same, GGCGC, but the total time of flight is much longer (92 days), since the first GC leg performs 1 full revolution and the CG leg 2 full revolutions, and a 10 m/s DSM is added to the second leg.
The 100 runs returned a number of feasible solutions. Four different solutions with a low value of v  are described in Table   5, and respectively plotted in Fig. 9. Note that we deliberately decided not to pick the best four solutions according to lowest v  , but rather to show the variety of different good solutions that were found through this approach, and demonstrate that this tool could be useful for finding valid alternative options. In the same way, Table 6 and Fig. 10 describe four different solutions among those with short total time of flight. Note that all the solutions with low v  follow the same planetary sequence GGCGC, while the short solutions use also the sequence GGGCC, but at the cost of a higher v  .
Finally, the best solution associated to each feasible sequence obtained by ACO-MGA is represented in the bar plot in Table   7.      Table 5). The first leg is not plotted.   Table 6). The first leg is not plotted.

Cassini case study
Cassini is the ESA-NASA mission to Saturn. The planetary sequence designed for the mission, EVVEJS is particularly long, allowing a substantial saving of propellant.
Since the launch date is not taken into account in the optimisation, in the following test it is considered fixed. In a real mission design case, where the launch date is to be determined, the entire launch window can be discretised with a given time step, and a systematic scan of several dates within the whole launch window should be run. The launch direction 0  is also kept fixed in these tests, although it is easy to find heuristics for determining the value of this parameter, or discretise it and include it in the optimisation process as an additional variable.
first bunch of runs and 6000 for the second bunch of runs. In particular, for 4000 function evaluations, the number of iterations allocated to the first step was 200 and the number of iterations allocated to the second step was 300. For 6000 function evaluations, instead, the number of iterations allocated to the first step was increased to 600 and the number of iterations allocated to the second step was left equal to 300. With these settings, a run of ACO-MGA requires, on average, 1900 function evaluations and 161 s, if the upper limit is 4000, and 3300 function evaluations and 273 s, if the upper limit is 6000 evaluations. This is considerably faster than the exhaustive scan of the solution space and the number of function evaluations is lower than in the case of GATBX and NSGA-II, however the computational time for each single run of GATBX, which is fully coded in MATLAB ® , is, on average, half of the one required to ACO-MGA as the access to the taboo and feasible lists is currently not optimal and coded in MATLAB ® . Current developments have addressed this bottleneck and the results of the new implementation will be presented in future papers. NSGA-II is even faster as the code is fully in C. Therefore, a fair comparison would allow GATBX and NSGA-II to perform a higher number of function evaluations.
The parameters used for GATBX and NSGA-II are reported in Table 8. The comparative results for the two sets of runs are shown in Table 10. It can be seen that, for 1900 function evaluations, ACO-MGA found feasible solutions in 91% of the runs, compared to 25% of GATBX and 26% of NSGA-II. The average ACO-MGA solution is also slightly better than GATBX, and considerably better than NSGA-II. The performances of ACO-MGA increase significantly by using 3300 evaluations: all the runs produce a feasible solution, and in 80% of the cases, the best solution found is below 16 km/s. The average value of the solution also decreases to 15.434 km/s. It is interesting to note that, for GATBX, the average best solution found with 6000 evaluations is higher than for 4000: this is partly balanced by the fact that it finds feasible solutions in 28% of the runs. Another thing worth noticing is that NSGA-II finds more often feasible solutions than GATBX, but their quality is in average worse.
The best solution found by ACO-MGA (sequence EVVEJS) has an objective value of 6.9686 km/s: The characteristics of this solution can be found in Table 9, compared to the best solution found for the Earth-Saturn transfer problem (see http:// www.esa.int/gsp/ACT/inf/op/globopt/edvdvdedjds.htm). The trajectory of the ACO-MGA solution is shown in Fig. 11 (a), while the 3D reference solution is in Fig. 11 (b).
It is worth to underline that the superior computing time of ACO-MGA is also due to the fact that for all runs, including the few that are not feasible, the algorithm finds partial solutions that are feasible up to the last transfer. As a consequence, the taboo lists of the last planetary encounter tend to grow significantly. A possible solution, on top of optimising the access to the list, would be to plan backward from the target planet in the case of transfers to the outer part of the solar system.   It is interesting to sort the feasible sequences found by ACO-MGA according to the best objective value that they can achieve. Table 11 shows the objective value, the relative velocity at arrival, the total DSM cost and the time of flight of the best solution found for each sequence. Note that every solution has a trajectory associated to it, modelled as shown before, and thus taking into account the phasing problem. This means that these solutions could be re-optimised with a more detailed model (in particular including the third dimension), leading to actual transfer solutions.
However, differently from the previous test case, the alternative solutions have a considerably longer time of flight (> 13 years) and higher v  , and therefore we consider them not interesting. Other promising solutions could possibly be found changing the definition of the problem, i.e. varying the launch angle and the sets , i j q .

Launch date analysis
As mentioned before, the algorithm, at the current state, does not perform any kind of search on the launch date 0 t . In fact, this variable is not even included in the solution vector s. Rather, if the launch date is not fixed, but a launch window is available, a systematic scan can be performed to find the best launch date, and the corresponding solutions. This procedure is not always applicable: in fact, if re-running the algorithm for a small change in the launch date, the solutions that ACO-MGA finds are substantially different, then the systematic scan along 0 t is not feasible, and this variable must be included in the optimization process. If, on the other hand, a small displacement along 0 t causes a small change in the best solution found (e.g. same planetary sequence, possibly different types of transfers, similar objective value), then the systematic scan over 0 t is applicable to the identification of promising launch dates.
In order to verify this assertion, a test was run using the BepiColombo mission as a reference case [16]. BepiColombo is a multiple gravity assist mission to Mercury, currently under study at ESA and JAXA. In Ref. [16] an optimal transfer solution is provided, using two swing-bys of Venus to reach Mercury (sequence EVVMe). The optimal launch date is the 15 August 2013, i.e. 0 4974.5 d, MJD2000 t   .
ACO-MGA was run, leaving the choice of the swing-by sequence and the other transfer parameters, like the number of revolutions, free. The objective was to minimise the relative velocity v  at Mercury.
Five different launch dates, in a window of 10 days around the one chosen by ESA, were considered, and for each one of them, 100 runs of ACO-MGA were performed. The best solutions for each launch date can be found in Table 12. Note that, the optimal launch date occurs 1 day before 0 t  , while earlier or later launch dates appear to be less convenient. In addition, all the solutions have the same planetary sequence.
The discrepancy between the value of v  found by ACO-MGA and the one in [16] has two reasons: the first is that ACO-MGA does not take into account the inclination of the planets, and the orbit of Mercury is highly inclined. The second is that the ESA solution was found as a part of a longer trajectory, and thus with a different objective.
The same reasons explain why, according to ACO-MGA, the ideal launch date is 1 day earlier. As a matter of fact, this is not a problem, and a subsequent local optimisation of the ACO-MGA solutions with a full model would tune the launch date.
This test demonstrates that the 2D model and the search process in ACO-MGA yield a distribution of solutions, along the departure time axis, that is consistent with the existing reference solution for this mission. The convergence of ACO-MGA is robust against variations of the launch date, as it consistently provides, solutions with unchanged planetary sequence that display a small variations in the cost function for small variations of the launch date. These properties together with the small computing time suggest that ACO-MGA can be used to systematically scan the launch dates in search for an optimal one. Note that, as the model in ACO-MGA is intended for the generation of first guess solutions, the scan of the launch dates is expected to provide only an estimated location of the optimal point along the launch date coordinate.

CONCLUSION
The paper introduced a novel formulation of the automatic complete trajectory planning problem and proposed a new algorithm (ACO-MGA), based on the ant colony paradigm, to generate optimal solutions to this problem. Each solution is a complete, scheduled plan. A specific trajectory model was developed to efficiently generate families of scheduled trajectories for multi-gravity assist transfers, once a plan is available. The 2D trajectory model proved to be accurate enough to closely reproduce known MGA transfers even with moderate inclinations. Furthermore, the scheduling of the trajectories is fast and reliable allowing for the evaluations of thousands of plans in a short time.
ACO-MGA operates an effective search in the finite space of possible plans. The algorithm demonstrated a remarkable ability to find good solutions with a very high success rate, outperforming known implementations of genetic algorithms. As ACO-MGA requires very little information on the MGA problem under investigation, it represents a valuable tool for the complete automatic design of future space missions. Future work aims at a more efficient handling of the lists, which is currently the major bottleneck of the ACO-MGA implementation.