International Journal of Aerospace Engineering
Volume 2018 (2018), Article ID 4560173, 17 pages
https://doi.org/10.1155/2018/4560173

Research Article

Enhanced Hybrid Differential Evolution for Earth-Moon Low-Energy Transfer Trajectory Optimization

Yanyun Zhang,¹ Lei Peng,1,2 Guangming Dai,1,2 and Maocai Wang1,2

¹School of Computer, China University of Geosciences, No. 388 LuMo Road, Hongshan District, Wuhan, China
²Hubei Key Laboratory of Intelligent Geo-Information Processing, China University of Geosciences, Wuhan, China

Correspondence should be addressed to Lei Peng; lei.peng@cug.edu.cn

Received 27 November 2017; Revised 4 February 2018; Accepted 5 March 2018; Published 3 May 2018

Academic Editor: Christian Circi

Copyright © 2018 Yanyun Zhang et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

It is known that the optimization of the Earth-Moon low-energy transfer trajectory is extremely sensitive with the initial condition chosen to search. In order to find the proper initial parameter values of Earth-Moon low-energy transfer trajectory faster and obtain more accurate solutions with high stability, in this paper, an efficient hybridized differential evolution (DE) algorithm with a mix reinitialization strategy (DEMR) is presented. The mix reinitialization strategy is implemented based on a set of archived superior solutions to ensure both the search efficiency and the reliability for the optimization problem. And by using DE as the global optimizer, DEMR can optimize the Earth-Moon low-energy transfer trajectory without knowing an exact initial condition. To further validate the performance of DEMR, experiments on benchmark functions have also been done. Compared with peer algorithms on both the Earth-Moon low-energy transfer problem and benchmark functions, DEMR can obtain relatively better results in terms of the quality of the final solutions, robustness, and convergence speed.

1. Introduction

Deep space exploration, as a very important hotspot in space exploration, has attracted attention of various countries, while as the first step to explore the deep space, the lunar exploration has been receiving renewed attention, such as the NASA’s LADEE, ARTEMIS, and GRAIL missions and CNSA’s Chang’e 3 mission. To carry out lunar exploration successfully, the first priority is the Earth-Moon transfer trajectory design. The traditional way to construct the transfer trajectory is Hohmann transfer [1], which utilizes the two-body system dynamics model, and it was later proven optimal analytically by Barrar in [2]. Short as the duration of the Hohmann transfer orbit is, the fuel consumption, however, is high. And for a long-term unmanned space mission, it is essential to consider the increase of scientific loads and the reduction of the energy consumptions at the same time. Low-cost alternative methods rather than conventional methods have been explored as the knowledge of the solution space is expanded, and algorithms that employ the dynamical relationships are developed.

With the discovery of new types of particular solutions in the three-body problem, there has been accelerated interest in missions utilizing trajectories near libration points, and a number of missions have already incorporated the periodic halo orbits and/or quasiperiodic Lissajous trajectories as a part of the trajectory design, such as ISEE-3 (1978 launch), WIND (1994 launch), SOHO (1995 launch), and ACE (1997 launch). And for the trajectory design of a liberation point mission, dynamical systems theory was applied by Barden and Howell to better understand the geometry of the phase space in the three-body problem via stable and unstable manifolds [3]. In later years, low-energy trajectory designing based on hyperbolic invariant manifolds has also been proposed in [4–6].

Besides, in 1987, as a different methodology, the notion of a weak stability boundary (WSB) was first introduced heuristically by Belbruno for designing fuel-efficient space missions [7]. After that, a low-energy Earth-Moon transfer trajectory from the Earth to the Moon leading to ballistic capture was constructed based on WSB by Belbruno and Miller in [8], which used approximately 18% less cost in total than the Hohmann transfer to insert a spacecraft into a circular orbit about the Moon. In their work, it is presented that in the framework of the circular restricted three-body problem (CR3BP), a WSB region exists in the phase space near the Moon in the Earth-Moon system, where the spacecraft will be captured by the Moon’s gravity into lunar orbit requiring little or no fuel (called ballistic capture). The use of this transfer was firstly demonstrated by Japan’s Hiten spacecraft [9], and thus, an important role of WSB is to design the similar ballistic capture for space missions. Thereafter, Belbruno and Carrico [10], Belló et al. [11], and Belbruno [12] did further research on WSB gradually, which enhanced the theoretical basis of this kind of low-energy transfer. What is more, it has also been suggested that the hyperbolic invariant manifold method can be used to explain the trajectories obtained through the WSB method in some cases, while WSB exists even in models where the hyperbolic invariant manifolds are no longer well defined. It seems possible that the WSB may turn out to provide a good substitute for the hyperbolic invariant manifolds in such models [13].

In the past few years, varieties of Earth-Moon low-energy trajectories incorporating new concepts have been developed, such as the low-energy trajectory only in the three-body Earth-Moon system, corresponding to the interior capture around the Moon [14], and the low-energy trajectory which considers the gravity of the Sun and is also based on the concept of the WSB theory, corresponding to the so-called exterior WSB transfer [8]. The Earth-Moon low-energy transfer trajectory discussed in this paper belongs to the exterior WSB transfer. This type of low-energy transfer trajectory is mainly achieved by taking advantage of the invariant manifold structures associated with the 3D halo orbits (or 2D Lyapunov orbits) in the vicinity of the libration points in the Sun-Earth system and Earth-Moon system and can be modelled as two coupled CR3BPs. In fact, the Japanese Hiten mission was a paradigmatic example for a class of low-energy Earth-to-Moon orbits obtained by considering the gravitational effects of the Earth, the Moon, and the Sun on the motion of the spacecraft simultaneously [8, 9, 14, 15] and was later found to be related to the hyperbolic invariant manifolds of the CR3BP by employing the patched three-body approximation [5, 16, 17]. Gómez et al. [18] had also presented the use of the restricted three-body problem’s invariant manifold structural to design the low-energy transfer orbit in detail. And de Sousa-Silva and Terra [19] presented a survey of preliminary Earth-to-Moon transfers in the patched three-body approach, using the planar case of the CR3BP.

In general, there are two kinds of numerical optimization methods to deal with the trajectory designing, namely, the deterministic methods and the stochastic methods [20]. For the deterministic methods, they are local methods in nature and a suitable initial solution in the region of convergence, which is difficult to obtain for the space trajectory problem in the real system, is also needed [21]. For example, for the design of the Earth-Moon low-energy transfer trajectory based on the patched three-body approximation, the setting of the initial integral point in the Poincaré section of the unstable manifold of the Lyapunov orbit around the Sun-Earth L2 and the stable manifold of the Lyapunov orbit around the Earth-Moon L2 is sensitive to the traditional deterministic methods, which may easily cause the inefficiency of the traditional deterministic methods. To overcome these difficulties, stochastic methods, which do not require initial solution and are global optimization methods, are usually adopted.

Recently, evolutionary algorithm- (EA-) based optimization techniques for space trajectory design have received significant attention. As a branch of EAs, genetic algorithm (GA) has been proposed by several authors [22–25] for trajectory optimization. However, drawbacks (e.g., low speed, premature convergence, and degradation for highly interactive fitness functions) have been observed for GAs. Topputo et al. [26] proposed a hybridization of EAs with sequential quadratic programming method to look for the best intersection between invariant manifolds associated to the departure and arrival Sun-Planet-Spacecraft systems. By extension of evolutionary neurocontrol (ENC) to the planetary case of an Earth-Moon transfer, Ohndorf et al. [27] showed that automatic optimization of low-thrust trajectories crossing the boundary of the spheres of influence between the involved celestial bodies was feasible with ENC. The particle swarm optimization (PSO) algorithm has also been applied in the field of space trajectory widely. In [28], PSO was applied to a variety of space trajectory optimization problem, including the determination of periodic orbits in the context of the CR3BP and the optimization of (impulsive and finite thrust) orbital transfers, which was proved to be quite effective in finding the optimal solution to all of the applications. While in [29], a further study about the optimization problem of impulsive orbital transfer was proposed. Zotes and Peñas [30] used the swarm algorithms to optimize the tuning parameters of the trajectory in both the single-criteria case and multicriteria cases. In the single-criteria case, only the fuel consumption was considered as a variable to be minimized, while in the multicriteria case, the fuel consumption and the total time of the mission were simultaneously taken into account. Although PSO is easy to implement, and the convergent speed is fast, it still easily traps in a local optimum [31].

Differential evolution (DE) [32] is a simple yet powerful evolutionary algorithm developed by Storn and Price for global continuous optimization and has been successfully used in a variety of domains [33, 34], as well as the trajectory optimization. In 2007, an improved DE algorithm was applied by Lei et al. [21] into two-impulse transfer trajectory problem. In 2011, Attia [35] presented an adaptive probability of crossover technique as a variation of the differential evolution algorithm, for optimal parameter estimation in the general curve-fitting problem. And the technique was successfully applied to the determination of orbital elements of a spectroscopic binary system (eta Bootis). In 2015, Lu et al. [36] used DE to optimize the parameters describing the patched manifolds of low-energy transfer orbit to Mars with multibody environment effectively. In 2016, Nath and Ramanan [37] developed alternative single-level schemes for the mission design to a halo orbit around the libration points from the Earth based on DE. And this DE-based scheme for the transfer trajectory identified the precise location on the halo orbit that needs minimum energy for insertion and avoided exploration of multiple points. In 2017, in [38], a multistart DE enhanced with a deflection strategy with strong global exploration and bypassing abilities was adopted as a search engine to find multiple globally optimal regions in which potential periodic orbits (POs) were located.

The motivations and contributions of this paper are as follows:

(1)The Earth-Moon low-energy transfer trajectory is extremely sensitive with the changes of the initial condition in the Poincaré section, and with a small change of the initial condition (such as a small change in velocity at a fixed point), the destination of an orbit can be changed dramatically. So proper choices of the initial parameter values are key to the optimization of the Earth-Moon low-energy transfer trajectory.(2)As a powerful global search algorithm, DE possesses several unique features such as simple implementation and low complexity, while its performance is limited because of the slow converging at the late stage of evolution. Given this reason, the proposed algorithm combined DE with a local search sequential quadratic programming (SQP) to enhance the local search ability of DE and finally improve the search accuracy.(3)An extra storage is used in DEMR to collect the superior solutions, and it is associated with the reinitialization of the population later. The reinitialization procedure is triggered when the local search fails to find a better solution than DE. And the procedure includes both the random reinitialization and the biased reinitialization based on the extra storage.

To evaluate the performance of the proposed approach in the optimization of Earth-Moon low-energy transfer and the performance of the proposed algorithm itself, experiments on Earth-Moon low-energy transfer problem both in the two-dimensional space and three-dimensional space, as well as the experiments on the benchmark CEC2005 are conducted, and the results are compared with several peer algorithms.

The rest of this study is organized as follows. In Section 2, the Earth-Moon low-energy transfer trajectory model and objective function used in this work are briefly introduced. Section 3 presents the proposed DEMR method in detail, followed by the experimental results and analysis in Section 4. Lastly, in Section 5, this work is concluded, and several possible avenues of future work are discussed.

2. Formula and Optimization Model

2.1. Formula about CR3BP

2.1.1. Formula in Two-Dimensional Space

The planar CR3BP describes the motion of a spacecraft under the gravitational attractions of two primaries and [39]. Choose a rotating coordinate system so that the origin is at the center of mass between and . See Figure 1.

Figure 1: The planar CR3BP in the rotating frame.

Let be the position of the spacecraft in the plane. The motion equations of the spacecraft in rotating frame with normalized coordinates are as follows, and for simplicity and convenience, the differential equation is nondimensional and the results are then dimensionalized to show in tables.

where and are the second derivative of and against time , respectively, while and denote the derivative of and against time , and and represent the partial derivatives of the effective potential against and , respectively. The effective potential is given by

and

where and are the masses of and , respectively, while and are the distances from the spacecraft to the two celestial bodies and , respectively. With normalization, the normalized variables, including the distance between and , the sum of their masses, and their angular velocity around the barycenter, are normalized to one. We assume that the mass of is greater than the mass of , namely, , then , known as the mass parameter of the planar CR3BP, is the dimensionless mass of , while the dimensionless mass of is . In the rotating frame, and are fixed at and , respectively.

2.1.2. Formula in Three-Dimensional Space

Extending the planar problem to a three-dimensional situation, the above equations can be expanded as follows:

(1)The motion equations of the spacecraft in three-dimensional case with normalized coordinates:

(2)The effective potential in three-dimensional case:

    and

2.2. Optimization Model

The design of the low-energy transition orbit on the two-dimensional plane is a problem in the four-dimensional phase space . The main idea is to take the point of the intersection of the unstable manifold associated with the Lyapunov orbits around the Sun-Earth L2 and the stable manifold associated with the Lyapunov orbits around the Earth-Moon L2 in the Poincaré section as the integral initial point, and then do the following steps:

(1)From the initial point, make the reverse integration along the unstable manifold associated with the Lyapunov orbits around the Sun-Earth L2 to the Earth parking orbit, completing the orbit design of the launch section.(2)Adjust the current energy of the initial point and finally the energy of the initial point is the same as that of the current lunar Moon system.(3)From the initial point, make the positive integration along the stable manifold associated with the Lyapunov orbits around the Earth-Moon L2 to the parking orbit of the Moon, completing the design of the capture section.

Figure 2 is a brief look at the associated manifolds and a simulated Earth-Moon low-energy transfer trajectory.

Figure 2: The manifolds and the low-energy transfer trajectory in the rotating frame.

To get the initial patched point, the position of Poincaré section is selected firstly. At present, the location of the Earth is usually selected as the location of the Poincaré section. This also means the -axis position of the four-dimensional phase space is fixed. With the calculation of the curve in both the Sun-Earth system and the Earth-Moon system in Poincaré section, the initial and can be selected according to their intersection status. Since the energy constants of the Sun-Earth system and the Earth-Moon system are also determined as the amplitudes of the periodic orbits are given, the speed of the initial point in the -axis direction can also be determined. Thus, the initial patched point is obtained.

Given both the location and speed in the -axis direction is determined, the decision variable of the optimization model is

where and are the position and velocity of the spacecraft in the -axis direction, respectively, while denotes the conversion angle which is needed during the conversion process from the Earth-Moon system to the Sun-Earth system in the Poincaré section, namely, the angle between the Sun-Earth connection line and the Earth-Moon connection line. Different conversion angles make the different states of the intersection between the associated Earth-Moon stable manifold and the Sun-Earth unstable manifold. It is obvious that the conversion angle that makes the two manifolds intersect should be selected.

Similarly, the decision variable can also be expanded in the three-dimensional Earth-Moon problem as follows:

From the above discussion and the illustration of Figure 2, three impulsive vectors may be required in total during the process of low-energy transfer from the Earth to the Moon.

(1)Impulsive vector is required to put the spacecraft from the Earth parking orbit into the Sun-Earth stable manifold.(2)Impulsive vector is required to connect the two segments of the transfer trajectory in the Poincaré section.(3)Impulsive vector is required from unstable invariant manifolds associated to Lyapunov orbits around the Earth-Moon L2 point to the lunar parking orbit.

The optimization objective is to optimize the total cost:

3. Optimization Method

3.1. Differential Evolution

Differential evolution is a population-based optimization algorithm. In general, DE is used to optimize certain properties of a system by pertinently choosing the system parameters. For convenience, these parameters are usually described as a vector .

where represents the number of parameters or the dimension of a problem. The parameter vector is also known as an individual, a solution or a target vector and individuals form a so-called population . To evaluate the performance of the system under a specific set of parameters, an evaluate function (or object function, fitness function) is defined. So, the optimization aims at finding out the object vector , which minimizes (or maximizes) such an objective function , that is, for all and the is the best solution.

DE works on a population of D-dimensional objective vectors and goes through the same cycle of stages with the traditional EAs: initialization, mutation, crossover, and selection. Each round of the circle is called a generation and only if at least one of the terminated conditions is achieved can the loop be stopped. Either the maximum of evaluation number or the best solution may trigger the termination of the algorithm. A general mutation operator may be denoted as , where stands for the type of base vector in the mutation strategy, specifies the number of difference vectors in the mutation strategy, and refers to the crossover strategy. The five most frequently used mutation strategies suggested by Storn and Price have been presented in [40, 41]. An example of is illustrated in Algorithm 1, which means in this algorithm, DE is conducted with a randomly chosen base vector and only one difference vector at the mutation stage and with the exponential crossover strategy at crossover stage.

1: Generate the initial population , define as the -th individual in , is the  population size, is the number of function evaluations in each run, MaxNFEs is the number of max function evaluation, is the number of decision variable, is the mutation factor, is crossover rate.

2:

3: Evaluate the fitness for the each individual in .

4:

5: whileNFEs ≤ MaxNFEsdo

6:      for to do

7:            Select and

8:            

9:            

10:            

11:            repeat

12:                  

13:                  

14:                  

15:            until or

16:            Evaluate the trial vector

17:            

18:            if is better than then

19:                  

20:            end if

21:      end for

22: end while

Algorithm 1: DE/rand/1/exp.

3.2. DEMR

The new algorithm proposed in this paper is called DEMR. DEMR employs standard DE as a global exploration heuristic in combination with SQP as a local search mechanism and a reinitialization scheme for enhancement.

3.2.1. Hybridization of DE and SQP

In DEMR, local search is triggered according to two contraction criteria. One is the roughness in objective space, defined as

where is the objective value of the ith individual among the current population, while refers to the average value for the objective values of the current populations, and is the population size.

The other is the maximum distance in decision space:

where is the ith individual of the current population , while is the th individual.

is a measure of the degree of contraction in objective space while is a measure of the degree of contraction in decision space. When the population converges to small region, these two indexes can reflect the situation more or less. Especially, for , it illustrates the contraction of the population directly and should describe the situation more accurately. Only if at least one of those two is smaller than the given threshold , local search is started, shown in Algorithm 2.

1: Initialize population , define as the -th vector in at generation , and is the best vector in the current generation, while is the best solution we have found so far. is the population size, refers to the number of function evaluations, MaxNFEs is the number of max function evaluation. denotes the mutation factor and is crossover rate. and are used for contraction criterion.

refers to the resultant new local optimum found by SQP. and are the lower and upper boundaries of the variables. denotes a set which collects the superior solutions at the late stage of evolution.

2: , ,

3: Evaluate the fitness for the each individual in .

4:

5: whileNFEs < MaxNFEsdo

6:      for to do

7:            Select with the select probabilities and .

8:            ifNFEs < 0.2 × MaxNFEsthen

9:                 Using DE/rand/1/bin to generate

10:            else

11:                 Using DE/best/1/exp to generate

12:            end if

13:            Evaluate the trial vector

14:            

15:            if is better than then

16:                 

17:            end if

18:            Update and

19:      end for

20:      Calculate the contraction criteria

21:      if then

22:            Pick up the as the initial point of the local search SQP, and find the local optimum

23:            Update and ,

24:            if then

25:                 Re-Initialize the whole population as Algorithm 3 shows.

26:            end if

27:      end if

28: end while

Algorithm 2: DEMR.

1: is the number of re-initialization while is the maximum of the number of randomly re-initialization. denotes a set which collects the superior solutions at the late stage of evolution.

2: if then

3:        Re-initialize the whole population randomly

4: else

5:        Calculate the average value and standard deviation of the elements in

6:        Re-initialize the population as:

7: end if

Algorithm 3: Mix Reinitialization.

3.2.2. SQP

SQP methods represent the state of the art in nonlinear programming methods for constrained optimization, which were developed mainly by Han [42, 43], and Powell [44], based on the initial work of Wilson [45]. The method closely mimic Newton’s method for constrained optimization just as is done for unconstrained optimization. At each major iteration, an approximation is made of the Hessian of the Lagrangian function using a quasi-Newton updating method. This is then used to generate a QP subproblem whose solution is used to form a search direction for a line search procedure. Schittkowski [46], for example, has implemented and tested a version that outperforms every other tested method in terms of efficiency, accuracy, and percentage of successful solutions, over a large number of test problems. And Boggs and Tolle [47] have also proved it outperforms every other nonlinear programming method in terms of efficiency, accuracy, and percentage of successful solutions over a large number of test problems. An overview of SQP is found in Fletcher [48], Gill et al. [49], Powell [50], and Hock and Schittkowski [51].

3.2.3. Mix Reinitialization

As DEMR advances, SQP will be triggered to do the further exploitation when the whole population is relatively concentrated in the decision space or the objective space. And DEMR enables the reinitialization of the population only when the local search does not improve the best individual in the population, which helps DEMR jump out from the attraction of local optimal solutions. In order to guarantee the diversity of the population as well as the convergence speed at the same time, a mixed reinitialization, shown in Algorithm 3, is applied. Firstly, to guarantee the diversity of the population, the population is reinitialized randomly in the first times, and thus the algorithm can distribute individuals in the whole search space uniformly to a large extent; secondly, after times’ reinitialization, the population is reinitialized as follows:

where denotes the ith individual in the population, and is the population size. and stand for the average value and standard deviation of the elements in , respectively, while denotes the set which collects the superior solutions at the late stage of evolution. Once the local search is triggered, the extra storage starts its work, which means each time the algorithm goes to the local search process, the best solution it has found so far will be stored in . And after several times’ random reinitialization, the algorithm has searched superior solutions in the whole space as far as possible, and the relatively better solutions should have been stored in . At this time, reinitializing the population based on cannot only find the optimal value with great probability but also speed up the convergence rate. And the size of will not be restricted here since the is active only at the late stage of evolution, and the storage capacity will not be great.

4. Experiments

4.1. Experiments for Earth-Moon Low-Energy Transfer Trajectory Optimization

In this section, two groups of experiments are conducted, including the experiments with the two-dimensional planar CR3BP model and experiments in the three-dimensional space. In order to clarify the performance of DEMR in the Earth-Moon low-energy transfer trajectory optimization, comparisons have been made with peer algorithms, such as CMA-ES [52–54], LBBO [55], GL-25 [56], SaDE [57], SHADE [58], and MPEDE [59]. The reasons to choose those six algorithms are not only the competitive performance of them but also they represent the state of art in ES, BBO, GA, and DE, respectively. And among these algorithms, CMA-ES is a well-known local optimizer, GL-25 runs a global real-coded GA (RCGA) during the 25% of the available evaluations and triggers a local RCGA afterwards, and LBBO also includes global and local search operators, along with reinitialization logic, to improve performance. For SaDE and SHADE, the former is a classic adaptive DE, while the latter is an recent efficient adaptive DE which ranks the 4th at the 2013 IEEE International Conference on Evolutionary Computation (CEC2013). And finally, MPEDE is referred to as a multipopulation-based approach to realize an ensemble of multiple strategies, which simultaneously consists of three mutation strategies. All of these algorithms are applied in the optimization of Earth-Moon low-energy transfer trajectory, respectively, and for each algorithm, 25 rounds of independent tests are conducted on this real-case problem. The performances of the algorithms are evaluated in terms of best value, worst value, median value, mean value, and standard deviation of the 25 rounds tests. For the tolerance level of the experiments, all of the results are accurate to the integer part and keep a decimal as the unit is m/s.

The low-energy transfer trajectory is from the 200 km Earth parking orbit to the 100 km Moon parking orbit in this work. The amplitude of Sun-Earth L2 Lyapunov orbit is 201,000 km and Earth-Moon L2 Lyapunov orbit is 15,000 km. The parameters of DEMR are set as follows (Table 1).

Table 1

Parameters in DEMR
Value

Crossover rate

Scale factor

Dimension

Population size

Contraction criteria

Reinitialization control

Maximum evaluations
MaxNFEs = 10,000

Storn and Price [32] have indicated that the performance of DE is sensitive to the value of control parameters, and good choices of and are 0.5 and 0.9, respectively. In order to avoid tuning the parameters, a self-adaptive parameter control technology is adopted according to the scheme introduced above: and , where is a normal distribution. And according to the study of Storn and Price [32], a reasonable value for the population size could be chosen between and ; here, is set as 30. For the setting of the contraction criteria, and are suggested, and the maximum number of randomly reinitialization is set as 5. When the maximum number of function evaluations (MaxNFEs), set to , is reached, the algorithm will be stopped.

According to the analysis in [60, 61], is a widely used classic DE strategy, of which the slower convergence speed and stronger exploration capability make it more suitable for the early evolution process, while belongs to the kind of greedy strategies and always converges faster; thus, it is more reasonable to adopt in the late stage. And the binomial crossover has consistently performed well on several benchmarks and are known to be explorative, while the exponential crossover is yet another form which has been successful in exploitation [60]. Based on these analyses, the mixed strategy of and is employed that the is applied when the number of evaluation is less than 0.2 × MaxNFEs, while is used for the rest evolving generations, which will help to distribute the individuals in the whole decision space in the early stage of evolution and speed up converging in the late stage. The setting of the mixing ratio refers to the corresponding setting in [60].

As for the other algorithms used in the comparison, the corresponding parameter settings are the same in the original papers describing the algorithms.

4.1.1. Earth-Moon Low-Energy Transfer Trajectory in Two-Dimensional Space

According to some of our previous work, the ranges of decision variables set as below are more conducive to the optimization (Table 2).

Table 2

Decision variables
Range

Position in direction

Velocity in direction

Conversion angle

Table 3 summarizes the costs of the Earth-Moon low-energy transfer optimized by different algorithms, and the optimal value among all algorithms is highlighted in boldface. According to the results, DEMR outperforms other algorithms in all of the evaluation indexes, including the best value, worst value, median value, mean value, and standard deviation. For the best value, DEMR has a total cost for the Earth-Moon low-energy transfer, and the corresponding for worst value, median value, and mean value is , , and , respectively.

Table 3: Results for Earth-Moon low-energy transfer trajectory optimization in two-dimensional space.

Algorithms
Best (m/s)
Rank
Worst (m/s)
Rank
Median (m/s)
Rank
Mean (m/s)
Rank
Std
Rank

SaDE
3910.7
5
3919.4
2
3914.9
2
3914.5
2
2.2993
2

CMA-ES
3908.4
2
4080.1
7
4013.9
7
4002.3
7
42.3711
7

GL-25
3908.9
3
3922.2
3
3915.2
3
3915.5
3
3.0661
3

LBBO
3908.4
2
4007.0
6
3931.8
6
3933.0
6
26.3087
6

MPEDE
3913.4
6
3934.8
5
3928.3
5
3926.3
5
5.7536
5

SHADE
3910.4
4
3931.4
4
3919.8
4
3921.0
4
5.3871
4

DEMR
3908.3
1
3912.8
1
3909.0
1
3909.2
1
0.9387
1

The best value illustrates the maximum search accuracy that the algorithm can obtain on the Earth-Moon low-energy transfer trajectory designing to a certain extent. From Table 3, we can sort the algorithms according to the best values easily, and the rankings of SaDE, CMA-ES, GL-25, LBBO, MPEDE, SHADE, and DEMR are 5, 2, 3, 2, 6, 4, and 1, respectively. It is obvious that the algorithms DEMR, LBBO, CMA-ES, and GL-25 perform better than the rest, and that is because these algorithms either contain local operators or have strong local search ability, which may contribute to the accuracy. The rankings for other evaluation indexes are all the same, namely, 2, 7, 3, 6, 5, 4, and 1 for SaDE, CMA-ES, GL-25, LBBO, MPEDE, SHADE, and DEMR, respectively. That is to say the performances of all the algorithms are consistent on the worst value, median value, mean value, and standard deviation.

The worst value can briefly reflect the maximum deviation from the optimal solution. In this case, the local search algorithm CMA-ES shows great disadvantages since it is more likely to trap in the local optimal solution, while the algorithms SaDE, GL-25, SHADE, and DEMR have relatively better performance. It can be inferred that a simple local search algorithm does not necessarily have an advantage, even though the exploitation occupies a very important position for the Earth-Moon low-energy transfer optimizing. And the outstanding advantage of DEMR on the worst value also proves its reliability, when comparing with other algorithms. This is not only because DEMR adopts a powerful global optimizer, DE, just as SaDE and SHADE does, and more importantly, the mix reinitialization strategy employed by DEMR plays an important role in balancing the exploration and exploitation.

The good performances on both the mean value and standard deviation indicate DEMR can obtain a better value easily and owns a more stable optimization effect in Earth-Moon low-energy transfer trajectory optimization. And compared to the mean value, the median value excludes the influence of the extreme cases, such as the best value and the worst value, and shows the performance of an algorithm in the usual case. Normally, DEMR can achieve a greater advantage over other algorithms. According to test results for all evaluation indexes, DEMR achieves overall best optimization performance among seven algorithms. Hence, it is a good alternative for the optimizing of the Earth-Moon low-energy transfer trajectory.

To describe the convergence speeds of these 7 algorithms, 14 record points are set during the optimizing process for each algorithm, and when the evaluation number reaches any record point, the lowest total cost found so far will be recorded. The record points are . For each algorithm, 25-round independent tests are conducted on the Earth-Moon low-energy transfer trajectory optimization problem, and then the average value of the 25 recorded at each record point is calculated and the calculated 14 average values are finally used to illustrate the convergence graph as Figure 3 shows.

Figure 3: The convergence graph of the algorithms on the Earth-Moon low-energy transfer trajectory optimization problem.

From Figure 3, the general converging trends for all of the algorithms are shown in the big graph. The abscissa indicates the number of evaluations, and the coordinates represent the corresponding value of the total cost . While in order to get a clearer look at the details, a partially magnified graph is drawn in the upper left corner. And the scope of the magnified graph is . According to the results, SaDE converges the fastest in the early stage and DEMR is the second fastest, but DEMR can reach a better optimization accuracy than SaDE in the late stage. And as a local optimizer, it seems easier for CMA-ES to fall into preconvergence, resulting a worse optimization accuracy.

Figure 4 is a simple simulation result of the Earth-Moon low-energy transfer trajectory obtained by DEMR. The simulation result is drawn based on the best optimizing solution of DEMR, and the corresponding parameter setting is .

Figure 4: The best Earth-Moon low-energy transfer trajectory optimized by DEMR. .

According to the discussions above, DEMR can get a relatively better solution with a high stability. However, all of these results are obtained within a relatively small search space, which has high possibility to get a relatively good solution. So, to further argue the stability of the proposed algorithm on the design of Earth-Moon low-energy transfer trajectory, an extra group of experiments on a wider setting of the decision variables has been done. The widened ranges are as follows (Table 4):

Table 4

Decision variables
Range

Position in direction

Velocity in direction

Conversion angle

The results with the widened ranges are shown in Table 5. It is obvious that DEMR can also achieve a relatively better result in all of the evaluation indexes with a good stability. Comparing the standard deviations in Tables 3 and 5, the standard deviations are generally increased, and the range expansions have a big influence on almost all of the algorithms compared DEMR, except SaDE. This phenomenon can also been found in the column of the worst values. That means, with the expanding of the ranges, the possibility to find the bad solutions is increasing and the compared algorithms starts to become unstable, while DEMR can still hold a good solution stably.

Table 5: Results for Earth-Moon low-energy transfer trajectory optimization in two-dimensional space with widened ranges.

Algorithms
Best (m/s)
Rank
Worst (m/s)
Rank
Median (m/s)
Rank
Mean (m/s)
Rank
Std
Rank

SaDE
3912.0
3
3930.4
2
3917.8
2
3918.5
2
4.6163
2

CMA-ES
3908.1
1
4416.8
6
3961.3
4
3997.0
6
125.0366
6

GL-25
3908.1
1
4190.5
5
3918.6
3
3947.3
3
62.4762
4

LBBO
3908.6
2
4640.2
7
3984.5
6
4002.8
7
139.6852
7

MPEDE
3934.0
5
4075.9
3
3992.5
7
3993.2
5
33.2864
3

SHADE
3923.8
4
4185.6
4
3968.3
5
3990.8
4
63.9009
5

DEMR
3908.1
1
3913.8
1
3909.5
1
3909.7
1
1.2773
1

4.1.2. Earth-Moon Low-Energy Transfer Trajectory in Three-Dimensional Space

In this section, the optimization of the design for Earth-Moon low-energy transfer trajectory is extended to the three-dimensional space. To reflect the stability of different algorithms, the widened ranges of the decision variables in last section are used here and the ranges of the location and speed in the -axis direction will be supplied, while all of other parameter settings are the same as before. The ranges of the decision variables are as follows (Table 6):

Table 6

Decision variables
Range

Position in direction

Velocity in direction

Position in direction

Velocity in direction

Conversion angle

From Table 7, it is interesting that the best values of all the algorithms are good results and close to each other, which means all of the algorithms can achieve a relatively superior solution. But for the worst values, the mean values, and the standard derivations, the differences between the algorithms are a bit exaggerated. Actually, the exaggerated differences of the mean values and the standard derivations are caused by the abnormal values, such as those in the worst value column. From the results in the worst values column, it is obvious that the values of CMA-ES, GL-25, LBBO, and MPEDE are a bit overblown, which means these algorithms do not seem to converge in this case. In other words, CMA-ES, GL-25, LBBO, and MPEDE may be ineffective in the three-dimensional space. For other three algorithms, the performances of SaDE and DEMR are not only equally matched but also top-notch for all of the evaluation indexes, while the performance of SHADE is slightly inferior compared to these two, especially for the stability. Excluding the effects of extreme results, the median values help to show relatively normal performances of these algorithm. Among all of the algorithms, CMA-ES and DEMR can achieve the same optimization performance with the median values, followed by SaDE, GL-25, SHADE, LBBO, and MPEDE in order. In conclusion, DEMR can still remain stability of its superior performance in the three-dimensional space, while the performance of SaDE becomes better comparing to the two-dimensional case and is comparable to the performance of DEMR. Other algorithms, excluding SHADE, have the risk of failure in this three-dimensional case.

Table 7: Results for Earth-Moon low-energy transfer trajectory optimization in three-dimensional space with widened ranges.

Algorithms
Best (m/s)
Rank
Worst (m/s)
Rank
Median (m/s)
Rank
Mean (m/s)
Rank
Std
Rank

SaDE
3903.3
4
3905.5
1
3904.4
2
3904.3
2
0.5879
1

CMA-ES
3902.9
1
147694.8
6
3903.5
1
18581.0.0
7
4.0468E+ 3
6

GL-25
3903.6
5
12876.1
4
3905.2
3
4269.3
4
1.7931E+ 3
4

LBBO
3903.2
3
25788.9
6
3916.6
5
6679.8
6
4.9129E+ 3
7

MPEDE
3907.0
6
12606.4
3
3921.6
6
4469.8
5
1.9637E+ 3
3

SHADE
3903.6
5
3934.1
2
3911.0
4
3912.0
3
7.3546
5

DEMR
3903.0
2
3905.5
1
3903.5
1
3903.7
1
0.6650
2

According to the analysis above, although CMA-ES, GL-25, LBBO, and MPEDE may be ineffective, their performances with the best values and median values are acceptable. So, the success rate of each algorithm is calculated with the results of the 25 independent tests. Firstly, to determine if an optimization is successful, we set a threshold. If the result is below the threshold, the optimization process will be considered successful and vice versa. According to the total cost of Hohmann transfer, , in [14], the threshold is set as . Table 8 gives the number of successful results and the success rates among the 25 independent tests for all the algorithms. From the table, the “Success time(s)” refers to the number of successful optimization results among the 25 independent tests, while the “Success rate” denotes the ratio of “Success time(s)” to total number of tests, namely, 25 here. It can be found that SaDE, SHADE, and DEMR can be very successful in the optimization of the problem. LBBO and CMA-ES are more likely to get a bad solution, even their best values and median values are not that bad. While for GL-25 and MPEDE, their failures seem to be more by chance.

Table 8: Simple investigation of the success rates for all of the algorithms in three-dimensional space.

Algorithm
SaDE
CMA-ES
GL-25
LBBO
MPEDE
SHADE
DEMR

Success time(s)
25
19
24
16
23
25
25

Success rate
100%
76%
96%
64%
92%
100%
100%

4.2. Experiments on Benchmark Functions

To further confirm the effectiveness of DEMR, another group of experiments between DEMR and the same algorithms used in experiments for Earth-Moon low-energy transfer trajectory optimization (CMA-ES, LBBO, GL-25, SaDE, MPEDE, and SHADE) is conducted on 25 benchmark functions from CEC2005 special session [62] on real parameter optimization. Among the functions, F1–F5 are unimodal, while the rest are the multimodal, including basic functions from F6 to F12, expanded functions F13 and F14, and hybrid composition functions from F15 to F25. Each of the 7 algorithms under tests is executed 25 times with respect to each test function. Considering the optimization of Earth-Moon low-energy transfer trajectory is a low-dimensional problem, the validation experiment here on benchmark is only conducted under the dimension .

The performance of the algorithms is evaluated in terms of function error value (FEV), defined as , where is the global optimum of the test function, while denotes the actual optimal solution the algorithms find. The parameters in DEMR are set the same with those in the former experiment for Earth-Moon low-energy transfer trajectory optimization, except the maximum evaluations MaxNFEs and the maximum number of randomly reinitialization , which are set as MaxNFEs = 10,000D and . As for the other algorithms used in the comparison, the corresponding parameter settings are the same in the original papers describing the algorithms.

To provide a correct interpretation of the results, the nonparametric statistical tests, Wilcoxon rank-sum test and Friedman test, are employed in this section as similarly done in [63]. Wilcoxon signed-rank test at and is used to test the statistical significance of the experimental results between two algorithms in single-problem and multiple-problem. Friedman test is employed to obtain the average rankings of all the compared algorithms. The Wilcoxon rank-sum test in single-problem is calculated using MATLAB, while the multiple-problem Wilcoxon signed-rank test and Friedman test are calculated using the software KEEL [64].

The mean error and standard deviation of all the algorithms for solving each function over 25 independent runs when dimension is 10 are shown in Table 9. The symbols “−,” “+,” and “=”denote that the performance of the corresponding algorithm is, respectively, worse than, better than, and similar to DEMR, respectively, according to the Wilcoxon rank-sum test at . Generally speaking, DEMR performs better than all of the algorithms except SHADE. Among the 25 test functions, DEMR is significantly better than CMA-ES, LBBO, GL-25, SaDE, MPEDE, and SHADE in 15, 11, 16, 12, 10, and 9 functions, respectively, while the numbers of functions where CMA-ES, LBBO, GL-25, SaDE, MPEDE, and SHADE perform significantly better are 6, 7, 6, 10, 8, and 13. And for the remaining functions, DEMR provides error values comparable to those of the compared algorithms. According to the results, the advantages of DEMR over CMA-ES, LBBO, GL-25, SaDE, and MPEDE still persist, while the performance of DEMR is overtaken by SHADE. For SaDE and SHADE, the adaptive strategies play a big role as the dimension increases, and for MPEDE, its characteristics of multipopulation and multiple mutation strategies are also very efficient. According to the results of Table 9, the performances of the algorithms are compared on different types of benchmark functions in Table 10. And from Table 10, it can be seen that DEMR is more suitable for the multimodal functions rather than the unimodal functions when comparing to other algorithms, since the mix reinitialization not only enhances its ability to escape from the local optimal, but also helps to balance the diversity of the population and the convergence speed. Especially, for the hybrid composition functions, DEMR has a relatively obvious advantage.

Table 9: Comparisons of mean error and standard deviation between DEMR and other six algorithms on CEC2005 functions.

Func.
DEMR
CMA-ES
LBBO
GL-25
SaDE
MPEDE
SHADE

F1

⁺
⁺
⁺
⁻
⁻
⁻

F2

⁺
⁺
⁺
⁺
⁺
⁺

F3

⁺
⁺
⁻
⁻
3⁺
⁺

F4

⁻
⁻
⁺
⁺
⁺
⁺

F5

⁺
⁼
2⁻
⁺
⁺
0⁺

F6

⁼
⁼
⁻
⁻
⁻
⁺

F7

⁺
⁺
⁺
⁺
⁺
⁺

F8

⁼
⁼
⁻
⁻
⁻
⁻

F9

⁻
⁺
⁻
⁻
⁻
⁻

F10

⁻
⁻
⁼
⁼
⁼
⁺

F11

⁺
⁼
⁻
⁺
⁼
⁺

F12

⁻
⁻
⁻
⁻
⁻
⁻

F13

⁻
⁻
⁻
⁻
⁻
⁼

F14

⁻
⁼
⁻
⁺
⁼
⁺

F15

⁻
⁼
⁻
⁺
⁺
⁼

F16

⁻
⁻
⁼
⁼
⁼
9⁺

F17

4⁻
⁻
⁺
⁺
⁼
⁺

F18

⁻
⁻
⁻
⁻
⁻
⁻

F19

⁻
⁻
⁻
⁻
⁻
⁻

F20

⁻
⁻
⁻
⁻
⁻
⁻

F21

⁻
⁻
⁻
⁻
⁻
⁻

F22

⁼
⁼
⁼
⁼
⁼
⁻

F23

⁻
⁻
⁻
⁻
⁼
⁼

F24

⁼
⁺
⁻
⁺
⁺
⁺

F25

⁻
⁺
⁺
⁺
⁺
⁺

⁺
6
7
6
10
8
13

⁻
15
11
16
12
10
9

⁼
4
7
3
3
7
3

indicates that the intermediate results are reported at when the algorithms obtain the global optimum; “−,” “+,” and “=”denote that the performance of this algorithm is, respectively, worse than, better than, and similar to DEMR according to the Wilcoxon rank-sum test at .

Table 10: Comparison between DEMR and other algorithms on different types of functions.

(a) F1–F5: Unimodal functions

Algorithms
CMA-ES
LBBO
GL-25
SaDE
MPEDE
SHADE

4
3
3
3
4
4

1
1
2
2
1
1

=
0
1
0
0
0
0

(b) F6–F12: Multimodal-basic functions

Algorithms
CMA-ES
LBBO
GL-25
SaDE
MPEDE
SHADE

2
2
1
2
1
4

3
2
5
4
4
3

=
2
3
1
1
2
0

(c) F13-F14: Multimodal-expanded functions

Algorithms
CMA-ES
LBBO
GL-25
SaDE
MPEDE
SHADE

0
0
0
1
0
1

2
1
2
1
1
0

=
0
1
0
0
1
1

(d) F15–F25: Multimodal-hybrid composition functions

Algorithms
CMA-ES
LBBO
GL-25
SaDE
MPEDE
SHADE

0
2
2
4
3
4

9
7
7
5
4
5

=
2
2
2
2
4
2

The Wilcoxon results of the multiple-problem statistical analysis are listed in Table 11. The in Table 11 is a measure of the case where DEMR outperforms the compared algorithm and the greater the value is, the greater difference between DEMR and the compared algorithm is. Similarly, the is a measure of the case where DEMR is worse than the compared algorithm, and the size of also reflects difference between DEMR and the compared algorithm in this situation. According to Table 11, significant differences are obtained in three cases (DEMR versus CMA-ES, DEMR versus LBBO, and DEMR versus GL-25), which means that DEMR performs significantly better than CMA-ES and SaDE at both and . For MPEDE, it performs worse than DEMR when , while is comparable with DEMR at . And for other algorithms, it shows they have similar optimization effectiveness with DEMR, since the value goes beyond the at both and . The conclusions here based on multiple-problem statistical analysis are clearly consistent with the previous Wilcoxon rank-sum test analysis.

Table 11: Results obtained by the multiple-problem Wilcoxon test for DEMR and other six algorithms on CEC2005.

VS
R+
R−
value

CMA-ES
274
26
0.000375

LBBO
238.5
86.5
0.039554

GL-25
261
64
0.007727

SaDE
183
117
0.338495
=
=

MPEDE
216
84
0.057433
=

SHADE
183
117
0.338495
=
=

Eight functions are chosen from the 4 types of functions in CEC2005 benchmark to get a brief look at the convergence speed, including 2 unimodal functions, 2 basic functions (multimodal), 1 expanded functions (multimodal), and 3 hybrid composition functions (multimodal). The chosen functions are , and , and the convergence graphs are shown in Figure 5. In order to describe the convergences more clearly, the logarithmic coordinate is used in Figures 5(a), 5(b), 5(d) and 5(e). It can be observed that, DEMR can achieve both a fast convergence speed and a high search accuracy for multimodal functions, but for the unimodal functions and , DEMR converges a little early. What is more, it should be noted that, in the convergence Figure 5(a) for , the convergence of LBBO after about has not been shown, since LBBO has already found the optimal value by that time, which cannot be illustrated by the logarithmic coordinate. And the convergence of SHADE in Figure 5(b) is not shown after about for the same reason.

(a) The convergence graph on the unimodal function f1

(b) The convergence graph on the unimodal function f5

(c) The convergence graph on the basic function (multimodal) f9

(d) The convergence graph on the basic function (multimodal) f12

(e) The convergence graph on the expanded function (multimodal) f13

(f) The convergence graph on the hybrid composition function (multimodal) f18

(g) The convergence graph on the hybrid composition function (multimodal) f21

(h) The convergence graph on the hybrid composition function (multimodal) f23

Figure 5: The convergence graphs on different types of CEC2005 benchmark functions.

In order to get a more intuitive comparison for all of the algorithms from the overall performance, the results of the Friedman statistical test are shown in Figure 6. It shows that DEMR ranks the second with , behind SHADE with . And SaDE follows DEMR with . The results are consistent with previous analyses in Wilcoxon tests.

Figure 6: Average rankings for all algorithms on 25 benchmark functions.

4.3. Experiments on Extra Parameters

DEMR introduces its own parameters , , and , of which and determine the time to trigger the local search algorithm SQP, and controls the maximum number of random reinitialization. It is believed that these three parameters have a big influence on the performance because of the roles they play in DEMR. So an appropriate setting of these parameters is necessary.

As shown in Table 12 and Figure 7, the Friedman rankings for different combinations of and are illustrated: , , where is a constant . From Table 12, the setting is the best combination, followed by combinations and . Figure 7 is more intuitive for the performance of these settings. Generally speaking, the combination with both small and small is more reasonable. Table 13 gives a brief look at the test results on different with contraction criteria combination . And it can be seen that, for CEC2005 benchmark functions, is the most appropriate choice, while for the design of Earth-Moon low-energy transfer trajectory, the performance of DEMR ranks the first when .

Table 12: Results obtained by the Friedman test for contraction criteria combinations on CEC2005.

Friedman

1
1.5
2
2.5
3

1
10.34
11.08
10.12
13.28
12.32

1.5
11.14
11.74
13.98
12.24
13.90

2
15.04
14.66
13.08
15.46
16.84

2.5
11.14
12.78
14.10
14.48
14.66

3
9.12
12.34
13.90
14.94
13.32

Figure 7: Average Friedman rankings for contraction criteria combinations on 25 benchmark functions.

Table 13: Average Friedman rankings for the maximum number of random reinitialization.

(a) Average Friedman rankings for the maximum number of random reinitialization on 25 benchmark functions

Friedman rankings

5
3.16

10
3.52

15
2.72

20
2.58

25
3.02

(b) Performances for the maximum number of random reinitialization on the Earth-Moon low-energy transfer problem

Average (m/s)

5
3909.1

10
3909.8

15
3909.4

20
3909.3

25
3909.7

5. Conclusion

To optimize the Earth-Moon low-energy transfer trajectory, an hybridized differential evolution algorithm with an efficient reinitialization strategy, DEMR, is proposed in this paper. The DEMR method firstly hybridizes DE with local search algorithm SQP to balance the exploration and exploitation and then employs a mix reinitialization procedure to further avoid falling into local optimum. The effectiveness and efficiency of DEMR are validated by comparing the other six peer algorithms, which represent the state of art in ES, BBO, GA, and DE, respectively, on the Earth-Moon low-energy transfer problem in both the two-dimensional space and three-dimensional space. With an overall consideration for all the results, including the achieved best optimization performance, the success rate, and the stability, DEMR is a relatively better choice when comparing with those selected algorithms. The superiority of DEMR is also verified after comparing it with these six algorithms on benchmark problems. Although simple sensitive tests of the extra parameters have been done, a further more comprehensive analysis of the parameters’ features is necessary, and alternative contraction criteria will be developed based on the analysis.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

Acknowledgments

The authors sincerely thank the anonymous reviewers for their constructive and helpful comments and suggestions. This work is supported by the National Natural Science Foundation of China under Grant no. 61472375, the 13th Five-year Pre-research Project of Civil Aerospace in China, Joint Funds of Equipment Pre-Research, and Ministry of Education of China under Grant no. 6141A02022320, and Fundamental Research Funds for the Central Universities under Grant no. CUG160207 and no. CUG2017G01.

References

W. Hohmann, Die Erreichbarkeit der Himmelskorper (The Possibility of Reaching the Stars), Verlag R, R. Oldenbourg, Oldenburg, Munchen, 1925.

R. B. Barrar, “An analytical proof that the Hohmann-type transfer is the true optimal two-impulse maneuver,” Astronautica Acta, vol. 9, no. 1, pp. 1–11, 1963.

K. Howell, B. Barden, and M. Lo, “Applications of dynamical systems theory to trajectory design for libration point missions,” Journal of the Astronautical Sciences, vol. 45, no. 2, pp. 161–178, 1997.

E. M. Alessi, G. Gómez, and J. J. Masdemont, “Leaving the Moon by means of invariant manifolds of libration point orbits,” Communications in Nonlinear Science and Numerical Simulation, vol. 14, no. 12, pp. 4153–4167, 2009.

W. S. Koon, M. W. Lo, J. E. Marsden, and S. D. Ross, “Low energy transfer to the Moon,” Celestial Mechanics and Dynamical Astronomy, vol. 81, no. 1/2, pp. 63–73, 2001.

M. W. Lo, R. L. Anderson, G. Whiffen, and L. Romans, “The role of invariant manifolds in low thrust trajectory design,” in Advances in the Astronautical Sciences 119, AAS, pp. 2971–2990, Springfield, VA, USA, 2004.

E. A. Belbruno, “Lunar capture orbits, a method of constructing Earth Moon trajectories and the lunar GAS mission,” in 19th International Electric Propulsion Conference, Colorado Springs,CO,USA, May 1987.

E. A. Belbruno and J. K. Miller, “Sun-perturbed Earth-to-Moon transfers with ballistic capture,” Journal of Guidance, Control, and Dynamics, vol. 16, no. 4, pp. 770–775, 1993.

E. A. Belbruno and J. K. Miller, “A ballistic lunar capture trajectory for the Japanese spacecraft hiten,” in Internal document JPL IOM 312/90.4-1731, Jet Propulsion Laboratory, Pasadena, CA, USA, 1990.

E. A. Belbruno and J. Carrico, “Calculation of weak stability boundary ballistic lunar transfer trajectories,” in Astrodynamics Specialist Conference, Denver, CO, USA, August 2000.

M. Belló, F. Graziani, P. Teofilatto, C. Circi, M. Porfilio, and M. Hechler, “A systematic analysis on weak stability boundary transfers to the Moon,” in IAF-00-A. 6.03 51st International Astronautical Congress, Rio de Janeiro, Brazil, 2000.

E. Belbruno, “Analytic estimation of weak stability boundaries and low energy transfers,” Contemporary Mathematics, vol. 292, pp. 17–45, 2002.

E. Belbruno, M. Gidea, and F. Topputo, “Weak stability boundary and invariant manifolds,” SIAM Journal on Applied Dynamical Systems, vol. 9, no. 3, pp. 1061–1089, 2010.

E. A. Belbruno, Capture Dynamics and Chaotic Motions in Celestial Mechanics: with Applications to the Construction of Low Energy Transfers, Princeton University Press, Princeton, NJ, USA, 2004.

P. A. Sousa Silva and M. O. Terra, “Diversity and validity of stable-unstable transitions in the algorithmic weak stability boundary,” Celestial Mechanics and Dynamical Astronomy, vol. 113, no. 4, pp. 453–478, 2012.

W. S. Koon, M. W. Lo, J. E. Marsden, and S. D. Ross, “Shoot the Moon,” in Spaceflight Mechanics 2000. Advances in Astronautical Sciences, pp. 1017–1030, San Diego, CA, USA, 2000.

F. Topputo, M. Vasile, and F. Bernelli-Zazzera, “Earth-to-Moon low energy transfers targeting L₁ hyperbolic transit orbits,” Annals of the New York Academy of Sciences, vol. 1065, no. 1, pp. 55–76, 2005.

G. Gómez, W. S. Koon, M. W. Lo, J. E. Marsden, J. Masdemont, and S. D. Ross, “Connecting orbits and invariant manifolds in the spatial restricted three-body problem,” Nonlinearity, vol. 17, no. 5, pp. 1571–1606, 2004.

P. A. de Sousa-Silva and M. O. Terra, “A survey of different classes of Earth-to-Moon trajectories in the patched three-body approach,” Acta Astronautica, vol. 123, pp. 340–349, 2016.

J. T. Betts, “Survey of numerical methods for trajectory optimization,” Journal of Guidance, Control, and Dynamics, vol. 21, no. 2, pp. 193–207, 1998.

H. Lei, B. Xu, and Y. Sun, “Earth–Moon low energy trajectory optimization in the real system,” Advances in Space Research, vol. 51, no. 5, pp. 917–929, 2013.

P. J. Gage, R. D. Braun, and I. M. Kroo, “Interplanetary trajectory optimization using a genetic algorithm,” Journal of the Astronautical Sciences, vol. 43, no. 1, pp. 59–76, 1995.

G. A. Rauwolf and A. Friedlander, “Near-optimal solar sail trajectories generated by a genetic algorithm,” Astrodynamics, vol. 1999, pp. 489–505, 2000.

G. W. Hughes and C. R. McInnes, “Solar sail hybrid trajectory optimization for non-Keplerian orbit transfers,” Journal of Guidance, Control, and Dynamics, vol. 25, no. 3, pp. 602–604, 2002.

N. Yokoyama and S. Suzuki, “Modified genetic algorithm for constrained trajectory optimization,” Journal of Guidance, Control, and Dynamics, vol. 28, no. 1, pp. 139–144, 2005.

F. Topputo, M. Vasile, and F. B. Zazzera, “A hybrid optimization of the low energy interplanetary transfers associated to the invariant manifolds of the restricted three-body problem,” in 55th International Astronautical Congress, pp. 1–12, Vancouver, Canada, October 2004.

A. Ohndorf, B. Dachwald, and E. Gill, “Optimization of low-thrust Earth-Moon transfers using evolutionary neurocontrol,” in 2009. CEC’09. IEEE Congress on Evolutionary Computation, pp. 358–364, Trondheim, Norway, May 2009.

M. Pontani and B. A. Conway, “Particle swarm optimization applied to space trajectories,” Journal of Guidance, Control, and Dynamics, vol. 33, no. 5, pp. 1429–1441, 2010.

M. Pontani and B. A. Conway, “Particle swarm optimization applied to impulsive orbital transfers,” Acta Astronautica, vol. 74, pp. 141–155, 2012.

F. A. Zotes and M. S. Peñas, “Particle swarm optimisation of interplanetary trajectories from Earth to Jupiter and Saturn,” Engineering Applications of Artificial Intelligence, vol. 25, no. 1, pp. 189–199, 2012.

Y.-W. Leung and Y. Wang, “An orthogonal genetic algorithm with quantization for global numerical optimization,” IEEE Transactions on Evolutionary Computation, vol. 5, no. 1, pp. 41–53, 2001.

R. Storn and K. Price, “Differential evolution: a simple and efficient adaptive scheme for global optimization over continuous spaces,” Journal of Global Optimization, vol. 11, no. 4, pp. 341–359, 1997.

S. Das and P. N. Suganthan, “Differential evolution: a survey of the state-of-the-art,” IEEE Transactions on Evolutionary Computation, vol. 15, no. 1, pp. 4–31, 2011.

S. Das, S. S. Mullick, and P. N. Suganthan, “Recent advances in differential evolution – an updated survey,” Swarm and Evolutionary Computation, vol. 27, pp. 1–30, 2016.

A. F. Attia, “Precise determination of the orbital elements of binary stars using differential evolution,” Astrophysics and Space Science, vol. 334, no. 1, pp. 103–114, 2011.

Y. Lu, H. Li, J. Li et al., “Design and optimization of low-energy transfer orbit to Mars with multi-body environment,” Science China Technological Sciences, vol. 58, no. 10, pp. 1660–1671, 2015.

P. Nath and R. V. Ramanan, “Precise halo orbit design and optimal transfer to halo orbits from Earth using differential evolution,” Advances in Space Research, vol. 57, no. 1, pp. 202–217, 2016.

H. Shang, X. Wu, Y. Ren, and J. Shan, “An efficient algorithm for global periodic orbits generation near irregular-shaped asteroids,” Communications in Nonlinear Science and Numerical Simulation, vol. 48, pp. 550–568, 2017.

V. Szebehely, Theory of Orbits, Academic Press, New York, NY, USA, 1967.

K. V. Price, “An introduction to differential evolution,” in New Ideas in Optimization, D. Corne, M. Dorigo, and V. Glover, Eds., pp. 79–108, McGraw-Hill, London, U.K., 1999.

K. V. Price, R. M. Storn, and J. A. Lampinen, Differential Evolution—A Practical Approach to Global Optimization, Springer, Berlin, Germany, 2005.

S. P. Han, “Superlinearly convergent variable metric algorithms for general nonlinear programming problems,” Mathematical Programming, vol. 11, no. 1, pp. 263–282, 1976.

S. P. Han, “A globally convergent method for nonlinear programming,” Journal of Optimization Theory and Applications, vol. 22, no. 3, pp. 297–309, 1977.

M. J. Powell, “The convergence of variable metric methods for nonlinearly constrained optimization calculations,” in Nonlinear Programming 3, O. L. Mangasarian, R. R. Meyer, and S. M. Robinson, Eds., pp. 27–63, Academic Press, Madison, WI, USA, 1978.

R. B. Wilson, A Simplicial Algorithm for Concave Programming, Ph. D. Dissertation, Graduate School of Bussiness Administration, Boston, MA, USA, 1963.

K. Schittkowski, “NLPQL: a FORTRAN subroutine solving constrained nonlinear programming problems,” Annals of Operations Research, vol. 5, no. 1-4, pp. 485–500, 1986.

P. T. Boggs and J. W. Tolle, “Sequential quadratic programming for large-scale nonlinear optimization,” Journal of Computational and Applied Mathematics, vol. 124, no. 1-2, pp. 123–137, 2000.

R. Fletcher, Practical Methods of Optimization, John Wiley and Sons, New York, NY, USA, 1987.

P. E. Gill, W. Murray, and M. H. Wright, Practical Optimization, Academic Press, London, 1981.

M. J. Powell, “Variable metric methods for constrained optimization,” in Mathematical Programming the State of the Art, pp. 288–311, Springer, Berlin, Heidelberg, 1983.

W. Hock and K. Schittkowski, “A comparative performance evaluation of 27 nonlinear programming codes,” Computing, vol. 30, no. 4, pp. 335–358, 1983.

N. Hansen and A. Ostermeier, “Adapting arbitrary normal mutation distributions in evolution strategies: the covariance matrix adaptation,” in 1996 Proceedings of IEEE International Conference on Evolutionary Computation, pp. 312–317, Nagoya, Japan, May 1996.

N. Hansen and A. Ostermeier, “Completely derandomized self-adaptation in evolution strategies,” Evolutionary Computation, vol. 9, no. 2, pp. 159–195, 2001.

N. Hansen, S. D. Muller, and P. Koumoutsakos, “Reducing the time complexity of the derandomized evolution strategy with covariance matrix adaptation (CMA-ES),” Evolutionary Computation, vol. 11, no. 1, pp. 1–18, 2003.

D. Simon, M. G. H. Omran, and M. Clerc, “Linearized biogeography-based optimization with re-initialization and local search,” Information Sciences, vol. 267, pp. 140–157, 2014.

C. García-Martínez, M. Lozano, F. Herrera, D. Molina, and A. M. Sánchez, “Global and local real-coded genetic algorithms based on parent-centric crossover operators,” European Journal of Operational Research, vol. 185, no. 3, pp. 1088–1113, 2008.

A. K. Qin, V. L. Huang, and P. N. Suganthan, “Differential evolution algorithm with strategy adaptation for global numerical optimization,” IEEE Transactions on Evolutionary Computation, vol. 13, no. 2, pp. 398–417, 2009.

R. Tanabe and A. Fukunaga, “Success-history based parameter adaptation for differential evolution,” in 2013 IEEE Congress on Evolutionary Computation, pp. 71–78, Cancun, Mexico, June 2013.

G. Wu, R. Mallipeddi, P. N. Suganthan, R. Wang, and H. Chen, “Differential evolution with multi-population based ensemble of mutation strategies,” Information Sciences, vol. 329, pp. 329–345, 2016.

M. Asafuddoula, T. Ray, and R. Sarker, “An adaptive hybrid differential evolution algorithm for single objective optimization,” Applied Mathematics and Computation, vol. 231, no. 1, pp. 601–618, 2014.

W. Gong, Z. Cai, C. X. Ling, and H. Li, “Enhanced differential evolution with adaptive strategies for numerical optimization,” IEEE Transactions on Systems, Man, and Cybernetics, Part B (Cybernetics), vol. 41, no. 2, pp. 397–413, 2011.

P. N. Suganthan, N. Hansen, J. J. Liang et al., “Problem definitions and evaluation criteria for the CEC 2005 special session on real-parameter optimization,” Tech. Rep., Tech. Rep., Nanyang Technological University, Singapore, 2005.

S. García, D. Molina, M. Lozano, and F. Herrera, “A study on the use of non-parametric tests for analyzing the evolutionary algorithms behaviour: a case study on the CEC 2005 special session on real parameter optimization,” Journal of Heuristics, vol. 15, no. 6, pp. 617–644, 2009.

J. Alcalá-Fdez, L. Sánchez, S. García et al., “KEEL: a software tool to assess evolutionary algorithms for data mining problems,” Soft Computing, vol. 13, no. 3, pp. 307–318, 2009.