A hybrid optimization algorithm based on harmony search and differential evolution for continuous domain

Abstract

Harmony Search (HS), as a novel optimization algorithm, is interesting for researchers because of its simplicity and appropriate efficiency. This algorithm has been applied to different continuous optimization problems leading to successful results. The parameter setting significantly impacts the performance of Harmony Search, similar to other meta-heuristics. In order to enhance this algorithm and facilitate the task of parameter setting, a hybrid approach is introduced in this paper. The harmony algorithm is improved by combining with the Differential Evolution (DE), as a successful evolutionary optimization algorithm which have been yielded admissible results in solving difficult problems. The new algorithm is tested on some well-known benchmark functions. The results show that the new hybrid algorithm has superior performance in comparison to the original Harmony Search and a modified version which is called Improved Harmony Search.

Keywords

Continuous optimization harmony search differential evolution hybrid algorithms

1 Introduction

Nowadays, there is an increasing interest in the study of meta-heuristic algorithms, including evolutionary algorithms. Recently, Geem et al. proposed a new evolutionary algorithm called Harmony Search (HS) algorithm for continuous optimization problems which imitate the process of music improvisation by a musician [1, 2]. The structure of algorithm is simple and it is successfully applied to different optimization problems including function optimization, engineering optimization, NP-complete problems, vehicle routing, model checking, design of water distribution networks, groundwater modeling and structural design [2 , 14–19].

Similar to the most meta-heuristics algorithms, the effectiveness of HS is sensitive to parameter values. Finding appropriate values for different parameters is crucial in the application of these algorithms. Some works have been carried out to enhance the performance of HS by tuning its parameters or combining it with other search algorithms [4]. In this paper, a new hybrid version of HS is presented which combined with Differential Evolution (DE) as a powerful evolutionary optimization technique. DE is widely used in solving different continuous optimization problems. The proposed algorithm has been evaluated using eight well-known benchmark problems in continuous optimization. The results show that the performance of the new hybrid algorithm is improved in comparison to the original HS and another hybrid version.

The rest of this paper is organized as follows. In Section 2, the harmony search is briefly described. Another version of HS, namely Improved Harmony Search algorithm is presented in Section 3. The Differential Evolution algorithm is explained in the next section. In Section 5, the new algorithm is described and the experimental results is presented in Section 6. The concluding remarks are presented in the last section.

2 Harmony search algorithm

The steps of harmony search for the generalized orienteering problem are as follows: [4 –6].

Step 1. Initialize the parameters for problem and algorithm.

Step 2. Initialize the harmony memory (HM).

Step 3. Improvise new harmony.

Step 4. Update the harmony memory.

Step 5. Check the stopping criterion.

2.1 Parameters initialization

In Step 1, the optimization problem is specified as follows: $Minimize F (R)$ (1) $Subject to R_{i}, i = 1, 2, 3, \dots, N$ (2) where, F (R) is an objective function, R is the set of decision variable R _i and N is the number of decision variables. The HS algorithm parameters are also specified in this step: harmony memory size (HMS, number of solution vectors), harmony memory considering rate (HMCR), pitch adjusting rate (PAR) and the stopping criteria (number of improvisation). HMCR and PAR are the parameters of HS algorithm will be explained in Step 3.

2.2 Initialize harmony memory

In Step 2, the harmony memory (HM) matrix, as shown in Equation 3, is filled with as many randomly generated solution vectors as the size of the HM (HMS). $HM = [\begin{matrix} R_{1}^{1} & R_{2}^{1} & \dots & R_{N}^{1} \\ ⋮ & ⋮ & \dots & ⋮ \\ R_{1}^{HMS - 1} & R_{2}^{HMS - 1} & \dots & R_{N}^{HMS - 1} \\ R_{1}^{HMS} & R_{2}^{HMS} & \dots & R_{N}^{HMS} \end{matrix} | \begin{matrix} F (R^{1}) \\ ⋮ \\ F (R^{HMS - 1}) \\ F (R^{HMS}) \end{matrix}]$ (3)

2.3 Improvise new harmony

A new harmony vector, $R^{new} = (R_{1}^{new}, R_{2}^{new}, \dots, R_{N}^{new})$ is generated by the following three rules: HM consideration, pitch adjustment, pitch adjustment and totally random generation. For instance, the value of the first decision variable ( $R_{1}^{new}$ ) for the new vector can be chosen from values stored in HM ( $R_{1}^{1} \sim R_{1}^{HMS}$ ). Value of other variables ( $R_{i}^{new}$ ) can be chosen in the same style. There is also a possibility for totally random value to be chosen. HMCR parameter, which varies between 0 and 1, defines the rate whether a value stored in HM is chosen or a random value is chosen, as follows: $R_{1}^{new} \leftarrow {\begin{matrix} R_{1}^{new} \in {R_{1}^{1}, . . ., R_{1}^{HMS}} \\ \begin{matrix} with & probability & HMCR \end{matrix} \\ R_{1}^{new} \in R (t) \\ with probability (1 - HMCR) \end{matrix}$ (4)

The HMCR is the rate of choosing one value from historical values stored in HM while 1-HMCR is the rate of randomly choosing one value from the possible value range.

After choosing the new harmony vector $R^{new} = (R_{1}^{new}, R_{2}^{new}, \dots, R_{N}^{new})$ , pitch-adjusting decision is examined for each component of the new vector. This procedure uses the PAR parameter to set the rate of pitch adjustment as follows: $R_{1}^{new} \leftarrow {\begin{matrix} \begin{matrix} Yes & with & probability & PAR \end{matrix} \\ \begin{matrix} No & with & probability \end{matrix} 1 - PAR \end{matrix}$ (5)

The value of 1 - PAR sets the rate of doing nothing. If the pitch adjustment decision for $R_{i}^{new}$ is yes, $R_{i}^{new}$ is replaced as follows: $R_{i}^{new} \leftarrow R_{i}^{new} \pm rand (0, 1) * bw$ (6) where, bw is an arbitrary distance bandwidth, rand (0, 1) is a random number between 0 and 1.

2.4 Update harmony memory

If the new harmony vector, $R^{new} = (R_{1}^{new}, R_{2}^{new}, \dots, R_{N}^{new})$ is better than the worst harmony in the HM, judged in terms of the objective function value, the new harmony is included in the HM and the existing worst harmony is eliminated from the HM.

2.5 Checking the stopping criterion

If the stopping criterion (maximum number of improvisations) is satisfied, the algorithm is terminated. Otherwise, Steps 3 and 4 are repeated.

3 The improved harmony search algorithm

To address the shortcomings of the HS, Mahdavi et al. [7] proposed a new variant of the HS, called the improved harmony search (IHS). The IHS dynamically updates PAR according to the following equation [8]. $PAR (t) = {PAR}_{min} + \frac{{PAR}_{max} - {PAR}_{min}}{NI} \times t$ (7) where, PAR (t) is the pitch adjusting rate for generation t, PAR _max the minimum adjusting rate, PAR _min is the maximum adjusting rate and t is the generation number.

In addition, bw is dynamically updated as follows: $bw (t) = {bw}_{max} exp (\frac{ln (\frac{{bw}_{min}}{{bw}_{max}})}{NI} \times t)$ (8) where, bw (t) is the bandwidth for generation t, bw _min is the minimum bandwidth and bw _max is the maximum bandwidth.

A major drawback of the IHS is that the user needs to specify the values for bw _min and bw _max which are problem dependent and it is difficult to guess them accurately.

4 Differential evolution

Differential Evolution (DE) is a parallel direct search method [9] which utilizes NP parameter vectors by the size of D as a population for each generation G: $x_{i, G}, i : 1, 2, \dots, NP$ (9)

The parameters of NP does not change during the minimization process. The initial vector population is chosen randomly and should cover the entire parameter space. As a rule, we will assume that all random decisions to be uniformly distributed unless it was stated. In the case that a preliminary solution is available, the initial population might be generated by adding normally distributed random deviations to the nominal solution x _nom,0. DE generates new parameter vectors by adding the weighted difference between two population vectors to a third vector. This operation is called mutation. In the next step, the mutated vector’s parameters are then mixed with the parameters of another predetermined vector, the target vector, to yield the so-called trial vector. Parameter mixing is often referred to as “crossover” in the ES-community and will be explained later in more details. If the trial vector yields a lower cost function value than the target vector, the trial vector replaces the target vector in the next generation. This last operation is called selection. Each population vector has to be selected once as the target vector, so NP competitions take place in one generation. DE’s basic strategy will be described in the next subsections [10, 11, 10, 11].

4.1 Mutation

For each target vector; i = 1, 2, 3, …, NP, a mutant vector is generated according to $v_{i, G + 1} = x_{r_{1}, G} + F . (x_{r_{2}, G} - x_{r_{3}, G})$ (10) where, r ₁, r ₂ and r ₃ are mutually different random indexes and F > 0. The randomly chosen integers r ₁, r ₂ and r ₃ are also must be different from the running index i, so NP must be greater or equal to four to satisfy this condition. F is a real constant ∈ [0,2] which controls the amplification of the differential variation: (x _r₂,G - x _r₃,G).

4.2 Crossover

In order to increase the diversity of the perturbed parameter vectors, crossover is introduced. To this end, the trial vector is formed: $u_{i, G + 1} = (u_{1 i, G + 1}, u_{2 i, G + 1}, \dots, u_{Di, G + 1})$ (11) where,

$\begin{matrix} u_{ji, G + 1} \\ = {\begin{matrix} \begin{matrix} v_{ji, G + 1} & if & (rand (j) \leq CR) & \begin{matrix} or & j = rnbr (i) \end{matrix} \end{matrix} \\ \begin{matrix} x_{ji, G} & if & (rand (j) > CR) & \begin{matrix} or & j \neq rnbr (i) \end{matrix} \end{matrix} \end{matrix} \\ j = 1, 2, \dots, D \end{matrix}$ (12)

In Equation (12), rand (j) is the jth evaluation of a uniform random number generator with outcome∈ [0, 1]. CR is the crossover rate ∈[0, 1] which has to be determined by the user. rnbr (i) is a randomly chosen index ∈{1, 2, … , D} which ensures that u _i,G+1 gets at least one parameter from v _i,G+1.

4.3 Selection

To decide that the trial vector should become a member of generation G + 1 or not, u _i,G+1 is compared with the target vector x _i,G using the greedy criterion. If u _i,G+1 yields a smaller cost function value than x _i,G, then x _i,G+1 is set to u _i,G+1, otherwise the old value x _i,Gis retained.

4.4 Other variants of DE

The explained scheme is not the only variant of DE which has proven to be useful. In order to classify the different variants, the following notation is introduced [12, 13]: $DE / x / y / z$ where, x specifies the vector to be mutated which can be “rand” (a randomly chosen population vector) or “best” (the vector of lowest cost from the current population). y is the number of difference vectors used. z denotes the crossover scheme. The current variant is “bin”.

Using this notation, the basic DE-strategy can be written as: DE/rand/1/bin. This is the DE-variant which is used for all performance comparisons in this paper. Nevertheless, one highly beneficial method is DE/best/2/bin where: $v_{i, G + 1} = x_{r_{1}, G} + F . (x_{r_{1}, G} + x_{r_{2}, G} - x_{r_{3}, G} + x_{r_{4}, G})$ (13)

The usage of two difference vectors seems to improve the diversity of the population if the number of population vectors NP is great enough.

5 The new DE-HS algorithm

In this section, a novel hybrid algorithm, namely DE-HS, is introduced which combines characteristics of differential evolution algorithm based on the three following conditions:

To use HMCR control parameter.

To use BR control parameter.

To use vectors to create the difference vector.

As it is shown in Figure. 1, in the first line, some of the data structures and control parameters are initialized. HMS control parameter determines the number of solution and each solution is a D-dimensioned vector. HM is a two dimensional matrix in which the number of rows is HMS and the number of columns is D. In this case, UB is upper bound and LB is the lower bound. The application of HMCR control parameter is similar to HS algorithm and it determines the number of available solutions. Control parameters of F and CR are similar to DE algorithm, so F factor sets the magnitude of differential vector and CR is integration factor. The factor to select the best valuable solution is a new control parameter called BR which is defined in the DE-HS algorithm.

In line 2, a new solution is produced which is called sol1. In this algorithm, to determine the number of sol1 solutions, both the control parameter BR and HMCR are used.

First, the current dimension is determined using HMCR, we can use the values in the HM or assign a random value from the search space.

If the specified value of the current dimension of sol1 must be selected from HM, we can use BR again to choose the best available solution in HM or we may choose one of the available solutions randomly in HM.

One of the most important reasons to use the best existing solutions is to make the new solution of sol1 as similar as possible to the optimal solution.

In this way, we can also avoid searching other parts of the search space which are likely unnecessary to be explored.

However, to set the value of each dimension of sol1, the value of a random number is compared with a BR control parameter. If it was smaller, the value of the corresponding dimension is determined by the best solution available from HM in sol1. Otherwise, we will select one of the available solutions in the HM randomly and set the value for the corresponding dimension in sol1.

In line 3, the qualifications of new sol1 solutions are evaluated and maintained in f1. In line 4, the value of the new sol1 solution is set to another solution called sol2. In DE-HS algorithm, in each iteration, two new solutions are generated that the second one depends on the first. This means that the new sol2 solution is a new solution that takes its initial value from sol1.

To determine the new values for sol2 solution, in line 5, according to the 3rd condition, we use the DE algorithm.

First, we choose 4 random solution from HM. Then, we change the values of all sol2 dimensions according to the following equation: $\begin{matrix} Sol 2 [i] = Sol 2 [i] + \\ F . ((HM [j 2] [i] + HM [j 3] [i]) - \\ (HM [j 4] [i] + HM [j 5] [i])) \end{matrix}$ (14)

In the above equation, the f control parameter determines the magnitude of the differential vector. Finally using CR control parameter and k variable, we decide to accept or reject the new values of each dimension of sol2.

In the case of refusal, sol2 value is returned to its initial value which is the same as the final value of sol1.

In line 6, the value of qualification of the new sol2 solution is calculated and is set to f2.

Finally in line 7, one of the two solutions of sol2 and sol1 is selected based on f1 and f2 values. The selected solution is compared with the worst available solution in HM. If the selected solution is better than the worst available solution, the worst one will be replaced by the selected solution.

6 Experimental results

To evaluate the effectiveness of DE-HS algorithm and compare it with the original HS and IHS, some experiments have been conducted and the results are reported in this section. Some well-known benchmark optimization problems in the continuous domain are exploited to examine the performance of the algorithms. To run the algorithms and perform the experiments, the parameters of algorithms should be assigned by admissible values. The parameter setting which is used for HS, IHS and DE-HS in the experiments is shown in Table 1.

In Table 2, the employed benchmark problems are described. Some details including the mathematical representation, upper bound and lower bound ofvariables, multi-modality and the optimum value of the problems are also reported.

Each algorithm on a test function is executed 1000 and 5000 times independently. The criteria to compare algorithms are the mean and standard deviation of function values over 1000 and 5000 independent runs. Each experiment is repeated 20 times with different random values. For the first four benchmark functions, the experiments are done for dimensions of 2, 5 and 10. The other functions have some fixed dimensions. The results of applying three algorithms on the test functions are reported in Tables 3 to 8.

There are two important observations in the reported results. First of all, some experiments are performed on three different dimensions. By increasing the dimension, the size of search space is increased exponentially and the problem become more difficult. Therefore, the mean in the dimension of 10 is greater than the mean in the dimension of 2 (far from the optimum). Another important observation is that when there is more computational time (maximum number of iterations), the algorithm ends up with an improved result. These are the cases for all three algorithms.

Table 3 summarizes the result of all three algorithms on the Sphere function. However, by comparing results it is clear that the proposed algorithm, DE-HS, outperforms two other algorithms in all dimensions: the mean and standard deviation by this algorithm is less than others. Similar results are reported in Table 4 for the Rosenbrock function. It was expected that similar results are obtained for Sphere and Rosenbrock functions because they have similar characteristics.

Table 4 presents different results for Rastrigin function. The main difference is that this function is multi-modal while Sphere and Rosenbrock are single-modal. IHS outperforms HS by reaching improved mean value but DE-HS fails to perform superior than IHS. It seems that multi-modality of this function make it difficult for the proposed algorithm to find the optimums accurately.

The performance of algorithms on another multi-modal function of Greiwank is reported in Table 5. A glance on the obtained values of mean and standard deviation reveals that in this function DE-HS found the optimums accurately and gets the best performance among three algorithms. The range of input variables in this function is 40 times greater than Rastrigin. So, the search space is much greater (the exact size depends on the dimension of solution). This large search space impacts the performance of optimization in HS and IHS leading to poor results compared to DE-HS.

Table 7 shows the reported results for Schaffer function. Although this is a multi-modal function, the proposed algorithm outperforms the others with respect to the mean and standard deviation. One can say that the optimum is almost found by DE-HS while the result of HS and IHS are close to the optimum. Similar results are observed in Table 8 for the Goldstein and Price I. The performance of DE-HS is significantly better than other algorithms.

The results of optimizing the function of Goldstein and Price II are reported in Table 9. The performance of proposed algorithm is slightly better in this function. Impressive results are reported for the Powell function in Table 10. The optimum of this single-modal function is almost found by DE-HS. It is interesting that IHS is performed worse than the original HS in thisproblem.

7 Conclusion

In this section, concluding remarks about this research will be mentioned briefly.

The results show that the proposed approach (DE-HS) outperforms the original HS and the improved HS in near all experiments. For the number of iterations of 1000 and 5000, DE-HS reaches the global optima faster than the two other algorithms. However, in some functions and some dimensions the new algorithms does not work appropriately but one can say that the new algorithm outperforms the others. Using DE significantly improves the effectiveness of HS in the some well-known continuous optimization problems.

In addition to admissible performance, one of the main features of the new algorithm is that it has a number of reasonable control parameters and therefore, the new algorithm is easy to use and apply to different optimization problems.

References

Qian

Fan

2012

Hybrid harmony search and artificial bee colony algorithm for global optimization problems

Computers & Mathematics with Applications 64 8 2621 2634

Lee

Geem

2004

A new meta-heuristic algorithm for continuous engineering optimization: Harmony search theory and practice

Computer Methods in Applied Mechanics and Engineering 194 3902 3933

Chakraborty

Roy

Das

Jain

2009

An improved harmony search algorithm with differential mutation operator

Scientific Network for Innovation and Research Excellence 181 1 26

Geem

2010

Recent advances in harmony search algorithm

Studies in Computational Intelligence 270 1 10

Geem

2009

Music-inspired harmony search algorithm, theory and applications

Studies in Computational Intelligence 191 1 14

Omran

MGH

Mahdavi

2007

Global-best harmony search

Journal of Applied Mathematics and Computation 198 2 643 656

Mahdavi

Fesanghary

Damangir

2007

An improved harmony search algorithm for solving optimization problems

Applied Mathematics and Computation 188 2 1567 1579

Fesanghary

Mahdavi

2008

Hybridizing harmony search algorithm with sequential quadratic programming for engineering optimization problems

Computer Methods in Applied Mechanics and Engineering 197 3080 3091

Wang

Zhang

2007

A dynamic clustering based differential evolution algorithm for global optimization

European Journal of Operational Research 183 56 73

10.

Kaelo

Ali

2006

A numerical study of some modified differential evolution algorithms

European Journal of Operational Research 169 1176 1184

11.

Chakraborty

2008

Advances in Differential Evolution

Studies in Computational Intelligence 143 1 31

12.

Qin

Huang

Suganthan

2009

Differential evolution algorithm with strategy adaptation for global numerical optimization

IEEE Transaction on Evolutionary Computation 13 2 398 417

13.

Onwubolu

Davendra

2009 Differential Evolution: A Handbook for Global Permutation-Based Combinatorial Optimization Springer Berlin Heidelberg

14.

Geem

2008

Novel derivative of harmony search algorithm for discrete design variables

Applied Mathematics and Computation 199 1 223 230

15.

Das

Mukhopadhyay

Roy

Abraham

Panigrahi

2010

Exploratory power of the harmony search algorithm: Analysis and improvements for global numerical optimization

IEEE Transactions on Systems, Man, and Cybernetics, Part B: Cybernetics 99 1 18

16.

Yousefian

Rafe

Rahmani

2014

A heuristic solution for model checking graph transformation systems

Applied Soft Computing 24 169 180

17.

Rafe

Moradi

Yousefian

Nikanjam

2015

A meta-heuristic solution for automated refutation of complex software systems specified through graph transformations

Applied Soft Computing 33 136 149

18.

Rafe

Hajvali

2013

Designing an architectural style for pervasive healthcare systems

Medical Systems 37 1 13

19.

Mahdian

Rafe

Rafeh

Miralvand

MRZ

2009

Considering faults in service-oriented architecture: A graph transformation-based approach

International Conference on Computer Technology and Development