A differential evolution-based optimization tool for interplanetary transfer trajectory design

The extremely sensitive and highly nonlinear search space of interplanetary transfer trajectory design bring about big challenges on global optimization. As a representative, the current known best solution of the global trajectory optimization problem (GTOP) designed by the European space agency (ESA) is very hard to be found. To deal with this difficulty, a powerful differential evolution-based optimization tool named COoperative Differential Evolution (CODE) is proposed in this paper. CODE employs a two-stage evolutionary process, which concentrates on learning global structure in the earlier process, and tends to self-adaptively learn the structures of different local spaces. Besides, considering the spatial distribution of global optimum on different problems and the gradient information on different variables, a multiple boundary check technique has been employed. Also, Covariance Matrix Adaptation Evolutionary Strategies (CMA-ES) is used as a local optimizer. The previous studies have shown that a specific swarm intelligent optimization algorithm usually can solve only one or two GTOP problems. However, the experimental test results show that CODE can find the current known best solutions of Cassini1 and Sagas directly, and the cooperation with CMA-ES can solve Cassini2, GTOC1, Messenger (reduced) and Rosetta. For the most complicated Messenger (full) problem, even though CODE cannot find the current known best solution, the found best solution with objective function equaling to 3.38 km/s is still a level that other swarm intelligent algorithms cannot easily reach.


Introduction
Nowadays, heuristic search algorithms like Particle Swarm Optimization (PSO) [1] and Differential Evolution (DE) [2] have been widely employed to optimize the interplanetary trajectory design problem. The difficulty of optimizing interplanetary trajectory design problem is caused by extremely sensitive and highly nonlinear search space [3]. Figure 1 shows a highly nonlinear search space of Messenger (full) problem around the current known best solution. The real interplanetary exploration mission usually has to consider two basic objectives: minimum mission period and minimum fuel consumption. However, even though only one objective function is considered for mission optimization, the optimum is still very hard to be found. In 2005, the GTOP database published by ESA includes some interplanetary transfer trajectory design problems, which can be simply divided by multiple gravity-assist (MGA) problems and multiple gravity assists with deep space manoeuvre (MGA-1DSM) problems. In the GTOP database, MGA problem includes Cassini1 and GTOC1; MGA-1DSM problem includes Sagas, Rosetta, Cassini2, Messenger (reduced) and Messenger (full).
With the given search space and special boundaries on each variable, the GTOP problems have shown big challenges on global optimization. Three aspects of the challenge can be described as follows: • The attraction basin of global optimum is rather smaller than those of other local optima, and unluckily, this optimal basin is usually hidden in the "not promising" local spaces. During the evolutionary process, it is so difficult to locate these not promising local spaces when considering the population motivation based on greedy gradient information. This causes a result that the general swarm intelligent algorithms can only find sub-optimal results, and the difference between global optimum and found local optima is not acceptable.
• As a result of the above-mentioned reason, the optima of GTOP problems, like that of Messenger (full), are unknown, and the current known best results are still constantly updated.
Even though Messenger (full) is just a 26-D problem, its current known best result has been continuedly updated for 8 years from 2009 to 2017. However, we still do not know whether the current known best results are optima. • Even though the attraction basin of current known best result has been positioned, it does not mean that the current known best result has been found. It is also very difficult to define a "local space". For the extremely nonlinear search space, a "local search space" possibly contains smaller attraction basin. Hence, it will be a challenge to determine when performing the local optimizer. The most effective mode is to maximize the searchability of global optimizer and passing the found best solution to a powerful local optimizer.
Up to now, the known best solutions of these problems have been repeatedly updated, and some swarm intelligent algorithms like DE have shown promising performance. However, except MIDACO [4], which is a large-scale optimization platform based on tons of function evaluations, other swarm intelligent algorithms can only solve one or two benchmarks. In the previous studies, Zuo et al. [2][5] [6] designed a search framework named EP DE based on space exploration. it can try to explore promising regions from sensitive global search space and can find the global optimum of Cassini1 problem. Addis et al. [7] provided a remarkable solution for Rosetta problem when considering combining a multi-start technique and a two-phase local optimizer. Ampatzis et al. [8] found the optimum of Cassini1 problem based on an approximation method cooperating with a neural network. Biazzini et al. [9] proposed a distributed Hyper-Heuristics (HHs) for continuous search space and evaluated its performance on Cassini1 and Cassini2 problems. Biscani et al. [10] introduced an optimization software named PaGMO which is developed by Advanced Concepts Among all the related studies, MIDACO is a framework combining the supercomputing and ant colony optimization (ACO). As two crucial elements, the combination of swarm intelligence and high-performance computing has been proved to have the ability to solve the most complicated Messenger (full) problem. One of its original versions can reach an average result of 14.961 km/s using an average function evaluation of 232780 [13] among overall 1000 runs, and the found best result is 6.399 km/s. Compared with this result, the proposed CODE in this paper can reach an average result of 10.385 km/s, and the found best result is 5.7 km/s among 50 runs, which is a significant advantage. The improved version of MIDACO, named MXHPC, can reach a result of round 2.0km/s using a function evaluation of about 7 × 10 10 . With a multi-kernel gauss probability function, MIDACO can reach a refined result near to 1.959km/s. Recently, MIDACO was further improved to find a 2.1 km/s result within the lowest function evaluation of about 5.5 × 10 9 .
The success of MIDACO encourages us to concentrate on the development of one of the two crucial elements, namely heuristic search algorithms. The outstanding performance of DE on the interplanetary transfer trajectory optimization given in the previous studies is the motivation why it has been selected as the research target. Our previous research on the improved DE algorithm for interplanetary trajectory design, named CLDE, has been published recently [14]. In this research, a strictly limited case learning theory has been proposed. To test the performance, six GTOP problems were employed. In CLDE, two extended versions, including G-CLDE for global search and L-CLDE for local search are provided. However, L-CLDE is an iteratively called optimizer and usually generates a very high computational cost. For the most complicated problem Messenger (full), the found best result 6.487 km/s is far away from that of the known best result 1.9577 km/s. To further explore the potential of CLDE, an enhanced version CODE is now developed.
Compared with the CLDE algorithm, CODE needs lower function evaluations, and can still maintain an acceptable success rate of finding the current known best solutions. Compared with already published search algorithms on the GTOP database, CODE presents better universality, not only working on one or two problems. The main contributions of this paper can be summarized as follows: • A new two-stage DE has been proposed in this paper. At the first stage, a crossover operator with binary values and a mutation operator considering global exploration, local exploitation and random search are developed. The crossover operator considering binary values is motivated by the big valley and stripe structure of search space shown in Figure 1. At the second stage, the mutation factors of all individuals are adaptively generated according to the selected mutated individuals. Also, a multiple boundaries check technique has been developed and selectively employed based on the gradient information of different variables.
• A new evaluation criterion is proposed in this paper. Based on 150000 function evaluations, CODE can find the optima of Cassini1 and Sagas directly. Using 380000 function evaluations, CODE cooperated with CMA-ES can find the optima of Cassini2, Rosetta and GTOC1 with success rates of 2/500, 7/500 and 2/500 respectively. With 580000 function evaluations, CODE cooperated with CMA-ES can find the optimum of messenger (reduced) with a success rate of 1/500. The given maximum function evaluations and the achieved success rate on different GTOP problems can be a comparison criterion for the following researches.
The rest of this paper is organized as follows: Section 2 introduces the mathematical model of Lambert transfer, GA, MGA and MGA-1DSM problems. Section 3 introduces the DE algorithm. Section 4 presents the proposed CODE algorithm in details. The performance evaluation of CODE is given in Section 5. Finally, Section 6 concludes the paper and discusses future work.

Lambert's problem
A particle p moving around the center M can be regarded as a two-body motion. If denoting the position vectors of p at P 1 and P 2 as r 1 and r 2 respectively, the variation of true anomaly is where r 1 = √ r 1 ⋅ r 1 and r 2 = √ r 2 ⋅ r 2 . In Lambert's problem, the transfer time ∆t from P 1 to P 2 is a function of semi-major axis a, r 1 + r 2 and the chord length c between P 1 and P 2 . Assuming the transfer time ∆t from P 1 to P 2 is given, the Lambert's problem is transformed to determining the trajectory from P 1 to P 2 . When r 1 and v 1 are given, the motion status (r 2 and v 2 ) of any trajectory point can be obtained by where f and g stand for Lagrange coefficients, andḟ andġ stand for the corresponding derivatives.
Combing equations 2 and 3, The equation calculating v 2 is further represented as Considering the conservation of angular momentum, there is Therefore, equation 4 can be also expressed as It can be seen that the determination of orbital elements will be easy once r 1 and v 1 or r 2 and v 2 are already given, and the solution of Lambert's problem is transformed to calculating f , g,ḟ andġ. The methods of calculating f , g,ḟ andġ are introduced in [15], [16] and [17].

GA
The geometry simulation given in Figure 2 shows two kinds of GA: leading-side planetary GA and trailing-side planetary GA. The angles between gravity and velocity directions in these two cases are larger and smaller than 90 o , thus will generate velocity deceleration and acceleration respectively. The heliocentric velocity variation ∆v can be obtained by where V ∞1 and V ∞2 are hyperbolic velocities with the same scalar size, and δ is the angle between them. V is the velocity of gravity-assist planet. δ can be calculated by  where µ p is the gravitational constant, r p is the fly-by height

MGA
The geometry simulation given in Figure 3 is a MGA trajectory modelled with as conic arcs connecting three celestial bodies. In general, the encounter sequence of a MGA trajectory is given by P = {p 0 , ⋯, p i , ⋯, p n } , n > 1, 0 < i < n, where p 0 is the Earth, p i is the potential gravity-assist planet, and p n is the destination planet. Each conic arc is a Lambert's problem [18] considering the departure time at planet p i and the arrival time at planet p i+1 . Usually, the connection of two neighboring conic arcs has to consider a velocity change ∆v i at the pericenter of gravity assist hyperbola. Denoting the maneuver required to be captured by the destination planet as ∆v f , the objective function for evaluating a possible solution is where is departure fuel consumption.
B 0 and α 0 are declination and the right ascension of escape asymptote respectively. For a specific encounter sequence, the optimization of time sequence, namely optimal schedule of spacecraft, is the main work of this paper. It can be presented by where t 0 is the departure time, and T i is the transfer time length from p i−1 to p i .

MGA-1DSM problem
Compared with MGA problem, MGA-1DSM problem as shown in Figure 4 allows a chemical pulse at any time during each trajectory leg. The more uncertain variables will enlarge the problem dimension and make the optimization more difficult. A MGA-1DSM trajectory including n planets can be represented as where t 0 stands for departure time, [V 0 , α 0 , B 0 ] is the initial hyperbolic velocity, b incl k and r p k are related to the k th swing-by conditions, and T k is the transfer time length for each leg.
Each trajectory leg contains two conic arcs which are connected by the applied DSM i . Among them, the first conic arc of first trajectory leg is integrated with the initial velocity V 0 from t 0 to where ∆DSM i is the pulse implied by DSM in the i th leg and ∆v f is the manoeuvre required to be captured by the destination planet. More details about MGA and MGA-DSM can be learned from referrences [19], [20] and [21].

Introduction of DE
DE [22] is an evolutionary algorithm which tries to find the global optimum by repetitively generating a population and performing three operators (mutation operator, crossover operator and selection operator) on population in each generation, until the current generation G reaches the maximum generation number G max . The generated population in generation G can be denoted by where N P is population size, and x i,G is an individual. As shown in Algorithm 1, the operators mutation and crossover are controlled by the two parameters F and CR, respectively. The original version of these operators are introduced as follows: where G is generation number, N P is population size, F is mutation factor falling within [0, 1].
has exceeded the boundaries, it has to be re-generated by where LB j and U B j are the lower and upper boundaries on the j th variable respectively.
Crossover The individual is mixed with the mutated solution, using the following scheme, to yield the trial solution where x ji,G , if rand(0, 1) > CR and j ≠ r(D) where rand(0, 1) can generate a random number number between (0,1). r(D) can generate an integer less than D.
Selection A greedy selection scheme as follows is used to select the better individuals in the search process.
Only when the function value of u i,G+1 is less than that of

Introduction of CODE
In CODE, to improve the space exploration ability of global structure at the start and accelerate convergence at the last, the linear population reduction in L-SHADE [23] is employed. The initial population size and minimum population size are denoted as N P init and N P min respectively. If denoting the maximum function evaluations and already used function evaluations as M AX N F E and N F E respectively, the population size, denoted as N P can be obtained by Output: The found best solution x best,Gmax .
Perform the boundary check B1 on individuals; Perform crossover on v i,G+1 to generate u i,G+1 ; Print solution x best,Gmax ; includes three types: random re-generation (B1), interpolation re-generation (B2) and opposite re-generation (B3). Among them, B1 was employed in DE and can be presented as This randomness will enlarge the diversity of the population and can help to jump local spaces.
However, it also causes potential damage to the convergence. Hence it can be called more at an earlier stage.
For B2, it was employed in JADE algorithm [24], and can be presented as The feature of B2 is not in favour of enlarging the population diversity, while contributes to the convergence. Hence, it is expected to work at a later stage.
B3 inspired by [25] can be shown as and K2.
Output: The found best solution x best,Gmax .
Copy P to archive A; Update population N P according to equation 20; Update T according to equation 25; Perform the multiple boundary check B on individuals; Perform crossover on v i,G+1 to generate u i,G+1 ; Randomly replace a member of A with x i,G ; where pLB j and pU B j are the lower and upper boundaries of population distribution respectively. They can be updated by Under the motivation of the B3, the generated trial individuals will move farther. Therefore, it also can be employed at an earlier stage.
Two parameters K1 and K2 can be set to control the employment of these three boundary check methods. Besides K1 and K2, This paper also proposes a variable-related boundary check.
Here, the distribution of individuals on each variable is regarded as an estimation of linearization or nonlinearity. If a strong nonlinearity has shown on one variable, B1 and B3 will have high potential to improve the population diversity. Otherwise, B2 is more suitable for linearization and weak nonlinearity. The distribution T j of individuals at generation G on j th variable can be presented as The distribution information T indicates the stage of employing B. Given a threshold denoted as E, the calling of each boundary check method can be presented as Why this technique has been considered is that different boundary checks can motivate the individuals to different spatial positions, can also re-generate new elements at different variable locations.
For all N P individuals, three operators (two-stage mutation, crossover and selection) are employed to generate offsprings. Two control parameters F and CR are employed for mutation and crossover respectively. These operators are introduced as follows: • Two-stage Mutation (M). In each generation, mutation operator is firstly performed on the whole population. For {x i,G i = 1, 2, ..., N P }, its mutant solution generated in the first stage (M1) is as follows: where G stands for the generation number, N P stands for the population size. 1 ≤ r1 ≤ N P , 1 ≤ r2 ≤ N P × pbest, 1 ≤ r3 ≤ N P , and r1 ≠ r2.sx is the union of population set P and archive A storing the replaced parent individuals those have failed in comparison with their child individuals. To keep the stored information available, the size of A keeps fixed to N P .
That means when N P members have been already stored, the new member will replace a random one. x i,G .F is generated by where nx 0, , and F (x m,G ) means the objective function value of solution x m,G . Therefore, the advantage of M2 is a faster convergence speed. The mutation schedule can be presented as • Crossover. The target individual is mixed with the mutated individual, using the following scheme, to generate the trial solution x ji,G , if rand(0, 1) > x i,G .CR and j ≠ r(D) rand(0, 1) is a random number generator to generate number between (0, 1). r(D) is an integer generator to generate integer smaller than D. x i,G .CR is generated by If u ji,G+1 < L j or u ji,G+1 > U j , u ji,G+1 is re-generated using above mentioned boundary check technique.
• Selection. A greedy selection scheme is used to select the better individuals between parent and child vectors. As an example of minimization problem, this selection can be performed by If the fitness value of u i,G+1 is better than that of x i,G , x i,G+1 is replaced by u i,G+1 , and the replaced solutions will be stored in A. • Artificial Bee Colony optimization [34]: G max = 750, limit = 20;

Settings of K1 and K2
To explore the influences of  Table 1. The best mean results found by CODE with different combinations of K1 and K2 are also ranked, and the top three results are given in Table 2. From the experimental results, some interesting phenomenon can be observed as follows: • In the most cases, CODE has a very stable performance on Cassini1 problem, and the search result usually comes to 5.30 km/s. In some cases, CODE has successfully found the current known best solution with an objective function value of 4.93 km/s. In these cases, a surprise is that CODE can find the current known best result twice, when the combinations of K1 and • Although the dimension of Sagas is 12, higher than those of GTOC1 and Cassini1, CODE still can find the current known best result, and maintain a success rate higher than those of Cassini1 and GTOC1. In some cases, the success rates of finding the current known best result of Sagas reaches 5/50. These cases suggest that the good combinations of K1     • Having a lower dimension than that of Sagas, the current known best result of GTOC1 is more difficult to be found. Three combinations of K1 and K2 have helped CODE to find results lower than -1581700. Even though there is still a small gap to reach the current known best result, the found best result has shown the potential to be refined to the current known best result.
• Except for Messenger (full), Messenger (reduced) and Cassini2 are also difficult to be optimized. The found best result by CODE for Cassini2 has a difference of 0.1 in comparison with the current known best result. The found best result by CODE for Messenger (reduced) has a difference of 0.2 in comparison with the current known best result. The global search ability of CODE on these two problems should be further improved.
The found best result, mean result and worst result of GTOP problems by L-SHADE variants and PYGMO algorithms are given in Table 3 and Table 4 respectively. To provide a correct interpretation of the results and obtain the average rankings of CODE and other competitors, the non-parametric statistical Friedman test is employed in this section as similarly done in [39].
This test is calculated using the software KEEL [40] and the results in terms of the best and the mean results are shown in Tables 5 and 6 respectively. For the ranking results, smaller value means better ranking. CODE is proved to have the best overall performance in the comparison, and EIG L-SHADE SPS ranks the second with ranking values of 4.1667 and 3.6667 respectively.
As a relatively old algorithm, PSO also shows promising performance on the GTOP problems, comparable with L-SHADE algorithm.

Cooperation with CMAES algorithm
Since CODE can find very near results to the current known best ones within 150000 function evaluations, this section continues to explore whether a slightly larger function evaluation for CODE and the cooperation with CMA-ES algorithm can work. Here CMA-ES algorithm is regarded as a local optimizer, and the starting point is provided by CODE. For Cassini1 and Sagas problems, CODE has found the current known best result, here the local optimizer is not considered, and the maximum function evaluations are still fixed to 150000. For other problems, the objective is still to find the current known best results and try to improve the success rate. The maximum function evaluation is slightly improved from 150000 to 180000, and the maximum function evaluation for CMA-ES is 200000. The population size of CMA-ES algorithm is ⌊50 × log(D)⌋. For the Messenger (reduced) problem, within the given 200000 function evaluations, CMA-ES cannot help to find the current known best result, so the maximum function evaluation is improved to 400000. The success      rates of reaching the current known best result on each problem are given in Table 7. From the experimental results, some interesting phenomenon can be observed as follows: • Except for Messenger (full) problem, the cooperation of CODE and CMA-ES has found all the current known best results of Cassini1, Cassini2, GTOC1, Rosetta, Sagas and Messenger (reduced). The success rate of solving Messenger (reduced) is the lowest, while the highest success rate has been achieved on Sagas. Except for Messenger (full), Messenger (reduced) and Cassini2 are the most difficult problems to be optimized.
• The optimization performance of swarm intelligent algorithms is not so stale. CODE has achieved the success rates of 2/50 and 5/50 on Cassini1 and Sagas problems respectively among 50 independent runs. However, among the 500 independent runs, CODE does not achieve the expected success rates of 20/500 and 50/500, while realizes lower success rates of 6/500 and 20/500 respectively.
• CMA-ES has an important impact on optimizing the Messenger (full) problem. Within the 380000 function evaluations, the CODE cooperated with CMA-ES have found a result of 3.379885 km/s, which is a promising performance as a swarm intelligent algorithm.

The distribution of found local optima
To study the distribution of found local optima in section 5.   each interval. By comparing the success rates given in Table 7 and the distribution results, two interesting conclusions can be drawn as follows: • If the distribution of local optima is like a Gaussian function, the optimization difficulty is relatively higher. This type includes Messenger (full), Messenger (reduced) and GTOC1 problems. In this type, the individuals are more possibly stuck in sub-optimal local search spaces.
• If the distribution of local optima is like other shapes, it is easier to find the optima. This type including Cassini1, Sagas, Cassini2 and Rosetta has less sub-optimal local spaces and is easier to find the optimal solution.

Conclusions and future work
The effectiveness and efficiency of CODE have been illustrated by the experimental results in this paper. The influence of K1 and K2 combinations has been studied, and the optimal parameter settings are concluded for different GTOP problems. The cooperation of CODE with CMA-ES has helped to find the current known best results of Cassini1, Cassini2, GTOC1, Rosetta, Sagas and Messenger (reduced). The employment of "Linear population reduction" in L-SHADE does can improve the global exploration ability, and the comparisons with L-SHADE and its improved variants also illustrate the effectiveness of other CODE modules. In comparisons with other swarm intelligent algorithms in PYGMO, CODE achieves an obvious advantage. CODE has the best optimization performance on Sagas, and considering the relatively higher dimension of Rosetta, CODE also gives a surprise. The same surprise is also seen from the test result on Messenger (full) problem, where CODE can find a result of 5.407 km/s, and the CODE cooperated with CMA-ES can find a result of 3.38 km/s. This is a level not so easy for other optimization algorithms to reach.
However, the experimental results indicate that CODE still needs to be further developed from the following aspects: • The mutation in CODE is a two-stage operator, which concentrates on global exploration and fast convergence at different stages. Compared with the mutation operator "DE/current-topbest/1" in JADE and L-SHADE, the mutation operator in CODE just uses one difference vector, not two ones. This point should be re-considered in the following DE versions, at least at an earlier evolutionary stage for enlarging the randomness.