Advances in Engineering Software 69 (2014) 46–61

Contents lists available at ScienceDirect

Advances in Engineering Software

journal homepage: www.elsevier.com/locate/advengsoft

Grey Wolf Optimizer

Seyedali Mirjalili a,⇑, Seyed Mohammad Mirjalili b, Andrew Lewis a
a
School of Information and Communication Technology, Griffith University, Nathan Campus, Brisbane QLD 4111, Australia
b
Department of Electrical Engineering, Faculty of Electrical and Computer Engineering, Shahid Beheshti University, G.C. 1983963113, Tehran, Iran

a r t i c l e i n f o a b s t r a c t

Article history: This work proposes a new meta-heuristic called Grey Wolf Optimizer (GWO) inspired by grey wolves
Received 27 June 2013 (Canis lupus). The GWO algorithm mimics the leadership hierarchy and hunting mechanism of grey
Received in revised form 18 October 2013 wolves in nature. Four types of grey wolves such as alpha, beta, delta, and omega are employed for sim-
Accepted 11 December 2013
ulating the leadership hierarchy. In addition, the three main steps of hunting, searching for prey, encir-
Available online 21 January 2014
cling prey, and attacking prey, are implemented. The algorithm is then benchmarked on 29 well-known
test functions, and the results are verified by a comparative study with Particle Swarm Optimization
Keywords:
(PSO), Gravitational Search Algorithm (GSA), Differential Evolution (DE), Evolutionary Programming
Optimization
Optimization techniques
(EP), and Evolution Strategy (ES). The results show that the GWO algorithm is able to provide very com-
Heuristic algorithm petitive results compared to these well-known meta-heuristics. The paper also considers solving three
Metaheuristics classical engineering design problems (tension/compression spring, welded beam, and pressure vessel
Constrained optimization designs) and presents a real application of the proposed method in the field of optical engineering. The
GWO results of the classical engineering design problems and real application prove that the proposed algo-
rithm is applicable to challenging problems with unknown search spaces.
Ó 2013 Elsevier Ltd. All rights reserved.

1. Introduction problems since they mostly assume problems as black boxes. In

other words, only the input(s) and output(s) of a system are impor-
Meta-heuristic optimization techniques have become very pop- tant for a meta-heuristic. So, all a designer needs is to know how to
ular over the last two decades. Surprisingly, some of them such as represent his/her problem for meta-heuristics.
Genetic Algorithm (GA) [1], Ant Colony Optimization (ACO) [2], Third, the majority of meta-heuristics have derivation-free
and Particle Swarm Optimization (PSO) [3] are fairly well-known mechanisms. In contrast to gradient-based optimization ap-
among not only computer scientists but also scientists from differ- proaches, meta-heuristics optimize problems stochastically. The
ent fields. In addition to the huge number of theoretical works, optimization process starts with random solution(s), and there is
such optimization techniques have been applied in various fields no need to calculate the derivative of search spaces to find the opti-
of study. There is a question here as to why meta-heuristics have mum. This makes meta-heuristics highly suitable for real problems
become remarkably common. The answer to this question can be with expensive or unknown derivative information.
summarized into four main reasons: simplicity, flexibility, deriva- Finally, meta-heuristics have superior abilities to avoid local op-
tion-free mechanism, and local optima avoidance. tima compared to conventional optimization techniques. This is
First, meta-heuristics are fairly simple. They have been mostly due to the stochastic nature of meta-heuristics which allow them
inspired by very simple concepts. The inspirations are typically re- to avoid stagnation in local solutions and search the entire search
lated to physical phenomena, animals’ behaviors, or evolutionary space extensively. The search space of real problems is usually un-
concepts. The simplicity allows computer scientists to simulate dif- known and very complex with a massive number of local optima,
ferent natural concepts, propose new meta-heuristics, hybridize so meta-heuristics are good options for optimizing these challeng-
two or more meta-heuristics, or improve the current meta-heuris- ing real problems.
tics. Moreover, the simplicity assists other scientists to learn meta- The No Free Lunch (NFL) theorem [4] is worth mentioning here.
heuristics quickly and apply them to their problems. This theorem has logically proved that there is no meta-heuristic
Second, flexibility refers to the applicability of meta-heuristics best suited for solving all optimization problems. In other words,
to different problems without any special changes in the structure a particular meta-heuristic may show very promising results on a
of the algorithm. Meta-heuristics are readily applicable to different set of problems, but the same algorithm may show poor perfor-
mance on a different set of problems. Obviously, NFL makes this
⇑ Corresponding author. Tel.: +61 434555738. field of study highly active which results in enhancing current ap-
E-mail addresses: seyedali.mirjalili@griffithuni.edu.au (S. Mirjalili), proaches and proposing new meta-heuristics every year. This also
mohammad.smm@gmail.com (S.M. Mirjalili), a.lewis@griffith.edu.au (A. Lewis).

0965-9978/$ - see front matter Ó 2013 Elsevier Ltd. All rights reserved.
http://dx.doi.org/10.1016/j.advengsoft.2013.12.007
 S. Mirjalili et al. / Advances in Engineering Software 69 (2014) 46–61 47

motivates our attempts to develop a new meta-heuristic with inspired by the concepts of evolution in nature. The most popular
inspiration from grey wolves. algorithm in this branch is GA. This algorithm was proposed by
Generally speaking, meta-heuristics can be divided into two Holland in 1992 [13] and simulates Darwnian evolution concepts.
main classes: single-solution-based and population-based. In the The engineering applications of GA were extensively investigated
former class (Simulated Annealing [5] for instance) the search pro- by Goldberg [14]. Generally speaking, the optimization is done
cess starts with one candidate solution. This single candidate solu- by evolving an initial random solution in EAs. Each new population
tion is then improved over the course of iterations. Population- is created by the combination and mutation of the individuals in
based meta-heuristics, however, perform the optimization using the previous generation. Since the best individuals have higher
a set of solutions (population). In this case the search process starts probability of participating in generating the new population, the
with a random initial population (multiple solutions), and this new population is likely to be better than the previous genera-
population is enhanced over the course of iterations. Population- tion(s). This can guarantee that the initial random population is
based meta-heuristics have some advantages compared to single optimized over the course of generations. Some of the EAs are Dif-
solution-based algorithms: ferential Evolution (DE) [15], Evolutionary Programing (EP) [16,17],
and Evolution Strategy (ES) [18,19], Genetic Programming (GP)
 Multiple candidate solutions share information about the [20], and Biogeography-Based Optimizer (BBO) [21].
search space which results in sudden jumps toward the prom- As an example, the BBO algorithm was first proposed by Simon
ising part of search space. in 2008 [21]. The basic idea of this algorithm has been inspired by
 Multiple candidate solutions assist each other to avoid locally biogeography which refers to the study of biological organisms in
optimal solutions. terms of geographical distribution (over time and space). The case
 Population-based meta-heuristics generally have greater explo- studies might include different islands, lands, or even continents
ration compared to single solution-based algorithms. over decades, centuries, or millennia. In this field of study different
ecosystems (habitats or territories) are investigated for finding the
One of the interesting branches of the population-based meta- relations between different species (habitants) in terms of immi-
heuristics is Swarm Intelligence (SI). The concepts of SI was first gration, emigration, and mutation. The evolution of ecosystems
proposed in 1993 [6]. According to Bonabeau et al. [1], SI is ‘‘The (considering different kinds of species such as predator and prey)
emergent collective intelligence of groups of simple agents’’. The inspi- over migration and mutation to reach a stable situation was the
rations of SI techniques originate mostly from natural colonies, main inspiration of the BBO algorithm.
flock, herds, and schools. Some of the most popular SI techniques The second main branch of meta-heuristics is physics-based
are ACO [2], PSO [3], and Artificial Bee Colony (ABC) [7]. A compre- techniques. Such optimization algorithms typically mimic physical
hensive literature review of the SI algorithms is provided in the rules. Some of the most popular algorithms are Gravitational Local
next section. Some of the advantages of SI algorithms are: Search (GLSA) [22], Big-Bang Big-Crunch (BBBC) [23], Gravitational
Search Algorithm (GSA) [24], Charged System Search (CSS) [25],
 SI algorithms preserve information about the search space over Central Force Optimization (CFO) [26], Artificial Chemical Reaction
the course of iteration, whereas Evolutionary Algorithms (EA) Optimization Algorithm (ACROA) [27], Black Hole (BH) [28] algo-
discard the information of the previous generations. rithm, Ray Optimization (RO) [29] algorithm, Small-World Optimi-
 SI algorithms often utilize memory to save the best solution zation Algorithm (SWOA) [30], Galaxy-based Search Algorithm
obtained so far. (GbSA) [31], and Curved Space Optimization (CSO) [32]. The mech-
 SI algorithms usually have fewer parameters to adjust. anism of these algorithms is different from EAs, in that a random
 SI algorithms have less operators compared to evolutionary set of search agents communicate and move throughout search
approaches (crossover, mutation, elitism, and so on). space according to physical rules. This movement is implemented,
 SI algorithms are easy to implement. for example, using gravitational force, ray casting, electromagnetic
force, inertia force, weights, and so on.
Regardless of the differences between the meta-heuristics, a For example, the BBBC algorithm was inspired by the big bang
common feature is the division of the search process into two and big crunch theories. The search agents of BBBC are scattered
phases: exploration and exploitation [8–12]. The exploration phase from a point in random directions in a search space according to
refers to the process of investigating the promising area(s) of the the principles of the big bang theory. They search randomly and
search space as broadly as possible. An algorithm needs to have sto- then gather in a final point (the best point obtained so far) accord-
chastic operators to randomly and globally search the search space ing to the principles of the big crunch theory. GSA is another phys-
in order to support this phase. However, exploitation refers to the lo- ics-based algorithm. The basic physical theory from which GSA is
cal search capability around the promising regions obtained in the inspired is Newton’s law of universal gravitation. The GSA algo-
exploration phase. Finding a proper balance between these two rithm performs search by employing a collection of agents that
phases is considered a challenging task due to the stochastic nature have masses proportional to the value of a fitness function. During
of meta-heuristics. This work proposes a new SI technique with iteration, the masses are attracted to each other by the gravita-
inspiration from the social hierarchy and hunting behavior of grey tional forces between them. The heavier the mass, the bigger the
wolf packs. The rest of the paper is organized as follows: attractive force. Therefore, the heaviest mass, which is possibly
Section 2 presents a literature review of SI techniques. Section 3 close to the global optimum, attracts the other masses in propor-
outlines the proposed GWO algorithm. The results and discussion tion to their distances.
of benchmark functions, semi-real problems, and a real application The third subclass of meta-heuristics is the SI methods. These
are presented in Sections 4-6, respectively. Finally, Section 7 con- algorithms mostly mimic the social behavior of swarms, herds,
cludes the work and suggests some directions for future studies. flocks, or schools of creatures in nature. The mechanism is almost
similar to physics-based algorithm, but the search agents navigate
using the simulated collective and social intelligence of creatures.
2. Literature review The most popular SI technique is PSO. The PSO algorithm was pro-
posed by Kennedy and Eberhart [3] and inspired from the social
Meta-heuristics may be classified into three main classes: behavior of birds flocking. The PSO algorithm employs multiple
evolutionary, physics-based, and SI algorithms. EAs are usually particles that chase the position of the best particle and their
48 S. Mirjalili et al. / Advances in Engineering Software 69 (2014) 46–61

own best positions obtained so far. In other words, a particle is been observed, in which an alpha follows the other wolves in the
moved considering its own best solution as well as the best solu- pack. In gatherings, the entire pack acknowledges the alpha by
tion the swarm has obtained. holding their tails down. The alpha wolf is also called the dominant
Another popular SI algorithm is ACO, proposed by Dorigo et al. wolf since his/her orders should be followed by the pack [46]. The
in 2006 [2]. This algorithm was inspired by the social behavior of alpha wolves are only allowed to mate in the pack. Interestingly,
ants in an ant colony. In fact, the social intelligence of ants in find- the alpha is not necessarily the strongest member of the pack
ing the shortest path between the nest and a source of food is the but the best in terms of managing the pack. This shows that the
main inspiration of ACO. A pheromone matrix is evolved over the organization and discipline of a pack is much more important than
course of iteration by the candidate solutions. The ABC is another its strength.
popular algorithm, mimicking the collective behavior of bees in The second level in the hierarchy of grey wolves is beta. The be-
finding food sources. There are three types of bees in ABS: scout, tas are subordinate wolves that help the alpha in decision-making
onlooker, and employed bees. The scout bees are responsible for or other pack activities. The beta wolf can be either male or female,
exploring the search space, whereas onlooker and employed bees and he/she is probably the best candidate to be the alpha in case
exploit the promising solutions found by scout bees. Finally, the one of the alpha wolves passes away or becomes very old. The beta
Bat-inspired Algorithm (BA), inspired by the echolocation behavior wolf should respect the alpha, but commands the other lower-level
of bats, has been proposed recently [33]. There are many types of wolves as well. It plays the role of an advisor to the alpha and dis-
bats in the nature. They are different in terms of size and weight, cipliner for the pack. The beta reinforces the alpha’s commands
but they all have quite similar behaviors when navigating and throughout the pack and gives feedback to the alpha.
hunting. Bats utilize natural sonar in order to do this. The two main The lowest ranking grey wolf is omega. The omega plays the
characteristics of bats when finding prey have been adopted in role of scapegoat. Omega wolves always have to submit to all the
designing the BA algorithm. Bats tend to decrease the loudness other dominant wolves. They are the last wolves that are allowed
and increase the rate of emitted ultrasonic sound when they chase to eat. It may seem the omega is not an important individual in
prey. This behavior has been mathematically modeled for the BA the pack, but it has been observed that the whole pack face internal
algorithm. The rest of the SI techniques proposed so far are as fighting and problems in case of losing the omega. This is due to
follows: the venting of violence and frustration of all wolves by the ome-
ga(s). This assists satisfying the entire pack and maintaining the
 Marriage in Honey Bees Optimization Algorithm (MBO) in 2001 dominance structure. In some cases the omega is also the babysit-
[34]. ters in the pack.
 Artificial Fish-Swarm Algorithm (AFSA) in 2003 [35]. If a wolf is not an alpha, beta, or omega, he/she is called subor-
 Termite Algorithm in 2005 [36]. dinate (or delta in some references). Delta wolves have to submit
 Wasp Swarm Algorithm in 2007 [37]. to alphas and betas, but they dominate the omega. Scouts, senti-
 Monkey Search in 2007 [38]. nels, elders, hunters, and caretakers belong to this category. Scouts
 Bee Collecting Pollen Algorithm (BCPA) in 2008 [39]. are responsible for watching the boundaries of the territory and
 Cuckoo Search (CS) in 2009 [40]. warning the pack in case of any danger. Sentinels protect and guar-
 Dolphin Partner Optimization (DPO) in 2009 [41]. antee the safety of the pack. Elders are the experienced wolves who
 Firefly Algorithm (FA) in 2010 [42]. used to be alpha or beta. Hunters help the alphas and betas when
 Bird Mating Optimizer (BMO) in 2012 [43]. hunting prey and providing food for the pack. Finally, the caretak-
 Krill Herd (KH) in 2012 [44]. ers are responsible for caring for the weak, ill, and wounded wolves
 Fruit fly Optimization Algorithm (FOA) in 2012 [45]. in the pack.
In addition to the social hierarchy of wolves, group hunting is
This list shows that there are many SI techniques proposed so another interesting social behavior of grey wolves. According to
far, many of them inspired by hunting and search behaviors. To Muro et al. [47] the main phases of grey wolf hunting are as
the best of our knowledge, however, there is no SI technique in follows:
the literature mimicking the leadership hierarchy of grey wolves,
well known for their pack hunting. This motivated our attempt  Tracking, chasing, and approaching the prey.
to mathematically model the social behavior of grey wolves, pro-  Pursuing, encircling, and harassing the prey until it stops
pose a new SI algorithm inspired by grey wolves, and investigate moving.
its abilities in solving benchmark and real problems.  Attack towards the prey.

These steps are shown in Fig. 2.

3. Grey Wolf Optimizer (GWO)
In this work this hunting technique and the social hierarchy of
grey wolves are mathematically modeled in order to design GWO
In this section the inspiration of the proposed method is first
and perform optimization.
discussed. Then, the mathematical model is provided.

3.1. Inspiration

Grey wolf (Canis lupus) belongs to Canidae family. Grey wolves

are considered as apex predators, meaning that they are at the top
of the food chain. Grey wolves mostly prefer to live in a pack. The
group size is 5–12 on average. Of particular interest is that they
have a very strict social dominant hierarchy as shown in Fig. 1.
The leaders are a male and a female, called alphas. The alpha is
mostly responsible for making decisions about hunting, sleeping
place, time to wake, and so on. The alpha’s decisions are dictated
to the pack. However, some kind of democratic behavior has also Fig. 1. Hierarchy of grey wolf (dominance decreases from top down).
 S. Mirjalili et al. / Advances in Engineering Software 69 (2014) 46–61 49

3.2. Mathematical model and algorithm points illustrated in Fig. 3. So a grey wolf can update its position in-
side the space around the prey in any random location by using
In this subsection the mathematical models of the social hierar- Eqs. (3.1) and (3.2).
chy, tracking, encircling, and attacking prey are provided. Then the The same concept can be extended to a search space with n
GWO algorithm is outlined. dimensions, and the grey wolves will move in hyper-cubes (or hy-
per-spheres) around the best solution obtained so far.
3.2.1. Social hierarchy
In order to mathematically model the social hierarchy of wolves 3.2.3. Hunting
when designing GWO, we consider the fittest solution as the alpha Grey wolves have the ability to recognize the location of prey
(a). Consequently, the second and third best solutions are named and encircle them. The hunt is usually guided by the alpha. The
beta (b) and delta (d) respectively. The rest of the candidate solu- beta and delta might also participate in hunting occasionally. How-
tions are assumed to be omega (x). In the GWO algorithm the ever, in an abstract search space we have no idea about the loca-
hunting (optimization) is guided by a, b, and d. The x wolves fol- tion of the optimum (prey). In order to mathematically simulate
low these three wolves. the hunting behavior of grey wolves, we suppose that the alpha
(best candidate solution) beta, and delta have better knowledge
3.2.2. Encircling prey about the potential location of prey. Therefore, we save the first
As mentioned above, grey wolves encircle prey during the hunt. three best solutions obtained so far and oblige the other search
In order to mathematically model encircling behavior the follow- agents (including the omegas) to update their positions according
ing equations are proposed: to the position of the best search agents. The following formulas
are proposed in this regard.
D ¼ j~
~ C ~
X p ðtÞ  ~
XðtÞj ð3:1Þ
Da ¼ j~
~ C1  ~
Xa  ~ Db ¼ j~
Xj; ~ C2  ~
Xb  ~ Dd ¼ j~
Xj; ~ C3  ~
Xd  ~
Xj ð3:5Þ
~ X p ðtÞ  ~
Xðt þ 1Þ ¼ ~ A~
D ð3:2Þ
~ Xa  ~
X1 ¼ ~ A1  ð~ X2 ¼ ~
Da Þ; ~ Xb  ~
A2  ð~ Xd  ~
X3 ¼ ~
Db Þ; ~ A3  ð~
Dd Þ ð3:6Þ
where t indicates the current iteration, ~ A and ~C are coefficient vec-
tors, X p is the position vector of the prey, and ~
~ X indicates the posi-
~
X1 þ ~
X2 þ ~
X3
tion vector of a grey wolf. ~
Xðt þ 1Þ ¼ ð3:7Þ
The vectors ~ A and ~C are calculated as follows: 3
~
A ¼ 2~ r1  ~
a ~ a ð3:3Þ Fig. 4 shows how a search agent updates its position according to
alpha, beta, and delta in a 2D search space. It can be observed that
~ the final position would be in a random place within a circle which
C ¼ 2 ~
r2 ð3:4Þ is defined by the positions of alpha, beta, and delta in the search
where components of ~ a are linearly decreased from 2 to 0 over the space. In other words alpha, beta, and delta estimate the position
course of iterations and r1, r2 are random vectors in [0, 1]. of the prey, and other wolves updates their positions randomly
To see the effects of Eqs. (3.1) and (3.2), a two-dimensional po- around the prey.
sition vector and some of the possible neighbors are illustrated in
Fig. 3(a). As can be seen in this figure, a grey wolf in the position of 3.2.4. Attacking prey (exploitation)
(X, Y) can update its position according to the position of the prey As mentioned above the grey wolves finish the hunt by attack-
(X, Y). Different places around the best agent can be reached with ing the prey when it stops moving. In order to mathematically
respect to the current position by adjusting the value of ~ A and ~C model approaching the prey we decrease the value of ~ a. Note that
vectors. For instance, (X–X, Y) can be reached by setting the fluctuation range of ~ a. In other words ~
A is also decreased by ~ A is
~
A ¼ ð1; 0Þ and ~ C ¼ ð1; 1Þ. The possible updated positions of a grey a random value in the interval [2a, 2a] where a is decreased from
wolf in 3D space are depicted in Fig. 3(b). Note that the random 2 to 0 over the course of iterations. When random values of ~ A are in
vectors r1 and r2 allow wolves to reach any position between the [1, 1], the next position of a search agent can be in any position

Fig. 2. Hunting behavior of grey wolves: (A) chasing, approaching, and tracking prey (B–D) pursuiting, harassing, and encircling (E) stationary situation and attack [47].
50 S. Mirjalili et al. / Advances in Engineering Software 69 (2014) 46–61

X*-X (X*-X,Y,Z) (X*,Y,Z) (X,Y,Z)

(X*-X,Y) (X*,Y) (X,Y) (X*-X,Y,Z*) (X*,Y,Z*) (X,Y,Z*)

(X*-X,Y,Z*-Z) (X*,Y,Z*-Z) (X,Y,Z*-Z) (X,Y*,Z)

Y*-Y
(X*,Y*,Z*)
(X*,Y*) (X,Y*,Z*)
(X*-X,Y*)
(X*-X,Y*,Z*-Z) (X*,Y*,Z*-Z) (X,Y*,Z*-Z) (X,Y*-Y,Z)
(X,Y*)

(X,Y,Z*)

(X*,Y*-Y) (X*-X,Y*-Y,Z-Z*) (X*,Y*-Y,Z*-Z) (X,Y*-Y,Z*-Z)

(X*-X,Y*-Y) (X,Y*-Y)

(a) (b)
Fig. 3. 2D and 3D position vectors and their possible next locations.

a1
C1

a2
C2
R

Dalpha

Dbeta

Move

Ddelta

a3
C3 or any other hunters

Estimated position of the

prey

Fig. 4. Position updading in GWO.

between its current position and the position of the prey. Fig. 5(a)
1
shows that |A| < 1 forces the wolves to attack towards the prey. If
|A
|>

With the operators proposed so far, the GWO algorithm allows 1

|<
|A
its search agents to update their position based on the location of If

the alpha, beta, and delta; and attack towards the prey. However,
the GWO algorithm is prone to stagnation in local solutions with
these operators. It is true that the encircling mechanism proposed
shows exploration to some extent, but GWO needs more operators
to emphasize exploration.

3.2.5. Search for prey (exploration)

(a) (b)
Grey wolves mostly search according to the position of the al-
Fig. 5. Attacking prey versus searching for prey.
pha, beta, and delta. They diverge from each other to search for
prey and converge to attack prey. In order to mathematically mod-
el divergence, we utilize ~ A with random values greater than 1 or Another component of GWO that favors exploration is ~ C. As may
less than -1 to oblige the search agent to diverge from the prey. be seen in Eq. (3.4), the ~
C vector contains random values in [0, 2].
This emphasizes exploration and allows the GWO algorithm to This component provides random weights for prey in order to sto-
search globally. Fig. 5(b) also shows that |A| > 1 forces the grey chastically emphasize (C > 1) or deemphasize (C < 1) the effect of
wolves to diverge from the prey to hopefully find a fitter prey. prey in defining the distance in Eq. (3.1). This assists GWO to show
 S. Mirjalili et al. / Advances in Engineering Software 69 (2014) 46–61 51

a more random behavior throughout optimization, favoring explo- However, we have kept the GWO algorithm as simple as possible
ration and local optima avoidance. It is worth mentioning here that with the fewest operators to be adjusted. Such mechanisms are
C is not linearly decreased in contrast to A. We deliberately require recommended for future work. The source codes of this algorithm
C to provide random values at all times in order to emphasize can be found in http://www.alimirjalili.com/GWO.html and http://
exploration not only during initial iterations but also final itera- www.mathworks.com.au/matlabcentral/fileexchange/44974.
tions. This component is very helpful in case of local optima stag-
nation, especially in the final iterations. 4. Results and discussion
The C vector can be also considered as the effect of obstacles to
approaching prey in nature. Generally speaking, the obstacles in In this section the GWO algorithm is benchmarked on 29 bench-
nature appear in the hunting paths of wolves and in fact prevent mark functions. The first 23 benchmark functions are the classical
them from quickly and conveniently approaching prey. This is ex- functions utilized by many researchers [16,48–51,82]. Despite the
actly what the vector C does. Depending on the position of a wolf, it simplicity, we have chosen these test functions to be able to compare
can randomly give the prey a weight and make it harder and far- our results to those of the current meta-heuristics. These benchmark
ther to reach for wolves, or vice versa. functions are listed in Tables 1–3 where Dim indicates dimension of
To sum up, the search process starts with creating a random the function, Range is the boundary of the function’s search space,
population of grey wolves (candidate solutions) in the GWO algo- and fmin is the optimum. The other test beds that we have chosen
rithm. Over the course of iterations, alpha, beta, and delta wolves are six composite benchmark functions from a CEC 2005 special ses-
estimate the probable position of the prey. Each candidate solution sion [52]. These benchmark functions are the shifted, rotated, ex-
updates its distance from the prey. The parameter a is decreased panded, and combined variants of the classical functions which
from 2 to 0 in order to emphasize exploration and exploitation, offer the greatest complexity among the current benchmark func-
respectively. Candidate solutions tend to diverge from the prey tions [53]. Tables 4 lists the CEC 2005 test functions, where Dim indi-
when j~ Aj > 1 and converge towards the prey when j~ Aj < 1. Finally, cates dimension of the function, Range is the boundary of the
the GWO algorithm is terminated by the satisfaction of an end function’s search space, and fmin is the optimum. Figs. 7–10 illustrate
criterion. the 2D versions of the benchmark functions used.
The pseudo code of the GWO algorithm is presented in Fig. 6. Generally speaking, the benchmark functions used are minimi-
To see how GWO is theoretically able to solve optimization zation functions and can be divided into four groups: unimodal,
problems, some points may be noted: multimodal, fixed-dimension multimodal, and composite func-
tions. Note that a detailed descriptions of the composite bench-
 The proposed social hierarchy assists GWO to save the best mark functions are available in the CEC 2005 technical report [52].
solutions obtained so far over the course of iteration. The GWO algorithm was run 30 times on each benchmark func-
 The proposed encircling mechanism defines a circle-shaped tion. The statistical results (average and standard deviation) are re-
neighborhood around the solutions which can be extended to ported in Tables 5–8. For verifying the results, the GWO algorithm
higher dimensions as a hyper-sphere. is compared to PSO [3] as an SI-based technique and GSA [24] as a
 The random parameters A and C assist candidate solutions to physics-based algorithm. In addition, the GWO algorithm is com-
have hyper-spheres with different random radii. pared with three EAs: DE [15], Fast Evolutionary Programing
 The proposed hunting method allows candidate solutions to (FEP) [16], and Evolution Strategy with Covariance Matrix Adapta-
locate the probable position of the prey. tion (CMA-ES) [18].
 Exploration and exploitation are guaranteed by the adaptive
values of a and A. 4.1. Exploitation analysis
 The adaptive values of parameters a and A allow GWO to
smoothly transition between exploration and exploitation. According to the results of Table 5, GWO is able to provide very
 With decreasing A, half of the iterations are devoted to explora- competitive results. This algorithm outperforms all others in F1, F2,
tion (|A| P 1) and the other half are dedicated to exploitation and F7. It may be noted that the unimodal functions are suitable for
(|A| < 1). benchmarking exploitation. Therefore, these results show the
 The GWO has only two main parameters to be adjusted (a and superior performance of GWO in terms of exploiting the optimum.
C). This is due to the proposed exploitation operators previously
discussed.
There are possibilities to integrate mutation and other evolu-
tionary operators to mimic the whole life cycle of grey wolves. 4.2. Exploration analysis

In contrast to the unimodal functions, multimodal functions

have many local optima with the number increasing exponentially
with dimension. This makes them suitable for benchmarking the

Table 1
Unimodal benchmark functions.

Function Dim Range fmin

P
f1 ðxÞ ¼ ni¼1 x2i 30 [100, 100] 0
P Q
f2 ðxÞ ¼ ni¼1 jxi j þ ni¼1 jxi j 30 [10, 10] 0
Pn Pi 2 30 [100, 100] 0
f3 ðxÞ ¼ i¼1 ð j1 xj Þ
f4 ðxÞ ¼ maxi fjxi j; 1 6 i 6 ng 30 [100, 100] 0
P 2 2 30 [30, 30] 0
f5 ðxÞ ¼ n1 2
i¼1 ½100ðxiþ1  xi Þ þ ðxi  1Þ 
Pn
f6 ðxÞ ¼ i¼1 ð½xi þ 0:5Þ 2 30 [100, 100] 0
P 30 [1.28, 1.28] 0
f7 ðxÞ ¼ ni¼1 ix4i þ random½0; 1Þ
Fig. 6. Pseudo code of the GWO algorithm.
52 S. Mirjalili et al. / Advances in Engineering Software 69 (2014) 46–61

Table 2
Multimodal benchmark functions.

Function Dim Range fmin

P pffiffiffiffiffiffiffi
F 8 ðxÞ ¼ ni¼1  xi sinð jxi jÞ 30 [500, 500] 418.9829  5
Pn
F 9 ðxÞ ¼ i¼1 ½x2i  10 cosð2pxi Þ þ 10 30 [5.12, 5.12] 0
 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
P ffi  P 
F 10 ðxÞ ¼ 20 exp 0:2 1n ni¼1 x2i  exp 1n ni¼1 cosð2pxi Þ þ 20 þ e 30 [32, 32] 0
Pn 2 Q n  
1
F 11 ðxÞ ¼ 4000 xiffi 30 [600, 600] 0
i¼1 xi  i¼1 cos þ1
p
i
Pn1 P 30 [50, 50] 0
F 12 ðxÞ ¼ n f10 sinðpy1 Þ þ i¼1 ðyi  1Þ2 ½1 þ 10 sin2 ðpyiþ1 Þ þ ðyn  1Þ2 g þ ni¼1 uðxi ; 10; 100; 4Þ
p
xi þ1
yi ¼ 1 þ 48
< kðxi  aÞm xi > a
uðxi ; a; k; mÞ ¼ 0 a < xi < a
: m
kðxi  aÞ xi < a
P P 30 [50, 50] 0
F 13 ðxÞ ¼ 0:1fsin2 ð3px1 Þ þ ni¼1 ðxi  1Þ2 ½1 þ sin2 ð3pxi þ 1Þ þ ðxn  1Þ2 ½1 þ sin2 ð2pxn Þg þ ni¼1 uðxi ; 5; 100; 4Þ
P   2
 2m 30 [0, p] 4.687
i:x
F 14 ðxÞ ¼  ni¼1 sinðxi Þ  sin pi ; m ¼ 10
h Pn Pn 2 i Q
2m
F 15 ðxÞ ¼ e i¼1 ðxi =bÞ  2e i¼1 xi  ni¼1 cos2 xi ; m ¼ 5 30 [20, 20] 1
P 2 P P 2
pffiffiffiffiffiffiffi
F 16 ðxÞ ¼ f½ ni¼1 sin ðxi Þ  expð ni¼1 x2i Þg  exp½ ni¼1 sin jxi j 30 [10, 10] 1

Table 3
Fixed-dimension multimodal benchmark functions.

Function Dim Range fmin

 1 2 [65, 65] 1
1
P25 1
F 14 ðxÞ ¼ 500 þ j¼1 P2
jþ ðxi aij Þ6
i¼1

P11 x1 ðb2i þbi x2 Þ

2 4 [5, 5] 0.00030
F 15 ðxÞ ¼ i¼1 ai  2
bi þbi x3 þx4

F 16 ðxÞ ¼ 4x21  2:1x41 þ 13 x61 þ x1 x2  4x22 þ 4x42 2 [5, 5] 1.0316

 2   2 [5, 5] 0.398
F 17 ðxÞ ¼ x2  45:1 2 5
p2 x1 þ p x1  6 þ 10 1  81p cos x1 þ 10
F 18 ðxÞ ¼ ½1 þ ðx1 þ x2 þ 1Þ2 ð19  14x1 þ 3x21  14x2 þ 6x1 x2 þ 3x22 Þ  ½30 þ ð2x1  3x2 Þ2  ð18  32x1 þ 12x21 þ 48x2  36x1 x2 þ 27x22 Þ 2 [2, 2] 3
P P
F 19 ðxÞ ¼  4i¼1 ci expð 3j¼1 aij ðxj  pij Þ2 Þ 3 [1, 3] 3.86
P4 P6
F 20 ðxÞ ¼  i¼1 ci expð j¼1 aij ðxj  pij Þ2 Þ 6 [0, 1] 3.32
P 1 4 [0, 10] 10.1532
F 21 ðxÞ ¼  5i¼1 ½ðX  ai ÞðX  ai ÞT þ ci 
P7 T 1 4 [0, 10] 10.4028
F 22 ðxÞ ¼  i¼1 ½ðX  ai ÞðX  ai Þ þ ci 
P10 T 1 4 [0, 10] 10.5363
F 23 ðxÞ ¼  i¼1 ½ðX  ai ÞðX  ai Þ þ ci 

exploration ability of an algorithm. According to the results of 4.4. Convergence behavior analysis
Tables 6 and 7, GWO is able to provide very competitive results
on the multimodal benchmark functions as well. This algorithm In this subsection the convergence behavior of GWO is investi-
outperforms PSO and GSA on the majority of the multimodal func- gated. According to Berg et al. [54], there should be abrupt changes
tions. Moreover, GWO shows very competitive results compare to in the movement of search agents over the initial steps of optimi-
DE and FEP; and outperforms them occasionally. These results zation. This assists a meta-heuristic to explore the search space
show that the GWO algorithm has merit in terms of exploration. extensively. Then, these changes should be reduced to emphasize
exploitation at the end of optimization. In order to observe the con-
vergence behavior of the GWO algorithm, the search history and
4.3. Local minima avoidance trajectory of the first search agent in its first dimension are illus-
trated in Fig. 11. The animated versions of this figure can be found
The fourth class of benchmark functions employed includes in Supplementary Materials. Note that the benchmark functions
composite functions, generally very challenging test beds for are shifted in this section, and we used six search agents to find
meta-heuristic algorithms. So, exploration and exploitation can the optima.
be simultaneously benchmarked by the composite functions. The second column of Fig. 11 depicts the search history of the
Moreover, the local optima avoidance of an algorithm can be search agents. It may be observed that the search agents of GWO
examined due to the massive number of local optima in such test tend to extensively search promising regions of the search spaces
functions. According to Table 8, GWO outperforms all others on and exploit the best one. In addition, the fourth column of Fig. 11
half of the composite benchmark functions. The GWO algorithm shows the trajectory of the first particle, in which changes of the
also provides very competitive results on the remaining composite first search agent in its first dimension can be observed. It can be
benchmark functions. This demonstrates that GWO shows a good seen that there are abrupt changes in the initial steps of iterations
balance between exploration and exploitation that results in high which are decreased gradually over the course of iterations.
local optima avoidance. This superior capability is due to the adap- According to Berg et al. [54], this behavior can guarantee that a
tive value of A. As mentioned above, half of the iterations are de- SI algorithm eventually convergences to a point in search space.
voted to exploration (|A| P 1) and the rest to exploitation To sum up, the results verify the performance of the GWO algo-
(|A| < 1). This mechanism assists GWO to provide very good explo- rithm in solving various benchmark functions compared to well-
ration, local minima avoidance, and exploitation simultaneously. known meta-heuristics. To further investigate the performance of
 S. Mirjalili et al. / Advances in Engineering Software 69 (2014) 46–61 53

Table 4
Composite benchmark functions.

Function Dim Range fmin

F24(CF1):
f1, f2, f3, . . ., f10 = Sphere Function 10 [5, 5] 0
½,1 ; ,2 ; ,3 ; . . . ; ,10  ¼ ½1; 1; 1; . . . ; 1
[k1, k2, k3 . . ., k10] = [5/100, 5/100, 5/100, . . ., 5/100]
F25(CF2):
f1, f2, f3, . . ., f10 = Griewank’s Function 10 [5, 5] 0
½,1 ; ,2 ; ,3 ; . . . ; ,10  ¼ ½1; 1; 1; . . . ; 1
[k1, k2, k3, . . ., k10] = [5/100, 5/100, 5/100, . . ., 5/100]
F26(CF3):
f1, f2, f3, . . ., f10 = Griewank’s Function 10 [5, 5] 0
½,1 ; ,2 ; ,3 ; . . . ; ,10  ¼ ½1; 1; 1; . . . ; 1
[k1, k2, k3, . . ., k10] = [1, 1, 1, . . ., 1]
F27(CF4):
f1, f2 = Ackley’s Function 10 [5, 5] 0
f3, f4 = Rastrigin’s Function
f5, f6 = Weierstras’s Function
f7, f8 = Griewank’s Function
f9, f10 = Sphere Function
½,1 ; ,2 ; ,3 ; . . . ; ,10  ¼ ½1; 1; 1; . . . ; 1
[k1, k2, k3, . . ., k10] = [5/32, 5/32, 1, 1, 5/0.5, 5/0.5, 5/100, 5/100, 5/100, 5/100]
F28(CF5):
f1, f2 = Rastrigin’s Function 10 [5, 5] 0
f3, f4 = Weierstras’s Function
f5, f6 = Griewank’s Function
f7, f8 = Ackley’s Function
f9, f10 = Sphere Function
½,1 ; ,2 ; ,3 ; . . . ; ,10  ¼ ½1; 1; 1; . . . ; 1
[k1, k2, k3, . . ., k10] = [1/5, 1/5, 5/0.5, 5/0.5, 5/100, 5/100, 5/32, 5/32, 5/100, 5/100]
f29(CF6):
f1, f2 = Rastrigin’s Function 10 [5, 5] 0
f3, f4 = Weierstras’s Function
f5, f6 = Griewank’s Function
f7, f8 = Ackley’s Function
f9, f10 = Sphere Function
½,1 ; ,2 ; ,3 ; . . . ; ,10  ¼ ½0:1; 0:2; 0:3; 0:4; 0:5; 0:6; 0:7; 0:8; 0:9; 1
[k1, k2, k3, . . ., k10] = [0.1  1/5, 0.2  1/5, 0.3  5/0.5, 0.4  5/0.5, 0.5  5/100, 0.6  5/100, 0.7  5/32, 0.8  5/32, 0.9  5/100, 1  5/100]

(F1) (F2) (F3) (F4)

(F5) (F6) (F7)

Fig. 7. 2-D versions of unimodal benchmark functions.

the proposed algorithm, three classical engineering design prob- designs, are employed. These problems have several equality and
lems and a real problem in optical engineering are employed in inequality constraints, so the GWO should be equipped with a con-
the following sections. The GWO algorithm is also compared with straint handling method to be able to optimize constrained prob-
well-known techniques to confirm its results. lems as well. Generally speaking, constraint handling becomes
very challenging when the fitness function directly affects the posi-
tion updating of the search agents (GSA for instance). For the fitness
5. GWO for classical engineering problems independent algorithms, however, any kind of constraint handling
can be employed without the need to modify the mechanism of
In this section three constrained engineering design problems: the algorithm (GA and PSO for instance). Since the search agents of
tension/compression spring, welded beam, and pressure vessel the proposed GWO algorithm update their positions with respect
54 S. Mirjalili et al. / Advances in Engineering Software 69 (2014) 46–61

(F8) (F9) (F10) (F11)

(F12) (F13)

Fig. 8. 2-D versions of multimodal benchmark functions.

(F14) (F16) (F17) (F18)

Fig. 9. 2-D version of fixed-dimension multimodal benchmark functions.

(F24) (F25) (F26)

(F27) (F28) (F29)

Fig. 10. 2-D versions of composite benchmark functions.

Table 5
Results of unimodal benchmark functions.

F GWO PSO GSA DE FEP

Ave Std Ave Std Ave Std Ave Std Ave Std
F1 6.59E28 6.34E05 0.000136 0.000202 2.53E16 9.67E17 8.2E14 5.9E14 0.00057 0.00013
F2 7.18E17 0.029014 0.042144 0.045421 0.055655 0.194074 1.5E09 9.9E10 0.0081 0.00077
F3 3.29E06 79.14958 70.12562 22.11924 896.5347 318.9559 6.8E11 7.4E11 0.016 0.014
F4 5.61E07 1.315088 1.086481 0.317039 7.35487 1.741452 0 0 0.3 0.5
F5 26.81258 69.90499 96.71832 60.11559 67.54309 62.22534 0 0 5.06 5.87
F6 0.816579 0.000126 0.000102 8.28E05 2.5E16 1.74E16 0 0 0 0
F7 0.002213 0.100286 0.122854 0.044957 0.089441 0.04339 0.00463 0.0012 0.1415 0.3522

to the alpha, beta, and delta locations, there is no direct relation be- are assigned big objective function values if they violate any of the
tween the search agents and the fitness function. So the simplest constraints, can be employed effectively to handle constraints in
constraint handling method, penalty functions, where search agents GWO. In this case, if the alpha, beta, or delta violate constraints, they
 S. Mirjalili et al. / Advances in Engineering Software 69 (2014) 46–61 55

Table 6
Results of multimodal benchmark functions.

F GWO PSO GSA DE FEP

Ave Std Ave Std Ave Std Ave Std Ave Std
F8 6123.1 4087.44 4841.29 1152.814 2821.07 493.0375 11080.1 574.7 12554.5 52.6
F9 0.310521 47.35612 46.70423 11.62938 25.96841 7.470068 69.2 38.8 0.046 0.012
F10 1.06E13 0.077835 0.276015 0.50901 0.062087 0.23628 9.7E08 4.2E08 0.018 0.0021
F11 0.004485 0.006659 0.009215 0.007724 27.70154 5.040343 0 0 0.016 0.022
F12 0.053438 0.020734 0.006917 0.026301 1.799617 0.95114 7.9E15 8E15 9.2E06 3.6E06
F13 0.654464 0.004474 0.006675 0.008907 8.899084 7.126241 5.1E14 4.8E14 0.00016 0.000073

Table 7
Results of fixed-dimension multimodal benchmark functions.

F GWO PSO GSA DE FEP

Ave Std Ave Std Ave Std Ave Std Ave Std
F14 4.042493 4.252799 3.627168 2.560828 5.859838 3.831299 0.998004 3.3E16 1.22 0.56
F15 0.000337 0.000625 0.000577 0.000222 0.003673 0.001647 4.5E14 0.00033 0.0005 0.00032
F16 1.03163 1.03163 1.03163 6.25E16 1.03163 4.88E16 1.03163 3.1E13 1.03 4.9E07
F17 0.397889 0.397887 0.397887 0 0.397887 0 0.397887 9.9E09 0.398 1.5E07
F18 3.000028 3 3 1.33E15 3 4.17E15 3 2E15 3.02 0.11
F19 3.86263 3.86278 3.86278 2.58E15 3.86278 2.29E15 N/A N/A 3.86 0.000014
F20 3.28654 3.25056 3.26634 0.060516 3.31778 0.023081 N/A N/A 3.27 0.059
F21 10.1514 9.14015 6.8651 3.019644 5.95512 3.737079 10.1532 0.0000025 5.52 1.59
F22 10.4015 8.58441 8.45653 3.087094 9.68447 2.014088 10.4029 3.9E07 5.53 2.12
F23 10.5343 8.55899 9.95291 1.782786 10.5364 2.6E15 10.5364 1.9E07 6.57 3.14

are automatically replaced with a new search agent in the next iter- mathematical approaches that have been adopted to solve this
ation. Any kind of penalty function can readily be employed in order problem are the numerical optimization technique (constraints
to penalize search agents based on their level of violation. In this correction at constant cost) [55] and mathematical optimization
case, if the penalty makes the alpha, beta, or delta less fit than any technique [56]. The comparison of results of these techniques and
other wolves, it is automatically replaced with a new search agent GWO are provided in Table 9. Note that we use a similar penalty
in the next iteration. We used simple, scalar penalty functions for function for GWO to perform a fair comparison [63]. Table 9
the rest of problems except the tension/compression spring design suggests that GWO finds a design with the minimum weight for this
problem which uses a more complex penalty function. problem.

5.2. Welded beam design

5.1. Tension/compression spring design
The objective of this problem is to minimize the fabrication cost
The objective of this problem is to minimize the weight of a ten-
of a welded beam as shown in Fig. 13 [60]. The constraints are as
sion/compression spring as illustrated in Fig. 12 [55–57]. The min-
follows:
imization process is subject to some constraints such as shear
stress, surge frequency, and minimum deflection. There are three
 Shear stress (s).
variables in this problem: wire diameter (d), mean coil diameter
 Bending stress in the beam (h).
(D), and the number of active coils (N). The mathematical formula-
 Buckling load on the bar (Pc).
tion of this problem is as follows:
 End deflection of the beam (d).
Consider ~
x ¼ ½x1 x2 x3  ¼ ½dDN;  Side constraints.
Minimize xÞ ¼ ðx3 þ 2Þx2 x21 ;
f ð~
x3 x
This problem has four variables such as thickness of weld (h),
2 3
Subject to g 1 ð~
xÞ ¼ 1  71785x 4 6 0; length of attached part of bar (l), the height of the bar (t), and thick-
1

4x2 x1 x2 ness of the bar (b). The mathematical formulation is as follows:
1
g2 ð~
xÞ ¼ 12566ðx
2
3 4 þ 5108x2 6 0;
2 x1 x1 Þ 1 Consider ~
x ¼ ½x1 x2 x3 x4  ¼ ½hltb;
4x22 x1 x2 1 Minimize xÞ ¼ 1:10471x21 x2 þ 0:04811x3 x4 ð14:0 þ x2 Þ;
ðf ~
g2 ð~
xÞ ¼ 12566ðx2 x31 x41 Þ
þ 5108x21
6 0; ð5:1Þ
Subject to xÞ ¼ sð~
g 1 ð~ xÞ  smax 6 0;
xÞ ¼ 1  140:45x
g3 ð~ x2 x
1
6 0; xÞ ¼ rð~
g 2 ð~ xÞ  rmax 6 0;
2 3
g 3 ð~
xÞ ¼ dð~
xÞ  dmax 6 0; ð5:2Þ
xÞ ¼ x11:5
g4 ð~ þx2
 1 6 0;
g 4 ð~
xÞ ¼ x1  x4 6 0;
Variable range 0:05 6 x1 6 2:00;
g 5 ð~
xÞ ¼ P  P c ð~
xÞ 6 0;
0:25 6 x2 6 1:30; g 6 ð~
xÞ ¼ 0:125  x1 6 0
2:00 6 x3 6 15:0 xÞ ¼ 1:10471x21 þ 0:04811x3 x4 ð14:0 þ x2 Þ  5:0 6 0
g 7 ð~

This problem has been tackled by both mathematical and heuristic Variable range 0:1 6 x1 6 2;
approaches. Ha and Wang tried to solve this problem using PSO 0:1 6 x2 6 10;
[58]. The Evolution Strategy (ES) [59], GA [60], Harmony Search
0:1 6 x3 6 10;
(HS) [61], and Differential Evolution (DE) [62] algorithms have
also been employed as heuristic optimizers for this problem. The 0:1 6 x4 6 2
56 S. Mirjalili et al. / Advances in Engineering Software 69 (2014) 46–61

Table 8
Results of composite benchmark functions.

F GWO PSO GSA DE CMA-ES

Ave Std Ave Std Ave Std Ave Std Ave Std
F24 43.83544 69.86146 100 81.65 6.63E17 2.78E17 6.75E02 1.11E01 100 188.56
F25 91.80086 95.5518 155.91 13.176 200.6202 67.72087 28.759 8.6277 161.99 151
F26 61.43776 68.68816 172.03 32.769 180 91.89366 144.41 19.401 214.06 74.181
F27 123.1235 163.9937 314.3 20.066 170 82.32726 324.86 14.784 616.4 671.92
F28 102.1429 81.25536 83.45 101.11 200 47.14045 10.789 2.604 358.3 168.26
F29 43.14261 84.48573 861.42 125.81 142.0906 88.87141 490.94 39.461 900.26 8.32E02

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
where sð~xÞ ¼ ðs0 Þ2 þ 2s0 s00 2R
x2 2
þ ðs00 Þ ; in Fig. 14. Both ends of the vessel are capped, and the head has a
s ¼ pffiffi2x1 x2 ; s ¼ J ; M ¼ PðL þ x22 Þ;
0 P 00 MR hemi-spherical shape. There are four variables in this problem:
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
x22  2
R¼ þ x1 þx 3
;
4
npffiffiffi
2
h2  2 io
 Thickness of the shell (Ts).
x2
J ¼ 2 2x1 x2 4 þ x1 þx 2
3
;  Thickness of the head (Th).
rð~xÞ ¼ x6PL ~
2 ; dðxÞ ¼
6PL
Ex2 x4
3
 Inner radius (R).
4x 3 3
qffiffiffiffiffiffi  Length of the cylindrical section without considering the head
4:013E
x2 x6
3 4  qffiffiffiffi
Pc ð~
xÞ ¼ L2
36
1  x2L3 E
4G
; (L).
P ¼ 6000 lb; L ¼ 14 in:; dmax ¼ 0:25 in: ; E ¼ 30  16 psi;G ¼ 12  106 psi;
smax ¼ 13600 psi; rmax ¼ 30000 psi This problem is subject to four constraints. These constraints
and the problem are formulated as follows:
Coello [64] and Deb [65,66] employed GA, whereas Lee and Geem Consider ~ x ¼ ½x1 x2 x3 x4  ¼ ½T s T h RL;
[67] used HS to solve this problem. Richardson’s random method, Minimize f ð~ xÞ ¼ 0:6224x1 x3 x4 þ 1:7781x2 x23 þ 3:1661x21 x4 þ 19:84x21 x3 ;
Simplex method, Davidon-Fletcher-Powell, Griffith and Stewart’s Subject to g 1 ð~
xÞ ¼ x1 þ 0:0193x3 6 0;
successive linear approximation are the mathematical approaches g2 ð~xÞ ¼ x3 þ 0:00954x3 6 0;
that have been adopted by Ragsdell and Philips [68] for this prob- g3 ð~xÞ ¼ px23 x4  43 px33 þ 1296000 6 0;
lem. The comparison results are provided in Table 10. The results ~
g4 ðxÞ ¼ x4  240 6 0;
show that GWO finds a design with the minimum cost compared
ð5:3Þ
to others.
Variable range 0 6 x1 6 99;
5.3. Pressure vessel design 0 6 x2 6 99;
10 6 x3 6 200;
The objective of this problem is to minimize the total cost con- 10 6 x4 6 200;
sisting of material, forming, and welding of a cylindrical vessel as

Fig. 11. Search history and trajectory of the first particle in the first dimension.
 S. Mirjalili et al. / Advances in Engineering Software 69 (2014) 46–61 57

Fig. 11 (continued)

Fig. 12. Tension/compression spring: (a) shematic, (b) stress heatmap (c) displacement heatmap.

This problem has also been popular among researchers and In summary, the results on the three classical engineering prob-
optimized in various studies. The heuristic methods that have been lems demonstrate that GWO shows high performance in solving
adopted to optimize this problem are: PSO [58], GA [57,60,69], ES challenging problems. This is again due to the operators that are
[59], DE [62], and ACO [70]. Mathematical methods used are aug- designed to allow GWO to avoid local optima successfully and con-
mented Lagrangian Multiplier [71] and branch-and-bound [72]. The verge towards the optimum quickly. The next section probes the
results of this problem are provided in Table 11. According to this ta- performance of the GWO algorithm in solving a recent real prob-
ble, GWO is again able to find a design with the minimum cost. lem in the field of optical engineering.
58 S. Mirjalili et al. / Advances in Engineering Software 69 (2014) 46–61

Table 9 Table 10
Comparison of results for tension/compression spring design problem. Comparison results of the welded beam design problem.

Algorithm Optimum variables Optimum Algorithm Optimum variables Optimum

weight cost
d D N h l t b
GWO 0.05169 0.356737 11.28885 0.012666 GWO 0.205676 3.478377 9.03681 0.205778 1.72624
GSA 0.050276 0.323680 13.525410 0.0127022 GSA 0.182129 3.856979 10.00000 0.202376 1.879952
PSO (Ha and Wang) 0.051728 0.357644 11.244543 0.0126747 GA (Coello) N/A N/A N/A N/A 1.8245
ES (Coello and Montes) 0.051989 0.363965 10.890522 0.0126810 GA (Deb) N/A N/A N/A N/A 2.3800
GA (Coello) 0.051480 0.351661 11.632201 0.0127048 GA (Deb) 0.2489 6.1730 8.1789 0.2533 2.4331
HS (Mahdavi et al.) 0.051154 0.349871 12.076432 0.0126706 HS (Lee and 0.2442 6.2231 8.2915 0.2443 2.3807
DE (Huang et al.) 0.051609 0.354714 11.410831 0.0126702 Geem)
Mathematical 0.053396 0.399180 9.1854000 0.0127303 Random 0.4575 4.7313 5.0853 0.6600 4.1185
optimization Simplex 0.2792 5.6256 7.7512 0.2796 2.5307
(Belegundu) David 0.2434 6.2552 8.2915 0.2444 2.3841
Constraint correction 0.050000 0.315900 14.250000 0.0128334 APPROX 0.2444 6.2189 8.2915 0.2444 2.3815
(Arora)

6. Real application of GWO in optical engineering (optical buffer There are two metrics for comparing the performance of slow
design) light devices: Delay-Bandwidth Product (DBP) and Normalized
DBP (NDBP), which are defined as follows [74]:
The problem investigated in this section is called optical buffer
DBP ¼ Dt  Df ð6:1Þ
design. In fact, an optical buffer is one of the main components of
optical CPUs. The optical buffer slows the group velocity of light where Dt indicates the delay and Df is the bandwidth of the slow
and allows the optical CPUs to process optical packets or adjust light device.
its timing. The most popular device to do this is a Photonic Crystal In slow light devices the ultimate goal is to achieve maximum
Waveguide (PCW). PCWs mostly have a lattice-shaped structure transmission delay of an optical pulse with highest PCW band-
with a line defect in the middle. The radii of holes and shape of width. Obviously, Dt should be increased in order to increase
the line defect yield different slow light characteristics. Varying ra- DBP. This is achieved by increasing the length of the device (L).
dii and line defects provides different environments for refracting To compare devices with different lengths and operating frequen-
the light in the waveguide. The researchers in this field try to cies, NDBP is a better choice [75]:
manipulate the radii of holes and pins of line defect in order to
NDBP ¼ ng  Dx=x0 ð6:2Þ
achieve desirable optical buffering characteristics. There are also
different types of PCW that are suitable for specific applications. where ng is the average of the group index, Dx is the normalized
In this section the structure of a PCW called a Bragg Slot PCW bandwidth, and x0 is the normalized central frequency of light
(BSPCW) is optimized by the GWO algorithm. This problem has wave.
several constraints, so we utilize the simplest constraint handling Since NDBP has a direct relation to the group index (ng), can be
method for GWO in this section as well. formulated as follows [76]:
BSPCW structure was first proposed by Caer et al. in 2011 [73].
C dk
The structure of BSPCWs is illustrated in Fig. 15. The background ng ¼ ¼C ð6:3Þ
slab is silicon with a refractive index equal to 3.48. The slot and
vg dx
holes are filled by a material with a refractive index of 1.6. The where x is the dispersion, k indicates the wave vector, C is the
Bragg slot structure allows the BSPCW to have precise control of velocity of light in free space, and shows the group index. Since ng
dispersion and slow light properties. The first five holes adjacent is changing in the bandwidth range, it should be averaged as
to the slot have the highest impact on slow light properties, as dis- follows:
cussed in [73]. As may be seen in Fig. 15, l, wl, and wh define the Z xH
shape of the slot and have an impact on the final dispersion and dx
ng ¼ n g ð xÞ ð6:4Þ
slow light properties as well. So, various dispersion and slow light xL Dx
properties can be achieved by manipulating the radii of holes, l, wl, The bandwidth of a PCW refers to the region of the ng curve where
and wh. ng has an approximately constant value with a maximum fluctua-

Fig. 13. Structure of welded beam design (a) shematic (b) stress heatmap (c) displacement heatmap.
 S. Mirjalili et al. / Advances in Engineering Software 69 (2014) 46–61 59

Fig. 14. Pressure vessel (a) shematic (b) stress heatmap (c) displacement heatmap.

Table 11
Comparison results for pressure vessel design problem.

Algorithm Optimum variables Optimum cost

Ts Th R L
GWO 0.812500 0.434500 42.089181 176.758731 6051.5639
GSA 1.125000 0.625000 55.9886598 84.4542025 8538.8359
PSO (He and Wang) 0.812500 0.437500 42.091266 176.746500 6061.0777
GA (Coello) 0.812500 0.434500 40.323900 200.000000 6288.7445
GA (Coello and Montes) 0.812500 0.437500 42.097398 176.654050 6059.9463
GA (Deb and Gene) 0.937500 0.500000 48.329000 112.679000 6410.3811
ES (Montes and Coello) 0.812500 0.437500 42.098087 176.640518 6059.7456
DE (Huang et al.) 0.812500 0.437500 42.098411 176.637690 6059.7340
ACO (Kaveh and Talataheri) 0.812500 0.437500 42.103624 176.572656 6059.0888
Lagrangian Multiplier (Kannan) 1.125000 0.625000 58.291000 43.6900000 7198.0428
Branch-bound (Sandgren) 1.125000 0.625000 47.700000 117.701000 8129.1036

R1 R2 R3 R4 R5 l wh wl
Consider : ~
x ¼ ½x1 x2 x3 x4 x5 x6 x7 x8  ¼ a a a a a a a a
;
a Super cell n Dx
Maximize : xÞ ¼ NDBP ¼ gx0 ;
f ð~

2R5 Subject to : maxðjb2 ðxÞjÞ < 106 a=2pc2 ;

2R4 xH < minðxup band Þ;
Si
2R3 xL > maxðxdown band Þ;
2R2 Filled kn > knH ! xGuided mode > xH ;
2R1 kn < knL ! xGuided mode < xL ; ð6:5Þ
wl wh
where : xH ¼ xðknH Þ ¼ xð1:1ng0 Þ;
l xL ¼ xðknL Þ ¼ xð0:9ng0 Þ;
kn ¼ 2kap
Dx ¼ xH  xL ;
a ¼ x0  1550ðnmÞ;

Fig. 15. BSPCW structure with super cell, nbackground = 3.48 and nfilled = 1.6. Variable range : 0 6 x15 6 0:5;
0 6 x6 6 1;
tion rage of ±10% [75]. Detailed information about PCWs can be
0 6 x7;8 6 1;
found in [77–80].
Finally, the problem is mathematically formulated for GWO as Note that we consider five constraints for the GWO algorithm.
follows: The second to fifth constraints avoid band mixing. To handle feasi-
bility, we assign small negative objective function values (100) to
Table 12 those search agents that violate the constraints.
Structural parameters and calculation results. The GWO algorithm was run 20 times on this problem and the
Structural parameter Wu et al. [81] GWO
best results obtained are reported in Table 12. Note that the algo-
rithm was run by 24 CPUs on a Windows HPC cluster at Griffith
R1 – 0.33235a
University. This table shows that there is a substantial, 93% and
R2 – 0.24952a
R3 – 0.26837a 65% improvement in bandwidth (Dk) and NDBP utilizing the
R4 – 0.29498a GWO algorithm.
R5 – 0.34992a The photonic band structure of the BSPCW optimized is shown
l – 0.7437a
in Fig. 16(a). In addition, the corresponded group index and opti-
Wh – 0.2014a
Wl – 0.60073a
mized super cell are shown in Figs. 16(b) and 17. These figures
a(nm) 430 343 show that the optimized structure has a very good bandwidth
g
n 23 19.6 without band mixing as well. This again demonstrated the high
Dk(nm) 17.6 33.9 performance of the GWO algorithm in solving real problems.
Order of magnitude of b2 (a/2pc2) 103 103
This comprehensive study shows that the proposed GWO algo-
NDBP 0.26 0.43
rithm has merit among the current meta-heuristics. First, the re-
60 S. Mirjalili et al. / Advances in Engineering Software 69 (2014) 46–61

100
0.26

Normalized Frequency (ωa/2πc=a/λ)

0.24 Light line 80

The Group Index—ng

of SiO2

0.22 60
Guided mode

0.2 40 Δω

ωH
0.18 20
ω0
ωL
0.16 0
0.25 0.3 0.35 0.4 0.45 0.5 0.214 0.216 0.218 0.22 0.222 0.224
Wavevector--ka/2π Normalized Frequency (ω a/2π c=a/ )
(a) (b)
Fig. 16. (a) Photonic band structure of the optimized BSPCW structure (b) The group index (ng) of the optimized BSPCW structure.

nally, the convergence analysis of GWO confirmed the convergence

R5=120 nm of this algorithm.
Moreover, the results of the engineering design problems also
R4=101 nm
showed that the GWO algorithm has high performance in un-
R3=92 nm known, challenging search spaces. The GWO algorithm was finally
R2=84 nm applied to a real problem in optical engineering. The results on this
problem showed a substantial improvement of NDBP compared to
R1=114 nm current approaches, showing the applicability of the proposed
wl=69 nm wh=206 nm algorithm in solving real problems. It may be noted that the results
on semi-real and real problems also proved that GWO can show
high performance not only on unconstrained problems but also
on constrained problems.
l=255 nm For future work, we are going to develop binary and multi-
objective versions of the GWO algorithm.

Appendix A. Supplementary material

Fig. 17. Optimized super cell of BSPCW.

Supplementary data associated with this article can be found, in

the online version, at http://dx.doi.org/10.1016/j.advengsoft.2013.
sults of the unconstrained benchmark functions demonstrate the
12.007.
performance of the GWO algorithm in terms of exploration, exploi-
tation, local optima avoidance, and convergence. Second, the re-
sults of the classical engineering problems show the superior
References
performance of the proposed algorithm in solving semi-real con-
strained problems. Finally, the results of the optical buffer design [1] Bonabeau E, Dorigo M, Theraulaz G. Swarm intelligence: from natural to
problem show the ability of the GWO algorithm in solving the real artificial systems: OUP USA; 1999.
problems. [2] Dorigo M, Birattari M, Stutzle T. Ant colony optimization. Comput Intell Magaz,
IEEE 2006;1:28–39.
[3] Kennedy J, Eberhart R. Particle swarm optimization, in Neural Networks, 1995.
In: Proceedings, IEEE international conference on; 1995. p. 1942–1948.
[4] Wolpert DH, Macready WG. No free lunch theorems for optimization. Evolut
7. Conclusion Comput, IEEE Trans 1997;1:67–82.
[5] Kirkpatrick S, Jr. DG, Vecchi MP. Optimization by simulated annealing. Science,
vol. 220; 1983. p. 671–80.
This work proposed a novel SI optimization algorithm inspired [6] Beni G, Wang J. Swarm intelligence in cellular robotic systems. In: Robots and
by grey wolves. The proposed method mimicked the social hierar- biological systems: towards a new bionics?, ed. Springer; 1993. p. 703–12.
chy and hunting behavior of grey wolves. Twenty nine test func- [7] Basturk B, Karaboga D. An artificial bee colony (ABC) algorithm for numeric
function optimization. In: IEEE swarm intelligence symposium; 2006. p. 12–4.
tions were employed in order to benchmark the performance of [8] Olorunda O, Engelbrecht AP. Measuring exploration/exploitation in particle
the proposed algorithm in terms of exploration, exploitation, local swarms using swarm diversity. In: Evolutionary computation, 2008. CEC 2008
optima avoidance, and convergence. The results showed that GWO (IEEE World Congress on Computational Intelligence). IEEE Congress on; 2008.
p. 1128–34.
was able to provide highly competitive results compared to well-
[9] Alba E, Dorronsoro B. The exploration/exploitation tradeoff in dynamic cellular
known heuristics such as PSO, GSA, DE, EP, and ES. First, the results genetic algorithms. Evolut Comput, IEEE Trans 2005;9:126–42.
on the unimodal functions showed the superior exploitation of the [10] Lin L, Gen M. Auto-tuning strategy for evolutionary algorithms: balancing
GWO algorithm. Second, the exploration ability of GWO was con- between exploration and exploitation. Soft Comput 2009;13:157–68.
[11] Mirjalili S, Hashim SZM. A new hybrid PSOGSA algorithm for function
firmed by the results on multimodal functions. Third, the results optimization. In: Computer and information application (ICCIA), 2010
of the composite functions showed high local optima avoidance. Fi- international conference on; 2010. p. 374–77.
 S. Mirjalili et al. / Advances in Engineering Software 69 (2014) 46–61 61

[12] Mirjalili S, Mohd Hashim SZ, Moradian Sardroudi H. Training feedforward [48] Digalakis J, Margaritis K. On benchmarking functions for genetic algorithms.
neural networks using hybrid particle swarm optimization and gravitational Int J Comput Math 2001;77:481–506.
search algorithm. Appl Math Comput 2012;218:11125–37. [49] Molga M, Smutnicki C. Test functions for optimization needs. Test functions for
[13] Holland JH. Genetic algorithms. Sci Am 1992;267:66–72. optimization needs; 2005.
[14] Goldberg D. Genetic Algorithms in optimization, search and machine learning, [50] Yang X-S. Test problems in optimization, arXiv, preprint arXiv:1008.0549;
Addison Wesley, New York. In: Eiben AE, Smith JE, editors. 2003 Introduction 2010.
to evolutionary computing. Springer. Jacq J, Roux C (1995) Registration of non- [51] Mirjalili S, Lewis A. S-shaped versus V-shaped transfer functions for binary
segmented images using a genetic algorithm. Lecture notes in computer Particle Swarm Optimization. Swarm Evolut Comput 2013;9:1–14.
science, vol. 905; 1989. p. 205–11. [52] Liang J, Suganthan P, Deb K. Novel composition test functions for numerical
[15] Storn R, Price K. Differential evolution – a simple and efficient heuristic for global optimization. In: Swarm intelligence symposium, 2005. SIS 2005.
global optimization over continuous spaces. J Global Optim 1997;11:341–59. Proceedings 2005 IEEE; 2005. p. 68–75.
[16] Yao X, Liu Y, Lin G. Evolutionary programming made faster. Evolut Comput, [53] Suganthan PN, Hansen N, Liang JJ, Deb K, Chen YP, Auger A, et al. Problem
IEEE Trans 1999;3:82–102. definitions and evaluation criteria for the CEC 2005 special session on real-
[17] Fogel D. Artificial intelligence through simulated evolution. Wiley-IEEE Press; parameter optimization, Technical Report, Nanyang Technological University,
2009. Singapore, 2005, http://www.ntu.edu.sg/home/EPNSugan.
[18] Hansen N, Müller SD, Koumoutsakos P. Reducing the time complexity of the [54] van den Bergh F, Engelbrecht A. A study of particle swarm optimization
derandomized evolution strategy with covariance matrix adaptation (CMA- particle trajectories. Inf Sci 2006;176:937–71.
ES). Evolut Comput 2003;11:1–18. [55] Arora JS. Introduction to optimum design. Academic Press; 2004.
[19] Rechenberg I. Evolution strategy. Comput Intel Imitat Life 1994;1. [56] Belegundu AD, Arora JS. A Study of mathematical programming methods for
[20] Koza JR. Genetic programming; 1992. structural optimization. Part I: Theory. Int J Numer Meth Eng 1985;21:1583–99.
[21] Simon D. Biogeography-based optimization. Evolut Comput IEEE Trans [57] Coello Coello CA, Mezura Montes E. Constraint-handling in genetic algorithms
2008;12:702–13. through the use of dominance-based tournament selection. Adv Eng Inform
[22] Webster B, Bernhard PJ. A local search optimization algorithm based on 2002;16:193–203.
natural principles of gravitation. In: Proceedings of the 2003 international [58] He Q, Wang L. An effective co-evolutionary particle swarm optimization for
conference on information and knowledge engineering (IKE’03), Las Vegas, constrained engineering design problems. Eng Appl Artif Intell 2007;20:89–99.
Nevada, USA; 2003. p. 255–61. [59] Mezura-Montes E, Coello CAC. An empirical study about the usefulness of
[23] Erol OK, Eksin I. A new optimization method: big bang–big crunch. Adv Eng evolution strategies to solve constrained optimization problems. Int J Gen Syst
Softw 2006;37:106–11. 2008;37:443–73.
[24] Rashedi E, Nezamabadi-Pour H, Saryazdi S. GSA: a gravitational search [60] Coello Coello CA. Use of a self-adaptive penalty approach for engineering
algorithm. Inf Sci 2009;179:2232–48. optimization problems. Comput Ind 2000;41:113–27.
[25] Kaveh A, Talatahari S. A novel heuristic optimization method: charged system [61] Mahdavi M, Fesanghary M, Damangir E. An improved harmony search
search. Acta Mech 2010;213:267–89. algorithm for solving optimization problems. Appl Math Comput
[26] Formato RA. Central force optimization: a new metaheuristic with applications 2007;188:1567–79.
in applied electromagnetics. Prog Electromag Res 2007;77:425–91. [62] Huang F, Wang L, He Q. An effective co-evolutionary differential evolution for
[27] Alatas B. ACROA: artificial chemical reaction optimization algorithm for global constrained optimization. Appl Math Comput 2007;186:340–56.
optimization. Expert Syst Appl 2011;38:13170–80. [63] Yang XS. Nature-inspired metaheuristic algorithms. Luniver Press; 2011.
[28] Hatamlou A. Black hole: a new heuristic optimization approach for data [64] Carlos A, COELLO C. Constraint-handling using an evolutionary multiobjective
clustering. Inf Sci 2012. optimization technique. Civil Eng Syst 2000;17:319–46.
[29] Kaveh A, Khayatazad M. A new meta-heuristic method: ray optimization. [65] Deb K. Optimal design of a welded beam via genetic algorithms. AIAA J
Comput Struct 2012;112:283–94. 1991;29:2013–5.
[30] Du H, Wu X, Zhuang J. Small-world optimization algorithm for function [66] Deb K. An efficient constraint handling method for genetic algorithms. Comput
optimization. In: Advances in Natural Computation, ed.: Springer; 2006. p. Methods Appl Mech Eng 2000;186:311–38.
264–73. [67] Lee KS, Geem ZW. A new meta-heuristic algorithm for continuous engineering
[31] Shah-Hosseini H. Principal components analysis by the galaxy-based search optimization: harmony search theory and practice. Comput Methods Appl
algorithm: a novel metaheuristic for continuous optimisation. Int J Comput Sci Mech Eng 2005;194:3902–33.
Eng 2011;6:132–40. [68] Ragsdell K, Phillips D. Optimal design of a class of welded structures using
[32] Moghaddam FF, Moghaddam RF, Cheriet M. Curved space optimization: a geometric programming. ASME J Eng Indust 1976;98:1021–5.
random search based on general relativity theory. arXiv, preprint [69] Deb K, Gene AS. A robust optimal design technique for mechanical component
arXiv:1208.2214; 2012. design. In: Presented at the Dasgupta D, Michalewicz Z, editors. Evolutionary
[33] Yang X-S. A new metaheuristic bat-inspired algorithm. In: Nature inspired algorithms in engineering applications, Berlin; 1997.
cooperative strategies for optimization (NICSO 2010), ed.: Springer; 2010. p. [70] Kaveh A, Talatahari S. An improved ant colony optimization for constrained
65–74. engineering design problems. Eng Comput Int J Comput-Aided Eng
[34] Abbass HA. MBO: Marriage in honey bees optimization – a haplometrosis 2010;27:155–82.
polygynous swarming approach. In: Evolutionary computation, 2001. [71] Kannan B, Kramer SN. An augmented Lagrange multiplier based method for
Proceedings of the 2001 congress on; 2001. p. 207–214. mixed integer discrete continuous optimization and its applications to
[35] Li X. A new intelligent optimization-artificial fish swarm algorithm. Doctor mechanical design. J Mech Des 1994;116:405.
thesis, Zhejiang University of Zhejiang, China; 2003. [72] Sandgren E. Nonlinear integer and discrete programming in mechanical
[36] Roth M. Termite: a swarm intelligent routing algorithm for mobile wireless design; 1988. p. 95–105.
ad-hoc networks; 2005. [73] Caer C, Le Roux X, Marris-Morini D, Izard N, Vivien L, Gao D, et al. Dispersion
[37] Pinto PC, Runkler TA, Sousa JM. Wasp swarm algorithm for dynamic MAX-SAT engineering of wide slot photonic crystal waveguides by Bragg-like
problems. In: Adaptive and Natural Computing Algorithms, ed.: Springer; corrugation of the slot. Photonics Technol Lett, IEEE 2011;23:1298–300.
2007. p. 350–57. [74] Baba T. Slow light in photonic crystals. Nat Photonics 2008;2:465–73.
[38] Mucherino A, Seref O. Monkey search: a novel metaheuristic search for global [75] Zhai Y, Tian H, Ji Y. Slow light property improvement and optical buffer
optimization. In: AIP conference proceedings; 2007. p. 162. capability in ring-shape-hole photonic crystal waveguide. Light Technol J
[39] Lu X, Zhou Y. A novel global convergence algorithm: bee collecting pollen 2011;29:3083–90.
algorithm. In: Advanced intelligent computing theories and applications. With [76] Wang D, Zhang J, Yuan L, Lei J, Chen S, Han J, et al. Slow light engineering in
Aspects of Artificial Intelligence, ed.: Springer; 2008. p. 518–25. polyatomic photonic crystal waveguides based on square lattice. Optics
[40] Yang X-S, Deb S. Cuckoo search via Lévy flights. In: Nature & Biologically Commun 2011;284:5829–32.
Inspired Computing, 2009. NaBIC 2009. World Congress on; 2009. p. 210–14. [77] Mirjalili SM, Mirjalili S. Light property and optical buffer performance
[41] Shiqin Y, Jianjun J, Guangxing Y. A dolphin partner optimization. In: Intelligent enhancement using Particle Swarm Optimization in Oblique Ring-Shape-
systems, 2009. GCIS’09. WRI Global Congress on; 2009. p. 124–28. Hole Photonic Crystal Waveguide. In: Photonics global conference (PGC);
[42] Yang X-S. Firefly algorithm, stochastic test functions and design optimisation. 2012. p. 1–4 [2012].
Int J Bio-Inspired Comput 2010;2:78–84. [78] Mirjalili SM, Abedi K, Mirjalili S. Optical buffer performance enhancement
[43] Askarzadeh A, Rezazadeh A. A new heuristic optimization algorithm for using Particle Swarm Optimization in Ring-Shape-Hole Photonic Crystal
modeling of proton exchange membrane fuel cell: bird mating optimizer. Int J Waveguide. Optik – Int J Light Elect Optics 2013;124:5989–93.
Energy Res 2012. [79] Mirjalili SM, Mirjalili S, Lewis A. A novel multi-objective optimization
[44] Gandomi AH, Alavi AH. Krill Herd: a new bio-inspired optimization algorithm. framework for designing photonic crystal waveguides. Photonics Technol
Commun Nonlinear Sci Numer Simul 2012. Lett IEEE 2014;26:146–9.
[45] Pan W-T. A new fruit fly optimization algorithm: taking the financial distress [80] Mirjalili SM, Mirjalili S, Lewis A, Abedi K. A tri-objective particle swarm
model as an example. Knowl-Based Syst 2012;26:69–74. optimizer for designing line defect photonic crystal waveguides. Photonics and
[46] Mech LD. Alpha status, dominance, and division of labor in wolf packs. Can J Nanostructures – Fundamentals and Applications.
Zool 1999;77:1196–203. [81] Wu J, Li Y, Peng C, Wang Z. Wideband and low dispersion slow light in slotted
[47] Muro C, Escobedo R, Spector L, Coppinger R. Wolf-pack (Canis lupus) hunting photonic crystal waveguide. Optics Commun 2010;283:2815–9.
strategies emerge from simple rules in computational simulations. Behav [82] Mirjalili S, Mirjalili SM, Yang X. Binary bat algorithm. Neural Comput Appl, in
Process 2011;88:192–7. press, DOI: 10.1007/s00521-013-1525-5.