Diversity-Based Evolutionary Population Dynamics: A New Operator for Grey Wolf Optimizer

Abstract

Evolutionary Population Dynamics (EPD) refers to eliminating poor individuals in nature, which is the opposite of survival of the fittest. Although this method can improve the median of the whole population of the meta-heuristic algorithms, it suffers from poor exploration capability to handle high-dimensional problems. This paper proposes a novel EPD operator to improve the search process. In other words, as the primary EPD mainly improves the fitness of the worst individuals in the population, and hence we name it the Fitness-Based EPD (FB-EPD), our proposed EPD mainly improves the diversity of the best individuals, and hence we name it the Diversity-Based EPD (DB-EPD). The proposed method is applied to the Grey Wolf Optimizer (GWO) and named DB-GWO-EPD. In this algorithm, the three most diversified individuals are first identified at each iteration, and then half of the best-fitted individuals are forced to be eliminated and repositioned around these diversified agents with equal probability. This process can free the merged best individuals located in a closed populated region and transfer them to the diversified and, thus, less-densely populated regions in the search space. This approach is frequently employed to make the search agents explore the whole search space. The proposed DB-GWO-EPD is tested on 13 high-dimensional and shifted classical benchmark functions as well as 29 test problems included in the CEC2017 test suite, and four constrained engineering problems. The results obtained by the proposal upon implemented on the classical test problems are compared to GWO, FB-GWO-EPD, and four other popular and newly proposed optimization algorithms, including Aquila Optimizer (AO), Flow Direction Algorithm (FDA), Arithmetic Optimization Algorithm (AOA), and Gradient-based Optimizer (GBO). The experiments demonstrate the significant superiority of the proposed algorithm when applied to a majority of the test functions, recommending the application of the proposed EPD operator to any other meta-heuristic whenever decided to ameliorate their performance.

1. Introduction

Optimization algorithms in evolutionary computation can be divided into two classes individual-based and population-based. Individual-based algorithms start the optimization by a single randomly generated search agent in the search space attempting to seek and detect the global optimum of the optimization problem. This category of algorithms benefits from some advantages and suffers from some disadvantages. The main shortcoming of these algorithms is premature convergence emerging in two different types of optimization problems, called uni-modal and multi-modal functions. In uni-modal functions, premature convergence usually occurs when the algorithm’s convergence rate during the optimization process is too slow.

In contrast, the major cause for premature convergence in multi-modal optimization problems is known as the local optima entrapment. Individual-based algorithms begin the optimization from an initial random point in the search space. Usually, they can only search for the optimal solution in the proximity of that initial random point without escaping from that region to search for other potentially high-fitness regions in the problem domain. On the contrary, the algorithms of type population-based commence the optimization process by generating a series of random solutions being improved throughout iterations. This category of algorithms is less likely to be stuck in the local optima in multi-modal problems and also benefits from a more rapid search process in uni-modal problems. However, having more computational costs and needing much more objective function evaluations are two significant disadvantages; this category suffers compared to the individual-based algorithms.

Among the stochastic population-based optimization algorithms, the Genetic Algorithm (GA) [1], Particle Swarm Optimization (PSO) [2], and Differential Evolution (DE) [3] are the pioneers, based on which a vast number of modified and improved algorithms have been developed so far and applied to a wide range of the practical problems in the science and engineering. There are a lot of nature-inspired and physics-based stochastic optimization algorithms in the literature as well, including but not limited to Artificial Bee Colony (ABC) algorithm [4], Firefly Algorithm (FA) [5], Bat Algorithm (BA) [6], Gravitational Search Algorithm (GSA) [7], Krill Herd (KH) algorithm [8], Sine Cosine Algorithm (SCA) [9], Monarch Butterfly Optimization (MBO) algorithm [10], and Invasive Weed Optimization (IWO) algorithm [11].

1.1. Hybrid Meta-Heuristic Optimization Algorithms

An effective way to enhance the performance of the optimization algorithms is the combination or hybridization of them with different operators and search mechanisms employed in other optimizers. Due to the consolidation of benefits from different algorithms during the hybridization process, many improved versions of the algorithms have been proposed so far. Some of these algorithms are described in the following paragraphs.

Tuba and Bacanin et al. [12] hybridized the ABC with the FA to improve the convergence rate of the ABC and also to well adjust the exploration-exploitation balance in the ABC algorithm. Wang et al. [13] introduced the Pitch adjustment operator employed in HS into the CS. This operator can act as a mutation operator to increase the cuckoo population diversity and compensate for the weakness of the Lévy flight incorporated in the CS, as the Lévy flight can lead the CS not to be able to converge to the best-fitted solutions whenever the size of the steps taken is too large. As a result, this hybridization can impede premature convergence. In order to improve the DE algorithm when attempting to balance local and global search processes while impeding the time-consuming control parameter tuning procedure, a new hybrid DE algorithm is proposed by Yi et al. [14]. This method divides the population into two sections to apply the most appropriate control parameters and mutation operators. Tuba and Bacanin [15] improved the global search procedure of the BA by incorporating the onlooker mechanism from the ABC algorithm to enhance the exploitation process of BA and applied the hybridized algorithm to solve a multi-objective radio frequency identification network planning problem. A hybrid PSO and GSA abbreviated as PSO-GSA was introduced by Das et al. [16], benefiting from a cooperative contribution of PSO velocity and GSA acceleration. In this hybridized algorithm, the PSO uses memory to save the best solutions found so far, and this ability is added to the ability of GSA to adjust the acceleration using the fitness value. Furthermore, in this algorithm, the global best particle of PSO can lead the other particles to move towards the best solution and accomplish exploitation slowly. A hybrid strategy is proposed by Abualigah et al. [17], aimed at increasing the KH’s population’s diversity by incorporating the HS operator using a new parameter named distance factor for fine-tuning the positions of the search agents before updating them. Nenavath and Jatoth [18] propose the hybridization of the SCA and DE. This hybrid algorithm improves the capability of the search agents to escape from the local optima and enhances the convergence speed compared to SCA and DE when operating separately. The proposed algorithm was then employed to effectively solve an object tracking problem as a real-life optimization problem. Ghanem proposes a new hybrid approach, and Jantan [19] to alter the butterfly tuning operator used in MBO while utilizing this improved operator in the role of a mutation operator to be applied as an alternative to the employee stage of the ABC. Another hybridization process has been carried out between the IWO algorithm and the FA. As any hybridization scheme aims at highlighting the best properties of a couple of the other algorithms to compensate for their drawbacks, the hybrid IWO-FA is proposed by Panda et al. [20] to improve both algorithms. As there is no known mechanism to hinder the fast movements of the search agents in IWO, the FA can also face local optima entrapments due to its dependence on the light intensity for attraction, and the IWO-FA can be taken into account as an approach to ameliorate these shortcomings simultaneously. Teng et al. [21] benefited from the high diversification ability of the GWO and high intensification ability of the PSO in local search to propose the hybrid PSO_GWO algorithm. In this algorithm, the best individual position experienced is added to the best wolves to guide the search agents in the population. The Tent chaotic sequence was utilized better to diversify the wolf pack in the algorithm initialization step. The use of the entire MBO algorithm while not using the Lévy flight that typically acts as a crossover operator for DE, proposed by Ibrahim and Tawhid [22], is another hybridized scheme named DEMBO. A modified version of the KH was proposed by Abualigah et al. [23], in which the genetic operators are employed to improve the KH performance highly. In this modified algorithm, the crossover and mutation processes are invoked after updating the krill positions because the nature of the search space of most optimization problems is ragged and deep. A hybrid method combining the SCA and ABC, called SCABC, is proposed by Gupta and Deep [24]. In the SCABC, first, the major problem of the SCA is the intensive diversification of the solutions in the early iterations while leaving the solutions less diversified in the later iterations. Then, this problem was mitigated by modifying its search procedures by integrating memory-based information. Then, the modified SCA search procedures are merged to be applied to the employed bee phase to enhance the global and local search abilities of the ABC algorithm. A new algorithm hybridizing the PSO and GWO is also proposed by Şenel et al. [21]. The hybridized scheme can replace some particles in PSO with a small probability by some other improved search agents obtained by GWO. The proposed mechanism can hamper the PSO from being trapped in local optima utilizing the exploration capability of GWO. Gupta and Deep [25] hybridized the Lévy-flight mechanism and the original GWO, yielding GLF-GWO. The Lévy-flight strategy helps better local search around the three leading wolves in GWO to enable the wolves to avoid local optima. Furthermore, a greedy selection mechanism is employed to prevent the wolves from diverging from the promising solutions during the optimization. A new hybrid HSCA is proposed by Gupta and Deep [26], mainly to impede the solutions of the SCA to skip the fitted regions of the decision space, to prevent the solutions from being trapped in local optima, and to balance the exploration-exploitation transition in the classical SCA more desirably. In the HSCA, the SCA and simulated quenching algorithm are hybridized to impart the elitism mechanism to the algorithm while easing the transition from exploration to exploitation. Mohammed and Rashid [27] hybridized the WOA with GWO to enhance the exploitation capability of the WOA. They found the exploitation capability of the GWO much stronger than the mechanism dedicated to exploitation in the WOA. That was their motivation to substitute the two mechanisms applied in two algorithms. A hybridization of SMA and AOA was proposed by Zheng et al. [28]. This hybrid algorithm replaces the contraction formula employed to perform exploration in the SMA is replaced by the multiplication and division operators used in the exploration phase of the AOA, as the position of each search agent may tend to zero in the late iterations, which, in, turn, can highly debilitate the performance of the SMA when solving an optimization problem with an unknown search space. Rezaei et al. [29] added the velocity term existing in the position updating formula of the PSO or some other meta-heuristics into the updating procedure of the GWO and proposed the VAGWO algorithm. This hybridized scheme can highly improve the exploration capability of the GWO by keeping push the search agents to move forwards to detect the high-fitness positions as much as possible. In addition, a decreasing radius circle is made around the elite agents to impede the potential drifts to occur when the agents are tending to the high-fitness positions. A hybridized Whale and Moth-Flame Optimizer, named WMFO, was proposed by Nadimi-Shahraki et al. [30]. The goal of this hybridization is enhancing the exploitation capability of the MFO via enforcing the flames in the MFO to move towards an average of the best-so-far agents found by the algorithm in the exploitation phase of the hybrid algorithm to expedite the agents’ access to the high-fitness region in the search space.

No Free Lunch (NFL) theorem [31] allows new evolutionary algorithms to be proposed, as based on this theory, all algorithms have the same performance when applied to a vast range of optimization problems. As a result, an algorithm may show too effective performance when solving a series of problems but suffers from serious shortcomings when going over solving other problems. Consequently, creating new effective algorithms via hybridizing them with the other ones or equipping them with some efficient operators can highly improve the performance of such algorithms to better solve a series of optimization problems the other optimizers have difficulty handling.

Selection, combination, and mutation are usually considered the widely used evolutionary operators incorporated into the meta-heuristics to enhance their performance, as illustrated by Lewis et al. [32]. However, another evolutionary operator, Evolutionary Population Dynamics (EPD), has been found to manipulate the whole population. EPD is inspired by the theory of Self-Organizing Critically (SOC), first introduced by [33]. Based on this theory, the critical state in nature can be reached dynamically. For instance, small perturbations can make operations precisely balance different features of the population while not imposing any external force on them. Extremal Optimization (EO), proposed by Boettcher and Percus [34], is a meta-heuristic inspired by the self-organizing critical model utilizing EPD. In EO, the worst solutions in the population are eliminated and repositioned around the best ones. This process can evolve the worst solutions found during the optimization course, in contrast to the typical process in the Genetic Algorithm (GA), where the best solutions are always combined and evolved during the optimization.

1.2. The Contribution of This Study

In this paper, we propose a novel EPD operator addressing the other side of the goodness of an individual in the population of the EPD-based evolutionary algorithms. In the original version of EPD, the worst individuals are considered low-fitness individuals, required to be improved by removing and repositioning them around the best (high-fitness) individuals. However, the other side of the goodness/fitness of an individual can also be interpreted as the diversity of that individual in the search space. Having this another side for fitness, we can evolve and improve the best (high-fitness) individuals by eliminating a number of them and then repositioning them around the most diversified individuals in the search space. In this way, the best individuals are getting better, not in their fitness, but in their diversity as the other aspect of fitness that is very important to address, especially in solving complex and large-scale optimization problems. Our novel EPD-based mechanism, which we name Diversity-Based EPD (DB-EPD), is applied to the Grey Wolf Optimizer (GWO) in this paper, and the results are compared to the original GWO and the GWO-EPD proposed by Saremi et al. [35], which is addressing the application of the traditional EPD, which we name the Fitness-Based EPD (FB-EPD), to improve the GWO.

The organization of the rest of this paper is as follows. Section 2 presents the methodology, in which we first introduce the theory of GWO and then formulate the FB-EPD for GWO. Then we introduce our novel DB-EPD mechanism and its application to improve the GWO in Section 3. Section 4 presents the results and some analyses of them. Section 5 concludes the paper and summarizes the main contributions of this paper.

2. Methodology

This section introduces the traditional grey wolf optimization and Fitness-Based EPD for GWO.

2.1. Original Grey Wolf Optimization (GWO) Algorithm

The GWO algorithm was proposed by Mirjalili et al. in 2014 [36]. This algorithm imitates the hunting behavior and hierarchical leadership of grey wolves in nature. The GWO initializes the optimization process by randomly producing a population of the solutions (wolves). At each iteration, the three best-fitted wolves, named alpha, beta, and delta, are identified as the guide for the other wolves named omega. Afterward, the omega wolves encircle their guides to find the high-fitness regions in the search space. The wolves are named the search agents in the concept of swarm-intelligence optimization techniques. As every omega wolf encircles the three best-fitted agents, its final position can be obtained by simply averaging the updated positions of the alpha, beta, and delta wolves. This procedure may boost the exploratory ability of the algorithm due to adopting more than one search agent as the guide to the other agents. The mathematical formulation of updating the omega wolves is as follows:

$$\overrightarrow{D}=\left|\overrightarrow{C}.\overrightarrow{X_p\left(t\right) }-\overrightarrow{X\left(t\right)}\right|$$

(1)

$$\overrightarrow{X(t+1) }=\overrightarrow{X_p\left(t\right) }-\overrightarrow{A}.\overrightarrow{D}$$

(2)

where t is the current iteration; $\overrightarrow{A}=2\overrightarrow{r_1}.\overrightarrow{a}-\overrightarrow{a}$; $\overrightarrow{C}=2\overrightarrow{r_2}$; $\overrightarrow{X_p\left(t\right)}$ is the position vector of the prey; $\overrightarrow{X\left(t\right)}$ is the position vector of a grey wolf. Moreover, $\overrightarrow{r_1}$ and $\overrightarrow{r_2}$ are two random vectors in [0, 1]. Furthermore, $\overrightarrow{a}$ is linearly decreased from 2 to 0 by lapse of iterations. The factor $\overrightarrow{A}$ facilitates the optimization process by controlling a safe and reliable transition from the exploration to the exploitation phases of the algorithm. Moreover, $\overrightarrow{C}$ is multiplied by the guiding vector to help the exploration process more precisely.

The position of each omega wolf is updated through several equations formulated below:

$$\overrightarrow{D_{\alpha }}=\left|\overrightarrow{C_1}.\overrightarrow{X_{\alpha }}-\overrightarrow{X}\right|$$

(3)

$$\overrightarrow{D_{\beta }}=\left|\overrightarrow{C_2}.\overrightarrow{X_{\beta }}-\overrightarrow{X}\right|$$

(4)

$$\overrightarrow{D_{\delta }}=\left|\overrightarrow{C_3}.\overrightarrow{X_{\delta }}-\overrightarrow{X}\right|$$

(5)

$$\overrightarrow{X_1}=\overrightarrow{X_{\alpha }}-\overrightarrow{A_1}.\overrightarrow{D_{\alpha }}$$

(6)

$$\overrightarrow{X_2}=\overrightarrow{X_{\beta }}-\overrightarrow{A_2}.\overrightarrow{D_{\beta }}$$

(7)

$$\overrightarrow{X_3}=\overrightarrow{X_{\delta }}-\overrightarrow{A_3}.\overrightarrow{D_{\delta }}$$

(8)

$$\overrightarrow{X\left(t+1\right) }=\frac{\overrightarrow{X_1}+\overrightarrow{X_{2 }}+\overrightarrow{X_3}}{3}$$

(9)

As Mirjalili et al. (2014) expressed, in GWO, the first half of the iterations are assigned to the exploration process when $\left|\overrightarrow{A}\right| > 1$, and the second half is allocated to the exploitation phase of the optimization process, when $\left|\overrightarrow{A}\right| < 1$.

2.2. Fitness-Based EPD for GWO (FB-GWO-EPD)

In this approach proposed by Saremi et al. [35], half of the worst search agents are detected and eliminated from the population. Then, the omitted agents are re-initialized around the best agents in terms of fitness. This process evolves the low-fitness solutions instead of the best-fitted ones. The procedure of repositioning the worst agents around the four random positions with equal probability can be implemented via Equations (10)–(13).

$$\overrightarrow{X(t+1)}=\overrightarrow{X_{\alpha }\left(t\right)} \pm \left(\left(\overrightarrow{\mathit{ub}}-\overrightarrow{\mathit{lb}}\right).\overrightarrow{r_1}+\overrightarrow{\mathit{lb}}\right); \mathrm{if} 0 \le { r}_5 \le 1/4$$

(10)

$$\overrightarrow{X(t+1)}=\overrightarrow{X_{\beta }\left(t\right)} \pm \left(\left(\overrightarrow{\mathit{ub} }-\overrightarrow{\mathit{lb}}\right).\overrightarrow{r_{2 }}+\overrightarrow{ \mathit{lb}}\right){; \mathrm{if} 1/4 < r}_5 \le 1/2$$

(11)

$$\overrightarrow{X(t+1)}=\overrightarrow{X_{\delta }\left(t\right)} \pm \left(\left(\overrightarrow{\mathit{ub}}-\overrightarrow{\mathit{lb}}\right).\overrightarrow{r_{3 }}+\overrightarrow{\mathit{lb}}\right){; \mathrm{if} 1/2 < r}_5 \le 3/4$$

(12)

$$\overrightarrow{X(t+1)}=\left(\left(\overrightarrow{\mathit{ub}}-\overrightarrow{\mathit{lb}}\right).\overrightarrow{r_4}+\overrightarrow{\mathit{lb}}\right){; \mathrm{if} 3/4 < r}_5 \le 1$$

(13)

where, $\overrightarrow{X_{\alpha }\left(t\right)}$, $\overrightarrow{X_{\beta }\left(t\right)}$ and $\overrightarrow{X_{\delta }\left(t\right)}$ are the alpha, beta, and delta wolf positions; $\overrightarrow{\mathit{ub}}$ indicates the upper bound vector of the search space; $\overrightarrow{\mathit{lb}}$ denotes the lower bound vector of the search space, and $\overrightarrow{r_1}$ to $\overrightarrow{r_5}$ are five uniformly distributed random vectors in [0, 1]. Equations (10)–(12) formulate the repositioning process of the worst solutions around the best ones, and Equation (13) formulates repositioning the worst solutions around randomly positions around the search space to maintain the diversity among the diversified worst solutions to promote the exploration and avoid missing a large number of suitable solutions when the algorithm approaches the leading solutions.

3. Proposed Method

As discussed in the previous subsection, the fitness-based GWO-EPD aims at eliminating the worst (low-fitness) individuals instead of evolving the best ones. This process is well done when reproducing the worst individuals around the three elite ones in GWO: alpha, beta, and delta. However, repositioning the low-fitness individuals around the high-fitness areas in the search space can, in turn, intensify the risk of engagement in premature convergence. In other words, the low-fitness individuals may benefit from being located in the high-diversified positions in the search space, in contrast to those in the low-diversified regions in the search space. The process of repositioning the worst individuals around the best ones can further exacerbate the diversity of the solutions in the search space when knowing that the best (high-fitness) individuals are being rapidly gathered in a closed region, leading the search space to lose its diversity. Since population diversification is a key phase in any stochastic population-based optimization algorithm, such as the GWO, there is an urgent need to diversify the solutions produced by the algorithm.

In the proposed diversity-based evolutionary population dynamics for GWO, not only the best-fitted solutions are repositioned around the highly diversified areas in the search space, but also all the solutions continue to converge to the three (alpha, beta, and delta) optimal solutions found at each iteration. This way, the solutions are frequently gathered inside the high-fitness region and dispersed outside of this region so that the solutions can be circulated across the search space. In this process, the poor-fitted solutions move to the best-fitted ones, on the one hand, and the best-fitted solutions are moved to the diversified areas, on the other hand. This process can better maintain the balance between the exploration capability of the GWO, in which the solutions are to be diversified to fully cover the search space, and the exploitation capability, in which the solutions must intensify, moving to the best regions in the search space to converge to the optimum of the optimization problem finally.

In the proposed diversity-based version of the GWO-EPD algorithm, we first eliminate half of the best individuals in terms of fitness and then reposition them around the three highly diversified solutions identified in the solution space. It is noteworthy that the alpha, beta, and delta agents are maintained in the algorithm’s memory to guide search agents. However, they are repositioned around the highly diversified agents if involved in the current population. The diversity of the solutions can be numerically calculated through Equation (14).

$$d_{i }=\underset{j}{\min }\left|f^i\left(X\right) -{ f}^j\left(X\right)\right|$$

(14)

where i = 1, 2, …, N; j $\in \left\{1, \dots ,i - 1, i+1, \dots , N\right\}$; N is the population size; $f^i\left(X\right)$ = objective (fitness) function value for the ith solution, and $f^j\left(X\right)$ = objective (fitness) function value for the jth solution. In addition, $d_i$ is the diversity index of the ith solution. As can be seen in Equation (14), the differences between the fitness value of each solution from the fitness values of the other solutions in the population are all calculated, and then, the minimum of these differences is taken as the intersection of them to result in the $d_i$ index. On the one hand, the differences in the fitness value of a solution from the other solutions’ can properly reflect how different that solution is from the other solutions and so how diversified it is across the search space. On the other hand, the higher the intersection of several values, the higher all of these values will be, as their intersection can always be regarded as those values in common. Hence, the more the $d_i$ value, the higher the diversity of the ith search agent will be. The three individuals of the swarm denoted by the three highest values of $d_i$ are then identified as the most diversified individuals and considered to guide half of the best-fitted individuals in the search space. Suppose that the three most diversified individuals at the tth iteration are denoted by $\overrightarrow{X_{\alpha }^{\mathit{div}}\left(t\right)}$, $\overrightarrow{X_{\beta }^{\mathit{div}}\left(t\right)},$ and $\overrightarrow{X_{\delta }^{\mathit{div}}\left(t\right)}$. Then, half of the best-fitted solutions are to be repositioned around these solutions with equal probability based on Equations (15)–(17), as follows:

$$\overrightarrow{X(t+1)}=\overrightarrow{{ X}_{\alpha }^{\mathit{div}}\left(t\right)} \pm \left(\left(\overrightarrow{\mathit{div}}-\overrightarrow{\mathit{lb}}\right).\overrightarrow{r_1}+\overrightarrow{\mathit{lb}}\right); \mathrm{if} 0 \le { r}_4 \le 1/3$$

(15)

$$\overrightarrow{X(t+1) }=\overrightarrow{X_{\beta }^{\mathit{div}}\left(t\right)} \pm \left(\left(\overrightarrow{\mathit{div}}-\overrightarrow{\mathit{lb}}\right).\overrightarrow{r_{2 }}+\overrightarrow{\mathit{lb}}\right){; \mathrm{if} 1/3 < r}_{4 }\le 2/3$$

(16)

$$\overrightarrow{X(t+1) }=\overrightarrow{{ X}_{\delta }^{\mathit{div}}\left(t\right) }\pm \left(\left(\overrightarrow{\mathit{div}}-\overrightarrow{\mathit{lb}}\right)\overrightarrow{r_{3 }}+\overrightarrow{\mathit{lb}}\right){; \mathrm{if} 2/3 < r}_{4 }\le 1$$

(17)

where, $\overrightarrow{\mathit{div}}$ is the upper bound vector of the search space, $\overrightarrow{\mathit{lb}}$ denotes the lower bound vector of the search space, and $\overrightarrow{r_1}$ to $\overrightarrow{r_3}$ are three uniformly distributed random vectors in [0, 1]. $r_4$ is also a random number in [0, 1], delineating which equation to be used to generate $X(t+1)$ as the new position of a solution included in the first half of the best-fitted solutions. In generating $X(t+1)$, the Equations (15)–(17) are to be used. The random number $r_4$ is first generated for each solution e.g. $\overrightarrow{X\left(t\right)}$ that is included in the first half of the best-fitted solutions, and then only one of these three equations is utilized to reposition that solution. Each equation is composed of two terms, including one of the three most-diversified solutions in the search space as the first term and a randomly generated solution in the search space as the second term. Note that the probability of choosing one of these three equations to reposition the focused solution is just the same as the $r_4$ is a uniformly distributed random number and has just the same probability of being chosen when its ranges are equal in length, just as considered in Equations (15)–(17). It is noticeable that in the proposed DB-GWO-EPD algorithm, the critical factor $\overrightarrow{a}$, introduced in Section 2.1, is linearly decreased from 2 to 0, such as the original GWO algorithm. As the last modification applied to the proposed approach, the factor $\overrightarrow{C}$ is eliminated. It is worth noticing that the factor $\overrightarrow{C}$ is used in the updating formulae in the original GWO algorithm to multiply the positions of the prey (alpha, beta, or delta wolves) and redirect the search agents to the random positions away from the accurate positions of the prey. This can help the GWO not to consider the fitness of the prey positions so certain and precise that it can, in turn, contribute to the search agents getting more diversified across the search space and avoiding local optima entrapments. The elimination of such a factor in the proposed DB-GWO-EPD algorithm is because of the well-maintaining of diversity in this algorithm, not necessitating the algorithm to take random positions around the alpha, beta, and delta wolves at each step of the optimization process, unlike the original GWO and the fitness-based GWO-EPD (FB-GWO-EPD). The rest of the calculations and updating equations are similar to the original GWO algorithm. Furthermore, since the proposed approach intrinsically maintains diversity over the whole search space, a random re-initialization process performed on the solutions in the FB-GWO-EPD is unnecessary. Eliminating this re-initialization process in the proposed DB-GWO-EPD can reduce the proposed algorithm’s complexity and runtime, at least when handling global optimization problems with slightly expensive objective function evaluations. The pseudo-code of the proposed DB-GWO-EPD algorithm is described in Algorithm 1. Moreover, the flowchart of the DB-GWO-EPD is depicted in Figure 1.

Algorithm 1: Pseudo-code of the DB-GWO-EPD algorithm.
1:	Initialize the DB-GWO-EPD parameters, a = [2, 0], population size (N), and the maximum number of iterations (M)
2:	whilet < M do
3:	if t = 1 then
4:	for (j = 1:N) do
5:	Initialize random positions for the jth solution at the first iteration as follows:
6:	$\overrightarrow{X_j}=\left(\overrightarrow{\mathit{div}}-\overrightarrow{\mathit{lb}}\right).\overrightarrow{r_j}+\overrightarrow{\mathit{lb}; }\overrightarrow{\mathit{div}}$ and $\overrightarrow{\mathit{lb}}$ are the upper and lower bounds of each dimension and $\overrightarrow{r_j}$ is a random vector
7:	end for
8:	else
9:	for (j = 1:N) do
10:	Calculate $\overrightarrow{X_1}$, $\overrightarrow{X_2}$, and $\overrightarrow{X_3}$, as the guiding solutions of the jth solution, using Equations (6)–(8)
11:	Calculate $\overrightarrow{X}$ as the arithmetic average of $\overrightarrow{X_1}$, $\overrightarrow{X_2}$, and $\overrightarrow{X_3}$, using Equation (9)
12:	Adopt $\overrightarrow{X}$ as the newly updated position of the jth solution at the tth itertaion
13:	end for
14:	end if
15:	Calculate the fitness function value for each solution
16:	Sort the fitness function values and their corresponding solutions in an ascending order for minimization purpose
17:	Identify the first half of the solutions in the solutions sorted as the ones to be repositioned
18:	Identify three best-fitted solutions ($\overrightarrow{X_{\alpha }\left(t\right)}$, $\overrightarrow{X_{\beta }\left(t\right)}$ and $\overrightarrow{X_{\delta }\left(t\right)}$), and save them in the memory of the algorithm
19:	Save $\overrightarrow{X_{\alpha }\left(t\right)}$ as the best-fitted solution found so far and name it $\overrightarrow{X_{\mathit{best}}}$
20:	Calculate the diversity index ($d_i$) for each solution using Equation (14)
21:	Appoint three solution positions with the highest $d_i$ values, as the $\overrightarrow{X_{\alpha }^{\mathit{div}}\left(t\right)}$, $\overrightarrow{X_{\beta }^{\mathit{div}}\left(t\right)},$ and $\overrightarrow{X_{\delta }^{\mathit{div}}\left(t\right)}$
22:	Reposition the first half of the best-fitted solutions around $\overrightarrow{X_{\alpha }^{\mathit{div}}\left(t\right)}$, $\overrightarrow{X_{\beta }^{\mathit{div}}\left(t\right)},$ and $\overrightarrow{X_{\delta }^{\mathit{div}}\left(t\right)}$, randomly, using Equations (15)–(17)
23:	Adopt the repositioned solution positions as their new positions
24:	t = t + 1
25:	end while
27:	Return $\overrightarrow{X_{\mathit{best}}}$ as the final result of the optimization process

The Pros and Cons of the Proposed DB-GWO-EPD and Its Computational Complexity

As comprehensively explained in the previous sub-section, the main advantage of the DB-GWO-EPD against the FB-GWO-EPD and the original GWO may be hidden in its highly enhanced exploration capability, as the search agents are frequently pushed towards the high-fitness region and pulled out of this region to be diversified and be enabled to discover much more good candidate positions in the search space. Another advantage is the elimination of the factor $\overrightarrow{C}$ may be so interesting to be applied to this algorithm, unlike the FB-GWO-EPD and the original GWO. This elimination is of higher importance when the runtime of these algorithms is discussed, especially when implemented on global optimization problems. The proposal also suffers from some weaknesses, the most important of which is difficulty in the intensification of the search process, especially at the later iterations of the optimization, which may result from the unbalanced contribution of the random auxiliary positions added to the diversified solutions in the repositioning procedure. This problem may be mitigated by setting an emotional weight for these random positions that highly decreases over the course of iterations.

The computational complexity of the DB-GWO-EPD is the summation of the complexity of the algorithm in four items: (1) initialization process, (2) mechanism for sorting the elite solutions, (3) calculation of the diversity each solution has in the search space, and (4) position updating procedure of the solutions. The complexity of the initialization process is $O(N \times D)$, and the complexity of the elites’ sorting mechanism as well as the diversity calculation are ${O(M\times N}^2)$. Finally, the complexity of the updating procedure of the solutions is $O(M \times N \times D)$. Consequently, the total computational complexity of the proposed algorithm is $O\left({M \times N}^2)+O(M \times N \times D\right)+O(N \times D)=O\left(M \times N \times (N+D)\right)$. It is worth noticing that $M$ is the maximum number of iterations, $N$ is the population size, and $D$ is the variable decision size.

4. Results and Discussion

4.1. Benchmark Functions

To test the competence of the proposed DB-GWO-EPD algorithm in solving optimization problems, a set of 13 shifted classical benchmark functions, as well as all the test functions, included in the technical report of the 2017 IEEE Congress on Evolutionary Computation (CEC), called CEC2017 test suite are adopted, and the proposal and its competitors are all implemented on these problems [37,38]. The proposed method was deliberately applied to these two sets of high-dimensional benchmark functions with shifted global optima to examine how well it can handle large-scale and hard-to-solve optimization problems compared to its competitors. The CEC2017 test suite was also chosen as another test bed to assess the competence of the optimizers to solve such hard-to-solve problems that can toughly challenge any optimizer when tackling them.

The mathematical formulation of the uni-modal and multi-modal shifted classical benchmark functions, as well as the CEC 2017 test suite, along with their optimal fitness values, can be seen in

Table 1. The shifted uni-modal benchmark functions.

Benchmark Function	Range	Shifted Position	fmin
$F_1\left(X\right)=\sum _{d=1}^n{x}_d^2$	${\left[-100, 100\right]}^n$	$[-30, -30, \dots , -30]$	0
$F_2\left(X\right)=\sum _{d=1}^n\left|x_d\right|+\prod _{d=1}^n\left|x_d\right|$	${\left[-10, 10\right]}^n$	$[-3, -3, \dots , -3]$	0
$F_3\left(X\right)=\sum _{d=1}^n{\left(\sum _{j=1}^d{x}_j\right)}^2$	${\left[-100, 100\right]}^n$	$[-30, -30, \dots , -30]$	0
$F_4\left(X\right)= \underset{d}{\mathit{max}}\left\{\left|x_d\right|,1 \le d \le n\right\}$	${\left[-100, 100\right]}^n$	$[-30, -30, \dots , -30]$	0
$F_5\left(X\right)=\sum _{d=1}^{n-1}\left[{100(x}_{d+1}-x_d^2{)}^2{+(x}_d-{1)}^2\right]$	${\left[-30, 30\right]}^n$	$[-15, -15, \dots , -15]$	0
$F_6\left(X\right)=\sum _{d=1}^n{\left(x_d+0.5\right)}^2$	${\left[-100, 100\right]}^n$	$[-750, -750, \dots , -750]$	4.225${ \times 10}^5$× n
$F_7\left(X\right)=\sum _{d=1}^n{\mathrm{dx}}_d^4+\mathrm{random} [0,1)$	${\left[-1.28, 1.28\right]}^n$	$[-0.25, -0.25, \dots , -0.25]$	0

,

Table 2. The shifted multi-modal benchmark functions.

Benchmark Functions	Range	Shifted Position	fmin
$F_8\left(X\right)=\sum _{d=1}^n-{ x}_d\sin \sqrt{\left|x_d\right|}$	${\left[-500, 500\right]}^n$	$[-300, -300, \dots , -300]$	−418.9829 × n
$F_9\left(X\right)=\sum _{d=1}^n\left[x_d^2- 10\cos \left({2\pi x}_d\right)+10\right]$	${\left[-5.12, 5.12\right]}^n$	$[-2, -2, \dots , -2]$	0
$F_{10}\left(X\right)=-20\exp \left(-0.2\sqrt{\frac{1}{n}\sum _{d=1}^n{x}_d^2}\right)-\exp \left(\frac{1}{n}\sum _{d=1}^n{\cos (2\pi x}_d)\right)+20+e$	${\left[-32, 32\right]}^n$	$[-16, -16, \dots , -16]$	0
$F_{11}\left(X\right)=\frac{1}{4000}\sum _{d=1}^n{x}_d^2-\prod _{d=1}^n\cos (\frac{x_d}{\sqrt{d}})+1$	${\left[-600, 600\right]}^n$	$[-400, -400, \dots , -400]$	0
$F_{12}\left(X\right)=\frac{\pi }{n}\left\{{10\sin (\pi y}_1)+\sum _{d=1}^{n-1}{\left(y_d-1\right)}^2\left[{1+10\sin }^2{(\pi y}_{d+1})\right]{+(y}_n-{1)}^2\right\}+\sum _{d=1}^{n }{u(x}_d, 10, 100, 40)$	${\left[-50, 50\right]}^n$	$[-30, -30, \dots , -30]$	0
$y_{d }=\frac{x_d+1}{4}$ u($x_d,a, k, m)=\left\{\begin{array}{l}{k(x}_d-{a)}^m\hspace{1em}\hspace{1em}{ x}_{d }> a\\ 0\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}-{a < x}_d < a\\ k(-x_d-{a)}^m\hspace{1em}\hspace{1em}{ x}_d < -a\end{array}\right.$
$F_{13}\left(X\right)=0{.1\{\sin }^2\left({3\pi x}_1\right)+\sum _{d=1}^n{\left(x_d-1\right)}^2\left[{1+\sin }^2\left({3\pi x}_d+1\right)\right]+ {\left(x_n-1\right)}^2\left[{1+\sin }^2\left({2\pi x}_n\right)\right]+\sum _{d=1}^n{u(x}_d, 5, 100, 4)$	${\left[-50, 50\right]}^n$	$[-100, -100, \dots , -100]$	*

* The minimum of this function varies as the dimensionality (n) varies, however, if n = 100, then f_min = 4.1006 ${\times 10}^{10}.$

and

Table 3. Summary of the CEC2017 test functions.

Description	NO.	Function	fmin
Uni-modal Functions	1	Shifted and Rotated Bent Cigar Function	100
	2	Shifted and Rotated Sum of Different Power Function *	200
	3	Shifted and Rotated Zakharov Function	300
Simple Multi-modal Functions	4	Shifted and Rotated Rosenbrock’s Function	400
	5	Shifted and Rotated Rastrigin’s Function	500
	6	Shifted and Rotated Expanded Scaffer’s F6 Function	600
	7	Shifted and Rotated Lunacek Bi_Rastrigin Function	700
	8	Shifted and Rotated Non-Continuous Rastrigin’s Function	800
	9	Shifted and Rotated Levy Function	900
	10	Shifted and Rotated Schwefel’s Function	1000
Hybrid Functions	11	Hybrid Function 1 (N = 3)	1100
	12	Hybrid Function 2 (N = 3)	1200
	13	Hybrid Function 3 (N = 3)	1300
	14	Hybrid Function 4 (N = 4)	1400
	15	Hybrid Function 5 (N = 4)	1500
	16	Hybrid Function 6 (N = 4)	1600
	17	Hybrid Function 6 (N = 5)	1700
	18	Hybrid Function 6 (N = 5)	1800
	19	Hybrid Function 6 (N = 5)	1900
	20	Hybrid Function 6 (N = 6)	2000
Composition Functions	21	Composition Function 1 (N = 3)	2100
	22	Composition Function 2 (N = 3)	2200
	23	Composition Function 3 (N = 4)	2300
	24	Composition Function 4 (N = 4)	2400
	25	Composition Function 5 (N = 5)	2500
	26	Composition Function 6 (N = 5)	2600
	27	Composition Function 7 (N = 6)	2700
	28	Composition Function 8 (N = 6)	2800
	29	Composition Function 9 (N = 3)	2900
	30	Composition Function 10 (N = 3)	3000
Search Range: [−100, 100]D; D is the dimensionality of the test problems

* F2 has been excluded because it shows unstable behavior, especially for higher dimensions, and significant performance variations for the same algorithm implemented in Matlab, C.

. In these tables, n denotes the number of dimensions of the corresponding problem. The global optima of these functions are all shifted from the center of the domain to a position away from the center to make these problems more challenging and complex to solve and more realistic for any optimizer when solving them. The shifted positions adopted for the global optima of the classical functions are all set similar to the same in extensive experiments conducted in the literature [35,39].

All of the test functions are to be minimized. The shifted classical benchmark test problems are assumed to be high-dimensional, with the dimensionality set to 100. In addition, the CEC2017 test functions are adopted to be 50-dimensional to increase the difficulty the algorithms may face when solving such hard-to-handle and complex optimization problems. All of these settings are made to examine better the proposed algorithm’s real strength against the other comparative algorithms. The first test bed consists of two categories of uni-modal functions (F1–F7) and multi-modal functions (F8–F13). The first one challenges the exploitation capability of the algorithms and tests their convergence rate of them when finding out the global optimum of the problems, and the second one tests the exploration capability of the algorithms, investigating if the optimization algorithms can avoid local optima entrapments and detect all the possible good candidate solutions in the search space. The CEC2017 test suite encompasses a variety of composition, hybrid, uni-modal, and multi-modal functions. Thus, it can suitably benchmark the maximum ability of the examined algorithms to solve a wide range of problems providing great difficulties to an optimizer when attempting to search for the optimum of such problems.

4.2. Comparison with Other Well-Known Algorithms

The proposed DB-GWO-EPD was compared with the FB-GWO-EPD and the original GWO along with four other popular and newly-proposed algorithms, including Aquila Optimizer (AO) [40]; Flow Direction Algorithm (FDA) [41]; Arithmetic Optimization Algorithm (AOA) [42]; and Gradient-based Optimizer (GBO) [43]. The parameter settings of all the competitive algorithms are presented in

Table 4. Parameter settings of the DB-GWO-EPD and its competitive algorithms.

Algorithm	Parameter Settings
AO	$r_1\in \left[1, 20\right]; U=0{.00565; D}_1=D; \omega =0.005; \alpha =\delta =0{.1; G}_{2 }=\left[2, 0\right]$
FDA	$\beta =1$
AOA	$\alpha =5; \mu =0.5$
GBO	${\beta }_{\min }=0{.2; \beta }_{\mathit{max}}=1.2; \mathrm{pr}=0.5$
GWO	$a=\left[2, 0\right]$
FB-GWO-EPD	$a=\left[2, 0\right]$
DB-GWO-EPD	$a=\left[2, 0\right]$

.

All algorithms were run 30 independent times, and the average and standard deviation of the final results achieved during all runs, abbreviated as “Ave” and “Std”, respectively, are adopted as the performance criteria and reported in

Table 5. Results of the DB-GWO-EPD and its rivals on the shifted benchmark functions, with n = 100.

	Criteria	AO	FDA	AOA	GBO	GWO	FB-GWO-EPD	DB-GWO-EPD
F1	Ave	6.6579 ${\times 10}^2$	9.2465$\times$10	8.1744${ \times 10}^4$	4.2565	2.6986${ \times 10}^4$	5.0429${ \times 10}^2$	8.3539 × 10−1
	Std	3.7689${ \times 10}^2$	2.7576$\times$10	3.5185${ \times 10}^3$	1.8961	4.6165${ \times 10}^3$	1.4285${ \times 10}^2$	3.0210 × 10−1
F2	Ave	6.2467$\times$10	1.2532$\times$10	8.2193${ \times 10}^{46}$	7.1612	1.5454${ \times 10}^2$	1.0983${ \times 10}^2$	1.9838
	Std	1.0766$\times$10	1.2938$\times$10	1.4548${ \times 10}^{47}$	2.5304	1.7258$\times$10	2.2135$\times$10	1.7873
F3	Ave	6.2975${ \times 10}^5$	5.1412${ \times 10}^4$	5.2178${ \times 10}^7$	5.0863${ \times 10}^4$	6.5265${ \times 10}^4$	4.5943${ \times 10}^4$	2.8288 × 104
	Std	3.3353${\times 10}^5$	1.2203${ \times 10}^4$	3.0106${ \times 10}^7$	1.1385${ \times 10}^4$	1.2188${ \times 10}^4$	2.3826 × 103	4.5466${ \times 10}^3$
F4	Ave	6.4199	5.7145$\times$10	3.0010$\times$10	3.0000$\times$10	3.0000$\times$10	3.0001$\times$10	1.0698$\times$10
	Std	1.0959	3.2665	1.0620${ \times 10}^{-2}$	0	1.4235 ${\times 10}^{-5}$	1.1221${ \times 10}^{-3}$	2.7652
F5	Ave	7.1099${ \times 10}^5$	2.6441${ \times 10}^4$	4.2934${ \times 10}^8$	2.4430${ \times 10}^3$	8.4413$\times {10}^7$	1.7271${ \times 10}^3$	7.9457 × 102
	Std	2.3792${ \times 10}^5$	1.0748${ \times 10}^4$	4.8797${ \times 10}^6$	1.0962${ \times 10}^3$	1.3821${ \times 10}^7$	6.9398 × 102	9.8350${ \times 10}^2$
F6	Ave	4.7087${ \times 10}^7$	4.2455${ \times 10}^7$	5.4015${ \times 10}^7$	4.2315 × 107	4.5852${ \times 10}^7$	4.2609${ \times 10}^7$	4.3322${ \times 10}^7$
	Std	1.1035${ \times 10}^6$	2.0457${ \times 10}^5$	5.9905${ \times 10}^5$	0	3.4135${ \times 10}^5$	1.4011${ \times 10}^5$	1.1990${ \times 10}^5$
F7	Ave	5.9662 × 10−2	1.8381	1.9591$\times$10	5.2178${ \times 10}^{-1}$	7.6576	2.0676	2.3332${ \times 10}^{-1}$
	Std	8.1757${\times 10}^{-2}$	2.3114${ \times 10}^{-1}$	6.6384 × 10−2	1.1624${ \times 10}^{-1}$	9.3946${ \times 10}^{-1}$	4.0188${ \times 10}^{-1}$	6.7209${ \times 10}^{-2}$
F8	Ave	−4.0026${ \times 10}^4$	−4.0001${\times 10}^4$	−1.6152${ \times 10}^4$	−4.6366 × 104	−3.0043${ \times 10}^4$	−3.1018${ \times 10}^4$	−2.9074${ \times 10}^4$
	Std	5.0202${ \times 10}^3$	2.8268${ \times 10}^3$	1.4243 × 103	3.3289${ \times 10}^3$	2.0813${ \times 10}^3$	6.5483${ \times 10}^3$	1.0665${ \times 10}^4$
F9	Ave	6.4745 × 10	4.8795${ \times 10}^2$	3.9917${ \times 10}^2$	3.7281${ \times 10}^2$	3.5234${ \times 10}^2$	3.7875${ \times 10}^2$	4.2068${ \times 10}^2$
	Std	2.4065$\times$10	5.6871$\times$10	3.0287 × 10−2	7.5392	1.5079$\times$10	3.7176$\times$10	2.7046${ \times 10}^2$
F10	Ave	9.2096	1.9712$\times$10	1.9185$\times$10	1.0791$\times$10	1.7895$\times$10	1.1858$\times$10	1.1055
	Std	8.7246${ \times 10}^{-1}$	2.6169${ \times 10}^{-1}$	4.3534 × 10−6	2.2658	9.7602${ \times 10}^{-1}$	7.1181	3.0835${ \times 10}^{-1}$
F11	Ave	2.0904${ \times 10}^2$	1.9947$\times$10	2.1440${ \times 10}^3$	1.0766	5.7018${ \times 10}^2$	2.1362	4.9285 × 10−1
	Std	2.5256${ \times 10}^2$	8.7853	1.2316${ \times 10}^2$	7.3605${ \times 10}^{-2}$	1.3720${ \times 10}^2$	6.0558${ \times 10}^{-1}$	7.3081 × 10−2
F12	Ave	9.2327${ \times 10}^2$	1.2972${ \times 10}^3$	1.5442${ \times 10}^9$	1.4721$\times$10	2.0258${ \times 10}^8$	1.2536$\times$10	1.5566
	Std	2.7363${ \times 10}^3$	1.4015${ \times 10}^3$	2.8499${ \times 10}^7$	4.6852	3.6813${ \times 10}^7$	4.0517	6.3168 × 10−1
F13	Ave	4.7531${ \times 10}^{10}$	4.1006 × 1010	7.9630${ \times 10}^{11}$	4.1006 × 1010	1.8052${ \times 10}^{11}$	4.2906${ \times 10}^{10}$	5.3127${ \times 10}^{10}$
	Std	1.5675${ \times 10}^{10}$	0	1.2203${ \times 10}^{10}$	0	3.8026${ \times 10}^{10}$	1.6574${ \times 10}^8$	1.0299${ \times 10}^9$

and

Table 6. Results of the DB-GWO-EPD and its rivals on the 50-dimensional CEC2017 test functions.

	Criteria	AO	FDA	AOA	GBO	GWO	FB-GWO-EPD	DB-GWO-EPD
F1	Ave	1.6838${ \times 10}^9$	3.8271${ \times 10}^4$	1.0679${ \times 10}^{11}$	7.2920 × 103	8.2896 ${\times 10}^9$	5.8980${ \times 10}^7$	1.8071${ \times 10}^4$
	Std	6.2457${ \times 10}^8$	4.4200${ \times 10}^4$	1.0241${ \times 10}^{10}$	8.4348 × 103	4.2006${ \times 10}^9$	2.0039${ \times 10}^7$	9.2804${ \times 10}^3$
F3	Ave	2.0248${ \times 10}^5$	2.7961 × 104	1.6642${ \times 10}^5$	3.9352${ \times 10}^4$	1.0484${ \times 10}^5$	6.6347${ \times 10}^4$	5.0470${ \times 10}^4$
	Std	4.4288${ \times 10}^4$	8.9157 × 103	2.0180${ \times 10}^4$	9.7975${ \times 10}^3$	1.6965${ \times 10}^4$	1.3566${ \times 10}^4$	1.2134${ \times 10}^4$
F4	Ave	1.0630${ \times 10}^3$	5.4393 × 102	3.0754${ \times 10}^4$	5.5786${ \times 10}^2$	1.1238${ \times 10}^3$	6.6857${ \times 10}^2$	5.9875${ \times 10}^2$
	Std	1.9549${ \times 10}^2$	5.8808$\times 10$	7.7284${ \times 10}^3$	5.3041$\times 10$	2.3977${ \times 10}^2$	4.1381 × 10	5.2919$\times 10$
F5	Ave	8.6665${ \times 10}^2$	8.2013${ \times 10}^2$	1.1525${ \times 10}^3$	8.1004${ \times 10}^2$	7.2367${ \times 10}^2$	6.7674${ \times 10}^2$	6.5307 × 102
	Std	3.3186$\times 10$	6.2431$\times 10$	4.0010$\times 10$	4.9971$\times 10$	2.7447 × 10	4.7522$\times 10$	9.5217$\times 10$
F6	Ave	6.6537${ \times 10}^2$	6.5096${ \times 10}^2$	6.9011${ \times 10}^2$	6.3725${ \times 10}^2$	6.1746${ \times 10}^2$	6.1596${ \times 10}^2$	6.0456 × 102
	Std	5.3658$\times 10$	8.7769	5.1072	9.7671	5.4024	5.9317	2.9473
F7	Ave	1.5161${ \times 10}^3$	1.4326${ \times 10}^3$	1.9341${ \times 10}^3$	1.2644${ \times 10}^3$	1.0724${ \times 10}^3$	1.0822${ \times 10}^3$	9.7839 × 102
	Std	1.2247${ \times 10}^2$	1.1392${ \times 10}^2$	6.3770$\times 10$	9.4996$\times 10$	5.9392$\times 10$	4.4810 × 10	1.5911${ \times 10}^2$
F8	Ave	1.1804${ \times 10}^3$	1.1620${ \times 10}^3$	1.4717${ \times 10}^3$	1.1110${ \times 10}^3$	1.0259${ \times 10}^3$	1.0124${ \times 10}^3$	9.1755 × 102
	Std	3.4016$\times 10$	4.3580$\times 10$	4.7283$\times 10$	5.6324$\times 10$	5.7317$\times 10$	9.3255$\times 10$	2.7887 × 10
F9	Ave	2.1934${ \times 10}^4$	1.0603${ \times 10}^4$	2.8784${ \times 10}^4$	7.4002${ \times 10}^3$	8.0871${ \times 10}^3$	6.2582${ \times 10}^3$	1.5274 × 103
	Std	3.8379${ \times 10}^3$	2.4207${ \times 10}^3$	4.3133${ \times 10}^3$	2.4743${ \times 10}^3$	3.3907${ \times 10}^3$	2.5762${ \times 10}^3$	1.2381 × 103
F10	Ave	9.1459${ \times 10}^3$	8.4369${ \times 10}^3$	1.3462${ \times 10}^4$	7.8672${ \times 10}^3$	7.5159 × 103	9.8846${ \times 10}^3$	8.1947${ \times 10}^3$
	Std	9.7917${ \times 10}^2$	1.0053${ \times 10}^3$	7.8133 × 102	9.4913${ \times 10}^2$	1.7420${ \times 10}^3$	4.0438${ \times 10}^3$	3.5965${ \times 10}^3$
F11	Ave	2.2437${ \times 10}^3$	1.3379 × 103	2.2605${ \times 10}^4$	1.3879${ \times 10}^3$	4.3152${ \times 10}^3$	1.5887${ \times 10}^3$	1.4279${ \times 10}^3$
	Std	2.7895${ \times 10}^2$	7.0421 × 10	3.8152${ \times 10}^3$	8.3613$\times 10$	1.5839${ \times 10}^3$	1.0617${ \times 10}^2$	8.1144$\times 10$
F12	Ave	6.1213${ \times 10}^8$	2.2862 × 106	6.8321${ \times 10}^{10}$	2.6664${ \times 10}^6$	5.9699${ \times 10}^8$	1.5015${ \times 10}^8$	2.3364${ \times 10}^7$
	Std	4.0105${ \times 10}^8$	1.4679 × 106	1.5355${ \times 10}^{10}$	2.2857${ \times 10}^6$	6.4348${ \times 10}^8$	8.5096${ \times 10}^7$	1.7026${ \times 10}^7$
F13	Ave	2.3863${ \times 10}^7$	6.1931 × 103	3.9481${\times 10}^{10}$	1.2570${ \times 10}^4$	4.1420${ \times 10}^8$	8.5343${ \times 10}^5$	8.4775${ \times 10}^4$
	Std	3.8111${ \times 10}^7$	7.0134 × 103	1.2296${ \times 10}^{10}$	9.9819${ \times 10}^3$	8.8631${ \times 10}^8$	4.4079${ \times 10}^5$	4.8870${ \times 10}^4$
F14	Ave	4.8743${ \times 10}^6$	3.1987 × 104	5.2754${ \times 10}^7$	3.9591${ \times 10}^4$	1.0451${ \times 10}^6$	3.7271${ \times 10}^5$	1.9834${ \times 10}^5$
	Std	4.3827${ \times 10}^6$	3.1930 × 104	4.4620${ \times 10}^7$	3.9226${ \times 10}^4$	1.2344${ \times 10}^6$	2.3847${ \times 10}^5$	9.5941${ \times 10}^4$
F15	Ave	6.1620${ \times 10}^5$	1.0266 × 104	4.4614${ \times 10}^9$	1.2206${ \times 10}^4$	1.8441${ \times 10}^7$	1.4188${ \times 10}^5$	4.2204${ \times 10}^4$
	Std	3.6634${ \times 10}^5$	6.4333 × 103	2.6802${ \times 10}^9$	7.4912${ \times 10}^3$	3.2540${ \times 10}^7$	1.8065${ \times 10}^5$	2.4322${ \times 10}^4$
F16	Ave	4.3612${ \times 10}^3$	3.6774${ \times 10}^3$	7.9190${ \times 10}^3$	3.4870${ \times 10}^3$	3.2102${ \times 10}^3$	2.9626${ \times 10}^3$	2.8461 × 103
	Std	5.3030${ \times 10}^2$	5.1378${ \times 10}^2$	1.2332${ \times 10}^3$	4.9718${ \times 10}^2$	4.4746${ \times 10}^2$	4.1833${ \times 10}^2$	3.4984 × 102
F17	Ave	3.6432${ \times 10}^3$	3.3813${ \times 10}^3$	9.1957${ \times 10}^3$	3.0819${ \times 10}^3$	2.9329${ \times 10}^3$	2.8919${ \times 10}^3$	2.8425 × 103
	Std	3.8545${ \times 10}^2$	3.7707${ \times 10}^2$	2.4943${ \times 10}^3$	3.4783${ \times 10}^2$	3.4609 × 102	3.9359${ \times 10}^2$	4.1003${ \times 10}^2$
F18	Ave	9.1232${ \times 10}^6$	2.2535${ \times 10}^5$	1.0356${ \times 10}^8$	2.1870 × 105	4.4394${ \times 10}^6$	2.9262${ \times 10}^6$	2.3528${ \times 10}^6$
	Std	6.4256${ \times 10}^6$	1.5052${ \times 10}^5$	4.8975${ \times 10}^7$	1.3210 × 105	5.2595${ \times 10}^6$	1.9711${ \times 10}^6$	2.3575${ \times 10}^6$
F19	Ave	2.2803${ \times 10}^6$	1.8197 × 104	2.7898${ \times 10}^9$	1.8302${ \times 10}^4$	4.1561${ \times 10}^6$	1.4102${ \times 10}^6$	8.0552${ \times 10}^5$
	Std	2.1570${ \times 10}^6$	1.0490 × 104	1.4850${ \times 10}^9$	1.1745${ \times 10}^4$	8.5345${ \times 10}^6$	1.0476${ \times 10}^6$	5.9277${ \times 10}^5$
F20	Ave	3.2709${ \times 10}^3$	3.4256 $\times {10}^3$	3.5912${ \times 10}^3$	3.2025${ \times 10}^3$	2.9462 × 103	3.1085${ \times 10}^3$	2.9489${ \times 10}^3$
	Std	2.6175 × 102	3.3007 ${\times 10}^2$	2.6817 ${\times 10}^2$	4.0586${ \times 10}^2$	3.6343${ \times 10}^2$	5.4458 ${\times 10}^2$	4.7087${ \times 10}^2$
F21	Ave	2.7043${ \times 10}^3$	2.6240${ \times 10}^3$	3.0782${ \times 10}^3$	2.5662${ \times 10}^3$	2.5173${ \times 10}^3$	2.4638${ \times 10}^3$	2.4115 × 103
	Std	6.3075$\times 10$	6.1326$\times 10$	8.5153$\times 10$	5.1954$\times 10$	5.5484$\times 10$	5.0296$\times 10$	2.3648 × 10
F22	Ave	1.0965${ \times 10}^4$	9.9308${ \times 10}^3$	1.5926${ \times 10}^4$	9.4638${ \times 10}^3$	9.6836${ \times 10}^3$	1.0955${ \times 10}^4$	8.7319 × 103
	Std	1.6841${ \times 10}^3$	8.3742${ \times 10}^2$	7.0190 × 102	1.6273${ \times 10}^3$	2.0126${ \times 10}^3$	3.9900${ \times 10}^3$	3.0635${ \times 10}^3$
F23	Ave	3.4393${ \times 10}^3$	3.1004${ \times 10}^3$	4.4239${ \times 10}^3$	3.0558${ \times 10}^3$	2.9820${ \times 10}^3$	2.9274${ \times 10}^3$	2.8798 × 103
	Std	9.2377$\times 10$	8.1408$\times 10$	2.3112${ \times 10}^2$	7.8350$\times 10$	6.3290 × 10	8.7578$\times 10$	6.8065$\times 10$
F24	Ave	3.5218${ \times 10}^3$	3.2665${ \times 10}^3$	4.9199${ \times 10}^3$	3.1724${ \times 10}^3$	3.1932${ \times 10}^3$	3.0855${ \times 10}^3$	3.0343 × 103
	Std	1.2135${ \times 10}^2$	9.9971$\times 10$	3.1459${ \times 10}^2$	6.0068 × 10	1.1563${ \times 10}^2$	1.2804${ \times 10}^2$	7.5755 $\times 10$
F25	Ave	3.4327${ \times 10}^3$	3.0855${ \times 10}^3$	1.5444${ \times 10}^4$	3.0855${ \times 10}^3$	3.5262${ \times 10}^3$	3.2034${ \times 10}^3$	3.0668 × 103
	Std	8.8058$\times 10$	2.3492 × 10	1.3989${ \times 10}^3$	2.3553$\times 10$	2.1484${ \times 10}^2$	5.9425$\times 10$	2.6992$\times 10$
F26	Ave	8.6249${ \times 10}^3$	9.0556${ \times 10}^3$	1.6923${ \times 10}^4$	7.0615${ \times 10}^3$	6.3760${ \times 10}^3$	5.8852${ \times 10}^3$	5.0518 × 103
	Std	2.4721${ \times 10}^3$	1.8166${ \times 10}^3$	1.2165${ \times 10}^3$	2.4977${ \times 10}^3$	5.2971${ \times 10}^2$	8.0009${ \times 10}^2$	3.4836 × 102
F27	Ave	4.0141${ \times 10}^3$	3.6061${ \times 10}^3$	6.7731${ \times 10}^3$	3.5950${ \times 10}^3$	3.6165${ \times 10}^3$	3.4305${ \times 10}^3$	3.4253 × 103
	Std	1.9720${ \times 10}^2$	1.5115${ \times 10}^2$	7.3983${ \times 10}^2$	1.3302${ \times 10}^2$	1.2148${ \times 10}^2$	5.3921 × 10	6.6431$\times 10$
F28	Ave	4.2899${ \times 10}^3$	3.3315${ \times 10}^3$	1.2297${ \times 10}^4$	3.3326${ \times 10}^3$	4.2560${ \times 10}^3$	3.4948${ \times 10}^3$	3.3136 × 103
	Std	2.7492${ \times 10}^2$	2.9990$\times 10$	1.3466${ \times 10}^3$	2.7941$\times 10$	4.0559${ \times 10}^2$	8.1306$\times 10$	2.5638 × 10
F29	Ave	6.1689${ \times 10}^3$	4.7010${ \times 10}^3$	3.6559${ \times 10}^4$	4.6781${ \times 10}^3$	4.6299${ \times 10}^3$	4.3602${ \times 10}^3$	4.2874 × 103
	Std	6.9454${ \times 10}^2$	4.4253${ \times 10}^2$	2.9085${ \times 10}^4$	3.7906${ \times 10}^2$	2.8006 × 102	2.8032${ \times 10}^2$	3.0278${ \times 10}^2$
F30	Ave	1.2248${ \times 10}^8$	1.1529${ \times 10}^6$	5.8046${ \times 10}^9$	1.0718 × 106	1.1702${ \times 10}^8$	9.2911${ \times 10}^7$	3.9885${ \times 10}^7$
	Std	4.9406${ \times 10}^7$	2.9984${ \times 10}^5$	2.6383${ \times 10}^9$	2.1165 × 105	5.0559 ${\times 10}^7$	2.0213${ \times 10}^7$	7.6628${ \times 10}^6$

. The results corresponding to the best-performing criteria are emboldened in these tables. The convergence curves of the algorithms, when implemented on the shifted benchmark functions, are also plotted in Figure 2 and Figure 3. Furthermore, for a fair and unbiased comparison, the population size of all the algorithms is set to 30 for the first benchmark function set and adopted to be 50 for the CEC2017 test suite. In addition, all the algorithms were implemented for 1000 iterations at each runtime. The stopping criterion is also considered fulfilled after reaching the maximum number of iterations.

4.2.1. Comparing the Algorithms on the Uni-Modal Functions

The results of applying the algorithms to the uni-modal benchmark functions show that the DB-GWO-EPD highly outperforms the original GWO and the FB-GWO-EPD as the improved version of the GWO. While the FB-GWO-EPD is slightly better than the GWO on almost all uni-modal functions, the DB-GWO-EPD performs much better than both other algorithms on the entire test functions falling in this category except for F6.

The DB-GWO-EPD can reach the slightest standard deviations among its competitors on F1, F2, F3, and F7. This can guarantee the use of this algorithm without any need to run the algorithm frequently. In other words, the user of the proposed algorithm can rely much more on the results of the DB-GWO-EPD compared to the other version of GWO. The proposed algorithm is superior to all the examined algorithms, including the four newly-proposed algorithms on F1, F2, F3, and F5. The closest rival of the proposed DB-GWO-EPD is the AO algorithm, outperforming its competitors on F4 and F7. Overall, the proposal is superior to the other algorithms on 6 out of 14 (43%) of the entire performance criteria (average and standard deviations of the best-achieved solutions). In comparison, each of the other examinees is superior to the others on a maximum of only 2 out of 14 (14%) criteria.

The main reason behind the superior results of the proposed algorithm is hidden in having a high exploration capability which, in turn, can impede the stagnation of the solutions in the search space by frequently diversifying them to accelerate the process of finding the optimal solutions across the search space and speed up the convergence of the proposal when implemented on such a high-dimensional uni-modal function set.

4.2.2. Comparing the Algorithms on the Multi-Modal Functions

The results of the multi-modal functions indicate the superiority of the proposed DB-GWO-EPD algorithm to the FB-GWO-EPD and the original GWO on most of the benchmark functions. The high capability of the proposed algorithm in the diversification of the search agents takes effect, especially when solving F10, F11, and F12. Compared with the other examined algorithms, including the popular and newly-proposed ones, the DB-GWO-EPD shows its superiority on 5 out of 12 (42%) of the total criteria, demonstrating its high efficacy in solving this set of problems.

Knowing that the multi-modal functions with many dimensions, as examined in this study, have a lot of local optima in their search space, an algorithm should have a powerful exploration capability to overcome these natural difficulties the multi-modal functions offer. Moreover, in all of the benchmark functions employed in this study, the optimal solution is shifted in all dimensions, and this alteration applied to these functions can make them significantly harder to solve. It can be inferred from the results obtained that the proposed DB-GWO-EPD algorithm is versatile and robust, tackling the difficulties mentioned for the high-dimensional multi-modal functions as the most realistic representatives of real-world optimization problems.

4.2.3. Comparing the Algorithms on the CEC2017 Test Suite

The CEC2017 benchmark functions are rated as a test bed that can challenge any optimization algorithm comprehensively. The problems in this set seem more realistic and more similar to the real-world optimization problems as compared to those in the first 13-problem set used to examine the eligibility of the proposed algorithm and its rivals. This test suite comprises a variety of composition, hybrid, uni-modal, and multi-modal functions. Here, all the functions included in the CEC2017 test suite are set to be 50-dimensional to challenge further the capabilities of the proposal to handle such large-scale and complex optimization problems.

As the results suggest, the proposed DB-GWO-EPD significantly outperforms the other comparative algorithms on F21–F29 problems. These problems are known as the composition functions comprising a variety of shifted, rotated, and biased multi-modal functions, representing them as the most challenging category of the test problems included in the CEC2017 benchmark functions. As a result, the superiority of the proposal on these functions can affirm the power of this algorithm to solve hard-to-solve optimization problems. Although the proposed DB-GWO-EPD performs much better than the FB-GWO-EPD and the original GWO algorithms on 10 out of 10 (100%) of the F11–F20 problems, the proposal loses the competition in solving these problems to the FDA algorithm having 60% outperformance rate that overtakes the DB-GWO-EPD with an outperformance rate of 15% on the F11–F20. This problem set is called the hybrid functions, the second most challenging category in the CEC2017 test suite. The proposal also outperforms both the original GWO and the FB-GWO-EPD on the first ten functions of this test bed that can technically be divided into uni-modal and multi-modal functions.

The proposal is also superior to all rivals on 5 out of 9 (56%) of the F1–F10 problems. Note that none of the algorithms is tested on F2, as recommended by the creators of the CEC2017 test suite. FDA is rated as the second-best algorithm outperforming the other examinees on 2 out of 9 (22%) of these problems. Overall, the proposed algorithm can achieve the best average objective function values on 16 out of 29 (55%) criteria in the CEC2017 test bed. In addition, the DB-GWO-EPD tackles its rivals in 23 out of 58 (40%) cases, including the average and standard deviation of the best-reached results by the optimizers on this test suite.

The main reason for the high competence of the proposed DB-GWO-EPD algorithm may be the strong capability of the proposed algorithm to preserve the diversity of the solutions during the whole optimization process, which especially matches the nature of this test problems suite. This property incorporated into the proposal can make this algorithm tackle its competitors when facing many local optima in the search space of multi-modal problems. In this case, the proposed algorithm explores the search space as much as possible as a result of being enabled to escape from such a large number of local optima and finds the best-fitted region and, finally, the optimal solution to the problem, especially when handling complex optimization problems. Preserving diversity not only guarantees reaching the optimal solution of a multi-modal problem but could also accelerate the convergence in solving the uni-modal problems when their global optimum is shifted to the corner of the problem domain, similar to the cases examined in this paper.

4.3. Statistical Analysis

In statistical analysis, the Wilcoxon rank-sum test is used as a non-parametric test to determine whether a couple of solutions follow the same probability distribution function (García et al., 2008). This test finally renders an output named p-value. If the p-value < 0.05 for an algorithm, it is assumed to be statistically significantly different compared to its competitor. This difference can appear in both cases of superiority and non-superiority in the algorithms’ results. The test results are presented in

Table 7. p-values calculated for the shifted classical functions (p-values $\ge$ 0.05 have been underlined).

Algorithms	AO	FDA	AOA	GBO	GWO	FB-GWO-EPD	DB-GWO-EPD
AO	N/A
FDA	7.7641 × 10−1	N/A
AOA	1.0217 × 10−3	2.5245 × 10−3	N/A
GBO	4.6500 × 10−2	3.2649 × 10−2	3.9261 × 10−5	N/A
GWO	1.5966 × 10−2	1.4304 × 10−2	5.1522 × 10−2	8.1663 × 10−5	N/A
FB-GWO-EPD	8.4754 × 10−1	8.7630 × 10−1	9.8442 × 10−4	4.6674 × 10−2	1.4854 × 10−2	N/A
DB-GWO-EPD	1.2496 × 10−1	1.2324 × 10−1	3.3946 × 10−4	9.1839 × 10−1	9.5116 × 10−4	1.9343 × 10−1	N/A

and

Table 8. p-values calculated for the CEC2017 benchmark functions (p-values $\ge$ 0.05 have been underlined).

Algorithms	AO	FDA	AOA	GBO	GWO	FB-GWO-EPD	DB-GWO-EPD
AO	N/A
FDA	4.4360 × 10−3	N/A
AOA	5.4935 × 10−15	4.0441 × 10−17	N/A
GBO	2.5666 × 10−3	8.0087 × 10−1	3.2827 × 10−16	N/A
GWO	9.3737 × 10−2	2.2604 × 10−1	3.9835 × 10−17	1.1853 × 10−1	N/A
FB-GWO-EPD	4.7753 × 10−3	3.5684 × 10−1	1.6881 × 10−15	3.8885 × 10−1	1.5497 × 10−1	N/A
DB-GWO-EPD	4.5442 × 10−10	7.3072 × 10−3	4.8602 × 10−18	5.4212 × 10−4	8.7465 × 10−8	1.5154 × 10−6	N/A

. In these tables, the expression N/A mentioned in each row denotes that the algorithm ahead of that row is “Not Applicable” in the test, meaning this algorithm cannot be compared. On the shifted classical benchmark functions, the proposed DB-GWO-EPD is assigned a p-value more significant than the level of 5% significance versus AO and GBO, despite being superior to them on a vast majority of the functions. The critical point is that the DB-GWO-EPD cannot achieve a p-value less than 0.05 when facing the FB-GWO-EPD. However, the p-value that the test produces for DB-GWO-EPD versus the original GWO is much less than that the test offers for the FB-GWO-EPD versus GWO, implicitly demonstrating that the results of the proposed DB-GWO-EPD are more significantly superior to the FB-GWO-EPD. While on the CEC2017 test problems, the proposed DB-GWO-EPD is significantly superior to all its rivals, as the p-values generated pairwise for this algorithm against the others clearly suggest. Since the CEC2017 test suite is a more appropriate bed for evaluating the overall performance of an algorithm, having significant superiority to a set of popular algorithms can properly delineate the high competence of the DB-GWO-EPD in solving optimization problems.

4.4. Comparative Results on Engineering Design Problems

An optimization algorithm may face some optimization problems having different difficulties, which in turn, can make these problems too hard-to-solve. The DB-GWO-EPD is implemented on four constrained engineering problems to validate further the proposal to tackle these difficulties. For constraint handling in these problems, the penalty function approach was adopted. In addition, 50 wolves, and 1000 iterations, are set for the DB-GWO-EPD algorithm when solving this problem set. The algorithm is run 30 times, and the best-achieved results are reported versus those of the other comparative algorithms. It is worth mentioning that the results found by the best-performing algorithm are emboldened in the tables showing the results of the optimization conducted on these problems.

4.4.1. Welded Beam Design Problem

This problem aims to minimize the construction cost of a welded beam. This problem is a high-constrained one, imposing several restrictions on the beam properties such as bending stress ($\sigma$), shear stress ($\tau$), deflection at the end of the beam ($\delta$), buckling load on the bar ($P$), and several other constraints. Four design components considered as the decision variables of the problem are illustrated in Figure 4 as ${h(x}_1)$, ${l(x}_2)$, ${t(x}_3)$, and ${b(x}_4)$. The problem formulation is as follows.

$$\mathrm{Minimize} f\left(\overrightarrow{x}\right)=1{.10471x}_1^2{x}_2+0{.04811x}_3{x}_4\left({14+x}_2\right)$$

(18)

$\mathrm{Subject} \mathrm{to}:$

$$g_1\left(\overrightarrow{x}\right)=\tau \left(x\right)-{ \tau }_{\mathit{max} }\le 0$$

(19)

$$g_2\left(\overrightarrow{x}\right)=\sigma \left(x\right)-{ \sigma }_{\mathit{max} }\le 0$$

(20)

$$g_3\left(\overrightarrow{x}\right){=x}_1-x_4\le 0$$

(21)

$$g_4\left(\overrightarrow{x}\right)=0{.10471x}_1^2+0{.04811x}_3{x}_4\left({14+x}_2\right)- 5 \le 0$$

(22)

$$g_5\left(\overrightarrow{x}\right)=0.125 -{ x}_1\le 0$$

(23)

$$g_6\left(\overrightarrow{x}\right)=\delta \left(x\right)-{ \delta }_{\mathit{max} }\le 0$$

(24)

$$g_7\left(\overrightarrow{x}\right)=P -{ P}_c\left(x\right) \le 0$$

(25)

$$0.1 \le { x}_i \le 2, i=1, 4$$

(26)

$$0.1 \le { x}_{i }\le 10, i=2, 3$$

(27)

where

$$\tau \left(x\right)=\sqrt{{{(\tau }^{\prime })}^2{+2\tau }^{\prime }{\tau }^{″}\frac{x_2}{2R}{+(\tau }^{″}{)}^2 }{, \tau }^{\prime }=\frac{P}{\sqrt{2}{x}_1{x}_2}{, \tau }^{″}=\frac{\mathit{MR}}{J}$$

(28)

$$M=P\left(L+\frac{x_2}{2}\right), R=\sqrt{\frac{x_2^2}{4}+{\left(\frac{x_1{+x}_3}{2}\right)}^2}, J=2\left\{\sqrt{2}{x}_1{x}_2\left[\frac{x_2^2}{12}+{\left(\frac{x_1{+x}_3}{2}\right)}^2\right]\right\}$$

(29)

$$\sigma \left(x\right)=\frac{6\mathit{PL}}{x_4{x}_3^2}, \delta \left(x\right)=\frac{{4\mathit{PL}}^3}{{\mathit{Ex}}_3^3{x}_4}{, P}_c\left(x\right)=\frac{4.013E\sqrt{\left(x_3^2{x}_4^6/36\right)}}{L^2}\left(1-\frac{x_3}{2L}\sqrt{\frac{E}{4G}}\right)$$

(30)

$$P=26.698 \mathit{KN}, L=35.56 \mathrm{cm}, E=2{.07 \times 10}^4 \frac{\mathit{KN}}{{\mathrm{cm}}^2}, G=8270 \frac{\mathit{KN}}{{\mathrm{cm}}^2}$$

(31)

$${\tau }_{\mathit{max} }=9.38 \frac{\mathit{KN}}{{\mathrm{cm}}^2}{, \sigma }_{\mathit{max}}=20.7 \frac{\mathit{KN}}{{\mathrm{cm}}^2}{, \delta }_{\mathit{max}}=6.35 \mathrm{mm}$$

(32)

Table 9. Results of solving the welded beam design problem.

Algorithm	x1	x2	x3	x4	fmin
SIMPLEX [45]	0.279200	5.625600	7.751200	0.279600	2.530700
DAVID [45]	0.243400	6.255200	8.291500	0.244400	2.384100
APPROX [45]	0.244400	6.218900	8.291500	0.244400	2.381500
GA [46]	0.248900	6.173000	8.178900	0.253300	2.430000
HS [47]	0.244200	6.223100	8.291500	0.240000	2.380700
CSCA [48]	0.203137	3.542998	9.033498	0.206179	1.733461
CPSO [49]	0.202369	3.544214	9.048210	0.205723	1.728020
RO [50]	0.203687	3.528467	9.004233	0.207241	1.735344
WOA [51]	0.205396	3.484293	9.037426	0.206276	1.730499
GSA [7]	0.182129	3.856979	10.000000	0.202376	1.879950
MVO [52]	0.205463	3.473193	9.044502	0.205695	1.726450
OBSCA [53]	0.230824	3.069152	8.988479	0.208795	1.722315
AOA [42]	0.194475	2.570920	10.000000	0.201827	1.716400
GWO	0.205711	3.254161	9.035520	0.205794	1.695653
FB-GWO-EPD	0.205486	3.259422	9.036860	0.205771	1.696094
DB-GWO-EPD	0.205699	3.253667	9.036660	0.205729	1.695281

shows the design variables and optimal cost to spend for designing such a welded beam when calculated by the proposed algorithm and its competitors. The results reveal that the DB-GWO-EPD can find the least cost among the other algorithms whose results are from the literature.

4.4.2. Three-Bar Truss Design Problem

This problem aims to minimize the weight of a three-bar truss subject to stress, deflection, and buckling constraints. A scheme of this problem is depicted in Figure 5. The truss focused on this problem has three cross-section areas denoted by A1, A2, and A3, out of which A1 is just equal to A3. Thus, there are only two areas to be optimized, which are denoted by $x_1$ and $x_2$ as the decision variables in the problem, formulation are illustrated below.

$$\mathrm{Minimize} f\left(\overrightarrow{x}\right)=\left(2\sqrt{2}{x}_1{+x}_2\right) \times l$$

(33)

$\mathrm{Subject} \mathrm{to}:$

$$g_1\left(\overrightarrow{x}\right)=\frac{\sqrt{2}{x}_1{+x}_2}{\sqrt{2}{x}_1^2{+2x}_1{x}_2}P - \sigma \le 0$$

(34)

$$g_2\left(\overrightarrow{x}\right)=\frac{x_2}{\sqrt{2}{x}_1^2{+2x}_1{x}_2}P - \sigma \le 0$$

(35)

$$g_3\left(\overrightarrow{x}\right)=\frac{1}{\sqrt{2}{x}_{2 }{+x}_1}P - \sigma \le 0$$

(36)

$$0 \le { x}_1{, x}_2 \le 1$$

(37)

$\mathrm{where}$

$$l=100 \mathrm{cm}, P=2\frac{\mathit{KN}}{{\mathrm{cm}}^2}, \sigma =2\frac{\mathit{KN}}{{\mathrm{cm}}^2}$$

(38)

Table 10. Results of solving three-bar truss design problem.

Algorithm	$x_1$	$x_2$	$g_1$	$g_2$	$g_3$	fmin
DEDS [54]	0.788675	0.408248	1.777971 × 10 −8	−1.464102	−0.535898	263.895841
SSA [55]	0.788665	0.408276	1.247906 × 10 −8	−1.464070	−0.535930	263.895842
MBA [56]	0.788565	0.408560	1.418869 × 10 −7	−1.463748	−0.536252	263.895834
PSO-DE [57]	0.788675	0.408248	1.427175 × 10 −7	−1.464102	−0.535898	263.895825
Tsa [58]	0.788000	0.408000	1.636731 × 10 −3	−1.463566	−0.534798	263.680057
Rai and Saini [59]	0.795000	0.395000	−3.375515${\times 10}^{-3}$	−1.480901	−0.522474	264.359956
CS [60]	0.788670	0.409020	−5.733901${\times 10}^{-4}$	−1.463512	−0.537062	263.971562
MFO [61]	0.788245	0.409467	−1.190244${\times 10}^{-9}$	−1.462717	−0.537283	263.895980
AOA [42]	0.793690	0.394260	−1.103007${\times 10}^{-5}$	−1.480113	−0.519898	263.915432
GWO	0.788398	0.409034	−8.468448${\times 10}^{-7}$	−1.463210	−0.536791	263.896012
FB-GWO-EPD	0.788911	0.407583	−5.171678${\times 10}^{-7}$	−1.464858	−0.535143	263.895952
DB-GWO-EPD	0.788539	0.408634	−5.484502 × 10−14	−1.463663	−0.536337	263.895857

shows the results of numerous algorithms attempting to solve this problem. As can be seen in the table, the results found by the first five algorithms in the table violate the constraint $g_1$. Hence, their optimal objective function values are overridden. All the constraints violated are underlined in the table. Among the remaining algorithms, the proposed DB-GWO-EPD offers the best objective value.

4.4.3. Cantilever Beam Design Problem

The goal of this problem is the minimization of the weight of a cantilever beam made of five hollow blocks with constant thicknesses and variable heights, as depicted in Figure 6. Thus, the blocks’ heights are the decision variables to be optimized, denoted by $x_1$ to $x_5$. The problem formulation is as follows.

$$\mathrm{Minimize} f\left(\overrightarrow{x}\right)=0.06224\left(x_{1 }{+x}_2{+x}_3{+x}_4{+x}_5\right)$$

(39)

$\mathrm{Subject} \mathrm{to}:$

$$g\left(\overrightarrow{x}\right)=\frac{61}{x_1^3}+\frac{27}{x_2^3}+\frac{19}{x_3^3}+\frac{7}{x_4^3}+\frac{1}{x_5^3}- 1 \le 0$$

(40)

$$0.01 \le { x}_1{, x}_2{, x}_3{, x}_4{, x}_5 \le 100$$

(41)

As the results in

Table 11. Results of solving cantilever beam design problem.

Algorithm	x1	x2	x3	x4	x5	fmin
MFO [61]	5.98487	5.31673	4.49733	3.51362	2.16162	1.33999
SOS [62]	6.01878	5.30344	4.49587	3.49896	2.15564	1.33996
CS [60]	6.00890	5.30490	4.50230	3.50770	2.15040	1.33999
MMA [63]	6.01000	5.30000	4.49000	3.49000	2.15000	1.34000
GCA1 [63]	6.01000	5.30000	4.49000	3.49000	2.15000	1.34000
GCA2 [63]	6.01000	5.30000	4.49000	3.49000	2.15000	1.34000
GWO	6.02059	5.31527	4.48660	3.50436	2.14698	1.33653
FB-GWO-EPD	6.00974	5.31982	4.48908	3.51235	2.14305	1.33654
DB-GWO-EPD	6.00305	5.30003	4.50193	3.51318	2.15568	1.33653

suggest, both the DB-GWO-EPD and GWO reach the lowest weight of the cantilever beam compared to those of their competitors.

4.4.4. Gas Transmission Compressor Design Problem

This problem is designed to minimize the cost of delivering 100 million cft gas daily. The decision variables $x_1$, $x_2$, and $x_3$ are the gas transmission parameters, including the distance between the two compressors, compression ratio value, and value of the gas pipe inside diameter. A scheme of this problem is shown in Figure 7. The problem formulation is presented as follows.

$$\mathrm{Minimize} f\left(\overrightarrow{x}\right)=3{.69 \times 10}^4{ \times x}_3+7{.72 \times 10}^8{x}_1^{-1}{x}_2^{0.219}- 765{.43 \times 10}^6{ \times x}_1^{-1}+8{.61 \times 10}^5{ \times x}_1^{\frac{1}{2}}{x}_2{\left(x_2^2- 1\right)}^{-\frac{1}{2}}{x}_3^{-\frac{2}{3}}$$

(42)

$\mathrm{Subject} \mathrm{to}:$

$$10 \le { x}_1 \le 55$$

(43)

$$1.1 \le { x}_2 \le 2$$

(44)

$$10 \le { x}_3 \le 40$$

(45)

As seen in

Table 12. Results of solving gas transmission compressor design problem.

Algorithm	x1	x2	x3	fmin
ECS-AGQPSO [64]	53.446716	1.190101	24.718579	2,964,375.495330
RCSOMGA [65]	53.446827	1.190101	24.718580	2,964,375.495330
SOMA [65]	53.347298	1.190142	24.737115	2,964,378.729000
RCGA [65]	53.520217	1.190361	24.723656	2,964,375.725000
PSO [66]	55.000000	1.195410	24.774900	2,964,460.000000
DE [66]	51.985700	1.183350	24.719500	2,964,480.000000
DE-PSO [66]	53.447400	1.190100	24.718600	2,964,375.503101
GP [67]	52.600000	1.187000	24.800000	2,964,419.625000
GWO	53.444416	1.190075	24.719897	2,964,375.499725
FB-GWO-EPD	53.440728	1.190088	24.717434	2,964,375.498014
DB-GWO-EPD	53.446720	1.190101	24.718583	2,964,375.495329

, the proposal can calculate the least value for the gas transmission cost, stressing the high capabilities of the DB-GWO-EPD to win such a highly close competition in solving this design problem.

5. Conclusions

This paper introduced a novel EPD operator for improving the performance of the Grey Wolf Optimization (GWO) algorithm. In the main version of EPD, half of the worst search agents are eliminated from the search space and repositioned around the best ones to evolve the worst individuals, while in our proposed EPD method, the best agents are eliminated and repositioned around the most diversified agents. The main reason justifying this new approach to applying EPD is the low-diversity nature of the best agents in the search space. Indeed, in any meta-heuristic, the best agents are gradually gathered around each other, such that the region the best individuals are placed in is becoming increasingly closed, and thus the diversity could be soon lost in the search space. Losing diversity could contribute to premature convergence. Hence, in the proposed Diversity-Based EPD (DB-EPD) mechanism, the enhanced diversity plays a major role in conducting the other search agents and accomplishing the exploration phase. It is worth mentioning that the three best-fitted alpha, beta, and delta wolves are always kept in the algorithm’s memory and appointed to guide the other wolves (agents), even if they are eliminated from the current population. In this paper, the proposed approach was applied to the GWO, and so the resulting algorithm was named DB-GWO-EPD.

The results of applying the proposed DB-GWO-EPD and its competitors to the high-dimensional shifted classical benchmark problems, as well as the CEC2017 challenging test bed, indicated that the proposed algorithm is highly superior to the original GWO and the GWO equipped with the EPD operator in its original definition, named Fitness-Based GWO-EPD (FB-GWO-EPD), as well as four other newly proposed optimization algorithms. The superior performance of the proposed DB-GWO-EPD algorithm is mainly due to its ability to desirably maintain diversity in the search space during the optimization process, helping the proposal to find as many candidate solutions in the search space as possible. This ability could help the proposed algorithm rapidly reach the single optimum in uni-modal functions while impeding the proposal from being trapped in local optima to explore as many optima as possible in the search space to find the fittest solution of the optimization problem finally.

In future work, we aim to apply the DB-EPD operator to other meta-heuristics to validate the proposed approach further to use such an EPD. The application of the DB-GWO-EPD in solving the more recently created and more complicated problems could be another field of interest through which the merits or shortcomings of the proposed algorithm may be better highlighted.