# Theoretical and Empirical Analyses of an Improved Harmony Search Algorithm Based on Differential Mutation Operator

Authors: Longquan Yong; Sanyang Liu; Jianke Zhang; Quanxi Feng
Journal: Journal of Applied Mathematics
 (2012)

---

## Abstract

Harmony search (HS) method is an emerging metaheuristic optimization algorithm. In this paper, an improved harmony search method based on differential mutation operator (IHSDE) is proposed to deal with the optimization problems. Since the population diversity plays an important role in the behavior of evolution algorithm, the aim of this paper is to calculate the expected population mean and variance of IHSDE from theoretical viewpoint. Numerical results, compared with the HSDE, NGHS, show that the IHSDE method has good convergence property over a test-suite of well-known benchmark functions.


---

## Body

## 1. Introduction

Most optimization algorithms are based on numerical linear and nonlinear programming methods that require substantial gradient information and usually seek to improve the solution in the neighborhood of an initial point. These algorithms, however, reveal a limited approach to complicated real-world optimization problems because gradient is often difficult to find out or does not exist. What is more, if there is more than one local optimum in the problem, the result may depend on the selection of the starting point, and the obtained optimal solution may not necessarily be the global optimum.Recently, a new class of metaheuristics, named harmony search (HS), has been developed. The HS algorithm proposed in [1] has been developed in an analogy with music improvisation process where musicians in an ensemble continue to polish their pitches in order to obtain better harmony. Jazz improvisation seeks to find musically pleasing harmony similar to the optimum design process which seeks to find optimum solution. The pitch of each musical instrument determines the aesthetic quality, just as the objective function value is determined by the set of values assigned to each decision variable [2]. In addition, HS uses a stochastic random search instead of a gradient search so that derivative information is unnecessary.HS may be viewed as a simple real-coded genetic algorithm (GA), since it incorporates many important features of GA like mutation, recombination, and selection. HS has been successfully applied to a wide variety of practical optimization problems like designing controller [3], economic dispatch problem [4, 5], optimal power flow problem [6], neural networks [7], medical diagnosis [8], broadcast scheduling in packet radio networks [9], university course timetabling [10], and other engineering optimization field [11].Similar to the GA and particle swarm algorithms, the HS method is a random search technique. It does not require any prior domain knowledge, such as the gradient information of the objective functions. Unfortunately, empirical study has shown that the original HS method sometimes suffers from a slow search speed, and it is not suitable for handling the multimodal problems [2].Recently, Omran and Mahdavi tried to improve the performance of HS by incorporating some techniques from swarm intelligence. The new variant called by them as global best harmony search (GHS) [12] reportedly outperformed three other HS variants over the benchmark problems. Chakraborty et al. proposed an improved harmony search algorithm with differential mutation operator [13]. Gao et al. proposed modified harmony search methods for unimodal and multimodal optimization [14]. Wang et al. proposed self-adaptive harmony search algorithm for optimization [15]. Pan et al. proposed a self-adaptive global best harmony search algorithm for continuous optimization problems [16]. Zou et al. proposed a novel global harmony search algorithm (NGHS, [17]). More latest HS algorithm can be found in [18, 19].To overcome the shortcoming of premature convergence and stagnation, in this paper, we replace the pitch adjustment operation in classical HS (CHS) with a mutation strategy borrowed from the realm of the differential evolution (DE) algorithms, and we use(1.1)xinew=xij+F×(xir1-xir2),
where F is chosen to be a uniformly distributed random variable between 0.6 and 1.0, F~U[0.6,1] instead of  F~U[0,1] in [13], and j,r1,r2 are randomly selected with uniform distribution from the set {1,2,…,HMS},  j≠r1≠r2. The new mutation strategy is inspired by Chakraborty’s [13], and S. Das’s [19] work, especially in Chakraborty’s work [13], where the author theoretically showed that the harmony search based on differential mutation scheme (HSDE) has greater explorative power than the classical HS with pitch adjustment operation.The new algorithm proposed in this paper, called IHSDE (improved harmony search methods based on differential mutation operator), has been extensively compared with the HSDE, and the classical HS. Mathematical analysis will show that the IHSDE, under certain conditions, possesses an increasing population variance (with generation) as compared to HSDE. The numerical experiments show that the proposed algorithm is effective in dealing with a test suite of well-known benchmark functions.The rest of the paper is organized in the following way. Section2 briefly outlines the classical HS. Section 3 presents IHSDE method based on differential mutation operator and analyzes the expected population mean and variance of IHSDE in terms of theory. Effectiveness of IHSDE method is demonstrated in Section 4 by solving well-known benchmarks. Section 5 concludes the paper.


## 2. Classical Harmony Search Algorithm

### 2.1. Harmony Search Algorithm Principles

Current metaheuristic algorithms imitate natural phenomena, and evolution in evolutionary algorithms. HS algorithm was conceptualized using the musical process of searching for a perfect state of harmony. In music improvisation, each player sounds any pitch within the possible range, together making one harmony vector. If all the pitches make a good harmony, that experience is stored in each player’s memory, and the possibility to make a good harmony is increased next time. Similarly, in engineering optimization, each decision variable initially chooses any value within the possible range, together making one solution vector. If all the values of decision variables make a good solution, that experience is stored in each variable’s memory, and the possibility to make a good solution is also increased next time. Figure1 shows the details of the analogy between music improvisation and engineering optimization.Figure 1
Analogy between music improvisation and engineering optimization.The HS algorithm does not require initial values for the decision variables. Furthermore, instead of a gradient search, the HS algorithm uses a stochastic random search that is based on the harmony memory considering rate and the pitch-adjusting rate so that derivative information is unnecessary. Compared to earlier metaheuristic optimization algorithms, the HS algorithm imposes fewer mathematical requirements and can be easily adopted for various types of engineering optimization problems.


### 2.2. Steps for Classical Harmony Search Algorithm

The steps in the procedure of standard harmony search algorithm are as follows:Step  1 (Initialize the problem and algorithm parameters). The optimization problem is specified as follows. 								(2.1)Minimizef(x)subjecttoxi∈Xi,i=1,2,…,N,
where f(x) is an objective function; x is the set of each decision variable xi; N is the number of decision variables; Xi is the set of the possible range of values for each decision variable,  Xi:xiL≤Xi≤xiU. The HS algorithm parameters are also specified in this step. These are the harmony memory size (HMS), or the number of solution vectors in the harmony memory; the harmony memory (HM); harmony memory considering rate (HMCR); pitch-adjusting rate (PAR); the number of improvisations (Tmax).The harmony memory (HM) is a memory location where all the solution vectors (sets of decision variables) are stored. HMCR and PAR are parameters that are used to improve the solution vector.Step  2 (Initialize the harmony memory). In Step 2, the HM matrix is filled with as many randomly generated solution vectors as the HMS:
(2.2)HM=[x1x2⋮xHMS|f(x1)f(x2)⋮f(xHMS)]=[x11x21…xN1x12x22…xN2⋮⋮…x1HMSx2HMS…xNHMS|f(x1)f(x2)⋮f(xHMS)].Step  3 (Improvise a new harmony). A new harmony vector, x′=(x1′,x2′,…,xN′), is generated based on three rules: (a) memory consideration, (b) pitch adjustment, and (c) random selection. Generating a new harmony is called “improvisation.” In the memory consideration, the value of the first decision variable x1′ for the new vector is chosen from any of the values in the specified HM with a probability HMCR. The HMCR is the rate of choosing one value from the historical values stored in the HM, while (1-HMCR) is the rate of randomly selecting one value from the possible range of values:
(2.3)xi′={xi′∈(xi1,xi2,…,xiHMS),ifrand<HMCR,xi′∈Xi,otherwise,
where rand is a uniformly distributed random variable between 0 and 1.Every component obtained by the memory consideration is examined to determine whether it should be pitch adjusted. This operation uses the PAR parameter, which is the rate of pitch adjustment as follows.Pitch adjusting decision forxi′:
(2.4)xi′={xi′±rand×bw,ifrand<PAR,xi′,otherwise,
where bw is an arbitrary distance bandwidth.In Step 3, HM consideration, pitch adjustment, or random selection is applied to each variable of the new harmony vector in turn.Step  4 (Update harmony memory). If the new harmony vector x′=(x1′,x2′,…,xN′) is better than the worst harmony in the HM, the new harmony is included in the HM and the existing worst harmony is excluded from the HM.Step  5 (Check stopping criterion). If the stopping criterion (Tmax) is satisfied, computation is terminated. Otherwise, Steps 3 and 4 are repeated.In the next section, we employ the differential mutation operator to improve the fitness of all the members in the HS memory so that the overall convergence speed of the original HS method can be accelerated.


## 2.1. Harmony Search Algorithm Principles

Current metaheuristic algorithms imitate natural phenomena, and evolution in evolutionary algorithms. HS algorithm was conceptualized using the musical process of searching for a perfect state of harmony. In music improvisation, each player sounds any pitch within the possible range, together making one harmony vector. If all the pitches make a good harmony, that experience is stored in each player’s memory, and the possibility to make a good harmony is increased next time. Similarly, in engineering optimization, each decision variable initially chooses any value within the possible range, together making one solution vector. If all the values of decision variables make a good solution, that experience is stored in each variable’s memory, and the possibility to make a good solution is also increased next time. Figure1 shows the details of the analogy between music improvisation and engineering optimization.Figure 1
Analogy between music improvisation and engineering optimization.The HS algorithm does not require initial values for the decision variables. Furthermore, instead of a gradient search, the HS algorithm uses a stochastic random search that is based on the harmony memory considering rate and the pitch-adjusting rate so that derivative information is unnecessary. Compared to earlier metaheuristic optimization algorithms, the HS algorithm imposes fewer mathematical requirements and can be easily adopted for various types of engineering optimization problems.


## 2.2. Steps for Classical Harmony Search Algorithm

The steps in the procedure of standard harmony search algorithm are as follows:Step  1 (Initialize the problem and algorithm parameters). The optimization problem is specified as follows. 								(2.1)Minimizef(x)subjecttoxi∈Xi,i=1,2,…,N,
where f(x) is an objective function; x is the set of each decision variable xi; N is the number of decision variables; Xi is the set of the possible range of values for each decision variable,  Xi:xiL≤Xi≤xiU. The HS algorithm parameters are also specified in this step. These are the harmony memory size (HMS), or the number of solution vectors in the harmony memory; the harmony memory (HM); harmony memory considering rate (HMCR); pitch-adjusting rate (PAR); the number of improvisations (Tmax).The harmony memory (HM) is a memory location where all the solution vectors (sets of decision variables) are stored. HMCR and PAR are parameters that are used to improve the solution vector.Step  2 (Initialize the harmony memory). In Step 2, the HM matrix is filled with as many randomly generated solution vectors as the HMS:
(2.2)HM=[x1x2⋮xHMS|f(x1)f(x2)⋮f(xHMS)]=[x11x21…xN1x12x22…xN2⋮⋮…x1HMSx2HMS…xNHMS|f(x1)f(x2)⋮f(xHMS)].Step  3 (Improvise a new harmony). A new harmony vector, x′=(x1′,x2′,…,xN′), is generated based on three rules: (a) memory consideration, (b) pitch adjustment, and (c) random selection. Generating a new harmony is called “improvisation.” In the memory consideration, the value of the first decision variable x1′ for the new vector is chosen from any of the values in the specified HM with a probability HMCR. The HMCR is the rate of choosing one value from the historical values stored in the HM, while (1-HMCR) is the rate of randomly selecting one value from the possible range of values:
(2.3)xi′={xi′∈(xi1,xi2,…,xiHMS),ifrand<HMCR,xi′∈Xi,otherwise,
where rand is a uniformly distributed random variable between 0 and 1.Every component obtained by the memory consideration is examined to determine whether it should be pitch adjusted. This operation uses the PAR parameter, which is the rate of pitch adjustment as follows.Pitch adjusting decision forxi′:
(2.4)xi′={xi′±rand×bw,ifrand<PAR,xi′,otherwise,
where bw is an arbitrary distance bandwidth.In Step 3, HM consideration, pitch adjustment, or random selection is applied to each variable of the new harmony vector in turn.Step  4 (Update harmony memory). If the new harmony vector x′=(x1′,x2′,…,xN′) is better than the worst harmony in the HM, the new harmony is included in the HM and the existing worst harmony is excluded from the HM.Step  5 (Check stopping criterion). If the stopping criterion (Tmax) is satisfied, computation is terminated. Otherwise, Steps 3 and 4 are repeated.In the next section, we employ the differential mutation operator to improve the fitness of all the members in the HS memory so that the overall convergence speed of the original HS method can be accelerated.


## 3. The Improved HS Based on Differential Mutation Operator

Experiments with the CHS algorithm over the standard numerical benchmarks show that the algorithm suffers from the problem of premature and/or false convergence, slow convergence especially over multimodal fitness landscape.To circumvent these problems of premature, Chakraborty et al. proposed harmony search algorithm with differential mutation operator (HSDE). They replaced the pitch adjustment operation (2.4) in classical HS with a mutation strategy borrowed from the realm of the DE algorithms [20]. The mutation strategy has been presented as
(3.1)xinew=xij+F×(xir1-xir2),
where j,r1,r2∈U{1,2,…,HMS},  j≠r1≠r2, and the scale factor F is chosen to be a uniformly distributed random variable between 0 and 1, that is,  F~U[0,1].In what follows, we resetF is chosen to be a uniformly distributed random variable between 0.6 and 1.0,  F~U[0.6,1]. This improved harmony search algorithm with differential mutation operator will be referred to as IHSDE. In the following, we will show that the IHSDE gives better performance than the HSDE algorithm in theory and experiment.


### 3.1. IHSDE Algorithm

The pseudocode of IHSDE is described in Algorithm1. Every variable in the HM needs to go through the above DE mutation refinement procedure. Thus, the resulting HM members are expected to have better fitness than that of the original ones. This strategy can also overcome the premature shortcoming of the regular HS method.Algorithm 1:Pseudocode of the IHSDE Algorithm.
Procedure IHSDE algorithm
Initiate_parameters
Initialize_HM
While (not_termination)
For  i=1 to Ndo    //N denotes the number of decision variables
if ( rand < HMCR) //memory consideration
Select one harmony from HM randomly: xnew=xj,j∈U{1,2,…,HMS};
Execute difference variation operation for the selected harmony:
xinew=xij+F×(xir1-xir2), where F~U[0.6,1],r1,r2∈U{1,2,…,HMS},j≠r1≠r2.
else
xinew=xiL+rand×(xiU-xiL)//random selection
end if
End for
Update harmony memory HM // if applicable
End while
End procedure


### 3.2. Theoretical Analyses of the Expected Population Variance

Theoretical analyses of the properties of HS algorithms are very important to understand their search behaviors and to develop more efficient algorithms [13, 19]. Compared to the plethora of works concerning the empirical study of parameter selection and tuning process in HS [14–18], not much research has so far been devoted to theoretically analyze the search mechanism and convergence properties of HS, and this area remains largely open to prospective future research.The evolution of the expected population variance over generations provides a measure of the explorative power of the algorithm. In the following, we will estimate the expected mean and variance of the population obtained by applying mutation operator.Our ideas are as follows, firstly we find an analytical expression for the population expected variance, and then we compare the expected population variance of IHSDE with HSDE to show that the IHSDE algorithm possesses greater explorative power.In HS type algorithms, since each dimension is perturbed independently, without loss of generality, we can make our analysis for single-dimensional population members.Let us consider a population of scalarsx={x1,x2,…,xm} with elements xk∈R, k=1,2,…,m. The variance of this population is given by
(3.2)Var(x)=1m∑i=1m(xi-x-)2=x2¯-x-2,
where x-=(1/m)∑i=1mxi is population mean and x2¯=(1/m)∑i=1mxi2 is quadratic population mean.If the elements of the population are perturbed with some random numbers or variables,Var(x) will be a random variable, and E[Var(x)] will be a measure of the explorative power. In the following discussion, we always suppose  Var(x)>0.Lemma 3.1.
Letx={x1,x2,…,xm}  be the current population of HS, y be an intermediate vector obtained after random selection with harmony memory consideration, and z the vector obtained after y by replacing the pitch adjustment operation with a mutation operator (3.1) borrowed from DE algorithms in classical HS. Let w={w1,w2,…,wm} be the final population after selection. If we let
(3.3)p=harmonicmemoryconsiderationprobability(HMCR),
the allowable range for the new values of x is {xmin,xmax} where xmin=-a,xmax=a, and the required random numbers are continuously uniformly distributed random variable between 0 and 1 except for random variables F, then
(3.4)E[Var(w)]=2mF2¯Var(x)+m-13m2a2(1-p)+1m2[(m-1)p+(m-1)2]x2¯-m-1m(1mpx-+m-1mx-)2,
where F2¯=E(F2).Proof.
Sincex={x1,x2,…,xm}= HM is the current population of HS and y is an intermediate vector obtained after random selection with harmony memory consideration, that is,
(3.5)y={xl,withprobabilitypfromcurrentpopulation{x1,x2,…,xm},xr,withprobability(1-p)allowablerange{xmin,xmax},
where l∈{1,2,…,m} and xr is a new random variable in allowable range {xmin,xmax}.
So,(3.6)E(xl)=x-,E(xl2)=x2¯.
According to [13, 19], we have
(3.7)E(xr)=0,E(xr2)=a23.
Using (3.6), we can obtain the following relations:
(3.8)E(y)=E(xl)p+E(xr)(1-p)=px-;E(y2)=E(xl2)p+E(xr2)(1-p)=px2¯+a23(1-p).
Now,z is the vector obtained after mutating y. It has the following structure:
(3.9)z=y+F×(xr1-xr2).
LetF-=E(F),F2¯=E(F2), xr1 and xr2, are two randomly chosen members from x={x1,x2,…,xm} such that r1,r2∈{1,2,…,m},r1≠r2.
Thus,(3.10)E(xr1)=E(xr2)=x-;E(xr12)=E(xr22)=x2¯.
According to [13], we have
(3.11)E(xr1xr2)=1m(m-1)[(mx-)2-mx2¯].
Therefore,(3.12)E(z)=E(y)+E(F)×E(xr1-xr2)=E(y)+F-×[E(xr1)-E(xr2)]=E(y)=px-,E(z2)=E(y+F×(xr1-xr2))2=E[y2+F2×(xr12+xr22-2xr1xr2)+2F×y(xr1-xr2)]=E[y2]+F2¯×(E[xr12]+E[xr22]-2E[xr1xr2])+2F-×E[y](E[xr1]-E[xr2])=E[y2]+2F2¯mm-1(x2¯-x-2).
Now,w={w1,w2,…,wm} is the final population. Each element wk of the final population may be represented by the following:
(3.13)wk={z,withprobability1m,xk,withprobability(1-1m).
Thus,(3.14)E(wk)=E(z)1m+E(xk)(1-1m)=1mpx-+m-1mx-,E(wk2)=E(z2)1m+E(xk2)(1-1m).
Letw-=1/m∑k=1mwk and w2¯=1/m∑k=1mwk2 represent the population mean and quadratic population mean of the final target population, respectively. Then,
(3.15)E(w-2)=E(1m∑k=1mwk)2=1m2[mE(wk2)+m(m-1)E2(wk)];E(w2¯)=E(1m∑k=1mwk2)=1m(∑k=1mE(wk2))=1m(mE(wk2))=E(wk2).
Therefore, the variance of final populationw={w1,w2,…,wm} is given as
(3.16)Var(w)=w2¯-w-2,
and its expected population variance
(3.17)E[Var(w)]=E[w2¯-w-2]=E[w2¯]-E(w-2)=E(wk2)-1m2[mE(wk2)+m(m-1)E2(wk)]=m-1m[E(wk2)-E2(wk)]=m-1mE(wk2)-m-1mE2(wk)=m-1m[E(z2)1m+E(xk2)(1-1m)]-m-1m(1mpx-+m-1mx-)2=m-1m[E[y2]1m+2F2¯mm-1(x2¯-x-2)1m+x2¯(1-1m)]-m-1m(1mpx-+m-1mx-)2=m-1m[1mpx2¯+a23(1-p)1m+2F2¯mm-1(x2¯-x-2)1m+x2¯(1-1m)]-m-1m(1mpx-+m-1mx-)2=2mF2¯(x2¯-x-2)+a23(1-p)m-1m2+1m2[(m-1)p+(m-1)2]x2¯-m-1m(1mpx-+m-1mx-)2=2mF2¯Var(x)+m-13m2a2(1-p)+1m2[(m-1)p+(m-1)2]x2¯-m-1m(1mpx-+m-1mx-)2.Remark 3.2.
The conclusion in [13] (Theorem 1) is incorrect. In this paper, we give the correct expression for the population expected variance.The main analytical result is expressed in the form of the following theorem.Theorem 3.3.
Letx={x1,x2,…,xm}, be the current population of HS. Giving the same values m,p, here m and p represent the harmony memory size (HMS) and harmony memory considering rate (HMCR), respectively. The intergeneration population expected variance of IHSDE algorithm is greater than HSDE.Proof.
In HSDE,F is chosen to be a uniformly distributed random variable between 0 and 1, and in IHSDE, F is chosen to be a uniformly distributed random variable between 0.6 and 1. For convenience, let F1~U[0,1], and the intergeneration expected variance of HSDE is  EHSDE[Var(w)]. Let F2~U[0.6,1], and the intergeneration expected variance of IHSDE is  EIHSDE[Var(w)].
Since(3.18)F12¯=E(F12)=F-12+Var(F1)=0.52+(1-0)212=13=0.3333,F22¯=E(F22)=F-22+Var(F2)=0.82+(1-0.6)212=0.6533.
Thus,(3.19)EIHSDE[Var(w)]-EHSDE[Var(w)]=2m(F22¯-F12¯)Var(x)>0.Remark 3.4.
Improper parametersF can lead to problem of premature and/or false convergence. Usually, F∈[0,2]. If F is chosen too small, it gets more difficult to escape local optima. A larger F (or F2¯) increases the probability (or the population expected variance of IHSDE) for escaping a local optimum. However, for F>1, the convergence speed decreases. It is more difficult to converge for a population when the perturbation is larger than the distance between two members. Thus, though standard Gaussian distribution N(0,1) has  much larger F2¯=E(F2)=1 than that of F2~U[0.6,1], it cannot get better convergence. In [21], Mperle et al. suggest that a good initial choice for the amplification factor should be F=0.6. If one suspects that with this setting only a local optimum is found, then F should be increased. This fully shows that our choosing F2~U[0.6,1] is appropriate.The above mathematical analysis show that the IHSDE possesses an increasing population variance as compared to HSDE. This ensures that the explorative power of IHSDE is on average greater than that of HSDE, which in turn results into better accuracy of the IHSDE algorithm.In the following section, we give some numerical experiments over standard test functions.


## 3.1. IHSDE Algorithm

The pseudocode of IHSDE is described in Algorithm1. Every variable in the HM needs to go through the above DE mutation refinement procedure. Thus, the resulting HM members are expected to have better fitness than that of the original ones. This strategy can also overcome the premature shortcoming of the regular HS method.Algorithm 1:Pseudocode of the IHSDE Algorithm.
Procedure IHSDE algorithm
Initiate_parameters
Initialize_HM
While (not_termination)
For  i=1 to Ndo    //N denotes the number of decision variables
if ( rand < HMCR) //memory consideration
Select one harmony from HM randomly: xnew=xj,j∈U{1,2,…,HMS};
Execute difference variation operation for the selected harmony:
xinew=xij+F×(xir1-xir2), where F~U[0.6,1],r1,r2∈U{1,2,…,HMS},j≠r1≠r2.
else
xinew=xiL+rand×(xiU-xiL)//random selection
end if
End for
Update harmony memory HM // if applicable
End while
End procedure


## 3.2. Theoretical Analyses of the Expected Population Variance

Theoretical analyses of the properties of HS algorithms are very important to understand their search behaviors and to develop more efficient algorithms [13, 19]. Compared to the plethora of works concerning the empirical study of parameter selection and tuning process in HS [14–18], not much research has so far been devoted to theoretically analyze the search mechanism and convergence properties of HS, and this area remains largely open to prospective future research.The evolution of the expected population variance over generations provides a measure of the explorative power of the algorithm. In the following, we will estimate the expected mean and variance of the population obtained by applying mutation operator.Our ideas are as follows, firstly we find an analytical expression for the population expected variance, and then we compare the expected population variance of IHSDE with HSDE to show that the IHSDE algorithm possesses greater explorative power.In HS type algorithms, since each dimension is perturbed independently, without loss of generality, we can make our analysis for single-dimensional population members.Let us consider a population of scalarsx={x1,x2,…,xm} with elements xk∈R, k=1,2,…,m. The variance of this population is given by
(3.2)Var(x)=1m∑i=1m(xi-x-)2=x2¯-x-2,
where x-=(1/m)∑i=1mxi is population mean and x2¯=(1/m)∑i=1mxi2 is quadratic population mean.If the elements of the population are perturbed with some random numbers or variables,Var(x) will be a random variable, and E[Var(x)] will be a measure of the explorative power. In the following discussion, we always suppose  Var(x)>0.Lemma 3.1.
Letx={x1,x2,…,xm}  be the current population of HS, y be an intermediate vector obtained after random selection with harmony memory consideration, and z the vector obtained after y by replacing the pitch adjustment operation with a mutation operator (3.1) borrowed from DE algorithms in classical HS. Let w={w1,w2,…,wm} be the final population after selection. If we let
(3.3)p=harmonicmemoryconsiderationprobability(HMCR),
the allowable range for the new values of x is {xmin,xmax} where xmin=-a,xmax=a, and the required random numbers are continuously uniformly distributed random variable between 0 and 1 except for random variables F, then
(3.4)E[Var(w)]=2mF2¯Var(x)+m-13m2a2(1-p)+1m2[(m-1)p+(m-1)2]x2¯-m-1m(1mpx-+m-1mx-)2,
where F2¯=E(F2).Proof.
Sincex={x1,x2,…,xm}= HM is the current population of HS and y is an intermediate vector obtained after random selection with harmony memory consideration, that is,
(3.5)y={xl,withprobabilitypfromcurrentpopulation{x1,x2,…,xm},xr,withprobability(1-p)allowablerange{xmin,xmax},
where l∈{1,2,…,m} and xr is a new random variable in allowable range {xmin,xmax}.
So,(3.6)E(xl)=x-,E(xl2)=x2¯.
According to [13, 19], we have
(3.7)E(xr)=0,E(xr2)=a23.
Using (3.6), we can obtain the following relations:
(3.8)E(y)=E(xl)p+E(xr)(1-p)=px-;E(y2)=E(xl2)p+E(xr2)(1-p)=px2¯+a23(1-p).
Now,z is the vector obtained after mutating y. It has the following structure:
(3.9)z=y+F×(xr1-xr2).
LetF-=E(F),F2¯=E(F2), xr1 and xr2, are two randomly chosen members from x={x1,x2,…,xm} such that r1,r2∈{1,2,…,m},r1≠r2.
Thus,(3.10)E(xr1)=E(xr2)=x-;E(xr12)=E(xr22)=x2¯.
According to [13], we have
(3.11)E(xr1xr2)=1m(m-1)[(mx-)2-mx2¯].
Therefore,(3.12)E(z)=E(y)+E(F)×E(xr1-xr2)=E(y)+F-×[E(xr1)-E(xr2)]=E(y)=px-,E(z2)=E(y+F×(xr1-xr2))2=E[y2+F2×(xr12+xr22-2xr1xr2)+2F×y(xr1-xr2)]=E[y2]+F2¯×(E[xr12]+E[xr22]-2E[xr1xr2])+2F-×E[y](E[xr1]-E[xr2])=E[y2]+2F2¯mm-1(x2¯-x-2).
Now,w={w1,w2,…,wm} is the final population. Each element wk of the final population may be represented by the following:
(3.13)wk={z,withprobability1m,xk,withprobability(1-1m).
Thus,(3.14)E(wk)=E(z)1m+E(xk)(1-1m)=1mpx-+m-1mx-,E(wk2)=E(z2)1m+E(xk2)(1-1m).
Letw-=1/m∑k=1mwk and w2¯=1/m∑k=1mwk2 represent the population mean and quadratic population mean of the final target population, respectively. Then,
(3.15)E(w-2)=E(1m∑k=1mwk)2=1m2[mE(wk2)+m(m-1)E2(wk)];E(w2¯)=E(1m∑k=1mwk2)=1m(∑k=1mE(wk2))=1m(mE(wk2))=E(wk2).
Therefore, the variance of final populationw={w1,w2,…,wm} is given as
(3.16)Var(w)=w2¯-w-2,
and its expected population variance
(3.17)E[Var(w)]=E[w2¯-w-2]=E[w2¯]-E(w-2)=E(wk2)-1m2[mE(wk2)+m(m-1)E2(wk)]=m-1m[E(wk2)-E2(wk)]=m-1mE(wk2)-m-1mE2(wk)=m-1m[E(z2)1m+E(xk2)(1-1m)]-m-1m(1mpx-+m-1mx-)2=m-1m[E[y2]1m+2F2¯mm-1(x2¯-x-2)1m+x2¯(1-1m)]-m-1m(1mpx-+m-1mx-)2=m-1m[1mpx2¯+a23(1-p)1m+2F2¯mm-1(x2¯-x-2)1m+x2¯(1-1m)]-m-1m(1mpx-+m-1mx-)2=2mF2¯(x2¯-x-2)+a23(1-p)m-1m2+1m2[(m-1)p+(m-1)2]x2¯-m-1m(1mpx-+m-1mx-)2=2mF2¯Var(x)+m-13m2a2(1-p)+1m2[(m-1)p+(m-1)2]x2¯-m-1m(1mpx-+m-1mx-)2.Remark 3.2.
The conclusion in [13] (Theorem 1) is incorrect. In this paper, we give the correct expression for the population expected variance.The main analytical result is expressed in the form of the following theorem.Theorem 3.3.
Letx={x1,x2,…,xm}, be the current population of HS. Giving the same values m,p, here m and p represent the harmony memory size (HMS) and harmony memory considering rate (HMCR), respectively. The intergeneration population expected variance of IHSDE algorithm is greater than HSDE.Proof.
In HSDE,F is chosen to be a uniformly distributed random variable between 0 and 1, and in IHSDE, F is chosen to be a uniformly distributed random variable between 0.6 and 1. For convenience, let F1~U[0,1], and the intergeneration expected variance of HSDE is  EHSDE[Var(w)]. Let F2~U[0.6,1], and the intergeneration expected variance of IHSDE is  EIHSDE[Var(w)].
Since(3.18)F12¯=E(F12)=F-12+Var(F1)=0.52+(1-0)212=13=0.3333,F22¯=E(F22)=F-22+Var(F2)=0.82+(1-0.6)212=0.6533.
Thus,(3.19)EIHSDE[Var(w)]-EHSDE[Var(w)]=2m(F22¯-F12¯)Var(x)>0.Remark 3.4.
Improper parametersF can lead to problem of premature and/or false convergence. Usually, F∈[0,2]. If F is chosen too small, it gets more difficult to escape local optima. A larger F (or F2¯) increases the probability (or the population expected variance of IHSDE) for escaping a local optimum. However, for F>1, the convergence speed decreases. It is more difficult to converge for a population when the perturbation is larger than the distance between two members. Thus, though standard Gaussian distribution N(0,1) has  much larger F2¯=E(F2)=1 than that of F2~U[0.6,1], it cannot get better convergence. In [21], Mperle et al. suggest that a good initial choice for the amplification factor should be F=0.6. If one suspects that with this setting only a local optimum is found, then F should be increased. This fully shows that our choosing F2~U[0.6,1] is appropriate.The above mathematical analysis show that the IHSDE possesses an increasing population variance as compared to HSDE. This ensures that the explorative power of IHSDE is on average greater than that of HSDE, which in turn results into better accuracy of the IHSDE algorithm.In the following section, we give some numerical experiments over standard test functions.


## 4. Computational Results

The effectiveness of the IHSDE algorithm has been evaluated on a test suite of well-known benchmark functions (Table1) [14–16, 22, 23]. In Table 1, n represents the number of dimensions. Here, n=5, except for function F12~F16, where n=2. The global minima of all the above functions are at F(x*)=0, except for F4, F10, F13 and F14. Table 1 summarizes the initialization and search ranges used for these functions.Table 1
Benchmark test functions.
Function name
Benchmark functions expression
Search range
Optimum value
F
1 Ackley function
F
1
(
x
)
=
-
20
e
(
1
/
5
)
(
1
/
n
)
∑
i
=
1
n
x
i
2
-
e
(
1
/
n
)
∑
i
=
1
n
cos
(
2
π
x
i
)
+
20
-
e
[
-
15,30
]
n
x
*
=
(
0,0
,
…
,
0
), F1(x*)=0
F
2 Dixon and Price function
F
2
(
x
)
=
(
x
1
+
1
)
2
+
∑
i
=
2
n
i
(
2
x
i
2
-
x
i
-
1
)
2
[
-
10,10
]
n
F
2
(
x
*
)
=
0
F
3 Levy function
F
3
(
x
)
=
sin
2
(
π
y
1
)
+
∑
i
=
1
n
-
1
[
(
y
i
-
1
)
2
(
1
+
10
sin
2
(
π
y
i
+
1
)
)
]
+
(
y
n
-
1
)
2
(
1
+
10
sin
2
(
2
π
y
n
)
)
y
i
=
1
+
x
i
-
1
4,i=1,2,…,n
[
-
10,10
]
n
x
*
=
(
1,1
,
…
,
1
), F3(x*)=0
F
4 Michalewics function
F
4
(
x
1
,
x
2
)
=
-
∑
j
=
1
2
sin
(
x
j
)
(
sin
(
j
x
j
2
/
π
)
)
2
m
;
m
=
10
[
0
,
π
]
n
n
=
5
;
F
4
(
x
*
)
=
-
4.687658
F
5 Perm function
F
5
(
x
)
=
∑
k
=
1
n
[
∑
i
=
1
n
(
i
k
+
β
)
(
(
x
i
/
k
)
k
-
1
)
]
2
[
-
n
,
n
]
x
*
=
(
1,2
,
…
,
n
), F5(x*)=0
F
6 Powell function
F
6
(
x
)
=
∑
i
=
1
n
/
4
[
(
x
4
i
-
3
+
10
x
4
i
-
2
)
2
+
5
(
x
4
i
-
1
-
x
4
i
)
2
+
(
x
4
i
-
2
-
x
4
i
-
1
)
2
+
10
(
x
4
i
-
3
-
x
4
i
)
4
]
[
-
4,5
]
n
x
*
=
(
3
,
-
1,0
,
1
,
…
,
3
,
-
1,0
,
1
), F6(x*)=0
F
7 Rastrigin function
F
7
(
x
)
=
10
n
+
∑
i
=
1
n
(
x
i
2
-
10
cos
(
2
π
x
i
)
)
[
-
5.12,5.12
]
n
x
*
=
(
0,0
,
…
,
0
), F7(x*)=0
F
8 Rosenbrock function
F
8
(
x
)
=
∑
i
=
1
n
-
1
[
100
(
x
i
2
-
x
i
-
1
)
2
+
(
x
i
-
1
)
2
]
[
-
5,10
]
n
x
*
=
(
1,1
,
…
,
1
), F8(x*)=0
F
9 Sphere function
F
9
(
x
)
=
∑
i
=
1
n
x
i
2
[
-
5.12,5.12
]
n
x
*
=
(
0,0
,
…
,
0
), F9(x*)=0
F
10
Trid function
F
10
(
x
)
=
∑
i
=
1
n
(
x
i
-
1
)
2
-
∑
i
=
2
n
x
i
x
i
-
1
[
-
n
2
,
n
2
]
n
n
=
5
;
F
10
(
x
*
)
=
-
30
F
11 Zakharov function
F
11
(
x
)
=
∑
i
=
1
n
x
i
2
+
(
∑
i
=
1
n
0.5
i
x
i
)
2
+
(
∑
i
=
1
n
0.5
i
x
i
)
4
[
-
5,10
]
n
x
*
=
(
0,0
,
…
,
0
), F11(x*)=0
F
12 Beale function
F
12
(
x
)
=
(
1.5
-
x
1
+
x
1
x
2
)
2
+
(
2.25
-
x
1
+
x
1
x
2
2
)
2
+
(
2.625
-
x
1
+
x
1
x
2
3
)
2
[
-
4.5,4.5
]
2
x
*
=
(
3,0.5
), F12(x*)=0
F
13 Easom function
F
13
(
x
)
=
-
cos
(
x
1
)
cos
(
x
2
)
exp
(
-
(
x
1
-
π
)
2
-
(
x
2
-
π
)
2
)
[
-
100,100
]
2
x
*
=
(
π
,
π
), F13(x*)=-1
F
14 Goldstein and Price function
F
14
(
x
)
=
[
1
+
(
x
1
+
x
2
+
1
)
2
(
19
-
14
x
1
+
3
x
1
2
-
14
x
2
+
6
x
1
x
2
+
3
x
2
2
)
]
×[30+(2x1-3x2)2(18-32x1+12x12+48x2-36x1x2+27x22)]
[
-
2,2
]
2
x
*
=
(
0
,
-
1
), F14(x*)=3
F
15  Hump function
F
15
(
x
)
=
4
x
1
2
-
2.1
x
1
4
+
(
1
/
3
)
x
1
6
+
x
1
x
2
-
4
x
2
2
+
4
x
2
4
[
-
5,5
]
2
x
*
=
(
0.0898
,
-
0.7126
) orx*=(-0.0898,0.7126),F15(x*)=0
F
16 My function
F
16
(
x
)
=
(
x
1
2
+
x
2
-
11
)
2
+
(
x
1
+
x
2
2
-
7
)
2
+
0.1
(
(
x
1
-
3
)
2
+
(
x
2
-
2
)
2
)
[
-
6,6
]
2
x
*
=
(
3,2
),F16(x*)=0All the experiments were performed on Windows XP 64 System running on an Hp desktop with Intel(R) Xeon(R)4×2.4 GHz and 6 GB RAM, and the codes were written in MATLAB R2008a. Simulations were carried out to compare the optimization (minimization) capabilities of the proposed method (IHSDE) with respect to (a) HSDE [13], (b) classical HS (CHS) [1, 2], (c) NGHS [17]. To make the comparison fair, the populations for all the competitor algorithms (for all problems tested) were initialized using the same random seeds. The HS-variant algorithm parameters were set the same parameters: harmony memory size HMS=10, harmony memory consideration rate HMCR=0.8, and the number of improvisations Tmax=10000. In HSDE, F is chosen to be a uniformly distributed random variable between 0 and 1, and in IHSDE F is chosen to be a uniformly distributed random variable between 0.6 and 1, F=0.6+0.4*rand, here  rand is a uniformly distributed random variable between 0 and 1. In classical HS, we set pitch-adjusting rate PAR=0.4.To judge the accuracy of different algorithms, 50 independent runs of each of the four algorithms were carried out and the best, the mean, the worst fitness values, and the standard deviation (Std) were recorded. Table2 compares the algorithms on the quality of the optimum solution.Table 2
The statistical results for 50 runs tested on sixteen benchmark functions.
Function name
Algorithm
Best
Mean
Worst
Std
F
1
Ackley function
CHS
2.8576
E
-
02
1.5581
E
-
01
3.9193
E
-
01
7.6256
E
-
02
HSDE
8.8818
E
-
16
2.5301
E
-
03
2.8793
E
-
02
6.2847
E
-
03
IHSDE
9.0683
E
-
13
6.7733
E
-
11
1.0303
E
-
09
1.5937
E
-
10
NGHS
4.4409
E
-
15
6.3594
E
-
15
1.5099
E
-
14
2.6087
E
-
15
F
2
Dixon and price function
CHS
8.8618
E
-
06
1.2597
E
-
01
6.1434
E
-
01
1.5982
E
-
01
HSDE
2.5471
E
-
08
2.3591
E
-
01
6.6820
E
-
01
2.7753
E
-
01
IHSDE
1.6000
E
-
13
1.8755
E
-
01
6.6667
E
-
01
3.0089
E
-
01
NGHS
5.0758
E
-
18
1.8658
E
-
01
6.6667
E
-
01
3.0222
E
-
01
F
3
Levy function
CHS
6.8353
E
-
07
3.2798
E
-
04
2.2443
E
-
03
4.5606
E
-
04
HSDE
1.4998
E
-
32
1.0481
E
-
06
2.3110
E
-
05
3.8915
E
-
06
IHSDE
2.2192
E
-
27
8.1646
E
-
22
1.2126
E
-
20
2.2009
E
-
21
NGHS
1.4998
E
-
32
1.4998
E
-
32
1.4998
E
-
32
1.3823
E
-
47
F
4
Michalewics function
CHS
-
4.6877
E
+
00
-
4.6859
E
+
00
-
4.6745
E
+
00
2.4486
E
-
03
HSDE
-
4.6877
E
+
00
-
4.6873
E
+
00
-
4.6778
E
+
00
1.4840
E
-
03
IHSDE
-
4.6877
E
+
00
-
4.6877
E
+
00
-
4.6876
E
+
00
3.1708
E
-
06
NGHS
-
4.6877
E
+
00
-
4.6615
E
+
00
-
4.4959
E
+
00
4.4690
E
-
02
F
5
Perm function
CHS
3.8305
E
-
01
3.0003
E
+
02
2.7696
E
+
03
5.6432
E
+
02
HSDE
1.1911
E
-
02
1.3767
E
+
02
4.9112
E
+
03
6.9678
E
+
02
IHSDE
7.7908
E
-
02
1.1573
E
+
02
1.1205
E
+
03
1.8330
E
+
02
NGHS
2.5956
E
-
02
3.0299
E
+
01
3.0618
E
+
02
6.7526
E
+
01
F
6
Powell function
CHS
4.9848
E
-
06
1.5997
E
-
02
1.1282
E
-
01
1.8387
E
-
02
HSDE
2.3029
E
-
07
5.6853
E
-
03
3.3194
E
-
02
7.9504
E
-
03
IHSDE
2.7878
E
-
07
9.6453
E
-
06
5.0805
E
-
05
1.2526
E
-
05
NGHS
3.4739
E
-
08
1.3536
E
-
06
3.0482
E
-
06
6.9536
E
-
07
F
7
Rastrigin function
CHS
1.2141
E
-
03
3.8574
E
-
02
3.1431
E
-
01
5.0927
E
-
02
HSDE
0.0000
E
+
00
4.0904
E
-
03
8.3710
E
-
02
1.3816
E
-
02
IHSDE
4.2633
E
-
14
2.5944
E
+
00
8.1800
E
+
00
2.6713
E
+
00
NGHS
0.0000
E
+
00
9.9496
E
-
02
9.9496
E
-
01
3.0152
E
-
01
F
8
Rosenbrock function
CHS
4.3234
E
-
03
1.7703
E
+
00
7.2601
E
+
00
1.9667
E
+
00
HSDE
2.2513
E
-
04
1.9541
E
+
00
3.3526
E
+
00
1.0636
E
+
00
IHSDE
9.5471
E
-
03
1.1625
E
+
00
2.1110
E
+
00
3.9876
E
-
01
NGHS
1.0784
E
-
03
9.3133
E
-
01
3.5246
E
+
00
6.5434
E
-
01
F
9
Sphere function
CHS
3.9500
E
-
06
2.2513
E
-
04
1.3357
E
-
03
2.6271
E
-
04
HSDE
1.2298
E
-
48
1.1425
E
-
06
1.7339
E
-
05
3.5860
E
-
06
IHSDE
4.9171
E
-
28
2.5756
E
-
23
5.8627
E
-
22
9.2448
E
-
23
NGHS
8.0887
E
-
67
1.2697
E
-
54
6.1078
E
-
53
8.6366
E
-
54
F
10
Trid Function
CHS
-
3.0000
E
+
01
-
2.9875
E
+
01
-
2.9156
E
+
01
2.3411
E
-
01
HSDE
-
3.0000
E
+
01
-
2.9962
E
+
01
-
2.9826
E
+
01
4.7768
E
-
02
IHSDE
-
3.0000
E
+
01
-
3.0000
E
+
01
-
3.0000
E
+
01
1.4565
E
-
10
NGHS
-
3.0000
E
+
01
-
3.0000
E
+
01
-
3.0000
E
+
01
9.1416
E
-
08
F
11
Zakharov Function
CHS
2.3830
E
-
05
1.9843
E
-
02
1.4567
E
-
01
3.0974
E
-
02
HSDE
2.2372
E
-
22
4.1294
E
-
03
1.2299
E
-
01
1.8304
E
-
02
IHSDE
2.1141
E
-
15
6.1236
E
-
10
2.9812
E
-
08
4.2138
E
-
09
NGHS
2.5239
E
-
29
1.0790
E
-
18
3.7173
E
-
17
5.3803
E
-
18
F
12
Beale function
CHS
4.3647
E
-
10
5.9406
E
-
03
9.3633
E
-
02
1.7594
E
-
02
HSDE
1.4369
E
-
12
3.6930
E
-
03
3.9410
E
-
02
6.6767
E
-
03
IHSDE
0.0000
E
+
00
9.5965
E
-
05
3.5321
E
-
03
5.0902
E
-
04
NGHS
0.0000
E
+
00
4.3172
E
-
01
9.3929
E
+
00
1.3416
E
+
00
F
13
Easom function
CHS
-
9.9999
E
-
01
-
6.3682
E
-
01
-
7.6526
E
-
05
4.7875
E
-
01
HSDE
-
1.0000
E
+
00
-
9.5984
E
-
01
-
7.9071
E
-
05
1.9790
E
-
01
IHSDE
-
1.0000
E
+
00
-
1.0000
E
+
00
-
1.0000
E
+
00
0.0000
E
+
00
NGHS
-
1.0000
E
+
00
-
7.8001
E
-
01
0.0000
E
+
00
4.1842
E
-
01
F
14
Goldstein and Price function
CHS
3.0000
E
+
00
3.0016
E
+
00
3.0277
E
+
00
4.1190
E
-
03
HSDE
3.0000
E
+
00
3.0003
E
+
00
3.0063
E
+
00
1.1281
E
-
03
IHSDE
3.0000
E
+
00
3.0000
E
+
00
3.0000
E
+
00
1.6324
E
-
07
NGHS
3.0000
E
+
00
1.5420
E
+
01
8.4000
E
+
01
2.5122
E
+
01
F
15
Hump function
CHS
4.6668
E
-
08
2.8159
E
-
05
3.5591
E
-
04
5.9106
E
-
05
HSDE
4.6510
E
-
08
1.7936
E
-
06
4.2616
E
-
05
8.1164
E
-
06
IHSDE
4.6510
E
-
08
3.3050
E
-
07
1.3125
E
-
05
1.8532
E
-
06
NGHS
4.6510
E
-
08
4.6510
E
-
08
4.6510
E
-
08
1.0318
E
-
16
F
16
My function
CHS
1.3297
E
-
08
4.7659
E
-
01
3.4872
E
+
00
7.9292
E
-
01
HSDE
0.0000
E
+
00
1.5107
E
-
01
1.5166
E
+
00
4.5650
E
-
01
IHSDE
0.0000
E
+
00
1.4701
E
-
07
7.1040
E
-
06
1.0045
E
-
06
NGHS
0.0000
E
+
00
2.1788
E
+
00
7.3673
E
+
00
2.5506
E
+
00Figure2 shows the convergence and its boxplot of the best fitness in the population for the different algorithms (CHS, HSDE, IHSDE, NGHS) for all benchmark functions. The values plotted for every generation are averaged over 50 independent runs.  The boxplot is the best fitness in the final population for the different algorithms (CHS, HSDE, IHSDE, NGHS) over 50 independent runs. We can see that the behavior of the two former algorithms (CHS and HSDE) is similar for all benchmark functions. From the graphical point of view, the classical HS algorithm is clearly the worst algorithm for most benchmark problems, while the IHSDE outperforms HSDE in most cases except for F2 and F7. Compared with the two former algorithms (CHS, HSDE), for most test functions, IHSDE can achieve much better optimization results within the same iterations, and IHSDE can reach or exceed the advanced level of the NGHS in most cases.The convergence and its boxplot.
(a)
(b)
(c)
(d)In order to accurately give search process of each algorithm, we set the same parameters, and run CHS, HSDE, IHSDE method, respectively, for multimodal function F16. Figure 3  shows iterations for three methods.Show search processes of CHS, HSDE, and IHSDE methods forF16.
(a)
Contour line and Optimal Solution
(b)
Iterations of CHS Method
(c)
Iterations of HSDE Method
(d)
Iterations of IHSDE MethodRemark 4.1.
We omitted plots for all the other functions (F12-F15) to save space and also in consideration of the fact that they display more or less the same trend.
Experimental results on benchmark functions show that the IHSDE method can outperform the other methods. From Figure3, we confirm that the IHSDE method has better ability of global search, and not easily fall in local optimum.


## 5. Conclusion

This paper has presented an improved harmony search algorithm by blending with it a different vector-based mutation operator borrowed from the DE algorithms. The HM members are fine tuned by the DE’s mutation operator to improve their affinities so that enhanced optimization performances can be achieved. Mathematical analysis indicates that the IHSDE posses an increasing population variance as compared to HSDE. This ensures that the explorative power of IHSDE is on average greater than that of HSDE, which in turn results in better accuracy of the IHSDE algorithm. Several simulation examples of the unimodal and multimodal functions have been used to verify the effectiveness of the proposed methods. Compared with the HSDE and CHS, better optimization results are obtained using IHSDE approaches in most cases. Checking the effect of variation of the scale factorF of the differential mutation operator may be a worthy issue for preventing the premature convergence in future investigations. Future works will also focus on studying the applications of IHSDE on engineering optimization problems.



---

*Source: 147950-2012-10-01.xml*