A hybridizing-enhanced differential evolution for optimization

Overview of DE and GWO algorithms

This section, at first, briefly presents the framework of DE and GWO algorithms, and then the proposed algorithm (HDE) is introduced by combining these two algorithms.

Differential evolution

The original version of the differential evolution algorithm was introduced by Storn & Price (1996). Das, Mullick & Suganthan (2016) has provided a comprehensive description of many applications of the DE algorithm and its improved versions, which appeared in recent years. Like other evolutionary algorithms, an initial population is generated at first. Then, some operators, including composition, mutation, and crossover, are applied to the population to form a new population. In the next step, i.e., the selection phase, the new population is compared to the current population based on objective function merits. Then, the best members enter the next stage as the next generation. This process continues until the desired results are achieved. This section describes the performance steps of the HDE algorithm.

From a mathematical point of view, a population can be represented using a vector of population members (i.e., search agents) $\left({\overrightarrow{X}}_1^{Iter},{\overrightarrow{X}}_2^{Iter},\dots ,{\overrightarrow{X}}_{N_{pop}}^{Iter}\right)$, where the $i$-th population member ${\overrightarrow{X}}_i^{Iter}$ has the form ${\overrightarrow{X}}_i^{Iter}=\left(x_{i,1}^{Iter},x_{i,2}^{Iter},\dots ,x_{i,D}^{Iter}\right)$. Here, $x_{i,j}^{Iter}$ is the $j$-th dimension of the $i$-th search agent ${\overrightarrow{X}}_i^{Iter}$ (i.e., the $j$-th decision variable), $N_{pop}$ is the number of search agents in a population, and $D$ is the number of decision variables (the dimension of the problem). The indices $i$ and $Iter$ are characterized as $Iter=1,2,\dots ,Ite{r}_{max}$ and $i=1,2,\dots ,N_{pop},$ where $Ite{r}_{max}$ is the maximal number of iterations.

1. Initial population generation

The minimal and maximal values for the $j$-th decision variable ( $j=1,2,\dots ,D$) of the problem are denoted by $x_{min,j}$ and $x_{max,j}$, respectively. The initial population (i.e., $Iter=1$) with the total number of members of the population $N_{pop}$ in $D$ dimensions can be generated as Eq. (1):

(1) $$x_{i,j}^{Iter}=ran{d}_j\left(0,1\right)\cdot \left(x_{max,j}-x_{min,j}\right),$$where $ran{d}_j\left(0,1\right)$ is a random function that generates a real number with a uniform probability distribution between 0 and 1 (i.e., $ran{d}_j\left(0,1\right)\in U\left[0,1\right]$).

2. Mutation

In the DE algorithm, different strategies can be used for mutating and creating a new population Mininno et al. (2010) and Ghasemi et al. (2016). The mutation operator selects random vectors for each possible range in the starting population named ${\overrightarrow{X}}_r^{Iter}$ (like ${\overrightarrow{X}}_{r_1}^{Iter}$, ${\overrightarrow{X}}_{r_2}^{Iter}$, ${\overrightarrow{X}}_{r_3}^{Iter}$, and ${\overrightarrow{X}}_{r_4}^{Iter}$) from starting population matrix and generates a new vector ${\overrightarrow{V}}_i^{Iter}=\left(v_{i,1}^{Iter},v_{i,2}^{Iter},\dots ,v_{i,D}^{Iter}\right)$ for the $Iter$-th iteration of the algorithm and the $i$-th member of the population based on the utilized mutation approach that creates a new population. Randomly chosen vectors ${\overrightarrow{X}}_r^{Iter}$must be dissimilar from each other. The resultant vector ${\overrightarrow{V}}_i^{Iter}$ produced by the mutation operator has one of the following configurations as presented in Eqs. (2) to (9) Mininno et al. (2010) and Ghasemi et al. (2016):

“DE/rand/1”:

(2) $${\overrightarrow{V}}_i^{Iter}={\overrightarrow{X}}_{r_1}^{Iter}+F\cdot \left({\overrightarrow{X}}_{r_2}^{Iter}-{\overrightarrow{X}}_{r_3}^{Iter}\right),$$

“DE/best/1”:

(3) $${\overrightarrow{V}}_i^{Iter}={\overrightarrow{X}}_{best}^{Iter}+F\cdot \left({\overrightarrow{X}}_{r_1}^{Iter}-{\overrightarrow{X}}_{r_2}^{Iter}\right),$$

“DE/current-to-best/1”:

(4) $${\overrightarrow{V}}_i^{Iter}={\overrightarrow{X}}_i^{Iter}+F\cdot \left({\overrightarrow{X}}_{best}^{Iter}-{\overrightarrow{X}}_i^{Iter}\right)+F\cdot \left({\overrightarrow{X}}_{r_1}^{Iter}-{\overrightarrow{X}}_{r_2}^{Iter}\right),$$

“DE/rand/2”:

(5) $${\overrightarrow{V}}_i^{Iter}={\overrightarrow{X}}_{r_1}^{Iter}+F\cdot \left({\overrightarrow{X}}_{r_2}^{Iter}-{\overrightarrow{X}}_{r_3}^{Iter}\right)+F\cdot \left({\overrightarrow{X}}_{r_4}^{Iter}-{\overrightarrow{X}}_{r_5}^{Iter}\right),$$

“DE/best/2”:

(6) $${\overrightarrow{V}}_i^{Iter}={\overrightarrow{X}}_{best}^{Iter}+F\cdot \left({\overrightarrow{X}}_{r_1}^{Iter}-{\overrightarrow{X}}_{r_2}^{Iter}\right)+F\cdot \left({\overrightarrow{X}}_{r_3}^{Iter}-{\overrightarrow{X}}_{r_4}^{Iter}\right),$$

“DE/rand-to-best/1”:

(7) $${\overrightarrow{V}}_i^{Iter}={\overrightarrow{X}}_{r_1}^{Iter}+F\cdot \left({\overrightarrow{X}}_{best}^{Iter}-{\overrightarrow{X}}_{r_1}^{Iter}\right)+F\cdot \left({\overrightarrow{X}}_{r_2}^{Iter}-{\overrightarrow{X}}_{r_3}^{Iter}\right),$$

“DE/rand-to-best/2”:

(8) $${\overrightarrow{V}}_i^{Iter}={\overrightarrow{X}}_{r_1}^{Iter}+F\cdot \left({\overrightarrow{X}}_{best}^{Iter}-{\overrightarrow{X}}_{r_1}^{Iter}\right)+F\cdot \left({\overrightarrow{X}}_{r_2}^{Iter}-{\overrightarrow{Y}}_{r_3}^{Iter}\right)+F\cdot \left({\overrightarrow{X}}_{r_4}^{Iter}-{\overrightarrow{X}}_{r_5}^{Iter}\right),$$

“DE/current-to-rand/1”:

(9) $${\overrightarrow{V}}_i^{Iter}={\overrightarrow{X}}_i^{Iter}+rand\left(0,1\right)\cdot \left({\overrightarrow{X}}_{r_1}^{Iter}-{\overrightarrow{X}}_i^{Iter}\right)+F\cdot rand\left(0,1\right)\cdot \left({\overrightarrow{X}}_{r_2}^{Iter}-{\overrightarrow{X}}_{r_3}^{Iter}\right),$$where $rand\left(0,1\right)\in U\left[0,1\right]$. The mutation operator has a parameter $F$, which is uniformly randomly picked in the range from 0 to 2 (i.e., $F\in U\left[0,2\right]$), and $\left({\overrightarrow{X}}_{r_1}^{Iter}-{\overrightarrow{X}}_{r_2}^{Iter}\right)$, $\left({\overrightarrow{X}}_{r_2}^{Iter}-{\overrightarrow{X}}_{r_3}^{Iter}\right),$ and $\left({\overrightarrow{X}}_{r_4}^{Iter}-{\overrightarrow{X}}_{r_5}^{Iter}\right)$ are diverse vectors that modify the base vector. Population members are recognized by ${\overrightarrow{X}}_r^{Iter}$, and the best individual vector with the best fitness merit in the present population at iteration $Iter$ is denoted by ${\overrightarrow{X}}_{best}^{Iter}$.

3. Crossover

This operator increases the algorithm’s strength and escapes from local optimal solutions. The crossover operation, with the constant of crossover $CR$ ( $0<CR<1$), is performed for each $j$-th decision variable of the $i$-th search agent of the population at iteration $Iter$, to compute the trial vector ${\overrightarrow{U}}_i^{Iter}=\left(u_{i,1}^{Iter},u_{i,2}^{Iter},\dots ,u_{i,D}^{Iter}\right)$ as Eq. (10):

(10) $$u_{i,j}^{Iter}=\{\begin{array}{cc}{v}_{i,j}^{Iter},& ran{d}_{i,j}\left(0,1\right)\le CR;\\ x_{i,j}^{Iter},& \mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e},\end{array}$$where $ran{d}_{i,j}\left(0,1\right)\in U\left[0,1\right]$, $i=1,2,\dots ,N_{pop}$, and $j=1,2,\dots ,D.$

4. Selection

In this step, the fitness function of a new population ${\overrightarrow{U}}_i^{Iter}$, i.e., $f\left({\overrightarrow{U}}_i^{Iter}\right)$ is evaluated and compared with the fitness of the current position $f\left({\overrightarrow{X}}_i^{Iter}\right)$, and if it has a better fitness function value, it replaces the current solution, otherwise, the $i$-th member keeps its current position as Eq. (11):

(11) $${\overrightarrow{X}}_i^{Iter+1}=\{\begin{array}{cc}{\overrightarrow{U}}_i^{Iter},& f\left({\overrightarrow{U}}_i^{Iter}\right)\le f\left({\overrightarrow{X}}_i^{Iter}\right);\\ {\overrightarrow{X}}_i^{Iter},& \mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}.\end{array}$$

GWO

This section presents the mathematical modeling of GWO (Mirjalili, Mirjalili & Lewis, 2014).

1. Social hierarchy

When developing GWO, the best solution based on the fitness function is considered as α wolf (the dominant wolf). Similarly, the second and third fittest solutions are denoted as $\beta$ and $\delta$ wolves. All the remaining candidates are denoted as ω wolves, which are directed by $\alpha$, $\beta$, and $\delta$ to pursue the hunting process.

2. Encircling prey

Eqs. (12) and (13) are proposed to mathematically model the encircling action of grey wolves.

(12) $${\overrightarrow{D}}_{}^{Iter}=\left|\overrightarrow{C}\ast {\overrightarrow{X}}_p^{Iter}-{\overrightarrow{X}}_{}^{Iter}\right|,$$

(13) $${\overrightarrow{X}}_{}^{Iter+1}={\overrightarrow{X}}_p^{Iter}-\overrightarrow{A}\ast {\overrightarrow{D}}_{}^{Iter},$$where $Iter$ denotes the present iteration, $\overrightarrow{A}$ and $\overrightarrow{C}$ represent coefficient vectors, ${\overrightarrow{X}}_p^{Iter}$ indicates the prey’s location vector in the iteration $Iter$, and vectors ${\overrightarrow{X}}_{}^{Iter}$ and ${\overrightarrow{X}}_{}^{Iter+1}$ are gray wolf position vectors in the iteration $Iter$ and $Iter+1$, respectively. The parameters $\overrightarrow{A}$ and $\overrightarrow{C}$ are computed by Eqs. (14) and (15).

(14) $$\overrightarrow{A}=2\overrightarrow{a}\ast \overrightarrow{r_1}-\overrightarrow{a},\mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}\overrightarrow{a}=2\left(1-{\displaystyle \frac{Iter}{Ite{r}_{max}}}\right),$$

(15) $$\overrightarrow{C}=2\overrightarrow{r_2},$$where the components of $\overrightarrow{a}$ are linearly decreasing from 2 to 0, concerning the growth of the number of iterations, and $\overrightarrow{r_1}$, $\overrightarrow{r_2}$ are uniform random vectors with components in the interval [0, 1].

3. Hunting

A pack of gray wolves can detect the presence of prey and surround and enclose it. The alpha is generally in charge of the hunt. It is likely that the beta and delta sometimes also involve in hunting. In an abstract solution area, however, we do not know where the optimum prey is located. In order to computationally recreate grey wolf hunting behavior, we assume that the alpha, beta, and delta wolves contain more information about possible prey locations. As a result, the first three best solutions found so far are kept, and the other search representatives (containing the omegas) update their locations in accordance with the best search agents’ placements. The following Eqs. (16) to (18) are proposed in this regard for $i=1,2,\dots ,N_{pop}$.

(16) $$\begin{array}{c}{\overrightarrow{D}}_{\alpha }^{Iter}=\left|\overrightarrow{C_1}\ast \overrightarrow{X_{\alpha }}-{\overrightarrow{X}}_i^{Iter}\right|,\\ {\overrightarrow{D}}_{\beta }^{Iter}=\left|\overrightarrow{C_2}\ast \overrightarrow{X_{\beta }}-{\overrightarrow{X}}_i^{Iter}\right|,\\ {\overrightarrow{D}}_{\delta }^{Iter}=\left|\overrightarrow{C_3}\ast \overrightarrow{Y_{\delta }}-{\overrightarrow{X}}_i^{Iter}\right|,\end{array}$$

(17) $$\begin{array}{c}\overrightarrow{X_{\mathrm{A}}}=\overrightarrow{X_{\alpha }}-\overrightarrow{A_1}\ast {\overrightarrow{D}}_{\alpha }^{Iter},\\ \overrightarrow{X_{\mathrm{B}}}=\overrightarrow{X_{\beta }}-\overrightarrow{A_2}\ast {\overrightarrow{D}}_{\beta }^{Iter},\\ \overrightarrow{X_{\mathrm{\Delta }}}=\overrightarrow{X_{\delta }}-\overrightarrow{A_3}\ast {\overrightarrow{D}}_{\delta }^{Iter},\end{array}$$

(18) $${\overrightarrow{X}}_i^{Iter+1}=\left(\overrightarrow{X_{\mathrm{A}}}+\overrightarrow{X_{\mathrm{B}}}+\overrightarrow{X_{\mathrm{\Delta }}}\right)/3$$

Greedy non-hierarchical GWO

During the update mechanism of the algorithm GWO, the best three options are steadily saved and denoted as wolves ( $\alpha$, $\beta$, and $\delta$) and these population members lead the other members in their upgrade calculation, similarly as it works in the natural social hierarchy of grey wolf packs. It was shown the GWO algorithm converges quickly to the optimal solution for most of the well-known benchmark functions. On the other hand, this mechanism has a crucial disadvantage in optimizing some real-world problems, i.e., the algorithm converges rapidly but to a locally optimal solution. To overcome these drawbacks in the GWO, the G-NHGWO algorithm (Akbari, Rahimnejad & Gadsden, 2021), the best personal optimal solution ${\overrightarrow{X}}_{i,best}^{Iter}$ for the $i$-th grey wolf, $i=1,2,\dots ,N_{pop},$ has been established and stored, like in the PSO algorithm. Then, three members, $r_1$, $r_2$, and $r_3$, with individual best positions, respectively ${\overrightarrow{X}}_{r_1,best}^{Iter}$, ${\overrightarrow{X}}_{r_2,best}^{Iter}$, and ${\overrightarrow{X}}_{r_3,best}^{Iter}$, are selected randomly and utilized to lead the population to update the new positions as Eqs. (19) to (21):

(19) $$\begin{array}{c}{\overrightarrow{D}}_{r_1}^{Iter}=\left|{\overrightarrow{C}}_1\ast {\overrightarrow{X}}_{r_1,best}^{Iter}-{\overrightarrow{X}}_{best}^{Iter}\right|,\\ {\overrightarrow{D}}_{r_2}^{Iter}=\left|{\overrightarrow{C}}_2\ast {\overrightarrow{X}}_{r_2,best}^{Iter}-{\overrightarrow{X}}_{best}^{Iter}\right|\\ {\overrightarrow{D}}_{r_3}^{Iter}=\left|{\overrightarrow{C}}_3\ast {\overrightarrow{X}}_{r_3,best}^{Iter}-{\overrightarrow{X}}_{best}^{Iter}\right|,\end{array}$$

(20) $$\begin{array}{c}{\overrightarrow{X}}_{r_1}={\overrightarrow{X}}_{r_1,best}^{Iter}-{\overrightarrow{A}}_1\ast {\overrightarrow{D}}_{r_1}^{Iter},\\ {\overrightarrow{X}}_{r_2}={\overrightarrow{X}}_{r_2,best}^{Iter}-{\overrightarrow{A}}_2\ast {\overrightarrow{D}}_{r_2}^{Iter},\\ {\overrightarrow{X}}_{r_3}={\overrightarrow{X}}_{r_3,best}^{Iter}-{\overrightarrow{A}}_3\ast {\overrightarrow{D}}_{r_3}^{Iter},\end{array}$$

(21) $${\overrightarrow{X}}_i^{Iter+1}=\left({\overrightarrow{X}}_{r_1}+{\overrightarrow{X}}_{r_2}+{\overrightarrow{X}}_{r_3}\right)/3.$$

HDE algorithm

In the GWO algorithm, $\alpha$, $\beta$, and $\delta$ wolves always participate in the hunting phase. Therefore, using the positions of the best wolves in the hunting phase increases the quality of the exploration and accelerates the convergence of the GWO method (Akbari, Rahimnejad & Gadsden, 2021). Since one of the main problems of the DE is the low convergence speed, in this study, the hunting phase of GWO is proposed as an auxiliary mutation vector to the DE to accelerate the convergence rate and achieve more optimal solutions. In the process of the proposed HDE algorithm, in each iteration, instead of the general vector in the DE algorithm, the GWO’s hunting vector is used as the mutation vector with the probability of H_m for each member. Then the crossover and selection operators are applied.

In the proposed HDE algorithm, instead of using ${\overrightarrow{X}}_{r_1,best}^{Iter}$, ${\overrightarrow{X}}_{r_2,best}^{Iter}$, and ${\overrightarrow{X}}_{r_3,best}^{Iter}$ in (19) and (20), three members of $\overrightarrow{X_{\alpha }}$, $\overrightarrow{X_{\beta }}$, and $\overrightarrow{X_{\delta }}$ are used. Therefore, the new position of the $i$-th member, $i=1,2,\dots ,N_{pop}$, can be calculated using Eqs. (22) to (24).

(22) $$\begin{array}{c}{\overrightarrow{D}}_1^{Iter}=\left|\overrightarrow{C_1}\ast \overrightarrow{X_{\alpha }}-{\overrightarrow{X}}_{best}^{Iter}\right|,\\ {\overrightarrow{D}}_2^{Iter}=\left|\overrightarrow{C_2}\ast \overrightarrow{X_{\beta }}-{\overrightarrow{X}}_{best}^{Iter}\right|\\ {\overrightarrow{D}}_3^{Iter}=\left|\overrightarrow{C_3}\ast \overrightarrow{X_{\delta }}-{\overrightarrow{X}}_{best}^{Iter}\right|,\end{array}$$

(23) $$\begin{array}{c}\overrightarrow{X_{\mathrm{A}}}=\overrightarrow{X_{\alpha }}-\overrightarrow{A_1}\ast {\overrightarrow{D}}_1^{Iter},\\ \overrightarrow{X_{\mathrm{B}}}=\overrightarrow{X_{\beta }}-\overrightarrow{A_2}\ast {\overrightarrow{D}}_2^{Iter},\\ \overrightarrow{X_{\mathrm{\Delta }}}=\overrightarrow{X_{\delta }}-\overrightarrow{A_3}\ast {\overrightarrow{D}}_3^{Iter},\end{array}$$

(24) $${\overrightarrow{X}}_i^{Iter+1}=\left(\overrightarrow{X_{\mathrm{A}}}+\overrightarrow{X_{\mathrm{B}}}+\overrightarrow{X_{\mathrm{\Delta }}}\right)/3.$$where the operation “ $\ast$” represents the Hadamard product of two vectors (i.e., the multiplication of vectors by multiplying their corresponding components), and the absolute value function $|\cdot |$ is calculated sequentially by the components of the vectors.

Afterward, ${\overrightarrow{X}}_i^{Iter+1}$ is compared with ${\overrightarrow{X}}_{best}^{Iter},$ and if it has a better objective function value, it will be considered the new personal best position (selection phase). This process is repeated until reaching the maximum number of fitness evaluations $NF{E}_{max}$ (i.e., $NF{E}_{max}=Ite{r}_{max}\cdot N_{pop}$).

The HDE algorithm can be written in the following pseudocode steps:

Step 0: Input the optimization problem information and values of parameters $N_{pop}$, $D$, $NF{E}_{max}$, $F$, $H_m$, and $CR$.

Step 1: Initialization of the population of search agents (grey wolfs), we set $Iter=1$: ${\overrightarrow{X}}_i^{Iter}$, $i=1,2,\dots ,N_{pop}$.

Step 2: Compute the fitness value for all search agents.

Step 3: Initialize parameters $\overrightarrow{a}$, $\overrightarrow{A}$, and $\overrightarrow{C}$ by (14) and (15).

Step 4: While $Iter\cdot N_{pop}<NF{E}_{max}$ Do

Step 5:       Evaluate the search agents $\overrightarrow{X_{\alpha }}$, $\overrightarrow{X_{\beta }}$, and $\overrightarrow{X_{\delta }}$ ( $\alpha ,\beta ,$ and $\delta$ wolves).

Step 6:       If $rand<H_m$ Then perform the hunting mutation using (22) to (24).

Step 7:       Else Perform the original mutation with a probability 1 – $H_m$ using (2) to (9).

Step 8:       Update parameters $\overrightarrow{a}$, $\overrightarrow{A}$, and $\overrightarrow{C}$ by (14) and (15).

Step 9:       Evaluate all search agents’ fitness values.

Step 10:      Update the value of $\overrightarrow{X_{\alpha }}$, $\overrightarrow{X_{\beta }}$, and $\overrightarrow{X_{\delta }}.$

Step 11:      Crossover operation by (10).

Step 12:      Evaluation of the fitness value of the population of all search agents.

Step 13:      Selection by (11).

Step 14:       $Iter=Iter+1$

Step 15: End While.

The computational complexity of HDE

It is essential to realize that three processes—initialization, fitness evaluation, and updating of the population—largely determine the HDE's computational complexity. First, the initialization process has a computational complexity equal to $O(N_{pop}$) for $N_{pop}$ individuals. During searching for the optimal location and updating the location vector of the entire population, the updating mechanism has a computational complexity equal to $O\left(Ite{r}_{max}\cdot N_{pop}\right)+O\left(Ite{r}_{max}\cdot N_{pop}\cdot D\right)$, where $Ite{r}_{max}$ is the maximal number of iterations and $D$ is the dimension of the problem. Finally, HDE has a total computational complexity of $O\left(N_{pop}\cdot \left(Ite{r}_{max}+Ite{r}_{max}\cdot D+1\right)\right).$

Experimental design

This section investigates the effectiveness and robustness of the proposed HDE algorithm on two standard benchmarks, including the CEC-2014 suite (Liang, Qu & Suganthan, 2013) and the CEC-2019 suite (Liang et al., 2019). Characteristics of the CEC-2014 suit of benchmark functions have been shown in Table 1. The CEC-2014 benchmarks have three unimodal functions ( $F_1$– $F_3$), thirteen simple multimodal functions ( $F_4$– $F_{16}$), six hybrid functions ( $F_{17}$– $F_{22}$), and eight composition functions ( $F_{23}$– $F_{30}$). The search range for each function is $\left[-100,100\right],$ while the minimum value for each function equals $100k$, where “ $k$” is the number of the benchmark function. The characteristics of CEC-2019 are defined in Table 2. In this scenario, different HDE algorithms are compared to clarify how different strategies can affect optimization performance. In contrast, the results of CEC-2014 benchmarks are compared with some other optimization algorithms. All experiments were conducted using MATLAB 2018a on a Windows 10 computer with a 2.2 GHz Core i7 processor and 16 GB of RAM.

Application of HDE for the CEC-2014 and CEC-2019 suit of benchmark functions

In this step of optimization, we require sophisticated functions to verify the capacity of the DE algorithms to optimize real functions, such as economic load dispatch, compared to other algorithms. Therefore, we chose the CEC-2014 test functions, which have already been applied in many previous articles. Here we have utilized the number of dimensions equal to 30, the number of function evaluations of 300,000, and the population's range of all algorithms was set to 30. The parameter $a$ decreases from 2 to 0. Each optimization method was run 30 times for each function. The obtained results were used to calculate the mean value, standard deviation, and a specified convergence graph. The set of control parameters of the algorithms is shown in Table 3. However, it should be noted that the value of the control parameter $F$ is not mentioned in this table since it is the same for all algorithms and is equal to $(0.1+0.8\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{d}())$.

The results characterized by the mean (the variable $\mathrm{M}\mathrm{e}\mathrm{a}\mathrm{n}$) and the standard deviation (the variable $\mathrm{S}\mathrm{t}\mathrm{d}.$) are summarized in Table 4 for the set of versions of both the proposed HDE and DE algorithms and the original GWO. In addition, the rank of all algorithms corresponding to each test function is shown in Table 4. The symbols plus sign ‘ $+$’ or minus sign ‘ $-$‘ and the symbol ‘ $=$’ help the readers determine the effectiveness of all proposed HDE compared to their original DE versions. The sign ‘ $+$’ indicates that our proposed algorithm performs better, the sign ‘ $-$’ indicates that it performs worse, and the symbol ‘ $=$’ demonstrates that it performs similarly.

Table 5 shows a global summary of the ranking sums (variable “ $\mathrm{S}\mathrm{R}\mathrm{N}$”) of all optimization algorithms from Table 4, then their qualitative ranking according to the value of $\mathrm{S}\mathrm{R}\mathrm{N}$, and the mean rank in the Friedman test. Furthermore, in Table 5, we find qualitative rankings by the partial summary of the ranking sums ‘ $\mathrm{S}\mathrm{R}\mathrm{N}-\mathrm{U}\mathrm{F}$’, ‘ $\mathrm{S}\mathrm{R}\mathrm{N}-\mathrm{S}\mathrm{M}\mathrm{F}$’, ‘ $\mathrm{S}\mathrm{R}\mathrm{N}-\mathrm{H}\mathrm{F}$’, and ‘ $\mathrm{S}\mathrm{R}\mathrm{N}-\mathrm{C}\mathrm{F}$’ for subgroups of functions unimodal function (UF), simple multimodal function (SMF), hybrid function (HF), and composition function (CF), respectively. Table 5 implies that the three methods that globally best optimize the CEC-2014 benchmark functions are “HDE/current-to-best/1”, “HDE/rand-to-best/2” and “HDE/rand/1”.

Further concretely, that:

• “HDE/best/2”, “HDE/best/1”, “HDE/rand-to-best/2”, and “HDE/rand/1” best optimize the subgroup UF,

• “HDE/rand-to-best/2”, "HDE/current-to-best/1”, and “HDE/current-to-rand/1” best optimize the subgroup SMF,

• “HDE/current-to-best/1”, “HDE/rand/1”, and “HDE/best/2” best optimize the subgroup HF,

• “HDE/rand/2”, “HDE/rand/1”, and “HDE/current-to-best/1” best optimize the subgroup CF.

The list of previous items shows that “HDE/current-to-best/1” and “HDE/rand/1” were the best optimizers three times placed top three positions in all four subgroups of CEC-2017 (similarly, “HDE/rand-to-best/2” two times placed top three positions).

In addition to the already mentioned great results of “HDE/current-to-best/1”, “HDE/rand-to-best/2”, and “HDE/rand/1” in some subsets of CEC-2017, it should be added that neither they did not do poorly on the other subsets, as they always ranked above the median, what indicates the robustness and dependability of the proposed HDE algorithms.

Moreover, finally, it can be seen that the proposed HDEs are generally much more successful than their original DE versions, proving their effectiveness and superiority over original DE versions. Hence, it is fair to say that the proposed HDE algorithms produce respectable solutions with the potential for further optimization toward becoming a more reliable algorithm. Figures 1–3 present good convergence properties of the HDEs on 30 test functions from the CEC-2014.

The outcomes of applying various proposed techniques to the CEC-2019 suite benchmark functions are shown in Table 6. The Friedman test demonstrates that the three best algorithms dealing with this set of benchmark functions are from the proposed HDEs, confirming the proposed method's effectiveness. The “HDE/current-to-rand/1” and the “HDE/rand-to-best/1” perform best on the CEC-2019 suite benchmark functions.

Also, to compare the proposed HDEs with other evolutionary algorithms, the results of “HDE/current-to-best/1” on the CEC-2014 suite with $D=30$ are compared to some modified versions of evolutionary algorithms in Table 7, such as m-SCA (Gupta & Deep, 2019a), MG-SCA (Gupta, Deep & Engelbrecht, 2020), CTLBO (Yang, Li & Guo, 2014), LJA (Iacca, dos Santos Junior & de Melo, 2021) and some modified version of GWO, e.g., Cauchy-GWO (Gupta & Deep, 2018), RW-GWO (Gupta & Deep, 2019b), and ERGWO (Long et al., 2020). This comparison confirms that the proposed HDE algorithms can be considered an effective evolutionary method dealing with a vast range of optimization problems.