Enhancing the search ability of a hybrid LSHADE for global optimization of interplanetary trajectory design

The global optimization of interplanetary trajectory design is a challenge and tough problem in deep space. To cope with the extreme nonlinearity of the search space and a large number of locally optimal solutions, a novel hybrid algorithm based on success-history based adaptive differential evolution with linear population size reduction (LSHADE), called HLSHADE, is proposed. In HLSHADE, after the global search finds a better solution, a new two-step local search with an adaptive control parameters strategy is designed to enhance the exploitation ability of HLSHADE. The presented method is tested on well-known global trajectory optimization problems (GTOPs) developed by the Advanced Concepts Team of the European Space Agency (ESA). Experimental results have demonstrated the competitive performance of HLSHADE with respect to the three sets of compared algorithms, including LSHADE and ten of its variants, six famous algorithms in the IEEE Competitions on Evolutionary Computation (CEC) and eleven algorithms in the software platform PyGMO developed by the ESA's Advanced Concepts Team.


Introduction
The design and optimization of interplanetary transfer trajectories is one of the most important tasks during the analysis and design of a deep space mission. In order to minimize the fuel consumption of the spacecraft or maximize the number of the visited asteroids, multi-gravity assist (MGA) manoeuvres that take advantage of the gravity of celestial bodies to change the momentum vector of a spacecraft without the use of any propulsion system are often used in the deep space exploration missions. The problem of optimizing multi-gravity assist trajectories with deep space manoeuvres has the characteristics of extreme nonlinearity, discontinuous, rugged landscapes.
In the last two decades, evolutionary algorithms such as differential evolution (DE), particle swarm optimization (PSO), genetic algorithms (GAs), ant colony optimization (ACO), etc. were introduced for the optimization of aerospace engineering problems (Yang et al. 2019), especially for the design of spacecraft interplanetary trajectories. Cage, Kroo, and Braun (1994) successfully introduced a genetic algorithm into the design of Jupiter's exploration orbit by combining the concept of gradient. This made the GA an effective method for orbit design. The multi-objective evolutionary algorithm was used to design the lowthrust transfer trajectory by Lee et al. (2005). In 2008, Vinkó and Izzo of the Advanced Concepts Team (ACT) of the European Space Agency (ESA) proposed a series of representative examples of space orbit design through analysis and classification of existing space tasks, which provided a unified test platform for the subsequent study of global optimization algorithms. Ampatzis and Izzo (2009) introduced machine learning into the approximate objective function of orbit optimization. The ant colony algorithm was used to solve the problem of optimizing the multi-gravity assist orbit planetary sequence by Ceriotti and Vasile (2010). In addition, Vasile, Minisci, and Locatelli (2011) demonstrated that the proposed discrete dynamical system had fixed points towards which it converged with probability one for an infinite number of generations and proposed an inflationary differential evolution algorithm (IDEA) for some space trajectory optimization problems. The hybrid hierarchical cellular genetic algorithm (HH-cGA) was proposed based on the hierarchical structure of cells by Danoy, Dorronsoro, and Bouvry (2012). This algorithm could effectively solve the Cassini2 problem published by ACT, and the optimization results were better than those of the compared algorithms. Bonami et al. (2013) designed a multi-stage mixed integer method for solving the problem of aircraft orbit optimization. Costanza and Casalino (2014) proposed a hybrid evolution algorithm combining differential evolution, a genetic algorithm and particle swarm optimization for solving multiple gravity-assisted tasks. A hybrid intelligent optimization algorithm was proposed by Zuo, Dai, and Peng (2016) for problems of multiple gravity-assisted orbit design, and achieved better results on these problems. Gao, Liu, and Wang (2017) hybridized PSO and biogeographic-based optimization (BBO) algorithms to design and compare low thrust transfer orbits from Earth to Mars. A hidden genetic algorithm was designed for solving space orbit optimization problems by Darani and Abdelkhalik (2018). The algorithm used the hidden and active state of the gene to improve search efficiency. McCarty, Burke, and McGuire (2018) proposed the parallel monotonic basin hopping algorithm, which was applied to low thrust trajectory optimization. Yao et al. (2018) proposed an improved differential evolution algorithm for solving two case studies on transfer orbit design and target observation. Zuo, Dai, and Peng (2019) proposed a multi-agent genetic algorithm with controllable mutation probability utilizing a back propagation neural network for global optimization of trajectory design. Schlueter and Munetomo (2019) studied the application of a multi-objective global optimization algorithm in the design of interplanetary orbits. Song and Gong (2019) employed the method of deep neural networks for the solar sail orbit design of near-Earth asteroid exploration. Zuo et al. (2020) proposed a case learning-based differential evolution (CLDE) algorithm for the global optimization of interplanetary trajectory design. Li et al. (2020) used deep feed-forward neural networks to estimate optimal transfer costs for three types of optimization problem in multitarget mission design.
Among the evolutionary algorithms, DE was first proposed by Storn and Price (1997) for solving global numerical optimization problems over continuous search spaces. It has emerged as one of the most competitive and versatile families of evolutionary algorithms. It is a simple yet powerful optimizer and exhibits excellent capability on a variety of numerical and real-world optimization problems. For about two decades, researchers have been working on improvements to DE and have made significant progress on DE developments.
An adaptive differential evolution (JADE) algorithm with external archiving based on a differential evolution algorithm was proposed by Jingqiao Zhang and Sanderson (2009). JADE made great improvements in parameter adaptation strategy and mutation strategy. It is called one of a few 'important variants of DE'. Based on JADE, Tanabe and Fukunaga (2013) proposed success-history based adaptive differential evolution (SHADE). SHADE used additional storage space to store the parameter values of previous generations of excellent individuals to guide population evolution. Based on JADE and SHADE, there were some variants which showed especially encouraging performance. Furthermore, Tanabe and Fukunaga (2014) proposed LSHADE, which further extended SHADE with linear population size reduction (LPSR). LPSR continually decreased the population size according to a linear function. This strategy greatly improved the convergence of the algorithm in the later period.
Since 2014, LSHADE has gained much attention, partly because LSHADE and its variants often obtained the winners in IEEE Competitions on Evolutionary Computation (CEC). On the basis of LSHADE, Guo et al. (2015) proposed the LSHADE_SPS_EIG algorithm which mainly employed the successful parent-selecting framework (SPS) and eigenvector-based (EIG) crossover operator. It showed the good global optimization capability. Awad et al. (2016) proposed the LSHADE_EpSin algorithm. It joined an effective adaptive framework, which provided a successful alternative adaptation through the selection of control parameters. In 2017, Awad, Ali, and Suganthan further proposed the LSHADE_cnEpSin algorithm, which developed an ensemble of sinusoidal approaches based on performance adaptation and covariance matrix learning for the crossover operator. Brest, Maučec, and Bošković (2017) also further proposed a new algorithm, called JSO, which was an improved variant of iL-SHADE with a new weighted version of the mutation strategy. Mohamed et al. (2017) proposed the LSHADE_SPA_CMA algorithm, which added a semi-parametric adaptive (SPA) strategy and a modified version of the covariance matrix adaptation ES (CMA-ES) hybrid framework on the basis of LSHADE, where CMA-ES is the covariance matrix adaptation evolution strategy . Hadi, Mohamed, and Jambi (2018) integrated the adaptability of mutation strategy parameters and some directed mutation strategies into LSHADE_SPA_CMA, and then proposed an ELSHADE_SPA_CMA algorithm with better performance. Stanovov, Akhmedova, and Semenkin (2018) studied the rank-based selection pressure (RSP) strategy on the basis of the LSHADE algorithm and proposed the LSHADE_RSP algorithm. Piotrowski (2018) studied whether the performance of LSHADE variants could be improved by adding a populationwide inertia term (PWI) to the mutation strategies. Poláková, Tvrdík, and Bujok (2019) presented a new mechanism for the adaptation of the population size, which they called the diversitybased strategy. The size of the population depended on its diversity. Stanovov, Akhmedova, and Semenkin (2019) further proposed a modification of differential evolution mutation strategies with the introduction of selective pressure, which was implemented by applying proportional, rankbased and tournament selection. Based on the new mutation strategies, a new algorithm called LSHADE-SP was proposed. Viktorin et al. (2019) proposed a simple yet effective modification to the scaling factor and crossover rate adaptation strategies in SHADE. The proposed distance-based parameter adaptation was designed to address the premature convergence of SHADE-based algorithms in higher dimensional search spaces to maintain a longer exploration phase. Chakraborty et al. (2021) further enhanced success history-based adaptive differential evolution (SHADE) by hybridizing it with a modified whale optimization algorithm (WOA). Wang et al. (2021) proposed L-SHADE-E, an ensemble of L-SHADE-EpSin and L-SHADE-RSP in this article. In L-SHADE-E, L-SHADE-EpSin or L-SHADE-RSP was selected randomly to run at the beginning. If progress on fitness was low for generations, the other constituent algorithm took over the population immediately.
The global optimization of interplanetary trajectory design has the characteristics of an extremely nonlinear search space, a large number of locally optimal solutions, and sensitivity to initial conditions. Although some previous studies used DE, PSO, GA and ACO, etc. to solve these problems, in some cases they have been shown to converge to a fixed point, a level set (Vasile, Minisci, and Locatelli 2011) or a hyperplane (Locatelli and Vasile 2015) not containing the global optimum and have slow local convergence. For example, when DE was utilized for the Cassini1 mission on the global trajectory optimization problem (GTOP) benchmark, it always converged to a sub-optimal solution with objective function value J = 5.30 km/s, while the global optimum with objective function value J = 4.93 km/s was very difficult to reach (Izzo 2010). Figure 1 shows the search space characteristics of the Messenger-full problem on GTOP. There are 26 variables in this problem. The variables of x 17 and x 24 are randomly selected among them, and plotted the relationship between the two variables and the objective function (V). From Figure 1, it can be seen that the objective function is highly nonlinear in the range of small independent variables, and a large number of locally optimal solutions around the global solution. Therefore, it requires the algorithm to have better global and local search capabilities. The contributions of this article are as follows.
• Through preliminary experimental analysis, it is found that LSHADE_SPS_EIG shows good global optimization capability in interplanetary trajectory design problems. However, the mining capacity of LSHADE_SPS_EIG around the local optimal attraction basin is relatively weak. Therefore, to alleviate the difficulty of trajectory optimization discussed above, a hybrid LSHADE with twostep local search, named HLSHADE is proposed. Firstly, HLSHADE obtains a good solution by LSHADE_SPS_EIG in the search space, and then iteratively uses CMA-ES to exploit a better solution around the basin of attraction of the already-found promising solution. Finally, in order to improve the exploitation ability further at a later stage, the interior-point method is employed to search for the optimal solution. In consideration of the extremely rugged objective function, this novel two-step local search strategy can effectively enhance the exploitation capability. • In HLSHADE, two extra parameters are introduced to switch between the local and global search processes. In this article, an adaptive control parameters strategy is proposed to determine dynamically when to start the local search (LS). If LS successfully updates the optimal solution, the control parameters will be rewarded by increasing, allowing LS to run more iterations. If LS does not update the optimal solution, then the parameters will be reduced as the penalty.
The rest of this article is organized as follows. Section 2 describes the mathematical formalization of MGA problems. Section 3 mainly describes the overall framework of HLSHADE, the novel two-step local search, and the adaptive control parameters strategy. The experimental results and comparative analysis are presented in Section 4. Finally, Section 5 concludes the article and discusses future work.

MGA trajectory model
Multi-gravity assist is an important technology for providing thrust to a spacecraft, which can help the spacecraft reach the target planet. The MGA problem is often used to find the best orbit of the spacecraft from Earth to the target planet. In this article, the sequence of gravity-assist planets is usually already determined. What needs to be optimized is the launch date and the transfer time of each leg. The generic mathematical form of the MGA trajectory problem can be written as where x is the decision vector, f (x) is the objective function and g(x) are nonlinear constraint functions that depend on the problem considered and are defined case by case. For example, for the Cassini1 mission on GTOP, the constraints g(x) limit the various fly-by pericenter to the values: r p1 ≥ 6351.8 km, r p2 ≥ 6351.8 km, r p3 ≥ 6778.1 km, r p4 ≥ 600,000 km. More details of the constraints g(x) on the various interplanetary missions can be found in Vinkó and Izzo (2008).
As shown in Figure 2, the planet sequence can be described as P = {P 0 , P 1 , . . . , P i , . . . , P N }, where P 0 is the departure planet (usually Earth), P i (i = 1, 2, . . . , N − 1) is a gravity-assist planet and P N is the target planet.
The transfer trajectory between P i−1 and P i can be regarded as Lambert's problem (Battin 1999). If the transfer time T i is known, the velocity vectors v i−1 and v i at positions P i−1 and P i can be calculated. Therefore, the objective function of the MGA problem can be described as follows: where v 0 is departure velocity, v i denotes the velocity change needed at the pericenter of the gravity assist hyperbola to connect two conic arcs, v t is the manoeuvre required to arrive at the target orbit, and x is the decision vector, which can be presented by where t 0 is the departure time and T i (i = 1, 2, . . . , N) is the transfer time of the spacecraft from the gravity-assist planet P i−1 to P i . Because the initial velocity is not a free parameter, v 0 is defined as the modulus of the vector difference between the velocity of Earth and the velocity of the spacecraft at time t 0 .

MGA trajectory model with deep space manoeuvres (DSMs)
Compared with the MGA problem, the MGA-DSMs problem can change the speed of a spacecraft by DSMs. This means that a spacecraft is able to thrust its engine during each trajectory leg in Figure 3. The objective function of a MGA-DSMs problem can be described by where DSM j is the velocity change due to the DSM between P j−1 and P j ; the definitions of v 0 , v i and v t are the same as for Equation (2). Given a trajectory including n planets, the decision vector x can be shown by where t 0 is the departure time, V ∞ , u, v define the heliocentric direction of the departure hyperbolic velocity v ∞ , b incl i (angle measured in planet B's plane of the planet approach vector) and r p i (radius of flyby) provide the ith swing-by conditions, and T i is the transfer time of the spacecraft from the gravity-assist planet P i−1 to P i . Given the transfer time T i and the variable η i relative to each leg i, it can propagate the spacecraft trajectory along a keplerian orbit for a time η i T i .

Proposed algorithm
In this section, details about the two-step local search strategy and adaptive strategy of the control parameters are described and the framework of HLSHADE is given in Algorithm 1.

The two-step local search
Because the global optimization of interplanetary trajectory design has the characteristics of an extremely nonlinear search space, a large number of locally optimal solutions and sensitivity to initial conditions, it calls for the design of hybrid optimization techniques that effectively explore and exploit the possible solutions. In HLSHADE, LSHADE_SPS_EIG (more details can be found in Guo et al. 2015) is used as the global search method, and then a two-step local search strategy is proposed in order to enhance the exploitation ability. The two-step local search consists of two parts: CMA-ES and the interior-point method.

The CMA-ES strategy
The covariance matrix adaptation evolution strategy (CMA-ES) is one of the evolution strategies with good performance and wide application range. In particular, CMA-ES very effective on medium scale optimization problems. The main idea of the CMA-ES algorithm is to use the covariance matrix C t in the normal distribution N(m t , σ 2 t C t ) to adjust and guide the evolution of the algorithm. So the key factor of the algorithm is how to adjust these parameters, especially the step size and covariance matrix. The adjustment idea is to increase the probability that the algorithm produces a better solution.
CMA-ES is used for local search in HLSHADE. The search range of CMA-ES is determined according to the optimal solution obtained by LSHADE_SPS_EIG, which is shown by Algorithm 1 The framework of HLSHADE 1: x best DE denotes the best solution found by LSHADE_SPS_EIG, x best CMA denotes the best solution found by CMA-ES, x best IP denotes the best solution found by the interior-point method, x gmin denotes the global optimum, f (x) is the objective function. ls_eval is the function evaluations scale factor of the interior-point method. NP is the population size. FES CMA denotes the number of function evaluations during the iterative process in CMA-ES. FES ip denotes the number of function evaluations during the iterative process in the interior-point method. 2: ρ 1 = 0, ρ 1,max = 20, ρ 2 = 0, ρ 2,max = 10, Bound init = 0.5, Bound min = 0.1, ls_eval = 0.01. 3: while FES < MAX_FES do 4: 22: 23: while FES < MAX_FES and ρ 2 < ρ 2,max do 24:  The main steps of the algorithm are as follows.
(1) Initialize parameters by where E n is the identity matrix of order n. Note that n is the dimension of the problem. (2) Use the normal distribution N(m t , σ 2 t C t ) to generate new individuals of the population as shown by (3) Sort the individuals in the population according to their fitness, and take the first μ best individuals to update the average vector m as where μ i=1 w i = 1, w 1 ≥ w 2 ≥ · · · ≥ w μ and w i is determined by (4) Update the step size σ t+1 and the covariance matrix C t+1 .

The interior-point method
The interior-point method is a traditional method for linear and (convex) quadratic programming. The idea is to construct an unconstrained objective function, also called a barrier function, and find the extreme point of the barrier function within the feasible domain. The interior-point method has low-degree polynomial worst-case complexity and an unrivalled ability to deliver optimal solutions in an almost constant number of iterations that depend very little, if at all, on the problem dimension. The constraints of general forms are represented as The barrier function associated with the linear optimization problem (LOP) is shown by where μ k > 0 is called the barrier parameter, which is a decreasing sequence of positive parameters tending to zero as shown by The basic steps of the interior-point method are as follows.
(2) Take the initial point x 0 within the feasible domain, k = 1.
In order to improve the exploitation ability further at a later stage, HLSHADE stipulates that, when the number of function evaluations exceeds 0.75 × MAX_FES, the interior-point method is triggered to find a local optimal solution. This can enhance the local search ability of HLSHADE.
FES ip denotes the number of function evaluations during the iterative process in the interior-point method. The initial value of FES ip is equal to 0.01 × MAX_FES. If the optimal solution is successfully updated, the value of FES ip will be increased. This can give the interior-point method more opportunities to exploit the optimal solutions. If the optimal solution is not updated, then the value of FES ip will be decreased, which can reduce the consumption of function evaluations. More details are shown in Algorithm 1.

The adaptive control parameters strategy
For successful incorporation of CMA-ES and the interior-point method in HLSHADE, two triggering parameters are designed to decide when the local searches start. The adaptive strategy is proposed so as to control the two triggering parameters. In HLSHADE, LSHADE_SPS_EIG is used as the global optimizer. When LSHADE_SPS_EIG stagnates in updating the optimal solution, the variable ρ 1 = 0 is defined to count stagnant generations and the threshold is set as ρ 1,max = 20. If ρ 1 > ρ 1,max , CMA-ES starts to perform a local search. If LSHADE_SPS_EIG successfully finds a better solution than the current optimal solution, ρ 1,max will be rewarded by increasing. Otherwise, ρ 1,max will be reduced by way of punishment. More details are shown by ρ 1,max = min(ρ 1,max + 5, 30) optimal solution becomes better max(ρ 1,max − 5, 5) otherwise.

Experimental study
In this article, the seven well-known global trajectory optimization problems (GTOP) are employed to evaluate the performances of HLSHADE. These benchmark problems involve MGA or MGA-DSMs, which are typical nonlinear constrained global optimization problems. Table 1 introduces basic information on the benchmarks used in this article. The range of independent variables of each benchmark can refer to the ESA-ACT team's global trajectory optimization problem database website http://www.esa.int/gsp/ACT/projects/gtop/. HLSHADE is compared with three sets of state-of-the-art algorithms, including LSHADE and ten of its variants, six famous CEC competition algorithms, and eleven algorithms on the software platform PyGMO developed by ESA-ACT. Furthermore, experimental analysis of the adaptive parameter settings and each operator of HLSHADE have been carried out to determine the initial values of the parameters and to confirm further the effectiveness of each part of HLSHADE. Finally, experiments In the experiments on GTOP, each test function is independently run 25 times with 150,000 function evaluations as the termination condition. The test parameters are referred to the CEC2011 competition on real-world numerical optimization problems (Das and Suganthan 2011). The best, worst, mean and standard deviation are recorded over 25 independent runs for each test problem. All algorithms are executed on a computer with a 2.80 GHz Intel Core tm i7-7700 HQ processor and 16 GB RAM on 64-bit Windows 10 and MATLAB 2017b.
For the above twelve algorithms, the same parameter settings are used in that the maximum number of function evaluations is 150,000 and the initial population size is NP = 18 × D, where D is the dimension of each benchmark on GTOP. Table S1 in the online supplemental data for this article, which can be accessed at https://doi.org/10.1080/0305215X.2021.2019250, shows the statistical results for HLSHADE and the other eleven algorithms on GTOP benchmarks. The optimal solutions have been shown in boldface. It can be seen that LSHADE_SPA_CMA, LSHADE_SPS_EIG, LSHADE_EpSin, JSO, LSHADE_RSP, DISH and HLSHADE obtain the optimal solution 4.9307 on the Cassini1 problem. HLSHADE obtains a better mean value than the algorithms with which it is compared. On the Cassini2 problem, LSHADE_PWI obtains the best solution, i.e. 8.6089, and HLSHADE has the best mean value, i.e. 12.0239. On the Gtoc1 problem, HLSHADE finds a better solution than the others, i.e. −1,576,143.0, and the mean value of LSHADE_SPS_EIG is slightly better than that of HLSHADE. On the Rosetta problem, JSO obtains the best solution, i.e. 1.3584, and LSHADE_EpSin has a better average value than the other algorithms, i.e. 1.9145. On the Sagas problem, HLSHADE is superior to the other algorithms in terms of both best value and mean value. Moreover, ELSHADE_SPA_CMA, LSHADE_EpSin and LSHADE_cnEpSin obtain a better solution, i.e. 8.7016, and the mean value of JSO is better than the others on the Messenger problem. As for the most difficult Messengerfull problem on GTOP, it is worth noting that HLSHADE obtains much better results than the algorithms with which it is compared for both the best value (5.2185) and the mean value (10.5344).
Note that '+', '−' and ' = ' denote that the performance of HLSHADE is respectively better than, worse than or similar to the compared algorithms according to the Wilcoxon rank sum test at α = 0.05.
To provide a correct interpretation of the results of all optimization algorithms and obtain their average rankings, two non-parametric statistical tests, as similarly done in García et al. (2008), are used in the experiments: (i) the Wilcoxon tests are performed to analyse the statistical significance of the experimental results between two algorithms in both the single-problem and multi-problem cases; (ii) the Friedman test is employed to obtain the average rankings of all the compared algorithms. If the Friedman test score of an algorithm is smaller, its performance is better. The Wilcoxon rank sum test at α = 0.05 for the single-problem case is calculated by MATLAB, while the multi-problem Wilcoxon test and the Friedman test are carried out by the KEEL software in Alcalá-Fdez et al. (2009). Table 2 shows the results of the Wilcoxon rank sum test between HLSHADE and the compared algorithms. It can be seen that HLSHADE performs significantly better than both SHADE and ELSHADE_SPA_CMA on seven test problems. Compared with LSHADE, HLSHADE shows better and worse performance on six test problems and one test problem, respectively. HLSHADE wins on five test problems over LSHADE_SPA_CMA and LSHADE_SPS_EIG. HLSHADE is significantly better on five test problems and worse on one test problem than JSO, LSHADE_EpSin and DISH. Moreover, HLSHADE outperforms LSHADE_cnEpSin and LSHADE_RSP on four test problems, while it is worse than these two algorithms on one test problem. HLSHADE is better on four test problems than LSHADE_PWI. Furthermore, HLSHADE obtains one result that is similar to those of JSO, LSHADE_EpSin and DISH, two results that are similar to those of LSHADE_SPA_CMA, LSHADE_SPS_EIG, LSHADE_cnEpSin and LSHADE_RSP, and three results that are similar to that of LSHADE_PWI.
To detect further differences in behaviour between HLSHADE and the eleven competitors, the multiple-problem Wilcoxon test and the Friedman test are carried out based on the mean  Table 3 indicate that HLSHADE has the best ranking (2.4286) among the twelve algorithms. In general, the above comparison clearly demonstrates that HLSHADE performs significantly better than its competitors on GTOP.
The convergence graphs of the twelve algorithms on GTOP are shown in Figure S1 of the online supplemental data. The y-axis indicates the objective function value, and the x-axis indicates the number of functions evaluations (FES). The curve represents the fact that the objective function value changes with the number of function evaluations. It can be seen from Figure S1 of the online supplemental data that HLSHADE has a faster convergence speed than the other algorithms, especially on the Cassini1, Cassini2, Gtoc1, Sagas, Messenger and Messenger-full benchmarks. In addition, it is worth noting that the black curve of HLSHADE obviously has a rapid drop in the later stages of the evolution process, and that the time of this descending curve is related to the start time of local search. It can be seen that the two-step local search and adaptive control parameters strategies in this article are effective in enhancing the exploitation ability on GTOP benchmarks.

Comparison with CEC competition algorithms
HLSHADE is also compared with six other EAs: UMOEAsII (Elsayed, Hamza, and Sarker 2016), EBOwithCMAR (Kumar, Misra, and Singh 2017), MVMO (Rueda and Erlich 2015), HSES (Geng Zhang and Shi 2018), GA-MPC (Elsayed, Sarker, and Essam 2011) and PaDE (Meng, Pan, and Tseng 2019). There are two reasons for choosing these six algorithms: (1) they are the top-level representative algorithms in recent CEC competitions, as shown in Table 4, and their performances are very competitive. (2) they are the new powerful ones.
The parameter settings for each CEC competition algorithm are as follows:  From Table S3 in the online supplemental data, EBOwithCMAR and HLSHADE obtain the optimal solution 4.9307 on the Cassini1 benchmark, and HLSHADE performs better with respect to mean value than the other algorithms. As for the Cassini2 benchmark, GA-MPC obtains the best solution 8.6095, and the mean value 11.3303 is also the best among the seven algorithms. On Gtoc1, HLSHADE finds the best solution −1,576,143.0, and MVMO outperforms HLSHADE on mean value. Moreover, UMOEAsII finds the best solution 1.3796 and EBOwithCMAR with an average value of 1.9528 performs better than the compared algorithms on the Rosetta benchmark. HLSHADE is superior to the other algorithms with respect to both the best value 18.5826 and the mean value 281.7730 on the Sagas benchmark. EBOwithCMAR obtains the best solution 8.6529, and the mean value of EBOwithCMAR is also better than the others on the Messenger benchmark. As for the Messengerfull benchmark, it is evident that HLSHADE outperforms the compared algorithms with respect to both the best value 5.2185 and the mean value 10.5344. Table 5 presents the Wilcoxon rank sum test results between HLSHADE and the compared algorithms. Specifically, HLSHADE outperforms UMOEAsII on three test problems. Compared with EBOwithCMAR, HLSHADE shows better and worse performance on three test problems and one test problem, respectively. Moreover, HLSHADE surpasses MVMO on four, HSES on seven, GA-MPC on five, and PaDE on six test problems. In addition, HLSHADE obtains results similar to those of UMOEAsII, EBOwithCMAR, MVMO, GA-MPC and PaDE on four, three, two, two and one case, respectively. Table S4 of the online supplemental data and Table 6 also present the statistical analysis results according to the multiple-problem Wilcoxon test and the Friedman test, respectively. It can be seen from Table S4 in the online supplemental data that HLSHADE performs significantly better than HSES, GA-MPC and PaDE at α = 0.05. Furthermore, it is significantly better than UMOEAsII, EBOwithCMAR, HSES, GA-MPC and PaDE at α = 0.1. On the other hand, the Friedman test results in Table 6 indicate that HLSHADE has the best ranking (2.0000) among the seven algorithms. The score of UMOEAsII is 2.7143 which is the second ranking. The score of EBOwithCMAR is 3.1429, which is the third ranking. The convergence graphs of the seven algorithms on GTOP are shown in Figure S2 of see the online supplemental data. The convergence performance of HLSHADE is better than that of the compared algorithms on the Cassini1, Sagas and Messenger-full benchmarks. In summary, the above comparison demonstrates that HLSHADE is significantly better than its six competitors.
Experimental analysis and results for eleven algorithms in PyGMO and five algorithms on the CEC2017 benchmark suite are presented in the online supplemental data. Furthermore, experimental analysis of the adaptive parameter settings and the effectiveness of each part of HLSHADE are also presented in the online supplemental data. Note that '+', '−' and ' = ' denote that the performance of HLSHADE is respectively better than, worse than or similar to the compared algorithms according to the Wilcoxon rank sum test at α = 0.05.

Conclusion
In this article, a new hybrid LSHADE, called HLSHADE, is proposed for enhancing searchability on interplanetary trajectory design problems. In HLSHADE, a two-step local search is proposed for improving its exploitation ability. The adaptive control parameters strategy is applied to decide dynamically when the two-step local search starts. These improvements make HLSHADE a powerful alternative for interplanetary trajectory design problems. To analyse the effectiveness of HLSHADE, experiments were conducted on GTOP benchmark problems. The superiority of HLSHADE was also verified by comparing it with LSHADE and ten of its variants, six top-level CEC competition algorithms and eleven algorithms that had been successfully used in PyGMO.
In future work, HLSHADE will be tested on other real-world applications. Moreover, it can be seen from Figure S2 in the online supplemental data that HLSHADE has a slower convergence speed in the early and middle stages than the compared algorithms, especially on the Cassini2, Rosetta and Messenger-full benchmark problems. Therefore, a more powerful mutation strategy could also be designed to enhance the early and mid-term exploration abilities for interplanetary trajectory design. Finally, according to experimental analysis and results on the CEC2017 benchmark suite in the online supplemental data, DISH has good global search on the CEC2017 benchmark suite, especially on high-dimensional functions. Therefore, hybrid improvements based on DISH and its variants will be studied to solve interplanetary trajectory optimization problems.

Disclosure statement
No potential conflict of interest was reported by the author(s).

Data availability statement
The data that support the findings of this study are available from the corresponding author, [Lei Peng. Email: lei.peng@cug.edu.cn], upon reasonable request.