Cellular Harmony Search for Optimization Problems

Abstract

Structured population in evolutionary algorithms (EAs) is an important research track where an individual only interacts with its neighboring individuals in the breeding step. The main rationale behind this is to provide a high level of diversity to overcome the genetic drift. Cellular automata concepts have been embedded to the process of EA in order to provide a decentralized method in order to preserve the population structure. Harmony search (HS) is a recent EA that considers the whole individuals in the breeding step. In this paper, the cellular automata concepts are embedded into the HS algorithm to come up with a new version called cellular harmony search (cHS). In cHS, the population is arranged as a two-dimensional toroidal grid, where each individual in the grid is a cell and only interacts with its neighbors. The memory consideration and population update are modified according to cellular EA theory. The experimental results using benchmark functions show that embedding the cellular automata concepts with HS processes directly affects the performance. Finally, a parameter sensitivity analysis of the cHS variation is analyzed and a comparative evaluation shows the success of cHS.

1. Introduction

The optimization techniques have the utility of navigating the search space using effective operators driven by control parameters. The tricky point of the success of any optimization method is its ability to strike a suitable balance between exploration (diversification) and exploitation (intensification) of the problem search space [1]. Exploration is the optimization method capability of navigating a promising region of the search space, if necessary, while exploitation refers to the capability of fine-tuning the already-navigated regions to converge into the local optima [1].

Harmony search (HS) algorithm is a recent evolutionary algorithm (EA) proposed by Geem et al. [2] to imitate the musical improvisation process. Due to its advantages over other optimization methods, it stipulates fewer mathematical requirements in the initial search [3]. It has a novel stochastic derivative which reduces the number of iterations required to converge towards local minima [4], in addition to being simple, adaptable, general, and scalable [5]. Therefore, HS algorithm has been intensively tailored for several optimization problems such as timetabling [6, 7], nurse restoring [8], space allocation [9], and many others [10–13]. However, due to the complex nature of some optimization problems and the avoidance of a premature convergence situation, HS theories are modified [14]. Furthermore, the control parameter adaptation for HS is also studied [5, 15–17].

In a procedural context, HS, which is an iterative improvement algorithm, initiates with a population of random individuals stored in harmony memory (HM). At each iteration, a new individual is generated based on three operators: (i) memory consideration, which selects the variables of a new individual from whole HM individuals; (ii) pitch adjustment, which is responsible for local improvement, and (iii) random consideration, used to provide random elements for the new individual. The new individual is then evaluated and replaces the worst individual in HM, if it is better. This process is cycled until a desired condition is met (see Table 3).

As other EAs, HS algorithm interacts with whole individuals in the HM during each breeding step. The update process selects the worst individual from a single HM and replaces it with a new one, if better. The decentralized methods used in other structured EAs have shown that the performance of EA is improved. Examples include cellular genetic algorithm (cGA) [18], distributed EA (dEA) [19], cellular PSO [20], and others [21]. The main idea of these structured methods is to partition the population into several sets with common features [22].

Cellular genetic algorithm (cGA), in particular, is a decentralized method where the population is represented as a toroidal grid of two-dimensions, as shown in Figure 3 [23, 24]. The individuals are located in this toroidal in a predefined topology and solely interact with their nearest neighbors in the breeding step [22]. Note that all the neighborhoods have the same size and identical shape. This concept embedded in cGA provides useful advantages for the optimization domain [23] and parallel implementations [25, 26] because it assists in providing a high-level of diversity and yields a small diffusion of solutions through the search. The cGA provides a proper exploitation power inside each neighborhood of an individual by the operators of GA [22]. Theoretically, it has been shown that nonrandom mating keeps genetic diversity at higher level, thus preventing the algorithms from converging prematurely to the local optima [27].

The main objective of this paper is to embed the cellular automata concepts in the HS algorithm optimization framework, where the HM is arranged as a two-dimensional toroidal grid. The population diffusion will be expectably preserved, maintaining a high-level of diversity during the search, thus avoiding genetic drift. The improvisation process of HS algorithm is adjusted to interact with the neighborhoods of specific individuals. The updating process of HM is done within the neighborhoods of that individual. The results show that the new decentralized version of HS (i.e., cHS) algorithm improves the performance of HS using standard benchmark functions.

The rest of the paper is organized as follows. The basics of HS algorithm are described in Section 2. The proposed cellular harmony search (cHS) is discussed in Section 3. Results of the experiments are presented in Section 4. Finally, the conclusion and promising future research directions are provided in Section 5.

2. The Harmony Search Algorithm

The harmony search (HS) algorithm is a recent evolutionary approach proposed by Geem et al. [2]. It is initiated with a set of individuals stored in an augmented matrix called harmony memory (HM). It is a centralized algorithm where at each breeding step, it generates a new individual by interacting with the whole individuals in HM. HS follows three rules in the breeding step to generate a new individual: memory consideration, random consideration, and pitch adjustment. If the new individual is better than the worst individual in the whole HM, the replacement process is triggered. This process is repeated as many times as the HS is stagnated. The steps of HS algorithm are flowcharted in Figure 1, where each step is described below in more detail.

Figure 1 

The flowchart of the HS algorithm.

Step 1 (initialize parameters). The optimization problem is initially represented as min{𝑓(𝐱)∣𝐱∈𝐗}, where 𝑓(𝐱) is the objective function and 𝐱={𝑥_𝑖∣𝑖=1,…,𝑁} is the set of decision variables. 𝐗={𝐗_𝑖∣𝑖=1,…,𝑁} is the possible value range for each decision variable, where 𝐗_𝑖∈[𝐿𝐵_𝑖,𝑈𝐵_𝑖], and 𝐿𝐵_𝑖 and 𝑈𝐵_𝑖 are the lower and upper bounds for the decision variable 𝑥_𝑖, respectively, and 𝑁 is the number of decision variables. The parameters of the HS algorithm are also predefined in the following step.(a)The harmony memory consideration rate (HMCR), used in the breeding step to determine whether the value of a decision variable is to be selected from the individuals stored in the harmony memory (HM). (b)The harmony memory size (HMS) which determines the number of individuals in HM. (c)The pitch adjustment rate (PAR), which is used to decide the adjustments of some decision variables selected from memory. (d)The distance bandwidth (BW) which determines the distance of the adjustment that occurs to the individual in the pitch adjustment operator.(e)The number of improvisations (NI) which is similar to the number of generations.

Step 2 (initialize the harmony memory (HM)). The harmony memory (HM) is a matrix of size 𝑁×HMS which includes sets of individuals determined by HMS (see (1)). In this step, these individuals are randomly generated as follows: 𝑥^𝑗_𝑖=𝐿𝐵_𝑖+(𝑈𝐵_𝑖−𝐿𝐵_𝑖)×𝑈(0,1), for all 𝑖=1,2,…,𝑁 and for all 𝑗=1,2,…,HMS, and 𝑈(0,1) generate a random number between 0 and 1. Consider 𝑥HM=¹¹𝑥¹₂⋯𝑥¹_𝑁𝑥²₁𝑥²₂⋯𝑥²_𝑁𝑥⋮⋮⋯⋮_HMS1𝑥^HMS₂⋯𝑥^HMS_𝑁. (1)

Step 3 (improvise a new individual). The HS algorithm generates a new individual, 𝐱^′=(𝑥^′₁,𝑥^′₂,…,𝑥^′_𝑁), using three operators: (1) memory consideration, (2) random consideration, and (3) pitch adjustment.

Memory Consideration. In memory consideration, the value of the first decision variable in the new individual 𝑥^′₁ is randomly selected from the historical values, {𝑥¹₁,𝑥²₁,…,𝑥^HMS₁}, stored in whole HM individuals. Values of the other decision variables, (𝑥^′₂,𝑥^′₃,…,𝑥^′_𝑁), are sequentially assigned in a similar way with a probability of HMCR, where HMCR∈(0,1). It is worth mentioning that this process interacts with the whole individuals in HM which might lead to a premature convergence situation due to the genetic drift.

Random Consideration. Decision variables that are not assigned with values according to memory consideration are randomly assigned according to their possible range by random consideration with a probability of (1-HMCR) as follows: 𝑥^′_𝑖⟵𝑥^′_𝑖∈𝑥¹_𝑖,𝑥²_𝑖,…,𝑥^HMS_𝑖𝑥w.p.HMCR,_′𝑖∈𝐗_𝑖w.p.1−HMCR. (2)

Pitch Adjustment. Each decision variable 𝑥^′_𝑖 of a new individual, 𝐱^′=(𝑥^′₁,𝑥^′₂,𝑥^′₃,…,𝑥^′_𝑁), that has been assigned a value by memory consideration is pitch adjusted with the probability of PAR, where PAR∈(0,1) as follows: Pitchadjustingdecisionfor𝑥^′_𝑖⟵Yesw.p.PAR,Now.p.1−PAR. (3)

In pitch adjustment, if the decision for 𝑥^′_𝑖 is Yes, the value of 𝑥^′_𝑖 is modified to its neighboring value as follows: 𝑥^′_𝑖=𝑥^′_𝑖±𝑈(0,1)×BW.

Step 4 (update HM). If the new individual, 𝐱^′=(𝑥^′₁,𝑥^′₂,…,𝑥^′_𝑁), has better objective function value than that of the worst individual 𝐱^worst stored in HM (i.e., 𝐱^worst=𝐱^HMS in case HM is sorted), the worst individual will be replaced by the new one. Note that the worst individual is selected from the whole HM where the decentralized structure of the population is not observed.

Step 5 (check the termination rule). The HS algorithm will repeat Steps 3 and 4 of HS until a termination rule is met, which is normally decided by the value of NI parameter.

3. Cellular Harmony Search (cHS) Algorithm

There has been much interest from researchers and scientists from different fields exploiting the cellular automata (CA) in physics, biology, social science, computer science, and so on. The initial concepts of CA were developed by Neumann [28] and have been an effective research tool to be incorporated with a wide variety of disciplines.

The concepts of cellular automata (CA) are normally concerned with individual perspective. The main idea from CA is to provide a population of a particular structure formulated as a toroidal grid. The cell in the toroidal grid refers to an individual who communicates with his closest neighboring individuals so that all the individuals have exactly the same number of neighbors. This leads us to a kind of locality known as isolation by distance. Normally, the Manhattan distance is used to measure the distance between any individual and his neighbors. Note that the neighboring individuals have identical shapes and the same size [22].

There exist two-different kinds of cellular models based on how the breeding cycle is performed to the individuals. To put it differently, if the cycle is performed to the whole individuals at the same time, the cellular model is said to be synchronous, where the individuals of the next generation are simultaneously build. On the other hand, if the individuals of the population are sequentially updated with a particular order policy, an asynchronous cellular model is stated. For more discussion about the theory of cellular automata, relevant papers can be seen in [22, 29].

Cellular harmony search (cHS) algorithm can be considered as a new decentralized variation of HS which hinges on the structured HM. The individuals in the HM are arranged in the form of two-dimensional toroidal grid. This is meant to keep a high level of diversity during the breeding step and thus increases the chance to converge into global minima. This can be achieved by avoiding the genetic drift and providing a more suitable population diffusion during the search.

Some steps of cHS algorithm to the original version of HS algorithm presented in Section 2 have been adjusted. The adjustments are flowcharted in Figure 2. The breeding step of cHS solely interacts with the neighboring individuals of a randomly selected individual using “cellular memory consideration” operators. The update of HM step is adjusted to replace the worst individual amongst the neighboring individuals with a new individual, but not amongst the whole HM. Figure 3 shows the population on cellular structure (toroidal grid population), which inspires the concept of small neighborhood regions in the cellular automata. Particularly, the overlap of the neighborhood provides implicit technique migration into population. Therefore, the best solutions are diffused smoothly in the whole population, where the diversity of the cellular harmony search is preserved throughout the search. The cellular HS model includes several components.(1)Cell: the selected random individual in the population (the number of individual is HMS). (2)Cell space: the set of the whole individuals in HM. (3)Neighborhood: the set of potential mates of any individual.(4)Neighborhood shapes: the way of selecting the neighborhoods of the cell as seen in Figure 5.(5)Discrete time limit: the number of generations in HS algorithm which is normally determined by NI parameter.

Figure 2 

The flowchart of cHS.

Figure 3 

The population arrangement based on cellular structure where each cell is an individual in HM.

The detailed steps of cHS are discussed in the steps below.

Step 1 (initialize cHS parameters and optimization problem). It is clear that the successful search of any metaheuristic method is based on skillful parameter setting. The parameters have different effects on optimization solutions. The parameters of cHS are harmony memory size (HMS), harmony memory consideration rate (HMCR), pitch adjusting rate (PAR), number of improvisation (NI), and the size of neighborhood (NH) determined by cellular structure (see Figure 5), where if NH=L9, this means the neighbors structure is length with 8 neighbors.

Step 2 (initialize the harmony memory (HM)). This step is the same as Step 2 in the original version of HS. Note that the individuals of the HM are arranged as a two-dimensional toroidal grid as shown in Figure 4.

Figure 4 

The individuals of HM mapped into the toroidal mesh Y.

Figure 5 

Neighboring shapes [22].

Step 3 (initialization of neighborhood matrix (NM)). NM is a binary matrix of size HMS×HMS (see (4)). The binary values are assigned to each element NM_𝑖,𝑗 based on NH parameter. This matrix is used during the breeding step to determine the neighboring individuals of any randomly selected individual. The NM matrix is filled by binary value as in (5). Consider 𝐍𝐌=NM_1,1……NM_1,HMS⋮⋱⋱⋮NM_HMS,1……NM_HMS,HMS,NM_𝑖,𝑗⟵1𝐱^𝑗∈𝒩𝐱^𝑖∀𝑖∧𝑗∈{1,…,HMS}.0otherwise. (4)(5)

The 𝒩(𝐱^𝐢) is a set of all neighboring individuals of the individual 𝐱^𝑖 arranged in two-dimensional toroidal mesh. This set is determined based on a neighborhood shape as seen in Figure 5.

Figure 4 shows how the HM individuals are mapped to toroidal grid. Note that the element 𝑌_1,1 reflects the index of the individual 1 in HM, while the 𝑌_𝐾,𝐾 reflects the individual index HMS in HM. The HMS value is the square value of K (i.e., 𝐾²=HMS).

To map the element (𝑌_𝑖,𝑗) in toroidal grid Y and the individual index 𝑟 in HM, the following will be used: 𝑟=(𝑖−1)×𝐾+𝑗. (6)

To map the index of the individual 𝑟 in HM to the element (𝑌_𝑖,𝑗) in the matrix NM, the following will be used: ()𝑖=int𝑟−1𝐾))+1,𝑗=((𝑟−1mod𝐾+1. (7)(8)

This mapping mechanism between the individual index of HM and the elements in Y is very useful to determine the neighboring of any individual in HM. As shown in Figure 5, there are several neighborhood shapes to determine the neighbors of any individual. For example, the L5 neighborhood shape takes the nearest neighbors of a given cell axial direction. Therefore, to determine the individual indexes in HM that belong to the neighbors of the individual index 𝑟 using L5, the 𝑖 and 𝑗 should be calculated using (7) and (8). And the set of neighbors of individual 𝑟 called 𝜗={𝑌_{𝑖−1,𝑗},𝑌_𝑖+1,𝑗,𝑌_{𝑖,𝑗−1},𝑌_𝑖,𝑗+1}, then (8) is used to map the elements in Y to the corresponding individual indexes in HM. The same step is used if different neighborhood shapes (i.e., L9, C9,…) are used to determine the neighboring individuals of 𝐱^𝑟.

As can be noted, the individuals in the set 𝜗 will be assigned by 1 in NM for the individual 𝐱^𝑟 while others will be assigned by 0.

Step 4 (select random individual from HM). The selection of a random individual 𝐱^𝑟, 𝑟∈{1,…,HMS}, is done to determine the neighboring individuals taken from the matrix NM  (i.e., 𝜉_𝑟={𝐱¹,𝐱²,…,𝐱^𝑚}), where 𝑚 is the number of the neighboring individuals; that is, 𝑚≤HMS. The neighboring individuals interact together in order to generate the new individual in the next step.

Step 5 (generate a new individual). In this step, a new individual, 𝐱=(𝑥́₁,𝑥́₂,…,𝑥́_𝑁), is generated based on the three operators: (1) cellular memory consideration, (2) pitch adjustment, and (3) random consideration. The whole process is drawn in Figure 6.

Figure 6 

The processes fluency of cHS algorithm.

The pitch adjustment and random consideration operators in cHS algorithm are the same as those in the original version of HS algorithm. However, the memory consideration is modified to be inline with the concepts of cellular automata as follows.

Cellular Memory Consideration. As in cellular GA, the cellular memory consideration solely interacts with the close individuals of 𝐱^𝑟 taken from the set of 𝜉_𝑟. The other individuals are not used in the breeding step. The value of the first decision variable 𝑥́₁ is randomly selected from the historical values, {𝑥¹₁,𝑥²₁,…,𝑥^𝑚₁}, stored individuals in the set 𝜉_𝑟. The other decision variables, (𝑥^′₂,𝑥^′₃,…,𝑥^′_𝑁), are sequentially assigned in a similar way with a probability of HMCR, where HMCR∈(0,1).

It is worth mentioning that cellular memory consideration is able to control the diffusion between the individuals in HM, and thus, it is able to preserve the cHS diversity as long as the search process is iterated. By this strategy, the population is structured and it is possible to improve the numerical behavior of the cHS algorithm.

Step 6 (update the harmony memory). This step is modified in cHS algorithm. The worst individual from the set of neighbors 𝜉_𝑟, that is, (𝐱^worst∣𝐱^worst∈𝜉_𝑟∧𝑓(𝐱^worst)≥𝑓(𝐱^𝑖),∀𝑖∈{1,…,𝑚}), is replaced by the new individual 𝑥́, if better. Note that the replacement process is done taking into account the neighboring individuals in the set 𝜉_𝑟 only.

Step 7 (check stopping criterion). In this step, the cHS will stop if the maximum number of the iteration (i.e., NI) is reached; otherwise the algorithm repeats Steps 4 to 6. The pseudocode of cHS is demonstrated in Algorithm 1.

Set HMCR, PAR, NI, HMS, bw. 𝑥𝑗𝑖=𝐿𝐵𝑖+𝑈(0,1)×(𝑈𝐵𝑖−𝐿𝐵𝑖), ∀𝑖=1,2,…,𝑁, ∀𝑗=1,2,…,HMS, 𝑈(0,1); Calculate (𝑓(𝑥́𝑖)), ∀𝑗=1,2,…,HMS; Generate(NM) {generate the neighborhood matrix} While  !𝑆𝑡𝑜𝑝𝐶𝑜𝑛𝑑𝑖𝑡𝑖𝑜𝑛()  do  select random individual 𝐱𝑟  𝑥́=Ø;  for  𝑖=1,…,𝑁  do   if  𝑈(0,1)≤HMS  then    𝑥́𝑖∈𝜉𝑟; {Memory Consideration}    if  𝑈(0,1)≤PAR  then       𝑥́𝑖=𝑥́𝑖+𝑈 (−1,1)×𝑏𝑤;  {Pitch Adjustment}    end if   else    𝑥́𝑖=𝐿𝐵𝑖+𝑈 (0,1)×(𝑈𝐵𝑖−𝐿𝐵𝑖);  {Random Consideration}   end if  end for  find(𝑥worst) where 𝑥worst∈𝜉𝑟  if  𝑓(𝑥)́<𝑓(𝑥worst)  then   include 𝑥́ to HM;   exclude 𝑥́worst from HM;  end if end while

Algorithm 1 

The pseudocode of cellular harmony search (cHS).

4. Experimental and Comparative Evaluation

In this section, cHS algorithm is evaluated using benchmark functions circulated in the literature used to evaluate different variations of HS algorithm. The comparative evaluation is demonstrated while the sensitivity analysis of the control parameters of the proposed method is carried out.

4.1. Comparative Evaluation

In this section, a set of test functions designed for the special session on real-parameter optimization organized in the 2005 IEEE Congress on Evolutionary Computation (CEC 2005) [30], is used. The CEC 2005 comprises 25 test functions including 5 unimodals and 20 multimodals functions, as shown in Table 1. Note that a full discussion about these functions can be taken from Sunganthan et al. [30]. The CEC 2005 provides a suitable number of comparative methods in which the proposed method can be evaluated, as abbreviated in Table 2.

Table 1 

CEC’2005 functions.

Table 2 

Key to CEC’2005 comparative methods.

Table 3 

Average error rate obtained in CEC’2005 special session in dimension 10.

The experiments done followed the conditions of CEC 2005 [30], where the 25 repeated runs have been performed for each test function. The 25 runs have been summarized in terms of average of the error (AE) of the best individual (i.e., AE=|𝑓(𝐱^∗)−𝑓(𝐱^best)|). Note that 𝐱^∗ is a given optimal solution while the 𝐱^best is the average best solution obtained in 25 runs. The dimension 𝑁=10 and the cHS is iterated 100,000 evaluations of the fitness function.

Notably, most of the winner comparative methods are hybrid versions of a particular EA, where their results are very efficient to the tested functions. It is appeared that the AE of HS and cHS are very close or sometimes better than those achieved by the comparative methods. In particular, the HS algorithm is able to achieve very powerful results for most test functions and excels some of the best results reported by the comparative methods. For example, HS has achieved the smallest AE for 𝑓₁,𝑓₂,𝑓₉,𝑓₁₄, and 𝑓₂₃.

It is noted that the AE obtained by cHS is not the best in most cases. This is because the main idea of proposing the cHS is to ensure a high-level of diversity, and thus, the cHS required a higher number of evaluations to hit better results. However, the results obtained by cHS are very close to the winner results.

4.2. Sensitivity Analysis of cHS Parameters

All the experiments are run using a computer with 2.66 Intel Core 2 Quad with 4 GB of RAM. The operating system used is Microsoft windows Vista Enterprise Service Pack 1. The source code is implemented using MATLAB Version 7.6.0.324(R2008a).

The common parameters among all algorithms used in the experiments are set based on empirical guidelines [2]. For the sake of studying the effect of different parameter settings, in general, the parameters setting used for evaluating the cHS method are as follows: HMS=25, HMCR=0.98, PAR=0.3, 𝑁=30, NH=9, and NI=5×10⁴. The neighborhood shape used is C9, presented in Figure 5.

All functions are implemented in 30 dimensions (30D). For the scalability study in Section 4.2.6, the functions are implemented in 10 dimensions (10D), 50 dimensions (50D), and 100 dimensions (100D). The stopping criterion used is the maximum number of improvisation (NI).

All the results in this section are presented in Tables 4 to 8, demonstrating the mean and standard deviation for thirty independent runs. The best solution has been highlighted in bold font. A comparative analysis between cHS algorithm and the original variant of HS for the common parameters is also conducted.

Table 4 

The effect of varying the HMCR parameter on HS and cHS for ten functions (𝑓₁--𝑓₁₀).

4.2.1. Benchmark Functions

The global minimization benchmark functions are used to study the sensitivity analysis of the parameters of the proposed method (cHS) against the original version of HS algorithm. Five functions are defined by Whitley et al. [42] and the other five were described by Yao et al. [43]. These functions provide a balance between unimodal and multimodal functions. These functions are commonly used to evaluate the state-of-the-art variations of harmony search algorithms [5, 44, 45].

Most of the benchmark functions have standard solution space range of the objective function. Otherwise, unsymmetrical initialization ranges are used for these functions whose global optima are at the center of the solution space. These benchmark functions are as follows.

(1) Sphere function, defined as 𝑓₁(𝑥)=^𝑁𝑑_𝑖=1𝑥²_𝑖, (9) where global optimum 𝑥^∗=0 and 𝑓(𝑥^∗)=0 (search range −100≤𝑥_𝑖≤100) (Figure 7(a)).

(a)  2D sphere function

(b) 2D Schwefel's problem 2.22

(c) 2D step function

(d) 2D Rosenbrock function

(e) 2D Rotated hyper-ellipsoid function

(f) 2D Schwefel's problem 2.26 function

(g) 2D Rastrigin function

(h) 2D Ackley function

(i) 2D Griewank function

(j) Six-Hump Camel-Back function

Figure 7 

The benchmark functions landscape where the value of 𝑁=2 (2 dimensions).

(2) Schwefel’s problem 2.22, defined as 𝑓₂(𝑥)=^𝑁𝑑_𝑖=1𝑥_𝑖+^𝑁𝑑_𝑖=1𝑥_𝑖, (10) where global optimum 𝑥^∗=0 and 𝑓(𝑥^∗)=0 (search range −10≤𝑥_𝑖≤10) (Figure 7(b)).

(3) Step function, defined as 𝑓₃(𝑥)=^𝑁𝑑_𝑖=1𝑥_𝑖+0.5², (11) where global optimum 𝑥^∗=0 and 𝑓(𝑥^∗)=0 (search range −100≤𝑥_𝑖≤100) (Figure 7(c)).

(4) Rosenbrock function, defined as 𝑓₄(𝑥)=^𝑁𝑑−1_𝑖=1100𝑥_𝑖+1−𝑥²_𝑖²+𝑥_𝑖−1², (12) where global optimum 𝑥^∗=(1,1,…,1) and 𝑓(𝑥^∗)=0 (search range −100≤𝑥_𝑖≤100) (Figure 7(d)).

(5) Rotated hyper-ellipsoid function, defined as 𝑓₅(𝑥)=^𝑁𝑑_𝑖=1^𝑖_𝑗=1𝑥_𝑗², (13) where global optimum 𝑥^∗=0 and 𝑓(𝑥^∗)=0 (search range −100≤𝑥_𝑖≤100) (Figure 7(e)).

(6) Generalized Schwefel’s problem 2.26, defined as 𝑓₆(𝑥)=−^𝑁𝑑_𝑖=1𝑥_𝑖sin𝑥_𝑖, (14) where global optimum 𝑥^∗=(420.9687,420.9687,…,…,…,420.9687) and 𝑓(𝑥^∗)=−12569.5(search range −500≤𝑥_𝑖≤500) (Figure 7(f)).

(7) Rastrigin function, defined as 𝑓₇(𝑥)=^𝑁𝑑_𝑖=1𝑥²_𝑖−10cos2𝜋𝑥_𝑖+10, (15) where global optimum 𝑥^∗=0 and 𝑓(𝑥^∗)=0 (search range −5.12≤𝑥_𝑖≤5.12) (Figure 7(g)).

(8) Ackley’s function, defined as 𝑓₈(𝑥)=−20exp−0.2130^𝑁𝑑_𝑖=1𝑥²_𝑖1−exp30^𝑁𝑑_𝑖=1cos2𝜋𝑥_𝑖+20+𝑒, (16) where global optimum 𝑥^∗=0 and 𝑓(𝑥^∗)=0 (search range −32≤𝑥_𝑖≤32) (Figure 7(h)).

(9) Griewank function, defined as 𝑓₉(𝑥)=14000^𝑁𝑑_𝑖=1𝑥²_𝑖−^𝑁𝑑_𝑖=1𝑥cos^𝑖√𝑖+1, (17) where global optimum 𝑥^∗=0 and 𝑓(𝑥^∗)=0 (search range −600≤𝑥_𝑖≤600) (Figure 7(i)).

(10) Six-Hump Camel-Back function, defined as 𝑓₁₀(𝑥)=4𝑥²₁−2.1𝑥⁴₁+13𝑥⁶₁+𝑥₁𝑥₂−4𝑥²₂+4𝑥⁴₂, (18) where global optimum 𝑥=(−0.08983,0.7126), (−0.08983,0.7126) and 𝑓(𝑥^∗)=1.0316285 (search range −5≤𝑥_𝑖≤5) (Figure 7(j)).

Figure 7 depicts the corresponding landscape of the search space of each function [46].

4.2.2. The Effect of HMCR

Table 4 summarizes the effect of HMCR on the performance of cHS and the original version of HS using four values of HMCR (i.e., 0.7, 0.9, 0.94, and 0.98). The best solution obtained at each corresponding value is highlighted in bold font. Figure 8 is the box plot visualizing the effect of various values of HMCR on the behavior of cHS algorithm using the ten global minimization functions.

(a) 2D sphere function

(b) 2D Schwefel's problem 2.22

(c) 2D step function

(d) 2D Rosenbrock function

(e) 2D Rotated hyper-ellipsoid function

(f) 2D Schwefel's problem 2.26 function

(g) 2D Rastrigin function

(h) 2D Ackley function

(i) 2D Griewank function

(j) Six-Hump Camel-Back function

Figure 8 

The box plots for showing the effect varying HMCR values using the ten global optimization functions.

The results show that increasing the HMCR value improves the performance of the cHS for all functions, except six-hump where the opposite function is true. Where a small value is used, HMCR increases the diversity and hence prevents the cHS from convergence (i.e., it results in random search). Thus, it is generally better to use a large value for the HMCR (i.e., ≤0.98).

The high value of HMCR means high probability of using the harmony memory that leads to less exploration of search space. Using a probability of HMCR close to 1 (high value) might lead the algorithm to fall into the local minima. Using less probability of HMCR allows more randomly generated solutions. Therefore, the diversity increases in a way that prevents the convergence. The results reveal that cHS and HS have identical sensitivity to the different values of HMCR for all functions, and in 0.98 probabilities to use HMCR they get the best results for the majority of optimization functions. Furthermore, the results produced by cHS are better than those produced by HS in almost all tested functions.

4.2.3. The Effect of HMS

Table 5 summarizes the effect of HMS on the performance of cHS and basic HS using four values of HMS (i.e., 16, 25, 36, and 100). The best solution obtained at each corresponding value is highlighted in bold font. Figure 9 is the box plot visualizing the effect of various values of HMS on the behavior of cHS algorithm using the ten global minimization functions.

Table 5 

The effect of varying the HMS parameter on HS and cHS for ten functions (𝑓₁--𝑓₁₀).

(a) 2D sphere function

(b) 2D Schwefel's problem 2.22

(c) 2D step function

(d) 2D Rosenbrock function

(e) 2D Rotated hyper-ellipsoid function

(f) 2D Schwefel's problem 2.26 function

(g) 2D Rastrigin function

(h) 2D Ackley function

(i) 2D Griewank function

(j) Six-Hump Camel-Back function

Figure 9 

The box plots for showing the effect varying HMS values using the ten global optimization functions.

It is revealed that both cHS and the basic HS are not sensitive to the HMS. For (𝑓₄,𝑓₆,𝑓₁₀), cHS obtained the best results when the value of HMS is large (see Figure 9). Apparently, with different values of HMS, cHS performance is better for most of the functions because the overlap of the neighborhoods provides an implicit mechanism of migration to the cHS. Since the best solutions spread smoothly through the whole population, the cHS diversity in the structured population is preserved longer than in the classical version of HS algorithm.

HM is analogous to the short-term memory of a musician that is known to be small. A plausible interpretation may rely on the high number of similar harmonies within the HM when the HMS is large that leads to shortages of diversity and, hence, lead to falling into local minima. Therefore, cHS is likely to be capable of maintaining the diversity than HS with the cellular structure.

4.2.4. The Effect of PAR

Table 6 summarizes the effect of PAR on the performance of cHS and basic HS using four values of PAR (0.1, 0.3, 0.7, and 0.9). The best solution at each corresponding value is highlighted in bold font. Figure 10 is the box plot visualizing the effect of various values of PAR on the behavior of cHS algorithm using the ten global minimization functions.

Table 6 

The effect of varying the PAR parameter on HS and cHS for ten functions (𝑓₁--𝑓₁₀).

(a) 2D sphere function.

(b) 2D Schwefel's problem 2.22

(c) 2D step function

(d) 2D Rosenbrock function

(e) 2D Rotated hyper-ellipsoid function

(f) 2D Schwefel's problem 2.26 function

(g) 2D Rastrigin function

(h) 2D Ackley function

(i) 2D Griewank function

(j) Six-Hump Camel-Back function

Figure 10 

The box plots for showing the effect varying PAR values using the ten global optimization functions.

It seems that using a relatively small value of PAR (i.e., ≤0.5) improves the performance of the cHA and HS. Most results at large value of PAR can increase the convergence speed of HS algorithm, while a small value of PAR increases diversity in HM. On the other hand, the small value of PAR allows more exploration of the search space, and the large value of PAR leads to a lower rate of exploration, where the diversity is reduced and then the algorithm might be trapped into the local optima. It is observed that cHS is able to get the best results than the original version of HS algorithm for some benchmark functions (i.e., 𝑓₄, 𝑓₂, 𝑓₈, 𝑓₉), but not so for the others.

4.2.5. The Effect Number of Neighbors NH

Table 7 summarizes the effect number of neighbors (NH) according to cellular structure on the performance of cHS. Table 7 also exhibits the effect of NH for (𝑓₁ to 𝑓₁₀) using four values of NH (4,8,12), where the value of HMS=25. The best solution at each corresponding value is highlighted in bold font. The best results of comparing cHS with original version HS are highlighted in bold. Figure 11 is the box plot visualizing the effect of various NH values on the behavior of cHS algorithm using the ten global minimization functions.

Table 7 

The effect of varying the number of neighbors on cHS for ten functions (𝑓₁−𝑓₁₀).

Table 8 

Mean and standard deviation of ten functions (𝑓₁--𝑓₁₀) when 𝑁=(10,30,50,100).

(a) 2D sphere function

(b) 2D Schwefel's problem 2.22

(c) 2D step function

(d) 2D Rosenbrock function

(e) 2D Rotated hyper-ellipsoid function

(f) 2D Schwefel's problem 2.26 function

(g) 2D Rastrigin function

(h) 2D Ackley function

(i) 2D Griewank function

(j) Six-Hump Camel-Back function

Figure 11 

The box plots for showing the effect varying NH values using the ten global optimization functions.

The results show that cHS obtained the best result when the number of neighbors is large, except (𝑓₆) which obtained the best result when the number of neighbors is small. Generally, this leads to the conclusion that cHS often shows better performance at larger values of NH between 8 and 12. Figure 11 shows the box plots of the recorded results which reveal the distribution of the 30 tested runs against the various value of NH.

4.2.6. Scalability Study of 𝑁

In this section, the results produced by cHS and HS, when the dimension of the function is set to 𝑁=10, 𝑁=30, 𝑁=50, and 𝑁=100 as shown in Table 8, are recorded.

In general, decreasing the dimensionality leads to better results in cHS and HS. This comes in line with the previous theory. However, it is observed from results (Table 8) that the proposed cHS algorithm outperforms the HS when dimensionality is large for the majority of the benchmark optimization functions. This shows that increasing the dimensionality of the problem needs a better algorithm like cHS.

5. Conclusion and Future Work

In this paper, a new version of HS algorithm called cellular harmony search (cHS) algorithm is proposed. cHS is an HS algorithm embedded with cellular automata concepts. The main idea of proposing cHS algorithm is to provide a structured population that preserves a high level of diversity during the search. In cHS, the HM individuals are arranged as a two-dimensional toroidal grid, where each individual is generated and only interacts with its neighboring individuals. The operators of the original version of HS algorithm are adjusted to observe the cellular GA theory, where the concepts of cell and cell search space are employed.

Using ten global optimization functions circulated in the literature, the cHS is evaluated. The results support the theory of cellular automata, where in almost all cases the cHS outperforms HS algorithm. The sensitivity analysis of cHS parameters has suggested that the cHS is sensitive to the values of HMCR, PAR, 𝑁, and NH. The comparative evaluation is also conducted with two versions of HS algorithms proposed in [44, 45], where comparable results were obtained.

This is an initial investigation of using cellular automata concepts in the HS algorithm optimization framework. The future, indeed, is pregnant with several research directions, such as(i)analyzing the selection pressure and time complexity concepts of cHS algorithm,(ii)studying the effect of the neighborhood shapes on the performance of cHS,(iii)studying a new migration strategy to empower the interaction between the individuals and their neighbors.

Acknowledgment

The first author is grateful to be awarded Postdoctoral Research Fellowship from the School of Computer Sciences (USM). This paper is supported by Grant no. 203/PTS6728001 (Grant Type: LRGS).