A Novel Complex-Valued Encoding Grey Wolf Optimization Algorithm

: Grey wolf optimization (GWO) is one of the recently proposed heuristic algorithms imitating the leadership hierarchy and hunting mechanism of grey wolves in nature. The aim of these algorithms is to perform global optimization. This paper presents a modiﬁed GWO algorithm based on complex-valued encoding; namely the complex-valued encoding grey wolf optimization (CGWO). We use CGWO to test 16 unconstrained benchmark functions with seven different scales and inﬁnite impulse response (IIR) model identiﬁcation. Compared to the real-valued GWO algorithm and other optimization algorithms; the CGWO performs signiﬁcantly better in terms of accuracy; robustness; and convergence speed.


Introduction
Over the last few decades, various heuristic optimizations have been developed to solve complex computational problems.Some of the most popular are: particle swarm optimization (PSO) [1], genetic algorithm (GA) [2], differential evolution (DE) [3,4], simulated annealing (SA) [5], ant colony optimization (ACO) [6], artificial bee colony optimization (ABC) [7], gravitational search algorithm (GSA) [8,9], charged system search (CSS) [10], magnetic optimization algorithm (MOA) [11], biogeography-based optimization algorithm (BBO) [12,13], and teaching learning-based optimization (TLBO) [14].The swarm intelligence optimization algorithm can flexibly and effectively deal with different problems that conventional optimization techniques cannot solve.For this reason, optimization techniques have been widely used in many domains of study.In addition, it has been demonstrated that there is no heuristic optimization that can solve all optimization problems [15].In other words, existing algorithms give satisfactory results in solving some problems, but not all problems.Therefore, several new heuristic algorithms are proposed every year, and research in this field is active.
The real-value grey wolf optimization (GWO) is a new population-based meta-heuristic approach initially proposed by Mirjalili et al. (2014) [16].Due to the fact that the GWO algorithm is simple, flexible and efficient, it can be successfully applied to practical applications [17][18][19].There are many modified versions of the GWO algorithm.For example, some people study the parameters of the GWO algorithm or combine the GWO algorithm with other heuristic optimization algorithms.However, these methods all use binary or decimal encoders to position the grey wolf, the individual gene contains information to be limited.
The complex-valued encoding method has been applied in expressing neural network weights [20] and representing individual genes for evolutionary algorithms [21].It uses a diploid to express genes and greatly expand the information capacity of the individual.Starting from the individual encoding method, this article studies the improvement of the performance of the grey wolf with the complex-valued method.The value of the independent variables for an objective function is determined by modules and the signthe independent variables is determined by angles.For the GWO algorithm, the grey wolf individual is divided into two parts, namely: a real part gene and an imaginary gene.Additionally, the two parts represent a variable being able to enhance the information of a grey wolf individual and a diverse population, avoiding the local optimum.Based on the above advantages, we provide a new method for the GWO algorithm to solve practical problems.
The rest of the paper is organized as follows: Section 2 presents a brief introduction to GWO.Section 3 discusses the basic principles of a complex-valued version of GWO.The experimental results of test functions and infinite impulse response (IIR) system identification problem are shown in Sections 4 and 5 respectively.Finally, Section 6 concludes the work and suggests some directions for future study.

Grey Wolf Optimization (GWO)
The GWO algorithm simulates the hunting and social leadership of grey wolves in nature [16].The algorithm is simple, robust, and has been used in various complex problems.In the GWO algorithm, the colony of grey wolves is divided into four groups: alpha (α), beta (β), delta (δ), and omega (ω).In every iteration, the first three best candidate solutions are named α, β, and δ.The rest of the grey wolves are considered as ω, and are guided by α, β, and δ to find the better solutions.The mathematical model of the ω wolves' encircling process is as follows [16]: where t indicates the current iteration, X p is the position vector of the prey, Ñ X indicates the position vector of a grey wolf, Ñ α is gradually decreased from 2 to 0, and r 1 , r 2 are random numbers over range r0, 1s.
It should be noted here that, during optimization, the ω wolves update their positions around α, β, and δ.Therefore, ω wolves are able to reposition with respect to α, β, and δ.The mathematical model of ω wolves readjusting the positions is as follows [16]: where Ñ X α is the position of the alpha, C that are random and adaptive vectors which help the algorithm to escape from local optima.When |A ą 1|, half of the iterations are devoted to exploration.The range of C is 2 ď C ď 0. The vector C also improves exploration when C ą 1.In contrast, the other half of the iterations are dedicated to exploitation when |A| ă 1, and the exploitation is also emphasized when C ă 1.Note here that A decreases linearly over the course of iterations in order to emphasize exploitation, and C is provided with a random value at all times to emphasize exploration not only during initial iterations but also during final iterations.The pseudo code of GWO (Algorithm 1) is presented as follows: Algorithm 1. GWO algorithm 1. Initialize the grey wolf population X i " pi " 1, 2, ..., nq 2. Initialize a, A, and C 3. For all X i do 4.
Calculate fitness FpX i q 5. End for 6.Get the first three best wolves as X α , X β , and X δ 7.While (t ă Max number o f iterations) 8.
For each search agent 9.
Update the position of the current search agent by Equation (5) 10.
End for 11.
Update a, A, and C 12.
For all X i do 13.
Calculate fitness FpX i q 14.
End for 15.

The Complex-Valued Encoding Method
In order to contain the M´variable function optimization problem within the M complex, corresponding to the M complex, the grey wolf location is recorded as follows: The genes of grey wolves can be expressed as a diploid, which is recorded as (R p , I p ).Where R p indicates the real part of the variable, I p indicates the imaginary part of the variable.Therefore, the ith grey wolf can be seen in Table 1.

Gene 1
Gene 2 Gene i Gene M pR p1 , I p1 q pR p2 , I p2 q . . .pR pM , I pM q

Initializing the Complex-Valued Encoding Population
The variable interval of the function is set to rA L , B L s, L " 1, 2, ... , M. Generate M modulus and M phase angle randomly, and they form a complex-valued M. The relationship of the modulus and the phase is provided as follows: So, we obtain M real and imaginary parts.The real and imaginary parts are updated according to Section 3.1.2.

Update the Real Parts
Update the Imaginary Parts C are coefficient vectors, so it is not necessary for them to be converted to real and imaginary.

The Calculation Method of Fitness Value
The complex domain is composed of real and imaginary parts, so the complex number must be changed into a real number when solving the fitness function value.The real value of an objective function is determined by its modulus, and the sign is determined by its phase angle.The concrete practices are as follows: RV n " ρ n sgnpsinp X In where ρ n indicates the nth dimension module, X 2 R n , X 2 I n indicate the real part and imaginary part of the nth dimension, respectively, and RV n is converted into a real variable.

CGWO Algorithm
Based on the above analysis, complex-valued encoding, which can be considered as a valid global optimization strategy, is applied to the grey wolf optimizer.The real and imaginary parts of complex numbers are updated independently due to two-dimensional properties of complex numbers.The strategy greatly expands the amount of information contained in the individual gene and enhances the diversity of individual populations.In addition, to improve the local search ability of the algorithm, the differential evolution strategy "DE/best/2" [22] is introduced.Under the circumstances, CGWO has the potential to balance global and local searches.In the following section, various benchmark functions are employed to explore the effectiveness of CGWO.The pseudo code of CGWO (Algorithm 2) is as follows: Algorithm 2. CGWO algorithm 1. Initialize the grey wolf population: Let ρ L " Get the real and imaginary part of the complex number by Equation (8) 3.
For all X i do 5.
Calculate fitness FpX i q 6. End for 7.
Get the first three best wolves as X α , X β , and X δ 8.While (t ă Max number o f iterations) 9.
For each search agent 10.
Update the real part of the current search agent by Equations (9-11) 11.
Update the imaginary part of the current search agent by Equations (12-14) 12.
End for 13.
Modify the real and imaginary part by "DE/best/2" 14.
Select randomly Update a, A, and C 17.
For all X i do 19.
Calculate fitness FpX i q 20.
End for 21.

Simulation Platform
All algorithms are tested in Matlab R2012a (7.14), and experiments are executed on AMD Athlon (tm) II X4 640 Processor 3.00 GHz PC with 3G RAM.Operating system is Windows 7.

Benchmark Functions
To evaluate the performance of the proposed method, termed CGWO, 16 standard benchmark functions are used to demonstrate the performance of the algorithm [23][24][25].They can be divided into two groups: unimodal and multimodal functions.Table 1 lists the functions, where Dim shows dimension of the function, Range is the boundary of search space of the function, and f min is the minimum of the function.
We compare ABC [7], GGSA [9], CGWO, GWO [16] on benchmark functions to verify the effectiveness of these algorithms dealing with problems of different scales.Moreover, in order to prove the superiority of CGWO over GWO, we derive the functions from complex-valued encoding rather than the differential evolution strategy "DE/best/2".CGWO is also compared to the basic GWO with a differential evolution strategy called differential evolution GWO (DGWO).The five algorithms have some parameters that are shown in Table 2.The experimental initial parameters of the five algorithms are provided in Table 3, the experimental results are provided in Tables 4 and 5.The results are averaged over 20 consecutive experiments.The best results are demonstrated in bold type.The average (AVE) and standard deviation (STD) of the best solution obtained in the last iteration are reported in the form of AVE ˘STD.Note that the Matlab code of the GGSA algorithm is given in [26].
To improve the performance evaluation of evolutionary algorithms, statistical tests should be conducted [27].In order to determine whether the results of CGWO differ from the best results of the other algorithms in a statistical method, a nonparametric test, which is known as Wilcoxon's rank-sum test [28][29][30][31], is performed at 5% significance level.The calculated p values in Wilcoxon's rank-sum test comparing CGWO and the other algorithms over all the benchmark functions are given in Tables 6-12.Usually, p values < 0.05 can be considered as sufficient evidence against the null hypothesis.

Unimodal Benchmark Functions
The unimodal benchmarks have only one global solution without any local optima for them.Consequently, they are useful in examining the convergence capability of heuristic optimization.The results of unimodal benchmark functions are shown in Table 3.As shown in this table, for 30-D unimodal benchmark functions, CGWO has the best results for all benchmark functions.For 10-D, 50-D, 100-D, 200-D, 300-D, and 400-D, this table shows that the CGWO algorithm outperforms the other algorithms on the majority of the unimodal benchmark test functions.Therefore, the proposed algorithm has high performance in finding the global solution of unimodal benchmark test functions.According to the p values of F 1 " F 8 with different dimensions in Tables 5-11 CGWO achieves significant improvement across the majority of dimensions compared to the other algorithms.Hence, this proves that CGWO has better performance than the other algorithms in forging for a global optimum solution of unimodal benchmark functions.
Figure 1 shows the averaged convergence curves of unimodal benchmark functions with 30-D over 20 independent runs.As shown from these curves, the CGWO has the fastest convergence rate.Figure 2 indicates the anova tests of the global minimum on unimodal functions with 30-D, which shows that CGWO is the most robust method.For F 1 " F 8 functions, the other algorithms that we compared (namely, DGWO, GWO, GGSA, ABC) were only robust for a few functions and were not able to be robust for all functions.On the basis of Tables 3 and 5-11 and Figures 1 and 3 it can be concluded that the CGWO tends to find the best solution in unimodal functions with a very fast convergence rate and a more stable performance.F F with different dimensions in Tables 5-11, CGWO achieves significant improvement across the majority of dimensions compared to the other algorithms.Hence, this proves that CGWO has better performance than the other algorithms in forging for a global optimum solution of unimodal benchmark functions.
Figure 1 shows the averaged convergence curves of unimodal benchmark functions with 30-D over 20 independent runs.As shown from these curves, the CGWO has the fastest convergence rate.
Figure 2 indicates the anova tests of the global minimum on unimodal functions with 30-D, which shows that CGWO is the most robust method.For 1 8 F F functions, the other algorithms that we compared (namely, DGWO, GWO, GGSA, ABC) were only robust for a few functions and were not able to be robust for all functions.On the basis of Tables 3, 5-11 and Figures 1 and 3, it can be concluded that the CGWO tends to find the best solution in unimodal functions with a very fast convergence rate and a more stable performance.

Multimodal Benchmark Functions
The results of the multimodal benchmark functions are provided in Table 4.Note that multimodal functions have many local minima with the number increasing exponentially with dimension.Therefore, they are very suitable to test the exploration capability of algorithms to avoid local minima.As the results of the mean and STD show, CGWO performs better than the other algorithms on 30-D, 50-D multimodal benchmark functions.For 9 F on 200-D, 300-D, 400-D, the GWO has better value than the CGWO.For 12 F on 10-D, the GGSA has best result.For 14 F on 100-D, the GWO is slightly better than the CGWO.For 16 F on 100-D, CGWO and GGSA have very close results.This table also shows that CGWO provides better results than the other algorithms across the majority of dimensions on 10-D, 100-D, 200-D, 300-D, and 400-D.Moreover, the p values of 9 F F with different dimensions reported in Tables 5-11 are less than 0.05 across the majority of dimensions, which is strong evidence against null hypothesis.So, this evidence demonstrates that the results of CGWO are statistically significant and do not occur by coincidence.Figure 3 shows the averaged convergence curves of multimodal benchmark functions with 30-D over 20 independent runs.As shown in Figure 2, the convergence rate of CGWO in all cases is better than the other algorithms.Figure 4 depicts the anova tests between CGWO and the other algorithms on multimodal functions with 30-D.As shown, the more stable performance belongs to CGWO.According to Tables 4, 5-11 and Figures 3 and 4, it can be stated that CGWO is capable of avoiding local minima with good convergence speed and stable property when dealing with multimodal benchmark functions.
According to this comprehensive comparative study and discussion, we state that this algorithm has merit compared to the other algorithms.The next section inspects the performance of the CGWO in solving the IIR model identification problem.

Multimodal Benchmark Functions
The results of the multimodal benchmark functions are provided in Table 4.Note that multimodal functions have many local minima with the number increasing exponentially with dimension.Therefore, they are very suitable to test the exploration capability of algorithms to avoid local minima.As the results of the mean and STD show, CGWO performs better than the other algorithms on 30-D, 50-D multimodal benchmark functions.For F 9 on 200-D, 300-D, 400-D, the GWO has better value than the CGWO.For F 12 on 10-D, the GGSA has best result.For F 14 on 100-D, the GWO is slightly better than the CGWO.For F 16 on 100-D, CGWO and GGSA have very close results.This table also shows that CGWO provides better results than the other algorithms across the majority of dimensions on 10-D, 100-D, 200-D, 300-D, and 400-D.Moreover, the p values of F 9 " F 16 with different dimensions reported in  are less than 0.05 across the majority of dimensions, which is strong evidence against null hypothesis.So, this evidence demonstrates that the results of CGWO are statistically significant and do not occur by coincidence.
Figure 3 shows the averaged convergence curves of multimodal benchmark functions with 30-D over 20 independent runs.As shown in Figure 2, the convergence rate of CGWO in all cases is better than the other algorithms.Figure 4 depicts the anova tests between CGWO and the other algorithms on multimodal functions with 30-D.As shown, the more stable performance belongs to CGWO.According to Tables 4-11 and Figures 3 and 4 it can be stated that CGWO is capable of avoiding local minima with good convergence speed and stable property when dealing with multimodal benchmark functions.
According to this comprehensive comparative study and discussion, we state that this algorithm has merit compared to the other algorithms.The next section inspects the performance of the CGWO in solving the IIR model identification problem.

IIR Model Identification
System identification is an effective method to obtain the mathematical model of an unknown system by analyzing its input and output values.In system identification, the use of infinite impulse response (IIR) has been widely studied since many problems encountered in signal processing can be characterized as a system identification problem (Figure 5).In this case, an optimization algorithm tries to iteratively adjust the parameters until the error between the output of the candidate model and the actual output of the plant is minimized.The use of IIR models for identification is preferred over their equivalent FIR (finite impulse response) models since the former produce more accurate models of physical plants for real world applications [32].Moreover, IIR models are able to use fewer model parameters to meet the performance specifications.The fitness function, which is used in the article, is calculated as follows [33].

IIR Model Identification
System identification is an effective method to obtain the mathematical model of an unknown system by analyzing its input and output values.In system identification, the use of infinite impulse response (IIR) has been widely studied since many problems encountered in signal processing can be characterized as a system identification problem (Figure 5).In this case, an optimization algorithm tries to iteratively adjust the parameters until the error between the output of the candidate model and the actual output of the plant is minimized.The use of IIR models for identification is preferred over their equivalent FIR (finite impulse response) models since the former produce more accurate models of physical plants for real world applications [32].Moreover, IIR models are able to use fewer model parameters to meet the performance specifications.The fitness function, which is used in the article, is calculated as follows [33].

IIR Model Identification
System identification is an effective method to obtain the mathematical model of an unknown system by analyzing its input and output values.In system identification, the use of infinite impulse response (IIR) has been widely studied since many problems encountered in signal processing can be characterized as a system identification problem (Figure 5).In this case, an optimization algorithm tries to iteratively adjust the parameters until the error between the output of the candidate model and the actual output of the plant is minimized.The use of IIR models for identification is preferred over their equivalent FIR (finite impulse response) models since the former produce more accurate models of physical plants for real world applications [32].Moreover, IIR models are able to use fewer model parameters to meet the performance specifications.The fitness function, which is used in the article, is calculated as follows [33].
Block diagram of IIR system identification.In an IIR system, the transfer function is provided as follows: where m and n are the number of numerator and denominator coefficients of the transfer function, respectively, where m and n are the number of numerator and denominator coefficients of the transfer function, respectively, a i pi P r1, ... , nsq is the pole of the IIR model, b j pj P r1, ... , msq is the zero parameters of the IIR model.Equation ( 17) is expressed in another form as follows: where uptq indicates the tth input of the system, yptq indicates the tth output of the system.Hence, the set of unknown parameters that simulates the IIR system is indicated by θ " ta 1 , ..., a n , b 0 , ..., b m u.The number of unknown parameters of θ is pn `m `1q.Correspondingly, the search space S of IIR feasible values for θ is pn`m`1q .According to Figure 5, the output of the plant is dptq, while the output of the IIR filter is yptq.The error eptq indicates the difference between the actual system and its model yields.Here the error is eptq " dptq ´yptq.Therefore, the problem of IIR model identification can be solved as a minimization problem of function.Additionally, the function can be molded as follows: where W indicates the number of the samples used in the simulation, and the aim is to minimize the function f pθq through adjusting θ.When the error function f pθq reaches its minimum value, the optimal model θ ˚or solution is obtained.The θ ˚is shown as follows:

IIR Model Identification Result
The result is reported considering a superior-order plant through a high-order IIR model.The unknown plant H p and the IIR model H M hold the following transfer function [27]: Since the plant is a sixth-order system and the IIR model a fourth-order system, the error surface f pθq is multimodal.A white sequence with 100 samples has been used as input.The result of the experience over 20 consecutive experiments is shown in Tables 5 and 6.Table 13 shows the best parameter values (ABP).Table 14 indicates the average f pθq value (AVE) and the standard deviation (STD).Figure 6 provides a comparison between CGWO and the other algorithms for IIR system identification.According to the AVE and STD indexes in Table 6, the CGWO algorithm finds better results than the other algorithms.The results show that CGWO provides better precision (AVE value) and robustness (STD).As shown in Figure 6, the CGWO has the fastest convergence rate.
In conclusion, the results show that the proposed method successfully outperforms the other algorithms in a majority of benchmark functions.In addition, IIR model identification shows that CGWO is able to obtain very competitive results.Therefore, as shown in the results of this comparative study, the proposed method has merit in the field of evolutionary algorithms and According to the AVE and STD indexes in Table 6, the CGWO algorithm finds better results than the other algorithms.The results show that CGWO provides better precision (AVE value) and robustness (STD).As shown in Figure 6, the CGWO has the fastest convergence rate.
In conclusion, the results show that the proposed method successfully outperforms the other algorithms in a majority of benchmark functions.In addition, IIR model identification shows that CGWO is able to obtain very competitive results.Therefore, as shown in the results of this comparative study, the proposed method has merit in the field of evolutionary algorithms and optimization.

Conclusions
In this work, the idea of complex-valued encoding is introduced into GWO, and a complex-valued version of GWO called CGWO was proposed based on complex-valued encoding.The proposed functions justified its excellent property on 16 benchmark functions in terms of improved avoidance of local minima and increased convergence speed.Furthermore, a complex system identification problem been employed to test the performance of the proposed method.The results verified the superior performance of CGWO on IIR model identification when compared to the other algorithms because of its improved ability.For future studies, it is recommended to apply the CGWO to solve more real-world engineering problems.

ÑXX
β is the position of the beta, Ñ X δ is the position of the delta, is the position of the current solution, and t indicates the number of iterations.As mentioned before, ω wolves update their positions based on positions of α, β, and δ.There are two vectors: Ñ A and Ñ the real and imaginary of the approximate distance between the current solution, alpha, beta, and delta respectively.Ñ X R and Ñ X L indicate the real and imaginary parts of the final position of the current solution.It should be noted here that Ñ A and Ñ

Figure 5 .
Figure 5. Block diagram of IIR system identification.In an IIR system, the transfer function is provided as follows:

Figure 5 .
Figure 5. Block diagram of IIR system identification.

Figure 6 .
Figure 6.Comparison between CGWO and the other algorithms for IIR system identification.

Figure 6 .
Figure 6.Comparison between CGWO and the other algorithms for IIR system identification.

Table 3 .
Initial parameters of algorithms (where T indicates the maximum number of interactions; t is the current iteration).

Table 4 .
Minimization results of the unimodal benchmark functions.

Table 5 .
Minimization results of the multimodal benchmark functions.
The unimodal benchmarks have only one global solution without any local optima for them.Consequently, they are useful in examining the convergence capability of heuristic optimization.The results of unimodal benchmark functions are shown in Table3.As shown in this table, for 30-D unimodal benchmark functions, CGWO has the best results for all benchmark functions.For 10-D, 50-D, 100-D, 200-D, 300-D, and 400-D, this table shows that the CGWO algorithm outperforms the other algorithms on the majority of the unimodal benchmark test functions.Therefore, the proposed algorithm has high performance in finding the global solution of unimodal benchmark test functions.According to the p values of 1 8 ´0.4z ´2 ´0.65z ´4 `0.26zH M pz ´1q " b 0 `b1 z ´1 `b2 z ´2 `b3 z ´3 `b4 z ´41 `a1 z ´1 `a2 z ´2 `a3 z ´3 `a4 z ´4

Table 14 .
The average f pθq value (AVE) and the standard deviation (STD).