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

Authors: Yanyun Zhang; Lei Peng; Guangming Dai; Maocai Wang
Journal: International Journal of Aerospace Engineering
 (2018)

---

## 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.


---

## Body

## 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 Section2, 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 primariesP1 and P2 [39]. Choose a rotating coordinate system so that the origin is at the center of mass between P1 and P2. See Figure 1.Figure 1
The planar CR3BP in the rotating frame.Letx,y 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.
(1)x¨−2ẏ=Ωx,y¨+2ẋ=Ωy,where x¨ and y¨ are the second derivative of x and y against time t, respectively, while ẋ and ẏ denote the derivative of x and y against time t, and Ωx and Ωy represent the partial derivatives of the effective potential Ω against x and y, respectively. The effective potential is given by
(2)Ωx,y=x2+y22+1−μr1+μr2+μ1−μ2 and
(3)μ=m2m1+m2,m1>m2,r1=x+μ2+y2,r2=x−1+μ2+y2,where m1 and m2 are the masses of P1 and P2, respectively, while r1 and r2 are the distances from the spacecraft to the two celestial bodies P1 and P2, respectively. With normalization, the normalized variables, including the distance between r1 and r2, the sum of their masses, and their angular velocity around the barycenter, are normalized to one. We assume that the mass of P1 is greater than the mass of P2, namely, m1>m2, then μ, known as the mass parameter of the planar CR3BP, is the dimensionless mass of P2, while the dimensionless mass of P1 is 1−μ. In the rotating frame, P1 and P2 are fixed at μ,0 and 1−μ,0, 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:(4)x¨−2ẏ=Ωx,y¨+2ẋ=Ωy,z¨=Ωz.
(2)
The effective potential in three-dimensional case:(5)Ωx,y,z=x2+y22+1−μr1+μr2+μ1−μ2
and(6)μ=m2m1+m2,m1>m2,r1=x+μ2+y2+z2,r2=x−1+μ2+y2+z2.


### 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 spacex,y,ẋ,ẏ. 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.Figure2 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 Earthx=1−μ is usually selected as the location of the Poincaré section. This also means the x-axis position of the four-dimensional phase space is fixed. With the calculation of the y−ẏ curve in both the Sun-Earth system and the Earth-Moon system in Poincaré section, the initial y and ẏ 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 x-axis direction ẋ can also be determined. Thus, the initial patched point is obtained.Given both the location and speed in thex-axis direction is determined, the decision variable of the optimization model is
(7)X=y,ẏ,θ,where y and ẏ are the position and velocity of the spacecraft in the y-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:(8)X=y,ẏ,z,ż,θ.From the above discussion and the illustration of Figure2, three impulsive vectors may be required in total during the process of low-energy transfer from the Earth to the Moon.
(1)
Impulsive vectorΔV1 is required to put the spacecraft from the Earth parking orbit into the Sun-Earth stable manifold.
(2)
Impulsive vectorΔV2 is required to connect the two segments of the transfer trajectory in the Poincaré section.
(3)
Impulsive vectorΔV3 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:(9)MinimizefX=ΔV1+ΔV2+ΔV3.


## 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 primariesP1 and P2 [39]. Choose a rotating coordinate system so that the origin is at the center of mass between P1 and P2. See Figure 1.Figure 1
The planar CR3BP in the rotating frame.Letx,y 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.
(1)x¨−2ẏ=Ωx,y¨+2ẋ=Ωy,where x¨ and y¨ are the second derivative of x and y against time t, respectively, while ẋ and ẏ denote the derivative of x and y against time t, and Ωx and Ωy represent the partial derivatives of the effective potential Ω against x and y, respectively. The effective potential is given by
(2)Ωx,y=x2+y22+1−μr1+μr2+μ1−μ2 and
(3)μ=m2m1+m2,m1>m2,r1=x+μ2+y2,r2=x−1+μ2+y2,where m1 and m2 are the masses of P1 and P2, respectively, while r1 and r2 are the distances from the spacecraft to the two celestial bodies P1 and P2, respectively. With normalization, the normalized variables, including the distance between r1 and r2, the sum of their masses, and their angular velocity around the barycenter, are normalized to one. We assume that the mass of P1 is greater than the mass of P2, namely, m1>m2, then μ, known as the mass parameter of the planar CR3BP, is the dimensionless mass of P2, while the dimensionless mass of P1 is 1−μ. In the rotating frame, P1 and P2 are fixed at μ,0 and 1−μ,0, 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:(4)x¨−2ẏ=Ωx,y¨+2ẋ=Ωy,z¨=Ωz.
(2)
The effective potential in three-dimensional case:(5)Ωx,y,z=x2+y22+1−μr1+μr2+μ1−μ2
and(6)μ=m2m1+m2,m1>m2,r1=x+μ2+y2+z2,r2=x−1+μ2+y2+z2.


## 2.1.1. Formula in Two-Dimensional Space

The planar CR3BP describes the motion of a spacecraft under the gravitational attractions of two primariesP1 and P2 [39]. Choose a rotating coordinate system so that the origin is at the center of mass between P1 and P2. See Figure 1.Figure 1
The planar CR3BP in the rotating frame.Letx,y 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.
(1)x¨−2ẏ=Ωx,y¨+2ẋ=Ωy,where x¨ and y¨ are the second derivative of x and y against time t, respectively, while ẋ and ẏ denote the derivative of x and y against time t, and Ωx and Ωy represent the partial derivatives of the effective potential Ω against x and y, respectively. The effective potential is given by
(2)Ωx,y=x2+y22+1−μr1+μr2+μ1−μ2 and
(3)μ=m2m1+m2,m1>m2,r1=x+μ2+y2,r2=x−1+μ2+y2,where m1 and m2 are the masses of P1 and P2, respectively, while r1 and r2 are the distances from the spacecraft to the two celestial bodies P1 and P2, respectively. With normalization, the normalized variables, including the distance between r1 and r2, the sum of their masses, and their angular velocity around the barycenter, are normalized to one. We assume that the mass of P1 is greater than the mass of P2, namely, m1>m2, then μ, known as the mass parameter of the planar CR3BP, is the dimensionless mass of P2, while the dimensionless mass of P1 is 1−μ. In the rotating frame, P1 and P2 are fixed at μ,0 and 1−μ,0, 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:(4)x¨−2ẏ=Ωx,y¨+2ẋ=Ωy,z¨=Ωz.
(2)
The effective potential in three-dimensional case:(5)Ωx,y,z=x2+y22+1−μr1+μr2+μ1−μ2
and(6)μ=m2m1+m2,m1>m2,r1=x+μ2+y2+z2,r2=x−1+μ2+y2+z2.


## 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 spacex,y,ẋ,ẏ. 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.Figure2 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 Earthx=1−μ is usually selected as the location of the Poincaré section. This also means the x-axis position of the four-dimensional phase space is fixed. With the calculation of the y−ẏ curve in both the Sun-Earth system and the Earth-Moon system in Poincaré section, the initial y and ẏ 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 x-axis direction ẋ can also be determined. Thus, the initial patched point is obtained.Given both the location and speed in thex-axis direction is determined, the decision variable of the optimization model is
(7)X=y,ẏ,θ,where y and ẏ are the position and velocity of the spacecraft in the y-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:(8)X=y,ẏ,z,ż,θ.From the above discussion and the illustration of Figure2, three impulsive vectors may be required in total during the process of low-energy transfer from the Earth to the Moon.
(1)
Impulsive vectorΔV1 is required to put the spacecraft from the Earth parking orbit into the Sun-Earth stable manifold.
(2)
Impulsive vectorΔV2 is required to connect the two segments of the transfer trajectory in the Poincaré section.
(3)
Impulsive vectorΔV3 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:(9)MinimizefX=ΔV1+ΔV2+ΔV3.


## 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 vectorX.
(10)X=x1,x2,x3,…,xDwhere D 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 N individuals form a so-called population P. To evaluate the performance of the system under a specific set of parameters, an evaluate function (or object function, fitness function) fX is defined. So, the optimization aims at finding out the object vector X∗, which minimizes (or maximizes) such an objective function fX, that is, fX∗<fX for all X and the X∗ is the best solution.DE works on a population ofN 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 DE/a/b/c, where a stands for the type of base vector in the mutation strategy, b specifies the number of difference vectors in the mutation strategy, and c 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 DE/rand/1/exp 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.Algorithm 1: DE/rand/1/exp.
1: Generate the initial populationP, define Xi as the i-th individual in P, N is the  population size, NFEs is the number of function evaluations in each run, MaxNFEs is the number of max function evaluation, D is the number of decision variable, F is the mutation factor, Cr is crossover rate.
2:NFEs=0
3: Evaluate the fitnessfXi for the each individual in P.
4:NFEs=NFEs+N
5:whileNFEs ≤ MaxNFEsdo
6:fori=1 to Ndo
7:   SelectXr1,Xr2,Xr3 and r1≠r2≠r3≠i
8:j=randint1,D
9:L=0
10:Ui=Xi
11:repeat
12:Ui,j=Xr1,j+F∗Xr2,j−Xr3,j
13:j=j+1%D
14:L=L+1
15:untilrandreal0,1>CrorL>D
16:   Evaluate the trial vectorUi
17:NFEs=NFEs+1
18:ifUi is better than Xithen
19:Xi=Ui
20:end if
21:end for
22:end while


### 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(11)ρ1=∑i=1NfiX−favgX2N−11/2,where fiX is the objective value of the ith individual among the current population, while favgX refers to the average value for the objective values of the current populations, and N is the population size.The other is the maximum distance in decision space:(12)ρ2=maxXi−Xj,∀Xi,Xj∈P,where Xi is the ith individual of the current population P, while Xj is the jth individual.ρ1 is a measure of the degree of contraction in objective space while ρ2 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 ρ2, 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 ρ1,max,ρ2,max, local search is started, shown in Algorithm 2.Algorithm 2: DEMR.
1: Initialize populationP, define Xi as the i-th vector in P at generation g, and Xbest is the best vector in the current generation, while Xoptimal is the best solution we have found so far. N is the population size, NFEs refers to the number of function evaluations, MaxNFEs is the number of max function evaluation. F denotes the mutation factor and Cr is crossover rate. ρ1 and ρ2 are used for contraction criterion.
Xlocal refers to the resultant new local optimum found by SQP. Xmin and Xmax are the lower and upper boundaries of the variables. Archive denotes a set which collects the superior solutions at the late stage of evolution.
2:NFEs=0, R=0, Archive=∅
3: Evaluate the fitnessfXi for the each individual in P.
4:NFEs=NFEs+N
5:whileNFEs < MaxNFEsdo
6:fori=1 to Ndo
7:      SelectXr1,Xr2,Xr3 with the select probabilities and r1≠r2≠r3≠i.
8:ifNFEs < 0.2 × MaxNFEsthen
9:     Using DE/rand/1/bin to generateUi
10:else
11:     Using DE/best/1/exp to generateUi
12:end if
13:      Evaluate the trial vectorUi
14:NFEs=NFEs+1
15:ifUi is better than Xithen
16:Xi=Ui
17:end if
18:      UpdateXbest and Xoptimal
19:end for
20:   Calculate the contraction criteria
21:ifρ1<ρ1,maxorρ2<ρ2,maxthen
22:      Pick up theXbest as the initial point of the local search SQP, and find the local optimum Xlocal
23:      UpdateXbest and Xoptimal, Xoptimal→Archive
24:ifXlocal>Xbestthen
25:     Re-Initialize the whole population as Algorithm3 shows.
26:end if
27:end if
28:end whileAlgorithm 3: Mix Reinitialization.
1:R is the number of re-initialization while Rmax is the maximum of the number of randomly re-initialization. Archive denotes a set which collects the superior solutions at the late stage of evolution.
2:ifR<Rmaxthen
3:  Re-initialize the whole population randomly
4:else
5:  Calculate the average valueXavg and standard deviation Xstd of the elements in Archive
6:  Re-initialize the population as:Xi=Xavg+rand0,1Xstd,i=1,…,N
7:end if


#### 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 Algorithm3, is applied. Firstly, to guarantee the diversity of the population, the population is reinitialized randomly in the first R times, and thus the algorithm can distribute individuals in the whole search space uniformly to a large extent; secondly, after R times’ reinitialization, the population is reinitialized as follows:
(13)Xi=Xavg+rand0,1Xstd,i=1,…,N,where Xi denotes the ith individual in the population, and N is the population size. Xavg and Xstd stand for the average value and standard deviation of the elements in Archive, respectively, while Archive denotes the set which collects the superior solutions at the late stage of evolution. Once the local search is triggered, the extra storage Archive 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 Archive. 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 Archive. At this time, reinitializing the population based on Archive cannot only find the optimal value with great probability but also speed up the convergence rate. And the size of Archive will not be restricted here since the Archive is active only at the late stage of evolution, and the storage capacity will not be great.


## 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 vectorX.
(10)X=x1,x2,x3,…,xDwhere D 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 N individuals form a so-called population P. To evaluate the performance of the system under a specific set of parameters, an evaluate function (or object function, fitness function) fX is defined. So, the optimization aims at finding out the object vector X∗, which minimizes (or maximizes) such an objective function fX, that is, fX∗<fX for all X and the X∗ is the best solution.DE works on a population ofN 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 DE/a/b/c, where a stands for the type of base vector in the mutation strategy, b specifies the number of difference vectors in the mutation strategy, and c 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 DE/rand/1/exp 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.Algorithm 1: DE/rand/1/exp.
1: Generate the initial populationP, define Xi as the i-th individual in P, N is the  population size, NFEs is the number of function evaluations in each run, MaxNFEs is the number of max function evaluation, D is the number of decision variable, F is the mutation factor, Cr is crossover rate.
2:NFEs=0
3: Evaluate the fitnessfXi for the each individual in P.
4:NFEs=NFEs+N
5:whileNFEs ≤ MaxNFEsdo
6:fori=1 to Ndo
7:   SelectXr1,Xr2,Xr3 and r1≠r2≠r3≠i
8:j=randint1,D
9:L=0
10:Ui=Xi
11:repeat
12:Ui,j=Xr1,j+F∗Xr2,j−Xr3,j
13:j=j+1%D
14:L=L+1
15:untilrandreal0,1>CrorL>D
16:   Evaluate the trial vectorUi
17:NFEs=NFEs+1
18:ifUi is better than Xithen
19:Xi=Ui
20:end if
21:end for
22:end while


## 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(11)ρ1=∑i=1NfiX−favgX2N−11/2,where fiX is the objective value of the ith individual among the current population, while favgX refers to the average value for the objective values of the current populations, and N is the population size.The other is the maximum distance in decision space:(12)ρ2=maxXi−Xj,∀Xi,Xj∈P,where Xi is the ith individual of the current population P, while Xj is the jth individual.ρ1 is a measure of the degree of contraction in objective space while ρ2 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 ρ2, 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 ρ1,max,ρ2,max, local search is started, shown in Algorithm 2.Algorithm 2: DEMR.
1: Initialize populationP, define Xi as the i-th vector in P at generation g, and Xbest is the best vector in the current generation, while Xoptimal is the best solution we have found so far. N is the population size, NFEs refers to the number of function evaluations, MaxNFEs is the number of max function evaluation. F denotes the mutation factor and Cr is crossover rate. ρ1 and ρ2 are used for contraction criterion.
Xlocal refers to the resultant new local optimum found by SQP. Xmin and Xmax are the lower and upper boundaries of the variables. Archive denotes a set which collects the superior solutions at the late stage of evolution.
2:NFEs=0, R=0, Archive=∅
3: Evaluate the fitnessfXi for the each individual in P.
4:NFEs=NFEs+N
5:whileNFEs < MaxNFEsdo
6:fori=1 to Ndo
7:      SelectXr1,Xr2,Xr3 with the select probabilities and r1≠r2≠r3≠i.
8:ifNFEs < 0.2 × MaxNFEsthen
9:     Using DE/rand/1/bin to generateUi
10:else
11:     Using DE/best/1/exp to generateUi
12:end if
13:      Evaluate the trial vectorUi
14:NFEs=NFEs+1
15:ifUi is better than Xithen
16:Xi=Ui
17:end if
18:      UpdateXbest and Xoptimal
19:end for
20:   Calculate the contraction criteria
21:ifρ1<ρ1,maxorρ2<ρ2,maxthen
22:      Pick up theXbest as the initial point of the local search SQP, and find the local optimum Xlocal
23:      UpdateXbest and Xoptimal, Xoptimal→Archive
24:ifXlocal>Xbestthen
25:     Re-Initialize the whole population as Algorithm3 shows.
26:end if
27:end if
28:end whileAlgorithm 3: Mix Reinitialization.
1:R is the number of re-initialization while Rmax is the maximum of the number of randomly re-initialization. Archive denotes a set which collects the superior solutions at the late stage of evolution.
2:ifR<Rmaxthen
3:  Re-initialize the whole population randomly
4:else
5:  Calculate the average valueXavg and standard deviation Xstd of the elements in Archive
6:  Re-initialize the population as:Xi=Xavg+rand0,1Xstd,i=1,…,N
7:end if


### 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 Algorithm3, is applied. Firstly, to guarantee the diversity of the population, the population is reinitialized randomly in the first R times, and thus the algorithm can distribute individuals in the whole search space uniformly to a large extent; secondly, after R times’ reinitialization, the population is reinitialized as follows:
(13)Xi=Xavg+rand0,1Xstd,i=1,…,N,where Xi denotes the ith individual in the population, and N is the population size. Xavg and Xstd stand for the average value and standard deviation of the elements in Archive, respectively, while Archive denotes the set which collects the superior solutions at the late stage of evolution. Once the local search is triggered, the extra storage Archive 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 Archive. 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 Archive. At this time, reinitializing the population based on Archive cannot only find the optimal value with great probability but also speed up the convergence rate. And the size of Archive will not be restricted here since the Archive is active only at the late stage of evolution, and the storage capacity will not be great.


## 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(11)ρ1=∑i=1NfiX−favgX2N−11/2,where fiX is the objective value of the ith individual among the current population, while favgX refers to the average value for the objective values of the current populations, and N is the population size.The other is the maximum distance in decision space:(12)ρ2=maxXi−Xj,∀Xi,Xj∈P,where Xi is the ith individual of the current population P, while Xj is the jth individual.ρ1 is a measure of the degree of contraction in objective space while ρ2 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 ρ2, 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 ρ1,max,ρ2,max, local search is started, shown in Algorithm 2.Algorithm 2: DEMR.
1: Initialize populationP, define Xi as the i-th vector in P at generation g, and Xbest is the best vector in the current generation, while Xoptimal is the best solution we have found so far. N is the population size, NFEs refers to the number of function evaluations, MaxNFEs is the number of max function evaluation. F denotes the mutation factor and Cr is crossover rate. ρ1 and ρ2 are used for contraction criterion.
Xlocal refers to the resultant new local optimum found by SQP. Xmin and Xmax are the lower and upper boundaries of the variables. Archive denotes a set which collects the superior solutions at the late stage of evolution.
2:NFEs=0, R=0, Archive=∅
3: Evaluate the fitnessfXi for the each individual in P.
4:NFEs=NFEs+N
5:whileNFEs < MaxNFEsdo
6:fori=1 to Ndo
7:      SelectXr1,Xr2,Xr3 with the select probabilities and r1≠r2≠r3≠i.
8:ifNFEs < 0.2 × MaxNFEsthen
9:     Using DE/rand/1/bin to generateUi
10:else
11:     Using DE/best/1/exp to generateUi
12:end if
13:      Evaluate the trial vectorUi
14:NFEs=NFEs+1
15:ifUi is better than Xithen
16:Xi=Ui
17:end if
18:      UpdateXbest and Xoptimal
19:end for
20:   Calculate the contraction criteria
21:ifρ1<ρ1,maxorρ2<ρ2,maxthen
22:      Pick up theXbest as the initial point of the local search SQP, and find the local optimum Xlocal
23:      UpdateXbest and Xoptimal, Xoptimal→Archive
24:ifXlocal>Xbestthen
25:     Re-Initialize the whole population as Algorithm3 shows.
26:end if
27:end if
28:end whileAlgorithm 3: Mix Reinitialization.
1:R is the number of re-initialization while Rmax is the maximum of the number of randomly re-initialization. Archive denotes a set which collects the superior solutions at the late stage of evolution.
2:ifR<Rmaxthen
3:  Re-initialize the whole population randomly
4:else
5:  Calculate the average valueXavg and standard deviation Xstd of the elements in Archive
6:  Re-initialize the population as:Xi=Xavg+rand0,1Xstd,i=1,…,N
7:end if


## 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 Algorithm3, is applied. Firstly, to guarantee the diversity of the population, the population is reinitialized randomly in the first R times, and thus the algorithm can distribute individuals in the whole search space uniformly to a large extent; secondly, after R times’ reinitialization, the population is reinitialized as follows:
(13)Xi=Xavg+rand0,1Xstd,i=1,…,N,where Xi denotes the ith individual in the population, and N is the population size. Xavg and Xstd stand for the average value and standard deviation of the elements in Archive, respectively, while Archive denotes the set which collects the superior solutions at the late stage of evolution. Once the local search is triggered, the extra storage Archive 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 Archive. 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 Archive. At this time, reinitializing the population based on Archive cannot only find the optimal value with great probability but also speed up the convergence rate. And the size of Archive will not be restricted here since the Archive 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 (Table1).Table 1
Parameters in DEMR
Value
Crossover rate
Cr=N0.9,0.1
Scale factor
F=N0.5,0.1
Dimension
D=10
Population size
N=30
Contraction criteria
ρ1,max=3.0,ρ2,max=1.0
Reinitialization control
Rmax=5
Maximum evaluations
MaxNFEs = 10,000Storn and Price [32] have indicated that the performance of DE is sensitive to the value of control parameters, and good choices of F and Cr 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: F=N0.5,0.1 and Cr=N0.9,0.1, where Nμ,σ is a normal distribution. And according to the study of Storn and Price [32], a reasonable value for the population size N could be chosen between 5D and 10D; here, N is set as 30. For the setting of the contraction criteria, ρ1,max=3.0 and ρ2,max=1.0 are suggested, and the maximum number of randomly reinitialization Rmax is set as 5. When the maximum number of function evaluations (MaxNFEs), set to 10,000, is reached, the algorithm will be stopped.According to the analysis in [60, 61], DE/rand/1 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 DE/best/1 belongs to the kind of greedy strategies and always converges faster; thus, it is more reasonable to adopt DE/best/1 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 DE/rand/1/bin and DE/best/1/exp is employed that the DE/rand/1/bin is applied when the number of evaluation NFEs is less than 0.2 × MaxNFEs, while DE/best/1/exp 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 (Table2).Table 2
Decision variables
Range
Position inY direction
0.0029≤y≤0.0065
Velocity inY direction
−0.022≤ẏ≤0.003
Conversion angle
60∘≤θ≤90∘Table3 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 ΔV=3908.3m/s for the Earth-Moon low-energy transfer, and the corresponding ΔV for worst value, median value, and mean value is 3912.8m/s, 3909.0m/s, and 3909.2m/s, 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
1The 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 Table3, 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ΔV found so far will be recorded. The record points are 0.01,0.02,0.03,0.05,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0∗MaxNFEs. 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 ΔV 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 Figure3, 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 ΔV. 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 NFEs∈0,10000,ΔV∈3900,4100. 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.Figure4 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 y,ẏ,θ=0.00305655131737,−0.00196665640171,76.39749140443999.Figure 4
The best Earth-Moon low-energy transfer trajectory optimized by DEMR.y,ẏ,θ=0.00305655131737,−0.00196665640171,76.39749140443999.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 (Table4):Table 4
Decision variables
Range
Position inY direction
0.001≤y≤0.008
Velocity inY direction
−0.03≤ẏ≤0.006
Conversion angle
0∘≤θ≤180∘The results with the widened ranges are shown in Table5. 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 thez-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 inY direction
0.001≤y≤0.008
Velocity inY direction
−0.03≤ẏ≤0.006
Position inZ direction
−0.0009≤z≤0.0009
Velocity inZ direction
−0.004≤ż≤0.004
Conversion angle
0∘≤θ≤180∘From Table7, 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
2According 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,ΔV=3990m/s, in [14], the threshold is set as 3990m/s. 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 D=10.The performance of the algorithms is evaluated in terms of function error value (FEV), defined asfX−fX∗, where X∗ is the global optimum of the test function, while X 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 Rmax, which are set as MaxNFEs = 10,000∗D and Rmax=20. 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 α=0.1 and α=0.05 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 Table9. 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 α=0.05. 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∗
1.57E−12±4.49E−15
3.12E−27±1.31E−27+
0.00E+00±0.00E+00+
2.30E−26±3.85E−26+
2.81E−07±2.19E−07−
2.19E+01±7.74E+00−
2.61E−02±2.02E−02−
F2
3.71E−11±2.37E−11
5.37E−27±1.86E−27+
0.00E+00±0.00E+00+
2.55E−25±1.26E−24+
1.62E−29±3.78E−29+
2.37E−15±4.82E−15+
0.00E+00±0.00E+00+
F3
3.44E−04±6.58E−05
3.41E−23±1.59E−23+
0.00E+00±0.00E+00+
5.76E+04±6.15E+04−
5.09E+03±6.33E+03−
3.06E−11±6.57E−11+
8.11E−26±2.62E−26+
F4
4.66E+00±2.96E+00
4.26E+04±9.37E+04−
8.03E+02±5.10E+02−
9.82E−05±3.17E−04+
1.06E−29±3.00E−29+
1.45E−13±2.47E−13+
0.00E+00±0.00E+00+
F5
4.57E−04±1.87E−04
4.69E−11±1.02E−11+
2.94E−03±4.74E−03=
1.54E+02±1.08E+02−
1.46E−13±7.28E−13+
3.19E−06±1.65E−06+
0.00E+00±0.00E+00+
F6
1.79E−09±9.18E−10
7.97E−01±1.63E+00=
1.59E−01±7.97E−01=
2.76E+00±2.00E+00−
1.40E+00±1.44E+00−
9.98E−02±1.16E−01−
0.00E+00±0.00E+00+
F7
1.15E+00±5.47E−01
1.39E−02±9.92E−03+
3.23E−01±2.58E−01+
1.36E−01±1.22E−01+
3.55E−02±2.69E−02+
2.06E−02±2.19E−02+
3.14E−03±5.84E−03+
F8
2.00E+01±4.31E−08
2.00E+01±2.05E−01=
2.00E+01±0.00E+00=
2.04E+01±6.12E−02−
2.03E+01±5.71E−02−
2.04E+01±6.03E−02−
2.02E+01±1.49E−01−
F9∗
4.46E+00±2.03E+00
1.26E+02±6.03E+01−
0.00E+00±0.00E+00+
9.65E+00±8.78E+00−
7.58E+00±1.44E+00−
2.57E+01±4.06E+00−
1.78E+01±3.65E+00−
F10
5.57E+00±2.57E+00
2.95E+01±4.76E+01−
1.56E+01±4.35E+00−
1.02E+01±9.72E+00=
6.21E+00±2.30E+00=
6.64E+00±2.50E+00=
3.25E+00±9.31E−01+
F11
5.33E+00±1.13E+00
2.30E+00±1.37E+00+
5.72E+00±1.54E+00=
9.21E+00±8.15E−01−
7.91E−01±6.03E−01+
5.73E+00±6.54E−01=
4.80E+00±5.41E−01+
F12
1.88E−10±1.46E−10
5.04E+03±9.51E+03−
6.30E+01±3.11E+02−
3.38E+03±3.57E+03−
6.36E+00±6.12E+00−
1.55E+00±4.54E+00−
1.64E+00±4.88E+00−
F13
3.32E−01±1.81E−01
1.19E+00±5.54E−01−
6.45E−01±1.87E−01−
1.18E+00±5.51E−01−
5.63E−01±9.31E−02−
5.25E−01±8.46E−02−
2.84E−01±4.59E−02=
F14
3.31E+00±3.96E−01
4.92E+00±7.83E−02−
3.34E+00±2.41E−01=
3.61E+00±4.01E−01−
2.98E+00±2.46E−01+
3.15E+00±2.49E−01=
2.54E+00±3.88E−01+
F15
1.06E+02±3.06E+01
5.63E+02±3.35E+02−
6.92E+01±1.62E+02=
3.00E+02±4.39E+01−
4.10E+01±1.10E+02+
3.64E+01±1.11E+02+
1.29E+02±1.77E+02=
F16
9.88E+01±2.21E+01
3.60E+02±4.11E+02−
1.41E+02±2.05E+01−
1.08E+02±1.87E+01=
9.70E+01±7.19E+00=
1.03E+02±5.90E+00=
9.49E+01±4.72E+00+
F17
1.22E+02±8.93E+00
4.59E+02±4.77E+02−
1.46E+02±1.55E+01−
1.03E+02±1.35E+01+
1.01E+02±6.33E+00+
1.19E+02±7.31E+00=
1.03E+02±3.83E+00+
F18
3.12E+02±2.59E+01
6.49E+02±3.11E+02−
9.30E+02±1.51E+02−
7.77E+02±1.34E+02−
7.49E+02±1.72E+02−
7.40E+02±1.66E+02−
6.29E+02±2.53E+02−
F19
3.41E+02±1.09E+02
6.86E+02±2.95E+02−
9.20E+02±1.29E+02−
8.17E+02±6.64E+01−
7.28E+02±1.93E+02−
6.40E+02±2.38E+02−
6.69E+02±2.36E+02−
F20
3.17E+02±4.09E+01
6.56E+02±2.83E+02−
9.05E+02±1.02E+02−
7.86E+02±1.33E+02−
7.87E+02±1.05E+02−
6.60E+02±2.29E+02−
6.22E+02±2.47E+02−
F21
3.69E+02±9.92E+01
8.63E+02±2.15E+02−
9.41E+02±3.02E+02−
6.04E+02±2.19E+02−
5.20E+02±2.16E+02−
4.60E+02±8.16E+01−
4.96E+02±1.79E+02−
F22
7.23E+02±1.24E+02
7.72E+02±2.49E+01=
6.84E+02±1.95E+02=
6.68E+02±1.03E+02=
7.21E+02±1.27E+02=
7.61E+02±4.63E+00=
7.47E+02±1.23E+01−
F23
5.64E+02±6.32E+01
1.07E+03±2.10E+02−
1.03E+03±2.54E+02−
5.76E+02±8.22E+01−
7.50E+02±1.88E+02−
6.18E+02±1.22E+02=
6.51E+02±1.23E+02=
F24
2.00E+02±9.22E−01
3.00E+02±1.80E+02=
2.00E+02±2.79E−03+
3.77E+02±3.69E+00−
2.00E+02±0.00E+00+
2.00E+02±0.00E+00+
2.00E+02±0.00E+00+
F25
4.01E+02±8.70E+01
4.42E+02±1.51E+02−
2.12E+02±6.00E+01+
3.79E+02±3.70E+00+
3.76E+02±2.99E+00+
3.79E+02±2.54E+01+
3.71E+02±2.65E+00+
+
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 NFEs=10000 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 α=0.05.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
2The Wilcoxon results of the multiple-problem statistical analysis are listed in Table11. The R+ in Table 11 is a measure of the case where DEMR outperforms the compared algorithm and the greater the R+ value is, the greater difference between DEMR and the compared algorithm is. Similarly, the R− is a measure of the case where DEMR is worse than the compared algorithm, and the size of R− 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 α=0.05 and α=0.1. For MPEDE, it performs worse than DEMR when α=0.1, while is comparable with DEMR at α=0.05. And for other algorithms, it shows they have similar optimization effectiveness with DEMR, since the P value goes beyond the α at both α=0.05 and α=0.1. 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−
P value
α=0.05
α=0.1
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 aref1,f5,f12,f13,f18,f21, and f23, 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 f1 and f5, DEMR converges a little early. What is more, it should be noted that, in the convergence Figure 5(a) for f1, the convergence of LBBO after about NFEs=1000 has not been shown, since LBBO has already found the optimal value 0 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 NFEs=50000 for the same reason.Figure 5
The convergence graphs on different types of CEC2005 benchmark functions.
(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) f23In 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 Figure6. It shows that DEMR ranks the second with Ranking=3.22, behind SHADE with Ranking=2.70. And SaDE follows DEMR with Ranking=3.56. 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ρ1,max, ρ2,max, and Rmax, of which ρ1,max and ρ2,max determine the time to trigger the local search algorithm SQP, and Rmax 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 Table12 and Figure 7, the Friedman rankings for different combinations of ρ1,max and ρ2,max are illustrated: ρ1,max∈1,1.5,2,2.5,3, ρ2,max∈1,1.5,2,2.5,3, where Rmax is a constant Rmax=5. From Table 12, the setting ρ1,max=3,ρ2,max=1 is the best combination, followed by combinations ρ1,max=1,ρ2,max=1 and ρ1,max=1,ρ2,max=2. Figure 7 is more intuitive for the performance of these settings. Generally speaking, the combination with both small ρ1,max and small ρ2,max is more reasonable. Table 13 gives a brief look at the test results on different Rmax with contraction criteria combination ρ1,max=3,ρ2,max=1. And it can be seen that, for CEC2005 benchmark functions, Rmax=20 is the most appropriate choice, while for the design of Earth-Moon low-energy transfer trajectory, the performance of DEMR ranks the first when Rmax=5.Table 12
Results obtained by the Friedman test for contraction criteria combinations on CEC2005.
Friedman
ρ1,max
1
1.5
2
2.5
3
ρ2,max
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.32Figure 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
Rmax
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
Rmax
AverageΔV (m/s)
5
3909.1
10
3909.8
15
3909.4
20
3909.3
25
3909.7


## 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 (Table1).Table 1
Parameters in DEMR
Value
Crossover rate
Cr=N0.9,0.1
Scale factor
F=N0.5,0.1
Dimension
D=10
Population size
N=30
Contraction criteria
ρ1,max=3.0,ρ2,max=1.0
Reinitialization control
Rmax=5
Maximum evaluations
MaxNFEs = 10,000Storn and Price [32] have indicated that the performance of DE is sensitive to the value of control parameters, and good choices of F and Cr 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: F=N0.5,0.1 and Cr=N0.9,0.1, where Nμ,σ is a normal distribution. And according to the study of Storn and Price [32], a reasonable value for the population size N could be chosen between 5D and 10D; here, N is set as 30. For the setting of the contraction criteria, ρ1,max=3.0 and ρ2,max=1.0 are suggested, and the maximum number of randomly reinitialization Rmax is set as 5. When the maximum number of function evaluations (MaxNFEs), set to 10,000, is reached, the algorithm will be stopped.According to the analysis in [60, 61], DE/rand/1 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 DE/best/1 belongs to the kind of greedy strategies and always converges faster; thus, it is more reasonable to adopt DE/best/1 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 DE/rand/1/bin and DE/best/1/exp is employed that the DE/rand/1/bin is applied when the number of evaluation NFEs is less than 0.2 × MaxNFEs, while DE/best/1/exp 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 (Table2).Table 2
Decision variables
Range
Position inY direction
0.0029≤y≤0.0065
Velocity inY direction
−0.022≤ẏ≤0.003
Conversion angle
60∘≤θ≤90∘Table3 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 ΔV=3908.3m/s for the Earth-Moon low-energy transfer, and the corresponding ΔV for worst value, median value, and mean value is 3912.8m/s, 3909.0m/s, and 3909.2m/s, 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
1The 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 Table3, 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ΔV found so far will be recorded. The record points are 0.01,0.02,0.03,0.05,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0∗MaxNFEs. 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 ΔV 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 Figure3, 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 ΔV. 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 NFEs∈0,10000,ΔV∈3900,4100. 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.Figure4 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 y,ẏ,θ=0.00305655131737,−0.00196665640171,76.39749140443999.Figure 4
The best Earth-Moon low-energy transfer trajectory optimized by DEMR.y,ẏ,θ=0.00305655131737,−0.00196665640171,76.39749140443999.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 (Table4):Table 4
Decision variables
Range
Position inY direction
0.001≤y≤0.008
Velocity inY direction
−0.03≤ẏ≤0.006
Conversion angle
0∘≤θ≤180∘The results with the widened ranges are shown in Table5. 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 thez-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 inY direction
0.001≤y≤0.008
Velocity inY direction
−0.03≤ẏ≤0.006
Position inZ direction
−0.0009≤z≤0.0009
Velocity inZ direction
−0.004≤ż≤0.004
Conversion angle
0∘≤θ≤180∘From Table7, 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
2According 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,ΔV=3990m/s, in [14], the threshold is set as 3990m/s. 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.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 (Table2).Table 2
Decision variables
Range
Position inY direction
0.0029≤y≤0.0065
Velocity inY direction
−0.022≤ẏ≤0.003
Conversion angle
60∘≤θ≤90∘Table3 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 ΔV=3908.3m/s for the Earth-Moon low-energy transfer, and the corresponding ΔV for worst value, median value, and mean value is 3912.8m/s, 3909.0m/s, and 3909.2m/s, 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
1The 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 Table3, 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ΔV found so far will be recorded. The record points are 0.01,0.02,0.03,0.05,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0∗MaxNFEs. 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 ΔV 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 Figure3, 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 ΔV. 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 NFEs∈0,10000,ΔV∈3900,4100. 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.Figure4 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 y,ẏ,θ=0.00305655131737,−0.00196665640171,76.39749140443999.Figure 4
The best Earth-Moon low-energy transfer trajectory optimized by DEMR.y,ẏ,θ=0.00305655131737,−0.00196665640171,76.39749140443999.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 (Table4):Table 4
Decision variables
Range
Position inY direction
0.001≤y≤0.008
Velocity inY direction
−0.03≤ẏ≤0.006
Conversion angle
0∘≤θ≤180∘The results with the widened ranges are shown in Table5. 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 thez-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 inY direction
0.001≤y≤0.008
Velocity inY direction
−0.03≤ẏ≤0.006
Position inZ direction
−0.0009≤z≤0.0009
Velocity inZ direction
−0.004≤ż≤0.004
Conversion angle
0∘≤θ≤180∘From Table7, 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
2According 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,ΔV=3990m/s, in [14], the threshold is set as 3990m/s. 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 D=10.The performance of the algorithms is evaluated in terms of function error value (FEV), defined asfX−fX∗, where X∗ is the global optimum of the test function, while X 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 Rmax, which are set as MaxNFEs = 10,000∗D and Rmax=20. 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 α=0.1 and α=0.05 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 Table9. 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 α=0.05. 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∗
1.57E−12±4.49E−15
3.12E−27±1.31E−27+
0.00E+00±0.00E+00+
2.30E−26±3.85E−26+
2.81E−07±2.19E−07−
2.19E+01±7.74E+00−
2.61E−02±2.02E−02−
F2
3.71E−11±2.37E−11
5.37E−27±1.86E−27+
0.00E+00±0.00E+00+
2.55E−25±1.26E−24+
1.62E−29±3.78E−29+
2.37E−15±4.82E−15+
0.00E+00±0.00E+00+
F3
3.44E−04±6.58E−05
3.41E−23±1.59E−23+
0.00E+00±0.00E+00+
5.76E+04±6.15E+04−
5.09E+03±6.33E+03−
3.06E−11±6.57E−11+
8.11E−26±2.62E−26+
F4
4.66E+00±2.96E+00
4.26E+04±9.37E+04−
8.03E+02±5.10E+02−
9.82E−05±3.17E−04+
1.06E−29±3.00E−29+
1.45E−13±2.47E−13+
0.00E+00±0.00E+00+
F5
4.57E−04±1.87E−04
4.69E−11±1.02E−11+
2.94E−03±4.74E−03=
1.54E+02±1.08E+02−
1.46E−13±7.28E−13+
3.19E−06±1.65E−06+
0.00E+00±0.00E+00+
F6
1.79E−09±9.18E−10
7.97E−01±1.63E+00=
1.59E−01±7.97E−01=
2.76E+00±2.00E+00−
1.40E+00±1.44E+00−
9.98E−02±1.16E−01−
0.00E+00±0.00E+00+
F7
1.15E+00±5.47E−01
1.39E−02±9.92E−03+
3.23E−01±2.58E−01+
1.36E−01±1.22E−01+
3.55E−02±2.69E−02+
2.06E−02±2.19E−02+
3.14E−03±5.84E−03+
F8
2.00E+01±4.31E−08
2.00E+01±2.05E−01=
2.00E+01±0.00E+00=
2.04E+01±6.12E−02−
2.03E+01±5.71E−02−
2.04E+01±6.03E−02−
2.02E+01±1.49E−01−
F9∗
4.46E+00±2.03E+00
1.26E+02±6.03E+01−
0.00E+00±0.00E+00+
9.65E+00±8.78E+00−
7.58E+00±1.44E+00−
2.57E+01±4.06E+00−
1.78E+01±3.65E+00−
F10
5.57E+00±2.57E+00
2.95E+01±4.76E+01−
1.56E+01±4.35E+00−
1.02E+01±9.72E+00=
6.21E+00±2.30E+00=
6.64E+00±2.50E+00=
3.25E+00±9.31E−01+
F11
5.33E+00±1.13E+00
2.30E+00±1.37E+00+
5.72E+00±1.54E+00=
9.21E+00±8.15E−01−
7.91E−01±6.03E−01+
5.73E+00±6.54E−01=
4.80E+00±5.41E−01+
F12
1.88E−10±1.46E−10
5.04E+03±9.51E+03−
6.30E+01±3.11E+02−
3.38E+03±3.57E+03−
6.36E+00±6.12E+00−
1.55E+00±4.54E+00−
1.64E+00±4.88E+00−
F13
3.32E−01±1.81E−01
1.19E+00±5.54E−01−
6.45E−01±1.87E−01−
1.18E+00±5.51E−01−
5.63E−01±9.31E−02−
5.25E−01±8.46E−02−
2.84E−01±4.59E−02=
F14
3.31E+00±3.96E−01
4.92E+00±7.83E−02−
3.34E+00±2.41E−01=
3.61E+00±4.01E−01−
2.98E+00±2.46E−01+
3.15E+00±2.49E−01=
2.54E+00±3.88E−01+
F15
1.06E+02±3.06E+01
5.63E+02±3.35E+02−
6.92E+01±1.62E+02=
3.00E+02±4.39E+01−
4.10E+01±1.10E+02+
3.64E+01±1.11E+02+
1.29E+02±1.77E+02=
F16
9.88E+01±2.21E+01
3.60E+02±4.11E+02−
1.41E+02±2.05E+01−
1.08E+02±1.87E+01=
9.70E+01±7.19E+00=
1.03E+02±5.90E+00=
9.49E+01±4.72E+00+
F17
1.22E+02±8.93E+00
4.59E+02±4.77E+02−
1.46E+02±1.55E+01−
1.03E+02±1.35E+01+
1.01E+02±6.33E+00+
1.19E+02±7.31E+00=
1.03E+02±3.83E+00+
F18
3.12E+02±2.59E+01
6.49E+02±3.11E+02−
9.30E+02±1.51E+02−
7.77E+02±1.34E+02−
7.49E+02±1.72E+02−
7.40E+02±1.66E+02−
6.29E+02±2.53E+02−
F19
3.41E+02±1.09E+02
6.86E+02±2.95E+02−
9.20E+02±1.29E+02−
8.17E+02±6.64E+01−
7.28E+02±1.93E+02−
6.40E+02±2.38E+02−
6.69E+02±2.36E+02−
F20
3.17E+02±4.09E+01
6.56E+02±2.83E+02−
9.05E+02±1.02E+02−
7.86E+02±1.33E+02−
7.87E+02±1.05E+02−
6.60E+02±2.29E+02−
6.22E+02±2.47E+02−
F21
3.69E+02±9.92E+01
8.63E+02±2.15E+02−
9.41E+02±3.02E+02−
6.04E+02±2.19E+02−
5.20E+02±2.16E+02−
4.60E+02±8.16E+01−
4.96E+02±1.79E+02−
F22
7.23E+02±1.24E+02
7.72E+02±2.49E+01=
6.84E+02±1.95E+02=
6.68E+02±1.03E+02=
7.21E+02±1.27E+02=
7.61E+02±4.63E+00=
7.47E+02±1.23E+01−
F23
5.64E+02±6.32E+01
1.07E+03±2.10E+02−
1.03E+03±2.54E+02−
5.76E+02±8.22E+01−
7.50E+02±1.88E+02−
6.18E+02±1.22E+02=
6.51E+02±1.23E+02=
F24
2.00E+02±9.22E−01
3.00E+02±1.80E+02=
2.00E+02±2.79E−03+
3.77E+02±3.69E+00−
2.00E+02±0.00E+00+
2.00E+02±0.00E+00+
2.00E+02±0.00E+00+
F25
4.01E+02±8.70E+01
4.42E+02±1.51E+02−
2.12E+02±6.00E+01+
3.79E+02±3.70E+00+
3.76E+02±2.99E+00+
3.79E+02±2.54E+01+
3.71E+02±2.65E+00+
+
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 NFEs=10000 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 α=0.05.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
2The Wilcoxon results of the multiple-problem statistical analysis are listed in Table11. The R+ in Table 11 is a measure of the case where DEMR outperforms the compared algorithm and the greater the R+ value is, the greater difference between DEMR and the compared algorithm is. Similarly, the R− is a measure of the case where DEMR is worse than the compared algorithm, and the size of R− 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 α=0.05 and α=0.1. For MPEDE, it performs worse than DEMR when α=0.1, while is comparable with DEMR at α=0.05. And for other algorithms, it shows they have similar optimization effectiveness with DEMR, since the P value goes beyond the α at both α=0.05 and α=0.1. 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−
P value
α=0.05
α=0.1
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 aref1,f5,f12,f13,f18,f21, and f23, 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 f1 and f5, DEMR converges a little early. What is more, it should be noted that, in the convergence Figure 5(a) for f1, the convergence of LBBO after about NFEs=1000 has not been shown, since LBBO has already found the optimal value 0 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 NFEs=50000 for the same reason.Figure 5
The convergence graphs on different types of CEC2005 benchmark functions.
(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) f23In 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 Figure6. It shows that DEMR ranks the second with Ranking=3.22, behind SHADE with Ranking=2.70. And SaDE follows DEMR with Ranking=3.56. 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ρ1,max, ρ2,max, and Rmax, of which ρ1,max and ρ2,max determine the time to trigger the local search algorithm SQP, and Rmax 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 Table12 and Figure 7, the Friedman rankings for different combinations of ρ1,max and ρ2,max are illustrated: ρ1,max∈1,1.5,2,2.5,3, ρ2,max∈1,1.5,2,2.5,3, where Rmax is a constant Rmax=5. From Table 12, the setting ρ1,max=3,ρ2,max=1 is the best combination, followed by combinations ρ1,max=1,ρ2,max=1 and ρ1,max=1,ρ2,max=2. Figure 7 is more intuitive for the performance of these settings. Generally speaking, the combination with both small ρ1,max and small ρ2,max is more reasonable. Table 13 gives a brief look at the test results on different Rmax with contraction criteria combination ρ1,max=3,ρ2,max=1. And it can be seen that, for CEC2005 benchmark functions, Rmax=20 is the most appropriate choice, while for the design of Earth-Moon low-energy transfer trajectory, the performance of DEMR ranks the first when Rmax=5.Table 12
Results obtained by the Friedman test for contraction criteria combinations on CEC2005.
Friedman
ρ1,max
1
1.5
2
2.5
3
ρ2,max
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.32Figure 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
Rmax
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
Rmax
AverageΔV (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.



---

*Source: 4560173-2018-05-03.xml*