An Efficient and Fast Hybrid GWO-JAYA Algorithm for Design Optimization

Optimization searches for the minimum or the maximum of an NDV–variables function $W\left(\overrightarrow{X}\right)$ subject to a set of N_CON inequality/equality constraint functions $G\left(\overrightarrow{X}\right)$ or $H\left(\overrightarrow{X}\right)$. Optimization problems may be iteratively solved with gradient based algorithms (GBAs) or metaheuristic algorithms (MHAs). GBAs formulate and solve a series of approximate sub-problems until the search process converges to the optimal solution; sub-problems are built by evaluating first-order and second-order derivatives of cost functions and constraints at the solution point found in the previous iteration. MHAs do not require gradients: local information provided by gradients in GBAs are replaced in MHAs by global information gathered from a population of candidate designs. Trial solutions of MHAs are randomly generated using a mathematical model inspired by evolutionary theory/processes, biology, physics, chemistry, mathematics, astronomy, astrophysics, herd behavior of animals, human behavior and activities, social sciences, etc. MHAs evaluate the new trial solutions randomly generated and then attempt to improve the design with respect to the previous iterations.

Metaheuristic algorithms can inherently deal with highly nonlinear and non-convex optimization problems. The random search allows MHAs to explore larger portions of search space than GBAs without being trapped in local optima. This ability together with their straightforward software implementation made MHAs become the standard approach to the solution of complex engineering optimization problems. Exploration and exploitation are the two typical phases carried out in metaheuristic optimization. Exploration is the leading search mechanism in the early optimization stages where large perturbations are given to optimization variables to quickly find the best regions of the search space. Exploitation governs the search process as the optimizer converges toward the global optimum: local search is carried out in properly selected neighborhoods of the most promising solutions. Exploration and exploitation must be appropriately combined in metaheuristic optimization to fully search the design space and find the globally optimal solution with low computational effort.

MHAs can be categorized in four groups: (i) evolutionary algorithms; (ii) science-based algorithms; (iii) human-based algorithms; (iv) swarm intelligence-based algorithms. Evolutionary algorithms such as genetic algorithms (GAs) [1,2], differential evolution (DE) [3,4], evolutionary programming (EP) [5], evolution strategies (ESs) [6], and biogeography-based optimization (BBO) [7] reproduce evolution theory and evolutionary processes.

Science-based MHAs mimic laws of physics, chemistry, astronomy and astrophysics, and mathematics. Simulated annealing (SA) [8,9], charged system search (CSS) [10], magnetic charged system search (MCSS) [11], colliding bodies optimization (CBO) [12], water evaporation optimization (WAO) [13], thermal exchange optimization (TEO) [14], equilibrium optimizer (EO) [15], gases Brownian motion optimization (GBMO) [16], and Henry gas solubility optimization (HGSO) [17] are physics/chemistry-based MHAs that reproduce equilibrium conditions of mechanical, electro-magnetic, physical/chemical, or thermal systems subject to external perturbations. Ray optimization (RO) [18] and light spectrum optimizer (LSO) [19] reproduce optics laws to define search directions containing high-quality trial solutions. Big bang–big crunch optimization (BB-BC) [20] and the gravitational search algorithm (GSA) [21] are inspired by astronomy/astrophysics phenomena such as the expansion (big bang)—contraction (big crunch) cycles leading to the formation of new star–planetary systems and gravitational interactions between masses. Mathematics inspired the sine cosine algorithm (SCA) [22], the Runge–Kutta optimizer (RUN) [23], and the arithmetic optimization algorithm (AOA) [24].

Tabu search (TS) [25], harmony search optimization (HS) [26], teaching–learning-based optimization (TLBO) [27], JAYA [28], the group teaching optimization algorithm (GTOA) [29], the mother optimization algorithm (MOA) [30], the preschool education optimization algorithm (PEOA) [31], the learning cooking algorithm (LCA) [32], the imperialist competitive algorithm (ICA) [33], and the political optimizer (PO) [34] are representative algorithms reproducing human activities, behaviors, learning/education processes, social sciences, international strategies, and politics.

Swarm intelligence-based algorithms represent the largest category in metaheuristic optimization. These algorithms mimic the social/individual behavior of animals (insects, terrestrial animals, birds and other flying animals, and aquatic animals) in reproduction, food search, hunting, migration, etc. Particle swarm optimization (PSO) [35], reproducing interactions between individuals of bird/fish swarms, is the most cited metaheuristic algorithm according to the Scopus database; in PSO, a population of candidate designs (the particles) are generated, and their positions and velocities are updated referring to the position of the leader(s) and the best positions of individual particles in each iteration. Insect behavior-inspired algorithms, for example, include ant colony optimization (ACO) [36], artificial bee colony (ABC) [37], the firefly algorithm (FFA) [38], and the ant lion optimizer (ALO) [39].

The grey wolf optimizer (GWO) [40], coyote optimization algorithm (COA) [41], snake optimizer (SO) [42], and snow leopard optimization algorithm (SLOA) [43] simulate the behavior of terrestrial animals. The GWO, reproducing the hunting behavior of grey wolves, is the 2nd most cited algorithm in terms of citations/year after PSO. However, the GWO is preferred to PSO in view of its simpler formulation that does not include internal parameters except population size and a limit to the number of iterations.

The bat algorithm (BA) [44], cuckoo search (CS) [45], crow search algorithm (CSA) [46], Harris hawks optimization (HHO) [47], and starling murmuration optimizer (SMO) [48] are inspired by flying animals like bats and birds. CS (reproducing the parasitic behavior of some cuckoo species that mix their eggs with those of other birds to guarantee the survival of their chicks) and HHO (reproducing the cooperative behavior of Harris hawks in nature, specifically their surprise attacks during the chase) are the most cited MHAs of this sub-category.

Last, the dolphin echolocation algorithm (DEA) [49], whale optimization algorithm (WOA) [50], salp swarm algorithm (SSA) [51], marine predators algorithm (MPA) [52], and giant trevally optimizer (GTO) [53] reproduce social behavior, hunting strategy, swarming, and foraging of aquatic animals; WOA, SSA, and MPA are the most cited MHAs in this sub-category.

Hybrid/improved/enhanced methods also were developed in the optimization literature by adding new equations in the original formulation or merging two or more MHAs. The goal always was finding the best balance between exploration and exploitation to minimize computational cost (i.e., the number of structural analyses or function evaluations required by the optimizer), improve robustness, and limit the number of internal parameters of the metaheuristic formulation. Hybrid metaheuristic algorithms may have a parallel or a serial architecture [54]. In the former case, the component algorithms are independently run on parallel computers, while in the latter case, they are sequentially executed on the same machine.

Metaheuristic optimization algorithms are widely employed in various engineering fields such as, for example, the mechanical characterization of materials and structural identification [55], static and dynamic structural optimization [56,57], damage identification [58], vehicle routing optimization [59], 3D printing process optimization [60], cancer classification [61], and image processing [62]. However, in spite of the large diffusion of MHAs in engineering practice, some issues remain open in metaheuristic optimization: (i) no metaheuristic algorithm can always outperform all other MHAs in all optimization problems; (ii) MHAs may require a very large number of function evaluations (analyses) for completing optimization process; (iii) sophisticated algorithmic variants and hybrid MHAs combining multiple methods often increase the computational complexity of the optimizer also because of the presence of additional internal parameters that are difficult to be tuned. As a matter of fact, newly developed MHAs often added very little to the optimization field and stopped attracting potential users even after a rather short time from their release.

The main objective of this study was to develop an efficient and robust hybrid metaheuristic algorithm for engineering optimization. Generally speaking, a hybrid optimizer should combine high-performance MHAs that complement each other in terms of exploration and exploitation. Furthermore, component algorithms should be versatile so as to successfully deal with as many different optimization problems as possible. Lastly, component algorithms should be simple enough in order to simplify the formulation of the hybrid optimizer and to limit the number of additional internal parameters governing the switch from one component algorithm to another. In view of this, the grey wolf optimizer (GWO) and the JAYA algorithms were selected in the present study and combined into the novel FHGWJA hybrid algorithm (the acronym stands for Fast Hybrid Grey Wolf JAYA) because they certainly satisfy the above-mentioned requisites.

The GWO mimics the leadership hierarchy and hunting mechanism of grey wolves [40]. The leadership hierarchy is simulated by defining four types of grey wolves (α, β, δ, and ω). The hunting mechanism is simulated with different algorithmic operators that represent searching for prey, encircling prey, and attacking prey. The α wolf is the group leader. The β wolf helps the α in decision-making or other pack activities. The lowest ranking is represented by ω wolves. The δ wolves rank between the α-β wolves and the ω wolves. Alpha, beta, and delta wolves estimate the position of the prey, and other wolves randomly update their positions around the prey.

As previously mentioned, the GWO is the 2nd most cited metaheuristic algorithm after PSO according to the Scopus database. It was successfully applied to many fields because of its simple formulation and computational efficiency (see, for example, the surveys on GWO applications presented in [63,64,65]). However, the GWO may suffer from lack of exploitation, risk of stagnation especially in complex problems, and limited adaptability to variations in the problem landscape during the optimization process. Enhanced GWO formulations including (i) new operators to capture other characteristic behaviors of wolves (i.e., gaze cues learning [64] and dimension learning-based hunting [66]), (ii) Cauchy-Gaussian mutation (increasing the search range of leader wolves when they tend to the local optimal solution) along with greedy selection (to maintain population diversity) and improved search strategy considering average position of all individuals [67], (iii) chaotic grouping and dynamic regrouping of individuals to increase population diversity [68], (iv) optimization of initial population along with a nonlinear control parameter convergence strategy and nonlinear tuning strategy of parameters [69], and (v) update of wolves’ positions with spiral movements [70], as well as (vi) hybrid algorithms combining GWO with other powerful MHAs (i.e., particle swarm optimization, biogeography-based optimization, differential evolution, Harris hawks optimization, and the whale optimization algorithm [71,72,73,74,75]) were proposed in order to overcome the above-mentioned limitations.

JAYA [28] utilizes the most straightforward search scheme amongst all MHAs consisting of only one equation to perturb optimization variables: to approach the population’s best solution and move away from the worst solution, thus achieving a significant convergence capability. JAYA is a very versatile algorithm with a large number of applications documented in the literature [76,77]. JAYA has an inherently hybrid nature combining basic features of evolutionary algorithms (the survival of fittest individual) and swarm-based algorithms where the swarm normally follows the leader in the search of the optimal solution. These characteristics make JAYA a very good candidate component of new hybrid metaheuristic algorithms.

Similar to the GWO, JAYA may present a rather weak exploitation phase. Furthermore, JAYA uses only one equation to perturb optimization variables involving only the best and the worst individuals of the population: this may produce stagnation and also limit the exploration phase. To overcome these issues, improved JAYA formulations using sub-populations [78], or involving also the current average solution and the historical solutions (a population of candidate solutions initially generated besides the standard population and permuted in the optimization process according to a probabilistic criterion) [79], were developed. Fuzzy clustering competitive learning, experience learning, and Cauchy mutation mechanisms were also implemented [80] to effectively utilize population information, speed up the convergence rate, improve the balance between exploiting the previously visited regions and exploring new regions of search space in the search process, and reduce the risk of being trapped into local optimum by fine-tuning the quality of the so-far best solution. The high computational cost is another issue in JAYA optimization. In [81,82], it was attempted to reduce the number of analyses and increase convergence speed by directly rejecting heavier designs than population individuals, but this strategy missed the global optimum in many structural optimization problems. Similarly, in [83], the population was updated not only if the new trial solutions improve the individuals currently stored in the population but also if they have the same values for the cost function or penalized cost function. In [84], generation of trial designs also relied on two randomly selected individuals $\overrightarrow{X_m}$ and $\overrightarrow{X_n}$, and the perturbation equation varied if m > n or m < n as well as if $W\left(\overrightarrow{X_m}\right)$ < $W\left(\overrightarrow{X_n}\right)$ or $W\left(\overrightarrow{X_m}\right)$ > $W\left(\overrightarrow{X_n}\right)$. JAYA was also hybridized with other efficient MHAs such as, for example, harmony search optimization [56], genetic algorithms [85], crow search [86], the Rao-1 optimizer [87], teaching–learning-based optimization [88], and differential evolution [89].

The main features of the novel FHGWJA algorithm developed in this study and the advancements with respect to the state-of-the-art can be summarized as follows:

The proposed FHGWJA algorithm was successfully tested in seven “real world” engineering problems. The selected test cases, including up to 29 optimization variables and 1200 nonlinear constraints, regarded, in particular, the following: (i–ii) 2D path planning (minimization of trajectory lengths, respectively, with 7 or 10 obstacles); (iii) shape optimization of a concrete gravity dam with additional earthquake load (volume minimization); (iv–v) calibration of nonlinear hydrologic models (the Muskingum problem solved with 3 or 25 unknown model parameters to be identified); (vi) optimal crashworthiness design of a vehicle subject to side impact; (vii) weight minimization of a planar 200-bar truss structure subject to three independent loading conditions. The present algorithm was compared with the best performing algorithms indicated in the literature for each test problem and advanced variants of state-of-the-art MHAs.

The rest of this manuscript is structured as follows. Section 2 describes the new hybrid optimization algorithm FHGWJA developed in this study. Test problems and optimization results are presented and discussed in Section 3. Section 4 summarizes the main findings of this study and outlines directions of future research.

The new hybrid metaheuristic algorithm FHGWJA developed in this study is now described in detail. The algorithm is composed of seven steps: (i) initialization; (ii) preliminary definition of trial solutions with the classical GWO; (iii) rank-based refinement of preliminary trial solutions defined with the classical GWO; (iv) evaluation/repair of trial solutions with elitist strategies and JAYA-based schemes; (v) population reordering to define the new best and worst individuals and update of the α-β-δ wolves trying to avoid stagnation; (vi) convergence check; (vii) end of optimization search and output.

FHGWJA randomly generates a population of N_POP candidate designs (i.e., wolves) as: where NDV is the number of optimization variables; $x_j^L$ and $x_j^U$, respectively, are the lower and upper bounds of the jth optimization variable; ${{\rho }_j}^i$ is a random number in the (0,1) interval.

In the minimization problem of the W($\overrightarrow{X}$) function (depending on NDV variables stored in the solution vector $\overrightarrow{X}$) subjected to N_CON inequality constraint functions G_k($\overrightarrow{X}$) ≤ 0, the penalized cost function W_p($\overrightarrow{\mathrm{X}}$) is defined as follows: where p is a penalty coefficient. The penalty function $\psi$ is defined as follows:

The W_p($\overrightarrow{X}$) penalized cost function coincides with the W($\overrightarrow{X}$) cost function if the trial solution $\overrightarrow{X}$ satisfies optimization constraints. No penalty function is defined in the case of unconstrained optimization problems. Any equality constraint H($\overrightarrow{X}$) = 0 included in the optimization problem is transformed into two inequality constraints G($\overrightarrow{X}$) ≤ 0 and −G_k($\overrightarrow{X}$) ≤ 0. Candidate solutions are sorted in terms of penalized cost function: the current best solution $\stackrel{⃑}{X_{best}}$ and the worst solution $\stackrel{⃑}{X_{worst}}$, respectively, correspond to the lowest and highest values of $\psi$.

FHGWJA utilizes the classical GWO scheme to preliminarily generate the new trial designs. The best three individuals stored in the population are ranked as wolves α, β, and δ. The other individuals of the population are ranked as wolves ω. Wolves encircle the prey during the hunt. Such a behavior is mathematically described as follows:

In Equations (4) and (5), it denotes the current iteration, $\overrightarrow{A}$ and $\overrightarrow{C}$ are coefficient vectors, $\overrightarrow{X_p}$ is the position vector of the prey, and $\overrightarrow{X_{it} }$ is the generic grey wolf (i.e., search agent) of the population updated in the current iteration to the new solution $\overrightarrow{X_{it+1}}$. The “·” notation denotes the term-by-term multiplication between vectors that leads to defining another vector.

Vectors $\overrightarrow{A}$ and $\overrightarrow{C}$ are defined as follows:

In Equations (6) and (7), the components of the $\overrightarrow{a}$ vector decrease linearly from 2 to 0 as the optimization progresses; $\overrightarrow{r_1}$ and $\overrightarrow{r_2}$ are random vectors in the [0,1] interval.

While grey wolves actually recognize the prey’s location in the hunting phase, the optimum position (prey) usually is not known a priori in the optimization search. In order to reproduce wolves’ behavior in the GWO algorithm, it is assumed that the three best solutions of the population (i.e., α, β, and δ wolves) have a better knowledge of the potential prey’s location (i.e., the target optimal solution). The rest of population must be updated using the $\overrightarrow{X_{\alpha }}$, $\overrightarrow{X_{\beta }}$, and $\overrightarrow{X_{\delta }}$ positions of the best three search agents. The generic solution $\overrightarrow{X_i}$ of population is updated to $\overrightarrow{X_{i,tr}}$ as follows:

In Equation (9), vectors $\stackrel{⃑}{X_1},\stackrel{⃑}{{ X}_2},\mathrm{a}\mathrm{n}\mathrm{d}\stackrel{⃑}{X_3}$ are defined for each candidate design $\stackrel{⃑}{X_i}$ of the population. Equations (8)–(10) update the positions of all search agents.

In the optimization process, the wolves α, β, and δ estimate the prey’s position (i.e., the position of the optimal solution) and each candidate solution in the population updates its distance from the prey. The parameter a is reduced from 2 to 0 to emphasize exploration and exploitation, respectively. In the exploration phase, candidate solutions tend to search for another prey when $‖\overrightarrow{A}‖$ > 1. Conversely, in the exploitation phase, candidate solutions converge towards the prey when $‖\overrightarrow{A}‖$ < 1.

In the classical GWO, the new population formed by the trial solutions $\overrightarrow{X_{i,tr}}$ obtained by perturbing the individuals $\overrightarrow{X_i}$ stored in the population in the previous iteration are sorted in terms of penalized cost function values to update positions $\overrightarrow{X_{\alpha }}$, $\overrightarrow{X_{\beta }}$, and $\overrightarrow{X_{\delta }}$. The process ends upon reaching the limit to the number of iterations or function evaluations defined by the user. However, the classical GWO search scheme has two inherent drawbacks: (i) since the position of the prey (optimal solution) is not known a priori, the search is directly biased toward the α wolf that has a null distance from the current best record X_best as it corresponds to the current best record itself in view of population sorting; (ii) there is no guarantee that each new trial solutions $\overrightarrow{X_{i,tr}}$ may improve the current best record $\stackrel{⃑}{X_{best}}$ or at least the corresponding individual $\overrightarrow{X_i}$ stored in the previous population.

In order to solve these issues, FHGWJA refines the preliminary trial solutions ${\left(\overrightarrow{X_{i,tr}}\right)}^{prel}$ obtained in Step 1 by means of a rank-based criterion. The population defined in the previous iteration is sorted with respect to the penalized cost function values W_p$\left(\overrightarrow{X_i}\right)$. The best solution is assigned the rank 1 while the worst solution is assigned the rank 1/N_POP: let rank⁽ⁱ⁾ = 1/i be the rank of the generic individual $\overrightarrow{X_i}$. The average rate of variation for the penalized cost function with respect to the current best record γ⁽ⁱ⁾ = ΔW_p⁽ⁱ⁾/ΔS⁽ⁱ⁾ is evaluated for each solution $\overrightarrow{X_i}$. The ΔW_p⁽ⁱ⁾ = [W_p$\left(\overrightarrow{X_i}\right)$ − W_p$\left(\overrightarrow{X_{best}}\right)$] numerator represents the increase in the penalized cost function that occurs when the design is perturbed by ΔS⁽ⁱ⁾ = $‖\overrightarrow{X_i}-\overrightarrow{X_{best}}‖$ from the current best record to the generic candidate solution $\overrightarrow{X_i}$. The average rate of variation ${\mathsf{\gamma }}^{\left(\mathrm{a}\mathrm{v}\mathrm{e}\mathrm{r}\right)}=\sum _{\mathrm{i}=1}^{{\mathrm{N}}_{\mathrm{P}\mathrm{O}\mathrm{P}}}{\mathsf{\gamma }}^{\left(\mathrm{i}\right)}/{\mathrm{N}}_{\mathrm{P}\mathrm{O}\mathrm{P}}$ is hence defined. The refinement step size ${\mu }^{\left(\mathrm{i}\right)}$ is defined for each individual $\overrightarrow{X_i}$ as follows:

The preliminary trial solution ${\left(\overrightarrow{X_{i,tr}}\right)}^{prel}$ is adjusted to the trial solution $\overrightarrow{X_{i,tr}}$ using the following equation:

The rationale behind Equations (11) and (12) is the following. All trial solutions generated via the GWO must be evaluated to see if they improve design. Individuals ranked as high-quality solutions in the previous iteration are very likely to preserve their rank in the current iteration. However, rank may decrease if the cost function changes sharply in the neighborhood of the currently perturbed solution: for example, this may occur if the perturbed solution turns infeasible, and the penalty term increases. For this reason, Equation (11) tends to preserve the overall ranks of population individuals by setting small refinement step sizes ${\mu }^{\left(\mathrm{i}\right)}$. Since the design improves by moving from $\overrightarrow{X_i}$ to $\overrightarrow{X_{best}}$, the direction $\overrightarrow{S_i}=\left(\overrightarrow{X_{best}}-\overrightarrow{X_i}\right)$ is a descent direction with respect to the $\overrightarrow{X_i}$ individual. Hence, FHGWJA attempts to perturb design along the $\overrightarrow{S_i}$ direction to further reduce the cost function at least with respect to $\overrightarrow{X_i}$.

The trial designs $\overrightarrow{X_{i,tr}}$ generated via FHGWJA in Steps 1 and 2 are then evaluated. In the classical GWO, if the new design $\overrightarrow{X_{i,tr}}$ improves the old design $\overrightarrow{X_i},$ it replaces $\overrightarrow{X_i}$ in the new population. However, this task entails a new constraint evaluation to compute the penalized cost function value for each new trial design.

In order to reduce the number of constraint evaluations required in the optimization search, FHGWJA utilizes an elitist strategy retaining only the trial solutions that are very likely to improve design. Hence, FHGWJA initially compares only the cost function W($\overrightarrow{X_{i,tr}}$) computed for the new trial design $\overrightarrow{X_{i,tr}}$ with the penalized cost function W_p$\left(\overrightarrow{X_{best}}\right)$ computed for the current best record. If W($\overrightarrow{X_{i,tr}}$) > W_p$\left(\overrightarrow{X_{best}}\right)$, the new trial design $\overrightarrow{X_{i,tr}}$ is directly rejected because it certainly cannot improve the current best record as the cost function value $\overrightarrow{X_{i,tr}}$ is by itself larger than the penalized cost of $\overrightarrow{X_{best}}$ which also accounted for any constraint violation.

If W($\overrightarrow{X_{i,tr}}$) < W_p$\left(\overrightarrow{X_{best}}\right)$, the cost function value W($\overrightarrow{X_{i,tr}}$) computed for the new trial solution $\overrightarrow{X_{i,tr}}$ is compared with 1.1 times the cost function value W($\overrightarrow{X_{best}}$) computed for the current best record. Ifit holds W($\overrightarrow{X_{i,tr}}$) ≤ 1.1W($\overrightarrow{X_{best}}$), the new trial solution $\overrightarrow{X_{i,tr}}$ is provisionally included in the new population. The 1.1W($\overrightarrow{X_{best}}$) threshold was selected upon the following rationale. Exploration is characterized by a high probability of improving design: hence, the W($\overrightarrow{X_{i,tr}}$) > 1.1W($\overrightarrow{X_{best}}$) condition is not likely to occur. In exploitation, local minima should be bypassed by the optimizer to converge to the global optimum. Similar to SA, FHGWJA is allowed to provisionally accept slightly worse candidate solutions than $\overrightarrow{X_{best}}$. The threshold level of acceptance probability in state-of-the-art SA formulations such as [81,90] is 0.9 if the ratio between the cost function increment [W($\overrightarrow{X_{i,tr}}$) − W($\overrightarrow{X_{best}}$)] and the annealing temperature T is 0.1. This means that there is a 90% probability to provisionally accept a worse design than $\overrightarrow{X_{best}}$ and make it improve in the next iterations. Since T is initially set equal to the expected optimum cost (or roughly corresponds to the cost function value W$\left(\overrightarrow{X_{best}}\right)$), it is reasonable to assume that solutions up to 10% worse than $\overrightarrow{X_{best}}$ may improve $\overrightarrow{X_{best}}$ itself in the next few iterations.

Since each new trial solution $\overrightarrow{X_{i,tr}}$ provisionally included in the new population was generated by FHGWJA using the classical GWO scheme based on the positions of wolves α, β, and δ, FHGWJA adopts a JAYA-based scheme to avoid stagnation. For that purpose, when W($\overrightarrow{X_{i,tr}}$) ≤ 1.1W($\overrightarrow{X_{best}}$), a new trial solution ${\left(\overrightarrow{X_{i,tr}}\right)}^{\prime }$ is defined as follows: where $\overrightarrow{{\lambda }_1}$ and $\overrightarrow{{\lambda }_2}$ are two vectors of NDV random numbers generated in the [0,1] interval. Equation (13) is relative to exploitation as it perturbs the good design $\overrightarrow{X_{i,tr}}$ that was provisionally included in the new population being close to the current best record or likely better than it. The δ wolf is temporarily selected as the worst individual of the population, and FHGWJA is forced to locally search in a region of design space containing only high-quality solutions.

The new trial solution $\left(\overrightarrow{X_{i,tr}}\right)\prime$ is compared with $\overrightarrow{X_{i,tr}}$ to select the solution that finally updates population. If W(($\overrightarrow{X_{i,tr}}$)′) < 1.1W($\overrightarrow{X_{best}}$) and W(($\overrightarrow{X_{i,tr}}$)′) < W($\overrightarrow{X_{i,tr}}$), $\left(\overrightarrow{X_{i,tr}}\right)\prime$ is finally retained in the new population, and $\overrightarrow{X_{i,tr}}$ is rejected. If W(($\overrightarrow{X_{i,tr}}$)′) < 1.1W($\overrightarrow{X_{best}}$) but it occurs that W(($\overrightarrow{X_{i,tr}}$)′) > W($\overrightarrow{X_{i,tr}}$), $\left(\overrightarrow{X_{i,tr}}\right)\prime$ is rejected, and $\overrightarrow{X_{i,tr}}$ is finally retained in the new population. If W(($\overrightarrow{X_{i,tr}}$)′) > 1.1W($\overrightarrow{X_{best}}$), $\overrightarrow{X_{i,tr}}$ is finally retained in the new population, and $\left(\overrightarrow{X_{i,tr}}\right)\prime$ is directly rejected.

FHGWJA adopts a repair strategy if the new trial solution $\overrightarrow{X_{i,tr}}$ is such that W($\overrightarrow{X_{i,tr}}$) > 1.1W($\overrightarrow{X_{best}}$). In this case, the trial solution $\overrightarrow{X_{i,tr}}$ was sufficiently good to avoid immediate rejection (i.e., it holds W($\overrightarrow{X_{i,tr}}$) < W_p$\left(\overrightarrow{X_{best}}\right)$) but it is not good enough to be provisionally included in the new population. Another JAYA-based scheme is adopted in FHGWJA to define the new trial solution ${\left(\overrightarrow{X_{i,tr}}\right)}^{\prime }$ as follows: where $\overrightarrow{{\lambda }_1}$ and $\overrightarrow{{\lambda }_2}$ are two vectors of NDV random numbers similar to those of Equation (13). Absolute values of optimization variables are taken for each component of vectors $\overrightarrow{X_i}$ or $\overrightarrow{X_{i,tr}}$, respectively, in Equations (13) and (14).

Equation (14) is related to exploration and resembles the classical JAYA solution updating equation where each individual is perturbed trying to approach the current best solution $\overrightarrow{X_{best}}$ and moving away from the worst solution of the population $\overrightarrow{X_{worst}}$. FHGWJA involves the whole population of N_POP search agents in the generation of new trial design ($\overrightarrow{X_{i,tr}}$)′ because wolves α, β, and δ that drive search in the GWO phase failed to move other individuals (i.e., $\overrightarrow{X_i}$) to good positions of the search space located near the prey (i.e., the current best record). Since population renewal aims to improve each individual $\overrightarrow{X_i}$, FHGWJA searches on the descent direction $\left(\overrightarrow{X_{best}}-\overrightarrow{X_i}\right)$ that leads to a better design than $\overrightarrow{X_i}$ and escape from the worst design of the population $\overrightarrow{X_{worst}}$, which certainly cannot improve $\overrightarrow{X_i}$.

Similar to the process implemented for Equation (13), FHGWJA compares the cost function value W(($\overrightarrow{X_{i,tr}}$)′) corresponding to the new trial design $\left(\overrightarrow{X_{i,tr}}\right)\prime$ defined with Equation (14) with 1.1 times the cost function value of the current best record W($\overrightarrow{X_{best}}$). The modified trial design ($\overrightarrow{X_{i,tr}}$)′ is included in the new population if it holds that W(($\overrightarrow{X_{i,tr}}$)′) ≤ 1.1W($\overrightarrow{X_{best}}$).

However, if the JAYA-based repair strategy fails and it still holds that W(($\overrightarrow{X_{i,tr}}$)′) > 1.1W($\overrightarrow{X_{best}}$), FHGWJA has to deal with three solutions, $\overrightarrow{X_i}$, $\overrightarrow{X_{i,tr}}$, and ${\left(\overrightarrow{X_{i,tr}}\right)}^{\prime }$, that do not improve the current best record $\overrightarrow{X_{best}}$. Figure 1 shows the mutual positions of the cost function gradient $\overrightarrow{\nabla W_{\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}}}$ evaluated at the current best record $\overrightarrow{X_{best}}$ and the non-descent directions $\overrightarrow{S_i}=\overrightarrow{X_i}-\overrightarrow{X_{best}}$, $\overrightarrow{S_{i,tr}}=\overrightarrow{X_{i,tr}}-\overrightarrow{X_{best}}$, and ${\left(\overrightarrow{S_{i,tr}}\right)}^{\prime }=\left(\overrightarrow{X_{i,tr}}\right)\prime -\overrightarrow{X_{best}}$; P_best is the point of the search space corresponding to the current best record $\overrightarrow{X_{best}}$ while the P_i, P_i,tr, and P_i,tr′ points, respectively, correspond to trial solutions $\overrightarrow{X_i}$, $\overrightarrow{X_{i,tr}}$, and ${ \left(\overrightarrow{X_{i,tr}}\right)}^{\prime }$. It can be seen from Figure 1 that non-descent directions make positive scalar products with the gradient $\overrightarrow{\nabla W_{\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}}}$. Design may be improved with respect to $\overrightarrow{X_{best}}$ by moving on the descent directions ${\overrightarrow{S_i}}^{\left(mirr\right)}=-\left(\overrightarrow{X_i}-\overrightarrow{X_{best}}\right)$, ${\overrightarrow{S_{i,tr}}}^{\left(mirr\right)}=-\left(\overrightarrow{X_{i,tr}}-\overrightarrow{X_{best}}\right)$, and ${\left(\overrightarrow{S_{i,tr}}\right)}^{\prime \left(mirr\right)}=\left(\overrightarrow{X_{i,tr}}\right)\prime -\overrightarrow{X_{best}}$ that are, respectively, opposite to the non-descent directions $\overrightarrow{S_i}$, $\overrightarrow{S_{i,tr}}$, and ${\left(\overrightarrow{S_{i,tr}}\right)}^{\prime }$ and hence make negative scalar products with respect to $\overrightarrow{\nabla {\mathrm{W}}_{\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}}}$ (see Figure 1). Hence, FHGWJA uses a “mirroring” strategy to define the new trial solution ${\left(\overrightarrow{X_{i,tr}}\right)}^{\prime \prime \left(mirr\right)}$ as follows: where η_mirr,i, η_mirr,i-tr, and η_{mirr,(i-tr)’} are three random numbers in the [0,1] interval. Basically, FHGWJA takes the descent direction $\left({\left(\overrightarrow{X_{i,tr}}\right)}^{\prime \prime \left(mirr\right)}-\overrightarrow{X_{best}}\right)$ obtained by combining three other descent directions that may improve design. Hence, the ${\left(\overrightarrow{X_{i,tr}}\right)}^{\prime \prime \left(mirr\right)}$ trial solution (corresponding to the point P_i,tr″^(mirr) of search space shown in Figure 1) is very likely to improve design. The new solution ${\left(\overrightarrow{X_{i,tr}}\right)}^{\prime \prime \left(mirr\right)}$ is evaluated as outlined above for the other trial solutions generated by FHGWJA.

Once FHGWJA updated all candidate solutions $\overrightarrow{X_i}$, the population is resorted and the N_POP individuals are ranked with′ respect to their penalty function values. Should W_p($\overrightarrow{X_{i,tr}}$) or W_p(($\overrightarrow{X_{i,tr}}$)′ or W_p(($\overrightarrow{X_{i,tr}}$)′) be greater than W_p($\overrightarrow{X_i}$), all new trial solutions generated/refined for $\overrightarrow{X_i}$ are rejected and the old design $\overrightarrow{X_i}$ is retained in the population. However, the elitist strategies based on the 1.1W($\overrightarrow{X_{best}}$) threshold, JAYA-based schemes of Equations (13) and (14), and the “mirroring” strategy of Equation (15) greatly increase the probability of updating the whole population with better designs than those stored in the previous population.

The best three individuals of the new population are reset as wolves α, β, and δ with the corresponding design vectors $\overrightarrow{X_{\alpha }}$, $\overrightarrow{X_{\beta }}$, and $\overrightarrow{X_{\delta }}$. The best and worst solutions are set as $\overrightarrow{X_{best}}$ and $\overrightarrow{X_{worst}}$, respectively.

In order to avoid stagnation, ranking of wolves α, β, and δ is checked with another elitist criterion if FHGWJA did not update their positions in the current iteration. The criterion is again based on the concept of descent direction. The $\overrightarrow{S_{\beta }}=\left(\overrightarrow{X_{best}}-\overrightarrow{X_{\beta }}\right)$ and $\overrightarrow{S_{\delta }}=\left(\overrightarrow{X_{best}}-\overrightarrow{X_{\delta }}\right)$ obviously are descent directions with respect to positions $\overrightarrow{X_{\beta }}$ and $\overrightarrow{X_{\delta }}$ of wolves β and δ, and the design improves by moving towards the α wolf. Hence, FHGWJA perturbs $\overrightarrow{X_{\beta }}$ and $\overrightarrow{X_{\delta }}$ along the descent directions $\overrightarrow{S_{\beta }}$ and $\overrightarrow{S_{\delta }}$. The positions of wolves β and δ are “mirrored” with respect to $\overrightarrow{X_{best}}$ as follows: where η_mirr,β and η_mirr,δ are two random numbers in the interval (0,1) that limit step sizes to reduce the probability of generating infeasible positions. The best three positions amongst $\overrightarrow{X_{\alpha }}$, $\overrightarrow{X_{\beta }}$, $\overrightarrow{X_{\delta }}$, ${\left(\overrightarrow{X_{\beta }}\right)}^{mirr}$, and ${\left(\overrightarrow{X_{\delta }}\right)}^{mirr}$ are taken as wolves α, β, and δ for the next iteration, while the two worst positions are compared with the rest of the population and may replace $\overrightarrow{X_{worst}}$ and the 2nd worst design of the old population. The latter operation also covers the case that ${\left(\overrightarrow{X_{\beta }}\right)}^{mirr}$ and ${\left(\overrightarrow{X_{\delta }}\right)}^{mirr}$ did not improve any of the wolves α, β, and δ.

The mirror strategy used in FHGWJA is illustrated in Figure 2. In particular, the figure shows the following: (i) the original positions of wolves α, β, and δ, respectively, correspond to points P_opt≡P_α, P_β, and P_δ of the search space; (ii) the positions of “mirror” wolves β^mirr and δ^mirr are defined with Equation (16) and, respectively, correspond to points P_β,mirr and P_δ,mirr of the search space; (iii) the cost function gradient vector $\overrightarrow{\nabla W_{\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}}}$ is evaluated at $\overrightarrow{X_{best}}$. It can be seen that the “mirror” wolves β^mirr and δ^mirr defined by Equation (16) lie on descent directions and may even improve $\overrightarrow{X_{best}}$, the position of wolf α. In fact, it holds that $\overrightarrow{\nabla W_{\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}}}\times \left({\left(\overrightarrow{X_{\beta }}\right)}^{mirr}-\overrightarrow{X_{opt}}\right)<0$ and $\overrightarrow{\nabla W_{\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}}}\times \left({\left(\overrightarrow{X_{\delta }}\right)}^{mirr}-\overrightarrow{X_{opt}}\right)<0$, where “×” denotes the scalar product between two vectors. Interestingly, since $\left({\left(\overrightarrow{X_{\delta }}\right)}^{mirr}-\overrightarrow{X_{opt}}\right)$ is in principle a steeper descent direction than $\left({\left(\overrightarrow{X_{\beta }}\right)}^{mirr}-\overrightarrow{X_{opt}}\right)$ as W_p($\overrightarrow{X_{\delta }}$) > W_p($\overrightarrow{X_{\beta }}$) > W_p($\overrightarrow{X_{best}}$), wolf δ may have a higher probability than wolf β to replace wolf α even though it originally occupied a worse position than wolf β in the search space. This elitist approach forces FHGWJA to carry out a new exploration of the search space rather than trying to exploit solutions that have not improved the design in the last iteration. Furthermore, by eliminating the two worst designs of the population, FHGWJA increases the average quality of candidate solutions and achieves a higher probability of generating high-quality trial solutions in the next iteration.

Standard deviations on optimization variables and penalty function values of candidate solutions included in the updated population should decrease as the optimization search approaches the global optimum. Hence, FHGWJA normalizes standard deviations with respect to average design $\overrightarrow{X_{aver}}=\left({\sum }_{\mathrm{i}=1}^{N_{\mathrm{P}\mathrm{O}\mathrm{P}}}\overrightarrow{X_i}\right)/N_{\mathrm{P}\mathrm{O}\mathrm{P}}$ and average value of penalty function $W_{p,\mathrm{a}\mathrm{v}\mathrm{e}\mathrm{r}}=\left({\sum }_{\mathrm{i}=1}^{N_{\mathrm{P}\mathrm{O}\mathrm{P}}}{W}_p\left(\overrightarrow{X_i}\right)\right)/N_{\mathrm{P}\mathrm{O}\mathrm{P}}$. The termination criterion is as follows: where the convergence limit ${\mathsf{\epsilon }}_{\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{v}}$ is equal to 10⁻⁷. Since convergence occurs in the last stage of the optimization process when exploitation dominates, Equation (17) requires all search agents be actually located in the best region of the search space hosting the global optimum and contributing to exploitation. Population diversity must decrease to aggregate search agents in the neighborhood of $\overrightarrow{X_{best}}$. Hence, Equation (17) normalizes the standard deviation of the search agents’ positions to average the solution to quantify the level of population diversity. The same rationale is adopted for the penalty function values by normalizing their standard deviation with respect to the average penalty function value: hence, competitive solutions correspond to the optimum solution only when they are effectively close to it, that is when the search process is near to end. Penalty function values are considered in Equation (17) to account for the effect of some solution that may turn locally infeasible if the set of active constraints changes as design is perturbed in the neighborhood of the global optimum. Penalty function values are equal to cost function values for feasible solutions or in the case of unconstrained optimization problems.

Steps 2 through 6 are repeated until FHGWJA converges to the global optimum.

FHGWJA terminates the optimization process and writes the output data (i.e., optimum design and optimized cost function value) in the results file.

Figure 3 shows the flow chart of the novel FHGWJA algorithm developed in this study. The present algorithm actually is an advanced grey wolf optimizer, which updates population taking care that the α, β, and δ wolves effectively lead to improve the current best record in each iteration. The JAYA schemes and the elitist strategies implemented in FHGWJA enhance exploration and exploitation forcing the present algorithm to increase diversity and select high-quality trial designs without performing too many function evaluations. The definition of all descent directions utilized in FHGWJA in the optimization search does not entail any new constraint evaluation, and all decisions are made in FHGWJA by comparing only cost function values. Interestingly, FHGWJA does not need any new internal parameters with respect to the classical GWO and JAYA formulations that require users to specify population size N_POP and limit number of iterations N_itermax. This feature of FHGWJA is not very common as simple hybrid metaheuristic algorithms may adopt a set of new heuristic parameters to switch from one component optimizer to the other. However, FHGWJA implements a master–slave algorithmic architecture where the slave algorithm, JAYA, is utilized only to refine or correct the trial designs generated by the master algorithm, GWO.

The computational cost of most GWO/JAYA implementations documented in the literature was N_POP × N_itermax function evaluations (also indicated as function calls or analyses for general problems or structural analyses for structural optimization problems) including constraint evaluations. The FHGWJA’s ability of generating high-quality trial solutions throughout the optimization process significantly reduces the computational cost of the optimization. As far as population size N_POP, it is well known that increasing N_POP in metaheuristic optimization may enhance exploration capabilities but also results in a very large number of function evaluations. Indeed, it is not necessary to work with a large population if the search engine can always generate high-quality designs that improve the current best record or the currently perturbed search agents such as is accomplished with FHGWJA. It should be noted that grey wolves hunt in nature in groups comprising, at most, 10–20 individuals (the typical family pack is formed by 5–11 animals but in some cases the pack can be formed by 2–3 families). In view of this, FHGWJA optimizations carried out in this study utilized a population of 10 individuals. Sensitivity analysis confirmed the validity of this setting.

The FHGWJA algorithm proposed in this study was tested in seven engineering problems including up to 29 optimization variables and 1200 nonlinear constraints. Test cases regarded the following: (i–ii) 2D path planning (minimization of trajectory lengths with 7 or 10 obstacles); (iii) shape optimization of a concrete gravity dam with additional earthquake load (volume minimization); (iv–v) calibration of nonlinear hydrologic models (the Muskingum problem solved with 3 or 25 unknown model parameters to be identified); (vi) optimal crashworthiness design of a vehicle subject to side impact; (vii) weight minimization of a planar 200-bar truss structure under three independent loading conditions.

FHGWJA was implemented in MATLAB as a standalone optimization code. Optimization runs were carried out on a standard PC equipped with a 3.1 Mhz AMD processor and 16 MB of RAM memory. The Mersenne-Twister (MT19937) MATLAB default algorithm was utilized to generate random numbers uniformly distributed between 0 and 1; this algorithm was selected in view of its ability to generate high-quality numbers with (2¹⁹⁹³⁷–1) long period. In order to draw statistically significant conclusions on the computational efficiency of FHGWJA, 20 independent optimization runs were executed for each test problem starting from different initial populations. The population size and limit number of iterations of FHGWJA were, respectively, set equal to 10 and 5000. Initial populations included some designs with up to 1000% constraint violation to check whether FHGWJA can quickly approach the best region of search space containing the optimum solution. Remarkably, FHGWJA always completed all optimization runs within a much lower number of analyses than the theoretical computational budget of 10 × 5000 = 50,000 analyses of the classical GWO and JAYA. The actual computational cost of FHGWJA was almost entirely due (about 90%) to the number of analyses performed in the optimization search. For example, the optimization runs performed with FHGWJA for the largest test problem solved in this study (weight minimization of a planar 200-bar truss structure including 29 sizing variables) were completed on average within about 25 min of wall clock time.