1 Introduction
Metaheuristic algorithms have demonstrated their potential in a large number of studies (Glover and Kochenberger, 2003; Blum and Roli, 2003; Dokeroglu et al., 2019; Tsai and Chiang, 2023) and applications (Soler-Dominguez et al., 2017; Osaba et al., 2021; Elaziz et al., 2021; Sharma and Tripathi, 2022; Yaqoob et al., 2024). From the perspective of algorithm development, it can be seen that various metaheuristic algorithms, e.g. tabu search (TS) in Glover (1989), simulated annealing in Kirkpatrick et al. (1983), genetic algorithm (GA) in Holland (1975), ant colony optimization (ACO) in Dorigo et al. (1991), particle swarm optimization (PSO) in Kennedy and Eberhart (1995), and differential evolution (DE) in Storn and Price (1995), have been presented to solve optimization problems for years. Many new metaheuristic algorithms have also been presented in recent years (Dokeroglu et al., 2019; Abualigah et al., 2022; Velasco et al., 2024). From the perspective of applications, several studies show that metaheuristic algorithms can be applied to many real-world applications, such as the internet of things (IoT) device control (Sharma and Tripathi, 2022), stock investments (Soler-Dominguez et al., 2017), and even the automatic design of deep neural network architectures (Elaziz et al., 2021; Liu et al., 2023). In addition to using metaheuristic algorithms for real-world applications, an interesting reference workflow for a metaheuristic algorithm can be found in a recent study (Osaba et al., 2021) that explained how to use a metaheuristic algorithm to solve real-world optimization problems, ranging from problem modelling, algorithm design, and solution encoding, to performance assessment, comparison, and replicability.
Researchers in this field have focused on enhancing search performance by presenting a new metaheuristic algorithm, modifying the transition operator, and combining two (meta)heuristic algorithms (Tsai and Chiang, 2023). It is well established that no metaheuristic algorithm is optimal for all optimization problems (Wolpert and Macready, 1997). However, this does not mean that they are not useful; on the contrary, every metaheuristic algorithm has its pros and cons. For example, the tabu list of TS (Glover, 1989) can be used to avoid searching for the same solution too often. The pheromone table of ACO (Dorigo et al., 1991) is used to accumulate the search results of all searches from the very beginning, while the personal and global best solutions of PSO (Kennedy and Eberhart, 1995) are used to guide the search to promising regions in the solution space. Moreover, in addition to modifying the transition operator and combining two (meta)heuristic algorithms to improve results or accelerate the convergence speed, “analysing the landscape of solution space” (Pitzer and Affenzeller, 2012; Ochoa and Malan, 2019; Zou et al., 2022) and “using a machine learning (ML) algorithm” (Talbi, 2021) to determine a suitable search strategy during the convergence process are two promising research directions that have been discussed in recent studies.
Since DE is one of the most well-known and successful evolution-based metaheuristic algorithms for solving the single-objective bound-constrained problem (SOP), many studies have proposed different ways to enhance its performance, such as an improved version of DE (Tanabe and Fukunaga, 2013, 2014) or combining it with other methods (Lin et al., 2023; Zhang et al., 2025). The transition operator (e.g. mutation) of DE has a strong impact on its search performance, just like the transition operators of the other search algorithms in the family of metaheuristic algorithms. Consequently, using two or more transition operators in a DE represents a viable approach to improving its search performance. Recently, an effective way to dynamically determine a suitable transition operator has been presented for the DE in solving an optimization problem (Tan et al., 2021, 2022b). The basic idea is to adopt the fitness landscape analysis (FLA) to understand the landscape characteristics of the optimization problem (such as roughness, neutrality, or gradient) first and then use this information to generate a learning model for the metaheuristic algorithm to further select the suitable search operator for different kinds of optimization problems. In other words, the machine learning algorithm will be used to understand the possible search performance of DE mutation operators for the landscape characteristics of the problem, and the learning model generated by the machine learning algorithm will be used to select a suitable mutation operator for different optimization problems. For example, the learning model generated by the machine learning algorithm will be used to determine the mutation operator before the search process begins (Tan et al., 2021). In a later study, Tan et al. (2022b) presented a similar method to analyse the landscape characteristics and use a learning model to determine the mutation operator at a fixed number of generations during the convergence process. Of course, using reinforcement learning (RL) as the machine learning component of such an improved DE-based algorithm is also a viable approach (Tan and Li, 2021; Tan et al., 2022a).
Based on our observations, most DEs (Tan et al., 2021; Tan and Li, 2021; Tan et al., 2022a, 2022b) do not account for: (1) the fact that the search state of DE might change after a different number of generations during every state change in the convergence process and (2) the landscape characteristics around the current searched solutions to further select the suitable mutation operator for DE in later generations. Therefore, we believe that there is still plenty of room to enhance the search performance of DE in solving single-objective optimization problems. Motivated by these observations, an effective DE based on FLA will be presented in this study that considers adaptively adjusting the timing to (1) analyse the landscape characteristics around current searched solutions and (2) use this landscape information to further select the appropriate mutation operator during the search process. The main contributions are as follows.
-
• A new method is proposed to generate detailed landscape characteristics and recommend the suitable mutation operator for the current search state.
-
• The analysis region is adaptively reduced during the search process to analyse the landscape characteristics in detail, as well as using more random walks for complicated landscapes and fewer random walks for simpler landscapes.
-
• When the search region differs from the previous analysis region, the proposed algorithm will re-analyse the landscape characteristics to determine the mutation operator, which can be used to reduce the number of times FLA is used in RL-based DE.
The remainder of the paper is organized as follows. Section 2 begins with the problem definition of the single-objective bound-constrained problem (SOP), followed by a brief review of FLA and metaheuristic algorithms based on FLA for solving optimization problems. Section 3 starts with the basic idea of the proposed algorithm, followed by a detailed description of each operator. To evaluate the performance of the proposed algorithm, we first compare a set of DEs with the proposed algorithm in Section 4. Moreover, we also compare the proposed algorithm with other state-of-the-art FLA-based DEs. Section 5 gives the conclusions of this research and shows some possible research directions of the proposed algorithm.
2 Related Work
2.1 Problem Definition
As one of the most well-known optimization problems, the SOP can be regarded as a representative optimization problem in the continuous space. This kind of optimization problem can also be found in many real-world applications. An SOP can generally be defined as follows (Mohamed et al., 2020; Kumar et al., 2021):
where $f()$ is the objective function to be minimized or maximized, $s=({s_{1}},{s_{2}},\dots ,{s_{d}})$ is a feasible solution in the solution space, $\mathbb{A}$ is the solution space of the optimization problem, and d is the number of dimensions of s. Moreover, ${L_{i}}$ and ${U_{i}}$ are the lower and upper bounds of the i-th dimension of s, respectively, ${h_{j}^{a}}()$ is the j-th inequality constraint, while ${h_{k}^{b}}()$ is the k-th equality constraint. From this definition, it can be easily seen that the aim of SOP is to find a solution that has the minimal or maximal objective value within the limited space (lower and upper bounds of each dimension) and that satisfies the optional constraints ${h^{a}}$ and ${h^{b}}$. The finding process of a metaheuristic algorithm for solving the SOP is usually called the search (or convergence) process. Also, because most SOPs are simple and can be easily realized, many recent studies (Tan et al., 2021, 2022b; Tong et al., 2021; Rani et al., 2024) used them as the touchstone to understand the performance of a metaheuristic algorithm.
(1)
\[ \begin{aligned}{}& \underset{s\in \mathbb{A}}{\min }f(s),\hspace{2.5pt}\text{subject to}\hspace{2.5pt}s=({s_{1}},{s_{2}},\dots ,{s_{d}}),\hspace{2.5pt}\text{and}\\ {} & \left\{\begin{array}{l@{\hskip4.0pt}l}{L_{i}}\leqslant {s_{i}}\leqslant {U_{i}},\hspace{1em}& i\in \{1,2,\dots ,d\},\\ {} {h_{j}^{a}}(s)\leqslant 0,\hspace{1em}& j\in \{1,2,\dots ,p\},\\ {} {h_{k}^{b}}(s)=0,\hspace{1em}& k,\end{array}\right.\end{aligned}\]2.2 Fitness Landscape Analysis
How to describe and display the landscape of the solution space can be dated back to the 1930s or even earlier. Among them is a very early study (Wright, 1932), in which Wright used contour lines to show the fitness values of a set of chromosomes to visualize a two-dimensional fitness landscape metaphor for evolutionary biology. From this perspective, valleys, peaks, ridges, and plateaus can be used to visualize the kind of fitness landscape a search algorithm confronts, as shown in Fig. 1. More precisely, in this example, “+” and “−” represent peaks and valleys, while darker colours represent lower fitness values and vice versa. Several later studies (Grefenstette, 1999; Malan and Engelbrecht, 2013) attempted to visualize the relationships between feasible solutions, fitness values of solutions, and the landscape of the solution space to further understand the search performance of a metaheuristic algorithm for various optimization problems. A general and formal definition of the fitness landscape, resting on three basic assumptions, can be found in Stadler (2002):
With this definition and the three assumptions, we can easily visualize the landscape of a discrete solution space. For instance, a candidate solution v can be generated from the current solution s by a transition operator in such a way that the Hamming distance between these two solutions is one. This definition can also be used for the continuous landscape. For example, a candidate solution v can also be a neighbour of the current solution s so that the Euclidean distance between these two solutions is less than a predefined small value (Stadler, 2002; Malan and Engelbrecht, 2013). Other studies (Pitzer and Affenzeller, 2012; Ochoa and Malan, 2019; Zou et al., 2022) also focus on understanding the fitness landscape of the solution space of an optimization problem via different perspectives when we want to use a metaheuristic algorithm to solve it. For example, the information from FLA can be used to know the ruggedness (Weinberger, 1990; Malan and Engelbrecht, 2009) or difficulty (Jones and Forrest, 1995) of a search algorithm for solving a particular optimization problem. This means that we can use the characteristics of the solution space of an optimization problem extracted by FLA (e.g. local optima and modality, basins of attraction and funnels, ruggedness and neutrality, gradients and deception, and searchability) to further design a better search algorithm. Nowadays, several methods for characterizing fitness functions and landscapes have been discussed in various studies (Malan and Engelbrecht, 2013; Richter and Engelbrecht, 2014; Ochoa and Malan, 2019) to understand the ruggedness, neutrality, and difficulty of the fitness landscape, such as autocorrelation function (Weinberger, 1990), ruggedness of information entropy (RIE) (Vassilev et al., 2000), neutral walk (Reidys and Stadler, 2001), and fitness distance correlation (FDC) (Jones and Forrest, 1995).
Among them, the autocorrelation function, presented in Weinberger (1990), is an early method based on random walks to understand the “ruggedness” of the binary fitness landscape. Figure 2 shows the landscapes of Rastrigin and Schwefel 2.26 test functions, respectively, which explain that different test functions might have different kinds of ruggedness. The autocorrelation function assumes that the behaviours of a random walk on a landscape will be the same even if they started from different points. Based on this assumption, the autocorrelation function will then use a sequence of fitness values generated by the random walk process on the landscape to calculate the correlation with the same sequence of values a small distance away to further know the ruggedness of the fitness landscape. By using this method, the fitness landscape of a solution space is smooth when the autocorrelation is high and vice versa. Vassilev et al. (2000) presented the ruggedness of information entropy (RIE) to measure the ruggedness of the fitness landscape. This RIE metric, of course, also uses a random walk process to sample a set of solutions. A string $\textbf{s}(\epsilon )=\{{\textbf{s}_{1}},{\textbf{s}_{2}},\dots ,{\textbf{s}_{n}}\}$ in RIE will be used to represent the ruggedness of the fitness landscape. An element ${\textbf{s}_{i}}$ of $\textbf{s}(\epsilon )$ can be defined as:
where ϵ is a real number to determine the accuracy of the calculation of the string $\textbf{s}(\epsilon )$ which is the observation value of the fitness landscape, n is both the number of solutions and the number of sampling steps of the random walk process, and ${f_{i}}$ is the objective value of the i-th sampling solution by a random walk process.
(2)
\[ {\textbf{s}_{i}}=\left\{\begin{array}{l@{\hskip4.0pt}l}\bar{1},\hspace{1em}& \text{for}\hspace{2.5pt}{f_{i}}-{f_{i-1}}\lt -\epsilon ,\\ {} 0,\hspace{1em}& \text{for}\hspace{2.5pt}|{f_{i}}-{f_{i-1}}|\leqslant \epsilon ,\\ {} 1,\hspace{1em}& \text{for}\hspace{2.5pt}{f_{i}}-{f_{i-1}}\gt \epsilon ,\end{array}\right.\]Malan and Engelbrecht (2009) proposed a different approach for measuring the ruggedness of the fitness landscape of an optimization problem, in which strings are generated between predefined maximum and minimum values. The string $\textbf{s}(\epsilon )$, based on Eq. (2) with ϵ defined by
and ${\varepsilon ^{\ast }}$ by
is used to calculate the information entropy value for generating nine different types of strings to map different states of the fitness landscapes, as shown in Table 1. This means that whether the fitness landscape is rugged, neutral, or smooth can then be measured by three neighbour candidate solutions.
(3)
\[ \epsilon =\bigg\{0,\hspace{2.5pt}\frac{{\epsilon ^{\ast }}}{128},\hspace{2.5pt}\frac{{\epsilon ^{\ast }}}{64},\hspace{2.5pt}\frac{{\epsilon ^{\ast }}}{32},\hspace{2.5pt}\frac{{\epsilon ^{\ast }}}{16},\hspace{2.5pt}\frac{{\epsilon ^{\ast }}}{8},\hspace{2.5pt}\frac{{\epsilon ^{\ast }}}{4},\hspace{2.5pt}\frac{{\epsilon ^{\ast }}}{2},\hspace{2.5pt}{\epsilon ^{\ast }}\bigg\}\](4)
\[ {\epsilon ^{\ast }}=\max \big(|{f_{2}}-{f_{1}}|,|{f_{3}}-{f_{2}}|,\dots ,|{f_{n}}-{f_{n-1}}|\big),\]The entropic measure $H(\epsilon )$ is calculated as follows:
where a and b, elements from the set $\{\bar{1},0,1\}$, are adjacent elements in the string $\textbf{s}(\epsilon )$, and ${P_{[ab]}}={n_{[ab]}}/n$ with ${n_{[ab]}}$ representing the number of sub-blocks $[ab]$ appearing in the string $\textbf{s}(\epsilon )$. In summary, by using Eq. (5) to measure the fitness landscape, a higher value of $H(\epsilon )$ implies a more rugged landscape in the random walk process.
Another type of FLA measurement can be found in Reidys and Stadler (2001), in which it is called the neutral walk and can be used to understand the “neutrality” of the fitness landscape. This method will start with a random solution ${s_{0}}$ in the search space and use the neutral walk process to generate all neutral neighbour solutions around ${s_{0}}$. The fitness values of these neutral neighbour solutions (e.g. ${s_{1}}$) from a neutral walk process need to be the same or similar if and only if $f({s_{0}})=f({s_{1}})$. This neutral walk process will continue until it cannot find more neutral neighbours farther away than the current neutral neighbours. The number of steps in the neutral walk will then be used as a measure of the neutrality of the fitness landscape. The fitness distance correlation (Jones and Forrest, 1995) provides an alternative way to measure the “difficulty” of an optimization problem for the search algorithm. The basic idea is to define the correlation coefficient between the objective values of solutions and the distances of the solutions from their corresponding best-searched solutions as follows:
where $\mathcal{F}$ is the correlation coefficient value, $F=\{{f_{1}},{f_{2}},\dots ,{f_{n}}\}$ is a set of objective values of solutions, and $D=\{{d_{1}},{d_{2}},\dots ,{d_{n}}\}$ is the corresponding distances of those solutions to the optimal solution. Besides, $\text{cov}(F,D)=\frac{1}{n}{\textstyle\sum _{i=1}^{n}}({f_{i}}-\bar{f})({d_{i}}-\bar{d})$ is the covariance between the objective values F and distances D, where ${\sigma _{F}}$, ${\sigma _{D}}$, $\bar{f}$, and $\bar{d}$ are the standard deviations and means of F and D, respectively. The value of $\mathcal{F}$ is in the range $[-1,+1]$, with a lower value implying that the analysed problem is an easy maximization problem while a higher value implying that the analysed problem is an easy minimization problem. This means that the optimization problem can be easily solved if the objective (fitness) values of the solutions are similar and their distances are small.
Shen and He (2010) also presented an FLA method for measuring terrain fitness called the “number of optimal values,” which is computed from the fitness values of individuals in the population after sorting. Assume that the population is $s=({s_{1}},{s_{2}},\dots ,{s_{n}})$, where ${s_{i}}$ is the i-th solution and n is the number of solutions in s. Further assume that each solution in the population is d-dimensional (i.e. contains d subsolutions); that is, ${s_{i}}=({s_{i,1}},{s_{i,2}},\dots ,{s_{i,d}})$, where ${s_{i,j}}$ represents the j-th dimension of the i-th solution in the population. We can then compute the distances between all solutions and the best one. All these distances will then be sorted in ascending order, yielding the sequence $\{{k_{1}},{k_{2}},\dots ,{k_{n}}\}$. This means that each solution ${s_{i}}$ is associated with an order index ${k_{j}}$ with $i\ne j$, so that the value of ${k_{j}}$ is i. Subsequently, the fitness values of individuals are compared in their sorted order. In case that ${f_{({k_{i+1}})}}\leqslant {f_{({k_{i}})}}$, the counter χ will be incremented by one. After completing the comparison of all individuals, this count is normalized by the population size as follows:
where ϕ can be used as an approximation of terrain roughness. When $\phi =0$, the structure is closer to unimodal, which makes the problem relatively easy to solve; conversely, when $\phi =1$, there are many “good” solutions in the region, resulting in terrain with varying elevations and greater problem complexity.
2.3 Fitness Landscape Analysis for Metaheuristics
Many studies (Tsai and Chiang, 2023; Talbi, 2021) found that metaheuristic and machine learning algorithms can complement each other; thus, the demarcation lines between these two different learning algorithms are becoming increasingly obscure and ambiguous. For example, we can use a metaheuristic algorithm to enhance the performance of a machine learning algorithm, such as finding a set of suitable hyperparameters for the deep learning algorithm or generating a suitable deep learning architecture (Tsai and Chiang, 2023). On the contrary, a machine learning algorithm can also be used to improve the performance of a metaheuristic algorithm, such as finding a suitable search operator for the metaheuristic algorithm in different search states or optimization problems. Talbi (2021) divides the ways of using machine learning as a supplementary mechanism for a metaheuristic algorithm into two categories: “offline” and “online,” respectively. Most of them imply that we can first use machine learning to roughly know the data distribution, problem characteristics, or even landscape of solution space and suitable search operator for different search states and then use the trained model to help the metaheuristic algorithm automatically determine the appropriate search operator in solving an optimization problem.
2.3.1 Offline Learning for the Metaheuristic Algorithm
Using machine learning as an offline mechanism is one of the possible solutions for the metaheuristic algorithms (Merz and Freisleben, 1999; Shen and He, 2010; Wang et al., 2018; Tan et al., 2021; Talbi, 2021; Tan et al., 2022b; da Costa Oliveira et al., 2023; Zhou et al., 2023). Figure 3 provides a simple example to illustrate how a machine learning algorithm is used as an offline mechanism to enhance the performance of a metaheuristic algorithm. In this example, the training process contains the “Input1” as training functions and “TA1” and “TA2” as steps for constructing the trained model, while the test process contains the “Input2” as the unseen functions and “TR” as steps for using the metaheuristic algorithm based on the trained model to determine the suitable search strategy to solve a new optimization problem.
In this way, the learning algorithm will first train the model via the input training problems from TA1 to TA2. Then, the metaheuristic algorithm will use the trained model to determine the suitable search strategy (e.g. search operator) for solving the new test problem and output the final result. Using the trained model of a machine learning algorithm to generate the elite solutions as the initial seeds of the metaheuristic algorithm and to select suitable search operators for designing a good metaheuristic algorithm are two representative research directions of the offline category (Talbi, 2021). Among them, using the information of the fitness landscape of the solution space to enhance the ability of machine learning for improving the end results of a metaheuristic algorithm has been a promising research topic in recent years. In an early study (Merz and Freisleben, 1999), Merz and Freisleben discussed the information of FLA that can be used to design components of a memetic algorithm in solving combinatorial optimization problems, as well as other metaheuristic algorithms. In Shen and He (2010), Shen and He attempted to use a training model with fitness landscape information to let evolutionary programming understand the relationship between the local fitness landscape and mutation strategies to determine suitable mutation operators for different search states during the convergence process. Wang et al. (2018) used the two fitness landscape measurements to measure the evolutionary probability (i.e. the population can evolve to a fitter state) and evolutionary ability (i.e. degree of evolution, fitness improvements, and diversity) of the population to adjust the search strategy during the convergence process.
Tan et al. (2021) presented a DE with an adaptive mutation operator based on the information of fitness landscape, called FLDE. The algorithm uses the random forest algorithm to train a learning model to evaluate which mutation operator is suitable for using which fitness landscape characteristic. So FLDE will be based on the fitness landscape characteristic of an optimization problem to further determine a suitable mutation strategy before the search process begins. In a later study, Tan et al. (2022b) further used the k-nearest neighbours (kNN) algorithm for constructing the relationship between dynamic fitness landscape characteristics (i.e. distance for each individual to the optimal global solution and information entropy to measure ruggedness of solution space) and mutation operators for determining the suitable mutation operators during the search process. Zhou et al. (2023) also attempted to use the FLA information to enhance the searchability of the artificial bee colony called FLABC. In this algorithm, the population can be regarded as a sample of the fitness landscape while the dispersion metric (i.e. average pairwise Euclidean distance between the best searched solutions) will be used to identify the landscape characteristics. Based on the result of the dispersion metric, FLABC can then know whether the local landscape is smooth or rugged. If the observation of the local landscape is smooth, the search behaviour of FLABC will then move toward “exploitation.” In case the observation of the local landscape is rugged, the search behaviour will tend to “exploration.”
2.3.2 Online Learning for the Metaheuristic Algorithm
In addition to using machine learning as an offline mechanism to enhance metaheuristic algorithm performance, it can also serve as an online mechanism to provide supplemental information for a metaheuristic algorithm to select a suitable transition operator and update the model incrementally during the convergence process, which is an alternative approach for metaheuristic algorithms (Zhang et al., 2023; Seyyedabbasi, 2023; Tessari and Iacca, 2022; Tan et al., 2022a). Figure 4 provides a simple example to illustrate how a machine learning algorithm is used as an online mechanism to enhance the performance of a metaheuristic algorithm. In this example, the “Input” is the unseen functions, and the TA step represents a learning model. When the metaheuristic algorithm is used to solve a new optimization problem, it will use the learning model to determine the suitable search strategy (i.e. the TR step). Then, the search results will be used to update the learning model and make the recommendation of the learning model more accurate in later iterations. In this example, of course, the learning model can also be trained by using the given training problems before it is used for the metaheuristic algorithm; however, the main focus is on incrementally updating the learning model with the current searched results from the metaheuristic algorithm.
In this online category, reinforcement learning (RL) is a representative learning algorithm (Zhang et al., 2023; Seyyedabbasi, 2023; Tessari and Iacca, 2022) because the RL model can be used to help the metaheuristic algorithm determine the appropriate search strategy (e.g. search operator) in different periods of the convergence process. This means that the RL model will be learned and updated based on the new searched results during the convergence process. In Seyyedabbasi (2023), Seyyedabbasi utilized a Q-learning algorithm (a variant of the reinforcement learning algorithm) to help sand cat swarm optimization algorithm (SCSO) tend to the “exploration” or “exploitation” on different search states of the convergence process. Tessari and Iacca (2022) presented a general-purpose framework by using the reinforcement learning algorithm for performing parameter adaptation in metaheuristic algorithms, such as step-size of covariance matrix adaptation evolution strategies (CMA-ES) and crossover mutation of differential evolution (DE).
2.4 Summary
The above discussion demonstrated that FLA methods provide several means to characterize optimization problems, such as whether the landscape is smooth or rugged. With the fitness landscape information, we can then provide a metaheuristic algorithm some additional “vision” to better understand the landscape of the solution space. Therefore, it might be able to enhance the searchability of the metaheuristic algorithm, and that is why some recent studies (Merz and Freisleben, 1999; Shen and He, 2010; Wang et al., 2018; Tan et al., 2021, 2022b; Zhou et al., 2023) attempted to combine the FLA with the metaheuristic algorithm for solving optimization problems. However, even though measurements of FLA might be able to provide us some possible ways to better understand the landscape of the solution space, they still have limitations because no FLA method can exactly depict the whole solution space. This means that most FLA methods can only roughly depict the whole/local solution space, but it can still provide some additional fitness landscape information to enhance the search performance of the metaheuristic algorithm compared with the metaheuristic algorithm without any fitness landscape information. That is why some studies (Wang et al., 2018; Tan et al., 2021) attempted to use two or more FLA methods to understand the solution space.
In addition to the accuracy of FLA, how to use this information to select the suitable search algorithm or operator for the metaheuristic algorithm is another critical research issue that has attracted the attention of researchers in recent years. This research topic can be regarded as knowing which algorithm or operator will suit the detected local fitness landscape. Of course, we can expect that a good way to map the search operator with the fitness landscape will enhance the metaheuristic algorithm performance if we have many search operators and many kinds of local fitness landscapes of the solution space. Among them, most combinations of metaheuristic algorithms with FLA determine the search strategy (or operator) before the search process begins (Tan et al., 2021), using the whole solution space information to further adjust the search strategy (Tan et al., 2022b) or using a predefined rule to determine the search strategy, as in Zhou et al. (2023). That is why this research focuses on developing an effective approach to use a machine learning algorithm based on the FLA information of partial solution space for the metaheuristic algorithm to adaptively determine the suitable search operator for different search states during the convergence process.
3 Proposed Algorithm
3.1 Basic Idea
This paper introduces an effective differential evolution algorithm with an adaptive mutation operator based on the information from the reinforcement learning model and FLA, called reinforcement learning based differential evolution with fitness landscape analysis (RLDEFLA). It is inspired by the studies of Tan et al. (2022b), Tan and Li (2021) because both use FLA to improve the performance of metaheuristic algorithms. In the studies of Tan et al. (2022b), Tan and Li (2021), an ML or RL algorithm is used to construct the relationship between landscape characteristics and mutation operators in the offline phase (training process). The trained model generated by the offline phase (training process) is then used to determine suitable mutation operators for DE before or during the convergence process. However, according to our observations, both studies (Tan et al., 2022b; Tan and Li, 2021) only use the random walk procedure to understand the fitness landscape of the whole solution space; the designs are not for the partial solution space around the current searched solution. We believe that the search strategy adjustment of DE during the convergence process needs to be based on the fitness landscape information around the current search solution to make such adjustments more accurate in solving an optimization problem. It is well known that the fitness landscape of the solution space might contain a couple of different landscapes in different regions. Hence, the basic idea of the proposed algorithm is to adjust the search strategy when the metaheuristic algorithm confronts different fitness landscape states (i.e. search state) in the “particular region of the solution space” instead of the “whole solution space.” This means that the search state information will be used to determine whether the landscape needs to re-analyse the landscape characteristic and change the search operator. Furthermore, the proposed algorithm will dynamically adjust the number of steps in the next random walk based on the FLA results obtained from the previous random walk, either to enable more random walks for comprehensive exploration when the landscape is complicated or to reduce the number of random walks when the landscape is simpler.
Figure 5 shows the system design of the proposed algorithm, RLDEFLA. It can be divided into two phases: offline and online. The first phase is the offline and can be regarded as a “supervised learning” process, whereas the proposed algorithm employs reinforcement learning during training. In this phase, RLDEFLA will use the random walk procedure from a set of training functions (i.e. Input 1) to understand which search strategy will be useful for which fitness landscape (i.e. TA1). More precisely, TA1 represents the landscape understanding procedure which contains the following modules: Initialization(), LandscapeAnalysis(), StrategySelection(), BoundaryReduction(), and ModelUpdate(). Subsequently, in TA2, the proposed algorithm will train a deep Q-learning (DQN) model with the results from TA1 via reinforcement learning. The goal of this step is to learn a better strategy for selecting search strategies under different landscape conditions.
The second phase is the online phase which can be regarded as an “unsupervised search” process. The new test function (i.e. Input 2) and learning model from phase 1 will be used as the inputs to the metaheuristic algorithm of the online phase. Note that the DE will be used as the metaheuristic algorithm in this study. Both inputs (i.e. Input 2 and learning model of TA2) will be used to select a suitable search operator (e.g. mutation operator of the DE) in the TR2. Moreover, the proposed algorithm then uses this mutation operator to update the population and uses a size reduction mechanism to enhance search performance during TR3 and TR4. At the end of the convergence process of this phase, the proposed algorithm will output a good result for the input test problem.
3.2 Phase 1: Offline Supervised Learning
Algorithm 1 shows the offline supervised learning phase of the proposed algorithm that contains the following operators: Initialization(), LandscapeAnalysis(), StrategySelection(), BoundaryReduction(), and ModelUpdate().
In addition to inputting the information of training problems to the Initialization() operator, a set of different operators is also randomly generated for the solution space by this operator. The LandscapeAnalysis() operator analyses the landscape characteristics based on the current search region. The StrategySelection() operator will use a learning model (e.g. DQN model) to select a suitable search operator, such as a mutation operator. The BoundaryReduction() operator is tasked to narrow down the analysis region based on the search state. The ModelUpdate() operator will then update the weights of the learning model using the reward from recent search results. All operators in the construction method of the landscape search strategy are described in detail below.
3.2.1 Initialization()
This operator randomly generates a set of solutions $s=({s_{1}},{s_{2}},\dots ,{s_{n}})$ in the solution space, where n is the number of individuals in the population. Then, the objective value of each individual is calculated using the objective function. The analysis upper bound ${U^{a}}$ and analysis lower bound ${L^{a}}$ are set equal to the upper bound U and lower bound L of the input training problems,1 respectively.
3.2.2 LandscapeAnalysis()
The main goal of this operator is to analyse the landscape characteristics of input training problems. It employs the random walk algorithm to sample the solution space to create a random walk path. The number of steps ${n_{s}}$ and the step size α need to be specified before using the random walk algorithm. The step size is defined as:
where ${\alpha _{j}}$ is the step size of the j-th dimension; $\bar{\alpha }$ is set equal to 10.0 according to Tan et al. (2021); ${U_{j}^{a}}$ and ${L_{j}^{a}}$ denote, respectively, the analysis upper bound and the analysis lower bound. A random walk path $P=\{{p_{1}},{p_{2}},\dots ,{p_{{n_{s}}}}\}$ is generated as follows:
where ${p_{i,j}}$ and ${p_{i+1,j}}$ are the values of the j-th dimension of the solution in the random walk path. If a solution is not in the analysis range during the generation of a random walk path, it will be adjusted as follows:
After generating the random walk path, the landscape characteristics are calculated. The encoding of the landscape characteristic vector ν is defined as:
where g represents the current generation, τ represents the autocorrelation value, ϕ is the number of optimal values, ρ is the fitness distance correlation value (Jones and Forrest, 1995), and ${R_{f}}$ denotes the maximum ruggedness of information entropy value (Malan and Engelbrecht, 2009). The landscape characteristics for each random walk path are computed, with each element in the landscape characteristic vector ν representing the average value of the corresponding metrics. These generated landscape characteristics serve as features for the dataset.
(10)
\[ {p_{i+1,j}}=\left\{\begin{array}{l@{\hskip4.0pt}l}{p_{i+1,j}}-({U_{j}^{a}}-{L_{j}^{a}}),\hspace{1em}& \text{if}\hspace{2.5pt}{p_{i+1,j}}\gt {U_{j}^{a}},\\ {} {p_{i+1,j}},\hspace{1em}& \text{otherwise}.\end{array}\right.\]3.2.3 StrategySelection()
In this study, one of the three mutation strategies will be selected using the DQN model, as defined in Eq. (12). The mutation strategy ${V^{\ast }}$ will be randomly selected from the three options if the random value φ is less than the predefined threshold ${t_{h}}$; otherwise, the mutation strategy will be selected by using DQN with the highest weight in its model.
where φ is a random value in $[0,1]$, ${t_{h}}$ is a predefined threshold (it is set equal to 0.05 in this study), and $\mathcal{V}\in \{{V^{1}},{V^{2}},{V^{3}}\}$. The three different mutation operators (i.e. DE/rand/1, DE/current-to-rand/1, and DE/current-to-pbest/1) are defined as
where ${X_{i,j}}$ is the j-th dimension of the i-th individual; ${X_{pbest}}$ is randomly chosen from the top $p\% $ of individuals; ${V_{i,j}}$ is the mutant vector of ${X_{i,j}}$; ${X_{{r_{1}},j}}$, ${X_{{r_{2}},j}}$, and ${X_{{r_{3}},j}}$ where ${r_{1}}\ne {r_{2}}\ne {r_{3}}$ are individuals randomly selected from the population. These three different mutation operators will search for a fixed number of generations I. Each mutation operator can obtain the R search results after being carried out for R runs.
(12)
\[ {V^{\ast }}=\left\{\begin{array}{l@{\hskip4.0pt}l}{V^{\text{rand}}},\hspace{1em}& \text{if}\hspace{2.5pt}\varphi \lt {t_{h}},\\ {} \underset{\mathcal{V}}{\text{arg\hspace{0.17em}max}}\hspace{0.1667em}Q(s,a),\hspace{1em}& \text{otherwise,}\end{array}\right.\](13)
\[\begin{aligned}{}{V_{i,j}^{1}}& ={X_{{r_{1}},j}}+{f_{i}}\cdot ({X_{{r_{2}},j}}-{X_{{r_{3}},j}}),\end{aligned}\]3.2.4 BoundaryReduction()
To imitate the converging process that the search region will tend to a smaller region gradually, this operator will reduce the analysis region to analyse the landscape characteristics of the particular region in detail. This operator updates the minimum and maximum values of each dimension in the entire population for the current analysis region. The maximum and minimum values will be, respectively, the analysis upper and lower bounds of the corresponding dimension defined as follows:
where ${L_{j}^{a}}$ and ${U_{j}^{a}}$ represent the analysed results of ${L_{j}}$ and ${U_{j}}$, respectively, and they will be used to control the lower and upper bounds of the j-th dimension of searches of the proposed algorithm, and ${X_{i,j}^{g}}$ represents the j-th dimension of the i-th individual at the g-th generation.
3.2.5 ModelUpdate()
The design of this operator is used to update the DQN weights. It means that after each DE generation, the evolutionary efficiency ${e_{p}}$ is computed from the fitness of new solutions, which will be used as a reward in reinforcement learning to update the DQN model. In this operator, ${\chi _{i}}$ is a survival counter for the i-th individual, as shown in Eq. (18) (Tan et al., 2022a). If the fitness value of a trial vector (a new individual) $f({V_{i}})$ is better than the fitness value of its parent vector2, its survival count ${\chi _{i}}$ will be incremented by one; otherwise, it is reset to zero.
As shown in Eq. (19), ${e_{p}}$ denotes the evolutionary efficiency of the population, defined as the ratio of the total survival count of all individuals to the current population size n.
Then, the reward is ${R_{t+1}}={e_{p}}$. With such information, we can then compute the reward to update the DQN model, as follows:
where $Q({S_{t}},{A_{t}})$ represents the behavioural value for the action ${A_{t}}$ in state ${S_{t}}$, α is the learning rate, γ is the discount factor to measure the importance of future rewards, the maximum value of $Q({S_{t+1}},a)$ represents the maximum expected value corresponding to all possible actions in the next state ${S_{t+1}}$.
(18)
\[ {\chi _{i}}=\left\{\begin{array}{l@{\hskip4.0pt}l}{\chi _{i}}+1,\hspace{1em}& \text{if}\hspace{2.5pt}f({V_{i}})\leqslant f({X_{i}}),\\ {} 0,\hspace{1em}& \text{otherwise.}\end{array}\right.\](20)
\[ Q({S_{t}},{A_{t}})\gets (1-\alpha )Q({S_{t}},{A_{t}})+\alpha \Big({R_{t+1}}+\gamma \underset{a}{\max }Q({S_{t+1}},a)\Big),\]3.3 Phase 2: Online Unsupervised Search
After training the machine learning model of phase 1, the proposed algorithm will use it to recommend the mutation operator of DE during the search process of phase 2. As shown in Algorithm 2, the search process of online unsupervised search algorithm starts with the Initialization() operator and contains three operators: StrategySelection(), PopulationUpdate(), and SizeReduction().
The Initialization() operator randomly generates a population of solutions in the solution space. The StrategySelection() operator is responsible for determining the recommended mutation operator. The trained model is used in this operator. The PopulationUpdate() operator is used to update the population using the recommended mutation operator, which is determined by the StrategySelection() operator. The SizeReduction() operator takes care of reducing the population size and improving the local search ability. Note that Algorithm 2 is built upon LSHADE, but instead of its original mutation operator, we employ three different mutation operators, as shown in Eqs. (13), (14), and (15), and remove the external archive. We refer to this new algorithm as the modified LSHADE.
3.3.1 Initialization()
In the initialization operator, the population $X=\{{X_{1}},{X_{2}},\dots ,{X_{n}}\}$ is randomly generated in the solution space, as follows:
where ${X_{i,j}}$ is the j-th dimension of the i-th individual; ${L_{j}}$ and ${U_{j}}$ are the lower and upper bounds, respectively. Moreover, each mutation operator maintains its historical table H to generate the control parameters. All the entries in the historical tables are initialized to 0.5 based on Tanabe and Fukunaga (2013).
3.3.2 StrategySelection()
This operator identifies recommended mutation operators for updating the population. To maintain the diversity of the algorithm, the proposed algorithm randomly selects a mutation operator from three mutation operators once after a fixed number of generations. When the current generation g equals a fixed number of generations, the trained model recommends a mutation operator. For analysing the landscape precisely, the analysis upper bound ${U^{a}}$ and the analysis lower bound ${L^{a}}$ are defined by the maximum and minimum values of the population. The random walk algorithm then generates a random walk path in the analysis range. Different from the original random walk algorithm, the starting point of the random walk path is set to the best solution so far, and the random walk step is calculated as follows:
where ${\alpha _{j}}$ is the j-th dimension of the step size, ${U_{j}^{a}}$ and ${L_{j}^{a}}$ are, respectively, the analysis upper and lower bounds, and ${\gamma _{c}}$ is a control value. Note that for Eqs. (21) and (22), ${L_{j}}$ and ${U_{j}}$ are the lower and upper bounds of the j-th dimension of the solution space of an optimization problem that are inputs by the problem descriptions, and both of them are fixed, while ${L_{j}^{a}}$ and ${U_{j}^{a}}$ represent the analysed results for ${L_{j}}$ and ${U_{j}}$ that are used to control the search scope of the proposed algorithm during the convergence process. Note also that in this study, the number of evaluations due to all random-walk sampling is counted and included in the overall function evaluation budget (MaxFES) of the proposed algorithm to ensure a fair comparison with other search algorithms. Landscape characteristics, calculated using a random walk path, are composed into a vector, as shown in Eq. (11), to be used as the input for the model. After inputting the landscape characteristic vector to the trained model, it will output the recommended mutation operator to update the solution. The best solution ${V^{\ast }}$ of random walk path ${p^{\prime }}$ will be the criterion to decide whether the landscape characteristics change, as shown in Eq. (23). If the worst individual ${X_{w}}$ in the population is better than ${p^{\prime }}$, the landscape needs to be analysed again. The analysis result will be input into the model to recommend the mutation operator.
where $\mathcal{V}\in \{{V^{1}},{V^{2}},{V^{3}}\}$ represents a set of mutation strategies and ${V^{\ast }}$ has the maximum weights of the DQN. With this model, the proposed algorithm will dynamically adjust the number of steps in the random walk based on the analysis results from the previous FLA. In other words, it will invest more random walk steps in the complex solution space to achieve more comprehensive exploration, and fewer steps in simpler solution-space landscapes (e.g. a flat landscape) to reduce computational cost.
3.3.3 PopulationUpdate()
This operator uses the recommended mutation operator to generate the mutant vector V and update the population. The historical memory (table) H generates the control parameters and can be defined as follows:
where ${f_{i}}$ is the scaling factor randomly generated from the Cauchy distribution with mean ${H_{{s_{r}}}^{f}}$ and variance 0.1; ${c_{i}}$ is the crossover rate randomly generated from the normal distribution with mean ${H_{{s_{r}}}^{c}}$ and variance 0.1; ${r_{i}}$ is a random number to randomly pick a value from the historical memory H for the random functions ${\text{randc}_{i}}$ and ${\text{randn}_{i}}$ to obtain the scaling factor and crossover rate for the i-th solution, and ${s_{r}}$ is a number that represents the recommended mutation operator. When the values of ${f_{i}}$ and ${c_{i}}$ are not in the range $[0,1]$, they have to be adjusted (Tanabe and Fukunaga, 2013). After generating the control parameters, the recommended mutation operator is used to generate the mutant vector. When the generated mutant vector is not in the solution space, it will be adjusted as follows:
After that, the crossover and selection operators are used to update the population. During the selection process, the control parameters ${f_{i}}$ and ${c_{i}}$ are recorded into the success parameter table as ${S^{f}}$ and ${S^{c}}$ if the trial vector is better than the corresponding individual. The success parameter table is then used to update the historical table as follows:
where $\Delta {I_{i}}$ is the improvement of the objective value of i-th individual; ${w_{i}}$, the weight of ${S_{i}^{f}}$ and ${S_{i}^{c}}$, is calculated by the improvement of the objective value; ${h_{{s_{r}}}}$ is the historical table counter of the ${s_{r}}$-th mutation operator. Moreover, in this paper, Eqs. (24) and (25) are used to generate the scaling factor ${f_{i}}$ and crossover rate ${c_{i}}$; so, both ${s_{r}}$ and ${r_{i}}$ are needed to represent the associated mutation operator and a random index to the historical memory H. Different from Eqs. (24) and (25), Eqs. (28) and (29) are used to update the historical memory; thus, they only need the index ${s_{r}}$ for the particular mutation operator, but not ${r_{i}}$ to randomly pick a value from H. The historical tables ${H_{{h_{{s_{r}}}}}^{f}}$ and ${H_{{h_{{s_{r}}}}}^{c}}$ are updated by calculating the weighted Lehmer mean value. Before the size reduction operator is invoked, the worst individual ${X_{w}}$ is recorded and used in the StrategySelection() operator.
(26)
\[ {V_{i,j}}=\left\{\begin{array}{l@{\hskip4.0pt}l}({X_{i,j}^{g}}+{L_{j}})/2,\hspace{1em}& \text{if}\hspace{2.5pt}{V_{i,j}}\lt {L_{j}},\\ {} ({X_{i,j}^{g}}+{U_{j}})/2,\hspace{1em}& \text{if}\hspace{2.5pt}{V_{i,j}}\gt {U_{j}}.\end{array}\right.\](27)
\[\begin{aligned}{}& {w_{i}}=\frac{\Delta {I_{i}}}{{\textstyle\textstyle\sum _{i=1}^{|S|}}\Delta {I_{i}}},\end{aligned}\]To prevent the proposed algorithm from stagnating at a local optimum for an extended period during the search, we design a population restart mechanism. This mechanism monitors the search status by counting the number of consecutive iterations in which no better solution is found, denoted as the number of stagnant generations ${G_{s}}$. When ${G_{s}}$ exceeds a predefined threshold ${G_{r}}$, the proposed algorithm replaces the κ worst-performing individuals in the current population with new ones generated by the same Initialization() method described earlier and resets ${G_{s}}$ to 0, thereby enabling the population to explore new search directions and enhancing its subsequent exploration capability. The value of κ is determined by the following equation.
where ${R_{r}}$ is the replacement ratio for the population.
3.3.4 SizeReduction()
The size reduction operator used in the proposed algorithm is based on the LSHADE (Tanabe and Fukunaga, 2014). In this operator, the population size is reduced throughout the search process to enhance the local search ability. The population size of the size reduction operator is defined as follows:
where $|X|$ is the population size after reduction; $|{X_{\min }}|$ is the minimum size of the population; $|{X_{\text{init}}}|$ is the initial population size; E is the number of evaluations; ${E_{\max }}$ is the maximum number of evaluations. Size reduction will delete the worst individual in the population until the size of the population is equal to $|X|$.
(31)
\[ |X|=\text{round}\bigg(\big(|{X_{\min }}|-|{X_{\text{init}}}|\big)\cdot \frac{E}{{E_{\max }}}+|{X_{\text{init}}}|\bigg),\]4 Experimental Results
4.1 Experimental Environment and Parameter Settings
All experiments were conducted on a PC with an Intel(R) i9-13900K processor (3.0 GHz, 32 MB cache, and 24 cores), 32GB of memory, and Mesa Intel(R) Graphics (RPL-S) running Linux Ubuntu 22.04.4 LTS. All the programs are written in C++ and compiled using g++ 11.4.0. The machine learning model is written in Python 3.10.12. The proposed algorithm is compared with seven algorithms: the adaptive differential evolution with optional external archive (JADE) (Zhang and Sanderson, 2009), SHADE (Tanabe and Fukunaga, 2013), LSHADE (Tanabe and Fukunaga, 2014), FLDE (Tan et al., 2021), DEDQN (Tan and Li, 2021), RLHPSDE (Tan et al., 2022a), and DFLDE (Tan et al., 2022b) on the CEC2021 and CEC2022 benchmark functions (Mohamed et al., 2020; Kumar et al., 2021). The DQN model used in the proposed algorithm is based on the deep Q-Network (Tan and Li, 2021; Mnih et al., 2015) with 2 fully connected hidden layers, 10 neurons per hidden layer, ReLU activation, 3 neurons for the output layer (each of which is a mutation strategy, i.e. DE/rand/1, DE/current-to-rand/1, and DE/current-to-pbest/1), Adam optimizer, $1\times {10^{-4}}$ learning rate, mean squared error as the loss function, and mini-batch size of 64. The CEC 2013 benchmark (Liang et al., 2013) is used as a function set to generate training samples as in Tan et al. (2021) and Tan et al. (2022b). All algorithms are carried out for 51 independent runs. The maximum numbers of evaluations ${E_{\max }}$ for 10 and 20 dimensions are 200 000 and 1 000 000 according to the settings of Mohamed et al. (2020) and Kumar et al. (2021). Note that CEC2021 contains 10 basic functions (categorized by their characteristics as unimodal, basic, hybrid, and composition functions), but they have been transformed via “shift,” “bias and shift,” “shift and rotation,” “bias, shift, and rotation” to generate another 40 functions. So there are a total of 50 test functions in CEC2021. Also note that CEC2022 contains 12 test functions.
Table 2
Parameter settings.
| Algorithm | Parameter setting |
| JADE | ${X_{\text{init}}}=100$, ${r_{\text{arc}}}=1$, $p=0.05$, $c=0.1$ |
| SHADE | ${X_{\text{init}}}=100$, ${r_{\text{arc}}}=1$, $H=100$ |
| LSHADE | ${X_{\text{init}}}=18\times $ d, ${X_{\min }}=4$, ${r_{\text{arc}}}=2.6$, $H=6$, $p=0.11$ |
| FLDE | ${X_{\text{init}}}=18\times $ d, ${X_{\min }}=4$, ${r_{\text{arc}}}=2.6$, $H=5$, ${n_{s}}=2$,000, $\alpha =10$ |
| DFLDE | ${X_{\text{init}}}=18\times $ d, ${X_{\min }}=4$, $I=100$, $H=5$, ${n_{s}}=1$,000, $\alpha =20$, $p=0.1$, $k=13$ |
| DEDQN | ${X_{\text{init}}}=20\times d$, ${X_{\min }}=7$, $I=1$, $H=8$, ${n_{s}}=200$, $\alpha =10$, $p=0.14$ |
| RLHPSDE | ${X_{\text{init}}}=18\times d$, ${X_{\min }}=4$, $I=1$, $H=5$, ${n_{s}}=200$, $\alpha =10$ |
| RLDEFLA | ${X_{\text{init}}}=22\times d$, ${X_{\min }}=10$, $I=110$, $H=15$, ${R_{n}}=3.1$, ${\gamma _{c}}={E_{\max }}/100000$, $p=0.15$, ${G_{r}}=70$, ${R_{r}}=0.05$ |
The parameters in the construction method of the landscape search strategy are set as follows: the number of runs R is set equal to 30; the maximum number of evaluations ${E_{\max }}$ is set equal to 100 × d, where d is the number of dimensions; and the step number ${n_{s}}$ is set equal to 1000. The parameter settings for all algorithms are shown in Table 2 where ${X_{\text{init}}}$ represents the initial population size, ${X_{\min }}$ the minimum population size, ${r_{\text{arc}}}$ the rate of the external archive, p the rate of the top individual, c the control parameter to adjust the crossover and mutation rates, H the size of the historical table, ${n_{s}}$ the number of steps, α the step size, I a fixed number of generations as an interval to reanalyse the fitness landscape and determine the suitable search strategy, ${R_{n}}$ the ratio of the step size, ${n_{s}}={E_{\max }}/2000\times {R_{n}}$, ${G_{r}}$ the restart threshold, and ${R_{r}}$ the replacement ratio for the population. Because the number of combinations of hyperparameters ${X_{\text{init}}}$, ${X_{\min }}$, H, p, I, ${R_{n}}$, ${G_{r}}$, and ${R_{r}}$ of RLDEFLA is excessively large, they are set by the search results of the tree-structured Parzen estimator (TPE) (Bergstra et al., 2011). The search range of these eight parameters is as given in Table 3.
Table 3
Search range of hyperparameter optimization for RLDEFLA.
| Parameters | Range |
| ${X_{\text{init}}}$ | $\{15,16,\dots ,24,25\}\cdot d$ |
| ${X_{\min }}$ | $\{5,6,\dots ,14,15\}$ |
| H | $\{5,6,\dots ,14,15\}$ |
| p | $\{0.11,0.12,\dots ,0.19,0.20\}$ |
| I | $\{50,60,\dots ,140,150\}$ |
| ${R_{n}}$ | $\{1.3,1.6,\dots ,3.7,4.0\}$ |
| ${G_{r}}$ | $\{40,50,\dots ,70,80\}$ |
| ${R_{r}}$ | $\{0.05,0.1,\dots ,0.2,0.25\}$ |
Using TPE for hyperparameter optimization typically requires a measurement to evaluate the quality of a set of parameter settings. In this study, the initial parameter settings for TPE are set to the midpoint values of the predefined search ranges for each parameter. The measurement of the improvement of a set of new parameter settings compared with initial parameter settings is defined as follows:
where ${P_{o}}$ and ${P_{t}}$ represent the average objective values by using initial parameter settings and parameter settings found by TPE for RLDEFLA, respectively.
Table 4
The initial and best hyperparameter settings of RLDEFLA.
| ${X_{\text{init}}}$ | ${X_{\min }}$ | I | H | p | ${R_{n}}$ | ${G_{r}}$ | ${R_{r}}$ | |
| Initial settings | $20\cdot d$ | 10 | 100 | 10 | 0.15 | 2.5 | 60 | 0.15 |
| Settings found by TPE | $22\cdot d$ | 10 | 110 | 15 | 0.15 | 3.1 | 70 | 0.05 |
Moreover, in this process, the TPE will be used for 500 different sets of hyperparameter settings. The initial hyperparameter settings and best hyperparameter settings found by TPE are shown in Table 4. To better understand the impact of each hyperparameter on the proposed algorithm, we further conduct single-parameter variation experiments for these eight parameters. We first vary one parameter while keeping the remaining parameters fixed, as determined by TPE. Figure 6 uses 0% as the baseline to represent the performance of the parameter combination determined by TPE; the other points show performance relative to this reference setting. These results show that the parameter values determined by TPE yield suitable performance in most cases, and that further improvement may be possible by varying a single parameter over a series of additional tests.
4.2 Ablation Study of Search Operators
This section presents the ablation study of the proposed algorithm based on the CEC2021 and CEC2022 benchmark functions. The ablation study focuses on the following mechanisms: timing to FLA (tFLA), step size (SS), and population restart (PR). As shown in Table 5, “Difference” represents the comparison difference between (1) the proposed algorithm using tFLA, SS, and PR mechanisms and (2) the proposed algorithm using some of these mechanisms, in terms of the objective value. Also, “✓” and “–” in this table represent whether the operator is used in the proposed algorithm or not, respectively. The ablation studies show the importance of the mechanisms. It can be easily seen that each of these mechanisms will affect the end result of the proposed algorithm regardless of whether they degrade the end result from 7.25% to 34.57% compared with using all of them in the proposed algorithm. The results further explain why we use these mechanisms in the proposed algorithm.
4.3 Comparisons with Other Algorithms
In this section, the CEC2021 and CEC2022 benchmarks, with 10D and 20D, are also used to test the performance of the proposed algorithm and other state-of-the-art search algorithms.
4.3.1 Evaluation Criteria for CEC2021
The evaluation criterion CEC2021 (Mohamed et al., 2020) is based on the error value and rank of a search algorithm and is defined by
where SNE denotes the average of all normalized error values over all functions and is defined as
where i and m represent the indices of the i-th function and the m-th transformation of a test function, $10d$ and $20d$ are the indices of the 10 and 20 dimensional test functions, respectively, $f({s_{e}})$ is the best result of a search algorithm out of 51 runs, $f({s_{b}})$ is the known optimal value of a function, and $f{({s_{e}})_{\max }}$ is the largest $f({s_{e}})$ among all algorithms. where r is the rank of an algorithm among all algorithms for a given test function.
4.3.2 Evaluation Criteria for CEC2022
The evaluation criterion for CEC2022 (Kumar et al., 2021) first considers accuracy and then MaxFES to determine the rank of an algorithm for all runs, which can be defined by
where ${S_{10}}$ and ${S_{20}}$ represent the sum of ranks for all their trials $\text{SR}$ for benchmarks of 10D and 20D, respectively. In case two algorithms cannot reach the minimum error value (${10^{-8}}$) before MaxFES is reached, the rank of an algorithm ${\mathcal{A}_{1}}$ will be better than that of algorithm ${\mathcal{A}_{2}}$ when the error value of ${\mathcal{A}_{1}}$ is smaller than that of ${\mathcal{A}_{2}}$. In case both algorithms can reach the minimum error value before MaxFES is reached, and the convergence of an algorithm ${\mathcal{A}_{1}}$ is faster than another algorithm ${\mathcal{A}_{2}}$, then the rank of ${\mathcal{A}_{1}}$ will be better than that of ${\mathcal{A}_{2}}$. Moreover, all the algorithms will use this approach to get the sum of ranks for all their trials, which is defined as follows
where R is the number of runs. The $\text{SR}$ can then be used to compare the performance of the search algorithms.
4.3.3 Score Results for CEC2021 and CEC2022
Tables 6 and 7 show the comparisons of the proposed algorithm (i.e. RLDEFLA) with other metaheuristic algorithms (i.e. JADE, SHADE, LSHADE), metaheuristic algorithms with FLA (i.e. FLDE and DFLDE), and metaheuristic algorithms with FLA and RL (i.e. DEDQN and RLHPSDE). The results of Table 6 show that the proposed algorithm can perform better than all the other metaheuristic algorithms compared in this study in terms of the score metric ${S_{21}}$ for CEC2021 benchmarks. SHADE achieves better results than RLDEFLA on ${S_{1}}$, which considers the average normalized error across all functions, and the difference between the best result of a search algorithm and the best-known solution. However, RLDEFLA can significantly outperform SHADE in terms of ${S_{2}}$, which is the rank of an algorithm among all algorithms for all test functions. These comparison results show that RLDEFLA can provide a balanced way for solving the SOP.
Table 6
The comparisons between the proposed algorithm and other algorithms for CEC2021.
| JADE | SHADE | LSHADE | FLDE | DFLDE | DEDQN | RLHPSDE | RLDEFLA | |
| ${S_{1}}$ | 32.89 | 50.00 | 41.53 | 44.62 | 42.15 | 34.44 | 41.92 | 41.80 |
| ${S_{2}}$ | 27.57 | 24.91 | 39.14 | 31.81 | 32.59 | 32.81 | 31.23 | 50.00 |
| ${S_{21}}$ | 60.46 | 74.91 | 80.67 | 76.43 | 74.74 | 67.25 | 73.15 | 91.80 |
Table 7
The comparisons between the proposed algorithm and other algorithms for CEC2022.
| JADE | SHADE | LSHADE | FLDE | DFLDE | DEDQN | RLHPSDE | RLDEFLA | |
| ${S_{10}}$ | 42636 | 152756 | 160329 | 105027 | 116481 | 104980 | 142668 | 123613 |
| ${S_{20}}$ | 58346 | 124541 | 163052 | 109676 | 129605 | 118809 | 103161 | 147997 |
| ${S_{22}}$ | 100982 | 277297 | 323381 | 214703 | 246086 | 223789 | 245829 | 271610 |
Table 7 further shows that the proposed algorithm can still obtain better results than metaheuristic algorithms with FLA (i.e. FLDE and DFLDE), and metaheuristic algorithms with FLA and RL (i.e. DEDQN and RLHPSDE) in terms of ${S_{22}}$, but it fails to surpass SHADE and LSHADE. Based on our observations, we attribute this to two reasons. First, FLA-based mechanisms require additional evaluations, which will degrade the score when convergence speed and results are considered together. Second, the parameter settings were determined by TPE over the 50 CEC2021 and 12 CEC2022 benchmarks, so they are strongly influenced by CEC2021. We therefore reran TPE on the CEC2022 benchmarks alone to find suitable parameter settings. The results of ${S_{10}}$, ${S_{20}}$, and ${S_{22}}$ for the proposed algorithm are 175142, 168445, and 343587, respectively, indicating that with suitable parameter settings, the proposed algorithm can achieve better results. These results are now better than SHADE and LSHADE, and significantly better than other FLA- and RL-based methods for metaheuristic algorithms, but it still has room for improvement to reduce the costs associated with the FLA.
4.3.4 Average Rank and Wilcoxon Test Results for CEC2021 and CEC2022
To compare the performance between RLDEFLA and other search algorithms, the average rank of each algorithm is shown in Fig. 7. Note that the average rank of an algorithm in this figure represents its average rank among all algorithms for all given test functions in all trials of CEC2021 or CEC2022, respectively. Also note that the “Rank” represents the average rank value of 51 independent runs of all test functions of CEC2021 (with 10D and 20D) and CEC2022 (with 10D and 20D), respectively. The results indicate that RLDEFLA achieves better performance than other algorithms for both CEC2021 and CEC2022 in most cases. More precisely, in these comparisons, the rank with a lower value means the search algorithm can get better results on average in terms of the objective value. That RLDEFLA can get lowest rank values of 3.06 and 2.78 for CEC2021 with 10D and 20D and rank values of 3.25 and 3.17 for CEC2022 with 10D and 20D than other search algorithms, respectively, explains that the proposed algorithm can find better results than other state-of-the-art search algorithms (including the search algorithms based on FLA) in most cases.
Fig. 8
The Wilcoxon test comparisons between RLDEFLA and the other algorithms for CEC2021 and CEC2022.
To further understand the performance of the proposed method, statistical analysis is used for the comparisons between the RLDEFLA and the other search algorithms. Figure 8 shows the results of the Wilcoxon test (Wilcoxon, 1945) for these algorithms. The terms “Better,” “Equal,” and “Worse” indicate that the result of RLDEFLA is significantly better than, not significantly better than (i.e. similar), and worse than the compared search algorithm, respectively. The results show that the proposed algorithm (RLDEFLA) outperforms other algorithms for both CEC2021 and CEC2022. As shown in Fig. 8, the proposed algorithm can find better results than the other seven search algorithms in most cases. For example, the results of CEC2021 with 10D show that the proposed algorithm can find better results than JADE in 28 test functions, similar results with JADE in 20 test functions, and worse results than JADE in only 2 test functions. Since the proposed algorithm is worse than JADE in 2 out of 50 test functions while the number of functions the proposed algorithm can get better results is more than the number of functions the proposed algorithm can get worse results (28 better vs. 2 worse), it can be easily seen that the proposed algorithm can find better results than JADE in most cases. Of course, similar situations can also be found when comparing the proposed algorithm with JADE, SHADE, LSHADE, FLDE, DFLDE, DEDQN, and RLHPSDE for CEC2021 with 10D and 20D.
Also shown in Fig. 8 are the comparisons between the proposed algorithm and the other seven search algorithms for the CEC2022 with 10D and 20D. Even though the number of test functions is not like that of CEC2021, which has 50, CEC2022 still has 12 test functions. From these comparisons, it can also be seen that the number of test functions the proposed algorithm can get better results than the other search algorithms is still more than the number of functions the proposed algorithm only gets worse results like the comparisons in CEC2021. For example, the proposed algorithm can get better results in 8 test functions with 10D, similar results in 3 test functions, and worse results in only 1 test function in the comparisons between the proposed and JADE. In the other comparisons of Fig. 8, the proposed algorithm still can get better or similar results in most cases. It is worth noting that the proposed algorithm does not achieve better results on CEC2022 with 10D than LSHADE, even though the results are very similar. We believe that it is because of the parameter setting of the proposed algorithm, as mentioned in Section 4.3.3. As an alternative approach, based on the same procedure as in Section 4.3.3, we can further apply TPE solely to the CEC2022 benchmarks to determine suitable parameter settings. As a result, the proposed algorithm outperforms LSHADE on 3 test functions, performs comparably on 7, and performs worse on 2.
4.3.5 Average Rank and Wilcoxon Test Results for CEC2020
In this study, we further compare the proposed algorithm with other metaheuristic algorithms on the CEC2020 benchmark (Yue et al., 2019) to assess its generalizability. Note that the maximum numbers of evaluations ${E_{\max }}$ for 5, 10, 15, and 20 dimensions are 50,000, 1,000,000, 3,000,000, and 10,000,000 for the CEC2020 according to the settings of Yue et al. (2019). As shown in Fig. 9, the proposed algorithm can still obtain better results than all the other metaheuristic algorithms for CEC2020, as measured by the metric from 5D to 20D.
Figure 10 further shows the comparison between the proposed algorithm and other metaheuristic algorithms for CEC2020 using the Wilcoxon test to determine whether the results of the proposed algorithm are better than those of the other metaheuristic algorithms. These results show that the proposed algorithm can find better results than all the other metaheuristic algorithms, demonstrating its generalizability for solving various SOPs. Note that from this comparison, although it can be easily seen that the proposed algorithm retains its advantage on the CEC2020 dataset without per-benchmark tuning, demonstrating that it has a good generalizability for CEC2020, CEC2021, and CEC2022, we cannot argue that the proposed algorithm will beat all the other algorithms compared in this study for all the other benchmarks to which it has not been applied. Instead, we can only argue that the proposed algorithm has a good generalizability for solving various SOPs, and it may provide better results than other search algorithms because the proposed algorithm embeds a good mechanism to automatically and dynamically pick a good transition operator from a set of operators to generate new candidate solutions during the convergence process. Based on these observations, we conclude that even if the proposed algorithm achieves good search performance on new optimization problems, it still needs to be modified to achieve optimal performance when the objectives, constraints, or input data distribution differ significantly from those of the test benchmarks used in this study.
5 Conclusions
This paper introduces a hybrid mutation operator for differential evolution based on the fitness landscape to improve search performance. The proposed algorithm represents an effective integration of supervised learning and unsupervised learning algorithms. The basic idea of the proposed algorithm is to first train a model using the known relationships between mutations and landscape characteristics. It can be viewed as a form of supervised learning. Then, this trained model is used for DE to adaptively determine the mutation operator for unseen single-objective optimization problems during the convergence process. It can be viewed as a form of unsupervised learning. The experimental results demonstrate that the proposed method improves the search performance on the CEC2020, CEC2021 and CEC2022 benchmarks, achieving better results than other DE-based algorithms for most functions. One of the important findings is that if we can have an accurate strategy detection model to dynamically choose the suitable search operator (e.g. mutation operator for DE) during the convergence process, the end results of the metaheuristic algorithm can then be significantly improved compared with those of a metaheuristic algorithm without such mechanism. Using a machine learning algorithm is a viable approach to generate such a strategy detection model for a metaheuristic algorithm. The comparison results of this study show that this can be a good solution to enhance the searchability of a metaheuristic algorithm. Further studies could investigate the potential of combining the proposed algorithm with other search algorithms. Also, future work could explore an effective alternative to the random walk for sampling fitness landscape characteristics.
