Adaptive harmony search algorithm utilizing differential evolution and opposition-based learning

1.   Introduction

Metaheuristic algorithms (the genetic algorithm (GA) ^[1], the particle swarm optimization (PSO) algorithm ^[2], the harmony search (HS) algorithm ^[3], etc.) are widely used to solve many optimization problems due to the mechanisms providing high-quality solutions, reasonable processing times, and other advantages. Among the series of classic metaheuristic algorithms, the HS algorithm has several advantages. For example, all individuals have a chance to influence the generation of new individuals, and the value of each dimension of the new individual is generated independently ^[4,5]. Therefore, the HS algorithm is used in various field optimization problems, such as structural design ^[6], object detection ^[7], economic dispatch ^[8], vehicle routing ^[9], the estimation of the optimum number of wind turbines ^[10], buffer allocation problems ^[11], weighted Fuzzy production rule extraction ^[12], and tool indexing ^[13]. Although research on the HS algorithm has gained fruitful and meaningful achievements, the development of the algorithm is still in the initial era because of the algorithm's shortcomings, such as its slow convergence and easy falling into a local optimum ^[4,5,14]. Currently, many improved HS algorithms have been suggested to obtain better performance, and they have focused on two aspects: 1) Automatically adjusting related parameters using dynamic mechanisms. E.g., ^[15] presents an improved harmony search (IHS) algorithm that includes dynamic adaptation for the pitch adjustment rate (PAR) and bandwidth (bw). ^[16] replaces the pitch adjustment step with a new strategy based on the optimal value in the harmony memory (HM). ^[17] presents a global dynamic harmony search algorithm (GDHS). In the GDHS, all parameters and the search domain can be adjusted adaptively. 2) Hybridization with other metaheuristic algorithms results in unique highlights. E.g., ^[18] presents an improved differential-based harmony search algorithm with linear dynamic domain (ID-HS-LDD), which has two characteristics. On one hand, the fine-tuning parameter is generated by utilizing differential evolution. On the other hand, an adaptively adjusted mechanism is used to control the search domain. ^[19] uses a new harmony selection mechanism inspired by a global best concept of PSO and the roulette wheel memory consideration. ^[20] presents an improved HS algorithm hybridized with differential evolution. ^[21] modifies the HS algorithm using the cooling strategy of the simulated annealing (SA) algorithm. The abovementioned HS variants all achieve better performance than the HS algorithm. However, some variants accelerate the convergence speed but easily fall into local optima, leading to premature convergence. Some variants do not consider the harmonies in HM when fine-tuning, leading to low precision. Some variants even adopt the forced convergence strategy, but this strategy is only effective for objective functions that can obtain the optimal solution at position zero. It makes sense to find a new variant of the HS algorithm to overcome the above shortcomings.

Inspired by the above references, an adaptive HS algorithm utilizing differential evolution and opposition-based learning (AHS-DE-OBL) is employed in this manuscript to solve the HS algorithm's existing shortcomings by hybridizing differential evolution (DE), opposition-based learning (OBL), and the adaptive adjustment of algorithm parameters. The DE algorithm is an optimization algorithm utilizing the theory of swarm intelligence to generate optimal values through cooperation and competition among the individuals in a group ^[22]. The mutation operation in DE considers the influence of the original individual on the new individual. In order to make the fine-tuning parameter change with the distribution of the harmonies in HM, AHS-DE-OBL generates a fine-tuning parameter utilizing DE in each iteration's improvisation stage. The OBL strategy can regularly generate new individuals and speed up convergence by comparing individuals' adaptability with their opposite individuals and retaining the fittest ^[23]. In order to increase the richness of HM and avoid easily becoming stuck in the local optimum, AHS-DE-OBL generates extra harmonies utilizing OBL. After updating, AHS-DE-OBL dynamically adjusts the search domain according to the distribution of harmony in HM. In the experiment, ten classic benchmark functions are used to compare AHS-DE-OBL with other HS variants. The experimental results show that the AHS-DE-OBL algorithm has better performance than the selected HS variants and the selected other metaheuristic algorithms.

The outline of the remainder of this article is summarized as follows. In Section 2, a simple introduction to the HS algorithm is given. In Section 3, the influences of the parameters on the HS algorithm are analyzed. In Section 4, the new algorithm, named AHS-DE-OBL, is described in detail. Section 5 shows the experiment and discussions. Finally, Section 6 presents a summary.

2.   Brief introduction of the HS algorithm

Assume $f(x)$ is the objective function, and the HS algorithm is used to search the minimum value in the interval [L, U], where L and U represent the lower and upper bounds of the search domain, respectively. The steps of the HS algorithm are as follows:

Step 1: Initialize the parameters. The parameters include the Harmony Memory Size (HMS), Harmony Memory Consideration Rate (HMCR), Pitch Adjustment Rate (PAR), and bandwidth (bw). The meaning of each parameter is shown in Table 1.

Step 2: Initialize the HM. As shown in Eq (1), randomly generate HMS harmonies and store the harmonies in the HM. In addition, a column is added to the HM to store the objective function value.

where $x_{i, j}$ is the jth dimension value of the ith harmony, and rand is a random number between zero and one.

Step 3: Improvisation. Generate a new harmony so that the HM can dynamically change. The specific process is shown below.

where D is the dimension of the harmony.

Step 4: Update the HM. If the objective function of the new harmony is better than the worst harmony in the HM, the worst harmony will be replaced by the new harmony.

Step 5: Check the termination conditions. If not satisfied, repeat Steps 3 to 5; otherwise, the algorithm ends.

3.   Parameter analysis

The elements in a harmony are affected by some random numbers; hence, the variance of the harmony vector in the HM is a random variable, and the expected variance will be a measure of the algorithm's exploration ability. The larger the expected variance is, the stronger the global search ability; conversely, the smaller the expected variance is, and the stronger the local search ability ^[18]. After the HM is updated, the new harmony generated during the improvisation stage may change the variance of the harmony in the HM. To avoid falling into the local optimum and to ensure that the searched area is as wide as possible, the parameters must be able to adjust the overall variance to achieve a proper balance between the global search and local search. In the HS algorithm, because the disturbance of each decision variable in the harmony is independent, the analysis of one-dimensional vectors does not lose generality and can be extended to multiple dimensions. Literature ^[24] analyzes the influence of the parameters in the HS algorithm on the exploration ability of the algorithm when L and U are symmetric about zero. To improve the universality of the conclusion in ^[24], the following theorem can be obtained.

Theorem 1: Assume that the harmony variables in the initial HM are $X = \{x_1, x_2, \ldots, x_{HMS}\}$. After several improvisation stages, the harmony variables in the HM are $Y = \{y_1\, y_2, \ldots, y_{HMS}\}$, and the expected variance of Y can be described as Eq (2):

where $Var(x)$ is the variance of the original harmony variable set, and $\bar{x}^{2}$ is the square of the mean of the original harmony variable set.

Proof: According to the literature ^[24], we have

According to the step of randomly generating the harmony vector, we have

where $R$ represents a uniformly distributed stochastic number between 0 and 1. $R$ is uniformly distributed in the interval $[0, 1]$. Then,

According to the continuous stochastic variable expectation calculation formula,

Therefore,

Then, from Eqs (5) and (6), we get Eqs (17) and (18) as

From Eqs (7) and (8), we get

Therefore, Eq (9) can be further described as Eq (21)

and the theorem is proved.

From Theorem 1, HMS has a limited impact on $E(Var(Y))$. For example, assuming $HMS = 10$, then $(HMS-1)/HMS = 0.9$; and assuming $HMS = 100$, then $(HMS-1)/HMS = 0.99$. Therefore, choosing a smaller suitable HMS helps to reduce the calculation but can obtain similar results with larger ones. The larger the HMCR is, the greater the influence of $Var(x)$ and $bw^{2}$ on $E(Var(Y))$; therefore, the process is most likely to select a harmony in the HM as the basis of new harmony that is generated. The smaller the HMCR is, the greater the influence of $(U-L)^{2}$on $E(Var(Y))$, and the new harmony is most likely to be generated randomly. The larger PAR is, the greater the influence of $bw^{2}$ on $E(Var(Y))$ and the greater the disturbance of the harmony vector. The influences of the parameters on the algorithm are shown in Table 2.

4.   Adaptive harmony search utilizing differential evolution and opposition-based learning

Based on the above analysis, this paper improves the HS algorithm from the following three aspects: 1) adaptively adjusts the HMCR, PAR and search domain with the iterations of the algorithm; 2) let the generated bw be related to the harmonies in the HM, which makes the new harmony approach the best harmony in the HM; and 3) to ensure that the algorithm does not fall into the local optimum when the search domain is reduced, the outside domain can also be searched. The corresponding improvement strategies are as follows: (i) dynamic parameter adjustment strategy, including the HMCR, PAR and search domain; (ii) bw generation strategy based on differential evolution; and (iii) new harmony generation strategy based on opposition-based learning. The three strategies are analyzed as follows. The main modifications are shown in Table 3.

4.1.   Analysis of the dynamic parameter adjustment strategy

In the early stage, the search domain should be as wide as possible to avoid becoming stuck in the local optimum, the HMCR needs to set to a small value, and the PAR needs to be as large as possible. However, in the later iteration process, the probability of fine-tuning should be increased to improve the search precision, and the HMCR should be as large as possible. Considering the stability, the PAR should be appropriately reduced to gradually make the harmonies in the HM approach the best harmony and increase the probability of fine-tuning the best harmony. The adjustment algorithm is as follows.

where gn is the current iteration number, and NI is the maximum number of iterations.

With the iterations, gradually narrowing the search boundary can make the algorithm focus on the optimal search domain. ^[18] proposes a linear dynamic domain mode that calculates the new boundary by utilizing the maximum and minimum values of the harmony to gradually change the boundary. For each dimension, the new boundary is calculated by the following equations.

The initial equations are:

The equations during iteration are:

where $x_{{maxbound, j}}$ and $x_{{minbound, j}}$ are the maximum and minimum values of the jth dimension in the HM, respectively, and $x_{U_{new, j}}$ and $x_{L_{new, j}}$ are the upper and lower boundary values of the new search domain of the jth dimension, respectively.

4.2.   Analysis of the bw generation strategy based on differential evolution

Differential evolution is an efficient population-based heuristic algorithm used to find the global optimal value through cooperation and competition between individuals in the population ^[22]. According to this idea, the adaptive ability of bw can be described by Eq (28) as follows.

where $x_{best, j}$ is the value of the jth dimension of the best harmony in the current HM, $x_{i, j}$ is the value of the jth dimension of a random harmony in the current HM, and $x_{worst, j}$ is the value of the jth dimension of the worst harmony in the current HM. Assuming each harmony is 2-dimensional, the geometric meaning of Eq (28) is shown in Figure 1.

where $x_{0}$ and $x_{1}$ respectively represent the 1st dimension and the 2nd dimension, $x_{best}$ is the best harmony in the HM, $x_{rand1}$ and $x_{rand2}$ are randomly selected harmonies from the HM, $x_{new} = x_{rand2}\pm bw\times rand$, and $rand\in(0, 1)$. $x_{new}$ has a 50 percent probability of $x_{rand2}+bw\times rand$ and a 50 percent probability of $x_{rand2}-bw\times rand$. In the case of Figure 1, $x_{new}$ has a certain probability of approaching $x_{best}$ or moving away from $x_{best}$. This strategy ensures that $x_{new}$ does not blindly approach $x_{best}$ and is helpful to increasing the diversity of the population.

4.3.   Analysis of new harmony generation strategy based on opposition-based learning

Although the above strategies can improve the local search ability of the algorithm, they may fall into the local optimum and cannot jump out. And the random generation of a solution from the current population often leads to revisiting the hopeless areas in the search space ^[25]. During the search process, generating a random solution and its opposition solution can search areas outside the new search domain and improve the efficiency of the algorithm ^[23]. In the improvisation stage, two extra opposition solutions are generated by Eqs (29) and (30).

where $x_{new, j}^{'}$ and $x_{new, j}^{''}$ express the jth dimension of the opposition solutions.

Assuming that the harmony has a dimension of two, according to $x_{best, j}$, the opposition solution is generated, as shown in Figure 2. When the new search area is the shaded part of Figure 2, the algorithm can still search for areas outside the new search domain.

4.4.   Description of the AHS-DE-OBL algorithm

In the early stage, the search domain should be as wide as possible The AHS-DE-OBL algorithm can be described by the following steps.

Step 1: Initialize the parameters. The pseudocode is as follows.

Step 2: Initialize the HM. Randomly generate HMS harmonies and store them in the HM.

Step 3: Improvisation. Generate three harmonies randomly. The pseudocode is as follows. In the code, $x_{U_{new}, j}$ and $x_{L_{new}, j}$ are the upper and lower bounds of the jth dimension, respectively. HMCR and PAR are dynamically adjusted with the number of iterations; and bw is generated by the best, the worst, and a random harmony.

Step 4: Update the HM. If the objective function value of the new harmony is better than that of the worst harmony, the worst harmony is replaced by the new harmony.

Step 5: Dynamically adjust the parameters. Update the parameters based on the current iteration number and harmony in the HM. The pseudocode is as follows.

Step 6: Check the termination conditions. If the condition is not satisfied, repeat Step 3 to Step 6; otherwise, the algorithm ends.

4.5.   Time complexity analysis

According to the steps of the algorithms, the time complexity of the HS algorithm and the AHS-DE-OBL algorithm is discussed below, where only one iteration is considered. In the HS algorithm, the initialization assignment of $HMS$ harmonies in the $N$-dimensional search space requires $N\times HMS$ operations, and $N$ operations are needed to calculate the fitness function of an individual. Thus, this stage's time complexity is O$(N\times HMS)$. The improvisation stage's time complexity is O$(N)$. The other stages' time complexity is O$(1)$. In the AHS-DE-OBL, Step 5's time complexity is O$(N\times HMS)$, the other stages are the same as the HS algorithm. Therefore, in T iterations, the time complexity of the HS algorithm is O$(T\times HMS\times N)$, and the AHS-DE-OBL is also. Theoretically, the improved algorithm does not increase the time complexity compared with the original algorithm.

5.   Experiments

5.1.   Test functions used in the experiment

The experiment is conducted using a computer with an Intel ® Core ™ i5 – 3470 CPU @ 3.20 GHz, 4 GB memory, the Windows 10 operating system, and the Python 3.7 programming environment. The ten benchmark functions used in the experiment are listed in Table 4.

5.2.   Experiment results and discussions

To verify the performance of AHS-DE-OBL, it was compared with IHS, GDHS, ID-HS-LDD using ten benchmark functions. Each algorithm was run 30 times independently with 7000 iterations per run. By referring to the relevant literature ^[15,17,18] of the selected algorithms, the parameter settings are shown in Table 5.

Figures 3–10 show the optimization curves of the ten test functions.

Table 6 shows the results of the four HS variants running 30 times for 10 functions. The values in bold font in Table 6 show the optimal values among all algorithms.

As seen from Figures 3–10 and Table 6, AHS-DE-OBL has better global search abilities and a faster convergence speed than the other three HS variants.

In AHS-DE-OBL, three strategies are adopted. The strategy of opposition-based learning can improve the search efficiency of the solution space, the bw generation strategy based on differential evolution, and the adaptive parameter adjustment strategy can improve the precision. Therefore, for most of the selected test functions, AHS-DE-OBL has a faster convergence speed and higher accuracy than the others. For unimodal functions F1 (a convex and bowl-shaped unimodal function), F2 (a tapering-shaped function), F3 (a bowl-shaped function with the optimal solution not at position zero), and F8 (a plate-shaped function), regardless of whether there are 2, 10 or 30 dimensions, AHS-DE-OBL has better performance than the other three methods. For multimodal functions F4 (a multimodal function with a regular distribution of the local minima), F6 (a nearly flat outer region and a large hole at the center multimodal function), F7 (a multimodal function with widely distributed local minima), F9 (a valley-shaped multimodal function), and F10 (a multimodal and highly complex function), which can detect the ability of the algorithm to jump out of the local optimum, AHS-DE-OBL obtains all of the best results. For F5, because ID-HS-LDD uses a forced convergence strategy, which forces the value of the harmony to approach zero in a later iteration stage, it can find the optimal value of F5. However, for F6, which is obtained by moving F5 one unit in the positive direction, the performance of ID-HS-LDD was significantly reduced, and AHS-DE-OBL still maintained fairly good performance. For the 30-dimensional F3, ID-HS-LDD also has a bad optimization result.

To verify the improved algorithm does not increase the time complexity compared with the original algorithm, the execution time of the four algorithms in the same experimental environment is recorded. As shown in Tables 6 and 7, although it takes less time to execute IHS and GDHS than the proposed algorithm, their optimization accuracy is lower than the proposed algorithm. The accuracy of ID-HS-LDD is the same as the proposed algorithm, however, it takes more time to execute than the proposed algorithm.

To further verify the performance of AHS-DE-OBL, it was compared with Sine Cosine Algorithm (SCA), self-adaptive DE (SaDE), and multi-scale cooperative mutatingly self-adaptive escape PSO (MAEPSO) using ten benchmark functions. Each algorithm was run 30 times independently with 7000 iterations per run. By referring to the relevant literature ^[18,26,27] of the selected algorithms, the parameter settings are shown in Table 8.

where $populaiton size$ is population size, $crm$ is the initial value of adaptive crossover rate, $learningPeriod$ is the update period of parameter $p$, $p$ is the selection probability of mutation strategy, $crPeriod$ controls the update frequency of the crossover rate for each individual, $crmUpdatePeriod$ controls the update frequency of the parameter $crm$, $M$ is the number of multi-scale Gaussian mutation operators, $subP$ is the number of individuals in a subgroup, $c1$ and $c2$ is the learning factor, $k1$ and $k2$ is the escape threshold. Table 9 shows the experimental reuslts.

As shown in Table 9, except for F3 and F5 of 10 dimensions, the proposed algorithm has achieved better results than other methods. About the SCA, it has obtained better values in seven cases, among which the F5 with 10-dimension is better than the proposed algorithm, but there is little difference in accuracy between them. For the SaDE, it has obtained better values in eight cases, among which the F3 is better than the proposed algorithm. However, it performs better in F7 with 30-dimension than in F7 with 10-dimension, showing its instability. While for the MAEPSO, it has obtained better values in three cases, among which the F3 with 10-dimension is better than the. But it performs relatively worse for the other cases.

To further analyze the experimental results, the Wilcoxon test is used, and Table 10 shows the results. Where the Wilcoxon column represents the comparison result of the selected algorithms and the proposed algorithm, '+' means the proposed algorithm is better than the selected algorithms in this case, '-' is the opposite. '$\approx$' means the two algorithms achieve the same results. the Rank column is the ranking of their mean solution accuracy. For the selected test function, it can be seen that the performance of IHS and GDHS is not good, therefore SCA, SaDE, MAEPOS, and ID-HS-LDD are chosen to compare with the proposed algorithm.

As shown in Table 10, for low-dimensional functions, the proposed algorithm can obtain the same results as some improved algorithms. As the dimensionality increases, in most cases, the performance of the proposed algorithm is significantly better than the selected improved algorithms.

6.   Conclusions

Aiming to address the inherent shortcomings of the HS algorithm, such as its slow convergence speed and low search precision, an improved HS algorithm, named AHS-DE-OBL, was proposed. AHS-DE-OBL is based on differential evolution, opposition-based learning, and a search domain adaptive adjustment strategy. In the improvisation stage, which uses the best and worst harmonies to affect bw, an opposition-based learning strategy is used to increase the diversity of the harmony in HM. After improvisation is completed, the range of the search domain is dynamically adjusted to increase the search efficiency. The experimental results also show that AHS-DE-OBL has better robustness and a better adaptive ability than the three selected HS variants and the selected metaheuristic algorithms.

Acknowledgments

This work was supported by the Natural Science Foundation of Hunan Province, China (Nos.2019JJ40234 and 2020JJ5458), the Research Foundation of Education Bureau of Hunan Province, China (No.19A414 and 18B317), the Special Project of Language and Writing Application Research of Hunan Provincial Language Commission (No. XYJ2019GB09), the National Natural Science Foundation of China (No.62066016), and Jishou University Graduate Research and Innovation Project (No. JGY202032).

Conflict of interests

The authors declare there is no conflict of interest.