A hybrid cuckoo–harmony search algorithm for optimal design of water distribution systems

Meta-heuristics have been widely used to handle water resources optimization problems over the past four decades (Maier et al. 2014). The relevant research is especially active in the water distribution system (WDS) design field as this problem type is typically difficult due to the large search space and the highly nonlinear constraints. Nowadays, a vast number of scientific works dealing with WDSs optimization are available in the literature (as, for example, Farmani et al. (2006), McClymont et al. (2013), etc.). The meta-heuristic optimization algorithms that have previously been applied to WDS design problems include genetic algorithms (GAs) (Simpson et al. 1994; Savic & Walters 1997), particle swarm optimization (PSO) (Suribabu & Neelakantan 2006), tabu search (TS) (Cunha & Ribeiro 2004), ant colony optimization (ACO) (Maier et al. 2003; Zecchin et al. 2006), simulated annealing (SA) (Cunha & Sousa 1999), differential evolution (DE) (Vasan & Simonovic 2010), harmony search (HS) (Geem et al. 2002), cross-entropy (Perelman & Ostfeld 2007), and charged system search (Sheikholeslami et al. 2014).

Despite the successes in finding optimal or near-optimal solutions, the practical applications of meta-heuristic algorithms are not without difficulties. One of the main issues is their slow convergence rate towards the optimum, resulting in high computational overheads that are typically beyond the availability for practical application (Maier et al. 2014). This is especially the case for the water resources problems as their simulation models used in the optimization process (function evaluation) are often computationally expensive.

A number of studies have been undertaken to improve the efficiency of the meta-heuristic algorithms applied to water resources problems. These studies can be divided into two main categories according to the means adopted to enhance the convergence speed.

The first category aims to reduce the computational time of the simulation model, which is often the most time-consuming part within the optimization process. By way of example, Broad et al. (2005) introduced a metamodeling approach to optimize a WDS design problem, in which the artificial neural networks (ANNs) were used to replace the hydraulic simulation model. More recently, Bi & Dandy (2014) also used ANNs to substitute the full hydraulic and water quality simulation model in their optimization framework, in order to expedite the convergence rate.

The second category attempts to reduce the computational overheads associated with the meta-heuristic algorithms by means of improving their searching efficiency. The improved efficiency is often achieved through the hybridization of the meta-heuristic algorithms and deterministic optimization techniques. For example, Cisty (2010) proposed a hybrid GA-linear programming (LP) method to seek the minimum-cost design of WDSs. In the GA-LP framework, a GA was used in the outer loop to decompose a complex looped network into a group of equivalent branched networks, followed by the implementation of a LP in the inner loop to optimize each branched network. Zheng et al. (2011) combined a nonlinear programming (NLP) and a DE algorithm for optimizing WDS design, where an NLP was used to approximate the optimal solutions, and a DE algorithm was subsequently employed to polish the identified approximate solutions. In a more recent study, Zheng et al. (2014) coupled a binary linear programming with a DE algorithm to improve the optimization efficiency of the WDS design problems.

In parallel with the development of the meta-heuristic-deterministic optimization models, research has also been carried out to hybridize multiple meta-heuristic algorithms for improving their searching ability. For instance, Keedwell & Khu (2005) proposed a hybrid method, in which a heuristic-based, local representative cellular automata approach was developed to provide a good initial population for GA runs. Geem (2009b) presented an improved HS method for the optimal design of WDSs by incorporating the searching mechanisms of the PSO. More recently, Sedki & Ouazar (2012) developed a hybrid PSO–DE model, taking advantage of the PSO's strong local searching strength and DE's great exploration ability.

The studies mentioned above have made great contributions in building knowledge for improving optimization efficiency of meta-heuristic algorithms. To further progress research in this field, this study aims to develop a new hybrid optimization framework in which a cuckoo search (CS) algorithm is combined with a HS method (denoted as CSHS). The utility of the proposed CSHS method in terms of tackling WDS design problems will be investigated in the present study. It is a well-established fact that the WDS optimization problem belongs to the class of non-deterministic polynomial-time hard problems known as NP-hard which means that for an N-pipe network, the computational time required for a rigorous algorithm is, at best, an exponential function of N and is thus enormous even for relatively small WDS. In particular, when large-scale problems are considered, meta-heuristics such as the proposed CSHS are one of the best alternatives that experts can often rely on, since exact algorithms take exponential time to find an optimal solution to NP-hard problems. Therefore, using a new meta-heuristic optimization technique is one of the interesting, if not the best, way of treating WDS optimization problems. Furthermore, the coupled CSHS–EPANET framework is a multi-purpose model and can be extended to tackle other water network management problems, which is appealing for its practical applications.

The CS algorithm is a relatively new optimization technique, which is designed according to the Lévy flight and brood parasitic behavior of some cuckoo species (Yang & Deb 2009). It has been proven to deliver good performance in dealing with a number of optimization problems, ranging from mathematical problem optimization, engineering design, neural network training, to the traveling salesman problem (Civicioglu & Besdok 2013; Gandomi et al. 2013; Yang & Deb 2014). Despite the great potential in efficiently finding optimal solutions for a range of problem types, the CS algorithm has received limited application in the WDS optimization field. To the authors' knowledge, the only work is Wang et al. (2012), who combined CS with the non-dominated sorting GA (NSGA-II) to optimize WDS design problems with multiple objectives.

This study focuses on the application of the CS algorithm to the WDS design problems in the single-objective domain, and more importantly, attempts to improve its performance through the incorporation of the HS mechanism. The specific contributions of this study include: (i) the investigation of the CS's ability in finding single-objective optimal solutions for WDS design problems; (ii) the proposal of a hybrid CSHS algorithm; and (iii) the demonstration of the CSHS's utility in optimizing WDS design.

Since the proposed CSHS method involves the standard CS and HS algorithms, it is necessary to briefly outline the main mechanisms of both methods, as shown in the following sub-sections.

The CS algorithm is inspired by the brood parasitism of some cuckoo species by laying their eggs in the nests of other host birds (Yang & Deb 2009). In this algorithm, each egg in a nest represents a candidate solution, and a cuckoo egg indicates a new solution. The strategy is to use the new, and potentially better, solutions (cuckoos) to replace a not-so-good solution in the nests, following the three rules below (Yang & Deb 2009):

• Rule 1 Each cuckoo lays one egg at a time, and dumps its egg in a randomly chosen nest, which corresponds to the initialization process of other algorithms such as GAs.

• Rule 2 The best eggs (i.e., solutions) are contained in a fraction of the nests and will carry over to the next generation. This is equivalent to the elitism scheme used in the GAs so that the best solutions are passed on to the next generation.

• Rule 3 The number of nests is fixed and a host bird can find an alien egg with a specified probability, p_a ∈ (0,1). In this case, the host bird can either throw the egg away or abandon the nest so as to build a completely new nest in a new location.

After producing the new solution, X_i^(t+1), it will be evaluated and compared to the X_i^(t). If the fitness of the new solution is better than the objective fitness of the previous solution, X_i^(t+1) is accepted as a new basic solution; otherwise X_i^(t) would be kept (Yang & Deb 2014).

The usage of the harmony memory (HM) as shown in Figure 2 is important because it ensures that good harmonies are considered as elements of new solution vectors. The initial HM consists of a certain number of randomly generated solutions for the optimization problem. For a D-dimension problem, the elements of HM can be represented as [x_i,j] (i = 1,…, HMS; j = 1,…, D). The harmony memory size (HMS) is typically set to be between 10 and 100 (Geem 2009a). After generating an initial HM, each component of a new solution (HM_n^(t)) is obtained based on the harmony memory considering rate (HMCR). The HMCR is defined as the probability of selecting a component from HM members, and 1–HMCR is, therefore, the probability of random generation. If a new solution comes from HM, it can be further mutated according to the pitch adjusting rate (PAR). The PAR determines the probability of a candidate from HM for mutation. In Figure 2, bw is a bandwidth, controlling the local range of pitch adjustment, which should be related to the scales of the problem of interest. In most cases, we can use bw = O (L/10), where L is the characteristic scale of the problem of interest, while in some cases bw = O (L/100) can be more effective and avoid flying too far.

Finally, the new solution is evaluated and if it yields a better fitness than that of the worst member in the HM, it will be replaced instead of that one. Otherwise, it is eliminated. The generation of new solutions and memory consideration steps are alternatively performed until a termination criterion is satisfied.

The standard CS algorithm is one of the most successful optimization techniques (Yang & Deb 2014). However, it has two main drawbacks, especially in the case of complex multi-modal problems. The first weakness is the slow convergence rate of the algorithm (Walton et al. 2011). In other words, given enough computation, the standard CS will always find the optimum solution, but as the search relies entirely on random walks a fast convergence cannot be guaranteed (Yang & Deb 2010). The second weakness of the algorithm is that, in the standard CS, there is no information exchange between individuals and, essentially, the searches are performed independently (Walton et al. 2011; Long & Jiao 2014). In order to overcome these limitations, the proposed hybrid CSHS algorithm combines the optimization capabilities of HS and CS. There are various hybrid models of HS with other meta-heuristics in the literature (Alia & Mandava 2011). The proposed hybrid CSHS integrates some HS components within the standard CS. In order to increase the exploitation capability of the algorithm so as to avoid being trapped into local optima and speeding up the convergence, the CSHS algorithm uses the memory (HM) and pitch adjusting operation from HS algorithm which serves as a mutation operator in the new hybrid method.

The proposed hybrid algorithm is a double stage optimization technique which employs the CS algorithm in the first stage and the HS algorithm in the second one. The first stage utilizes two search capabilities of the standard CS: local search and global search, controlled by a switching/discovery probability of CS (p_a). Also, its global search uses Lévy flights, rather than standard random walks. In the second stage, pitch adjustment is carried out by tuning the solution within a given bandwidth. A small random amount is added to or subtracted from an existing solution stored in HM. Essentially, pitch adjustment is a refinement process of local solutions. The use of a HM in CSHS allows the selection of the best solutions that may belong to different regions in the search space. Since the standard CS algorithm is memory-less, there is no information stored dynamically during the search, while the hybrid CSHS uses a memory that contains some information extracted online during the search. In other words, the history of the search stored in a memory can be used for generating new solution candidates. Therefore, incorporating the HM and PAR into the CSHS algorithm ensures that good local solutions are retained while the randomization makes the algorithm able to globally explore the search space more effectively. Such interactions between global search and local search components could be an important factor for the success of the CSHS algorithm over other algorithms because it guarantees a balanced combination of the exploration and exploitation. In the proposed hybrid method, exploration is represented in the form of randomization with a random component attached to a deterministic component in order to explore the search space effectively and efficiently (e.g., Equation (1)), while exploitation is the utilization of past solutions so as to select the potentially good solutions (Equation (2)) via elitism (Rule 2 in the CS) or use of memory (HM) or both.

In this article, the WDS design is formulated as a least-cost single objective optimization problem with a selection of pipe sizes as the decision variables, while the network layout and its connectivity, nodal demand, and minimum head requirements are imposed as design constraints.

Here, a design Φ is defined as a set of n_p decisions where n_p is the number of pipes to be sized, that is Φ = {D₁, D₂, …, D_np} where D_i is the selected diameter for pipe i. The optimization problem can be stated mathematically as follows:

The static penalty function method with constant δ has certain drawbacks, for example, trial-and-error process has to be performed to determine appropriate values for the penalty exponent. This repeated process is time-consuming and is even not allowed in some real-world applications. To cope with the above problems, in Equation (8), the dynamic penalty method is applied. Here, the exponent δ is selected so that it increases the penalties as the search progresses. The value of the penalty function exponent is important because it governs the rate of increase in the cost of infeasible designs, which directly affects the exploration of the CSHS. In the first iterations, if δ has a large value, the solutions tends to narrow the search space to designs that, while feasible, are more expensive than the optimal design and reduce the exploration of the solution space (Joines & Houck 1994). However, within the last iterations, if δ has a small value, the solutions have an undesirable tendency to converge to least-cost, but infeasible, designs that have a very small penalty. Setting a large value for δ in last iterations may help to prevent convergence to infeasible designs by increasing the applied penalty. Therefore, as the penalty increases it puts more and more selective pressure on the CSHS to find a feasible solution. Thus, in the first step of the search process, δ is set to 1.5 and ultimately increased to 2.5. (These two values were selected after a number of fine-tuning trials.)

Here the CSHS algorithm is coupled with the widely used water distribution network software, EPANET 2 (Rossman 2000). After analyzing a model by EPANET 2, the pressure of each node is obtained, then these values are compared to the allowable limits. EPANET programmer's toolkit was provided by the United States Environmental Protection Agency that is a dynamic link library of functions which allows developers to customize EPANET's computational engine for the user's specific needs. Here, the algorithms were implemented in MATLAB and optimization runs were carried out on a computer with an Intel Core i5 CPU, 2.53 GHz processor, and 3.00 GB RAM. A brief description of the steps in the proposed algorithms for pipe network optimization is given below:

1. Generate N population of nests randomly in the solution space. Each of the N populations represents a possible combination of pipe indices.

2. Compute the network cost C(Φ) for each of the N solutions after converting the randomly generated pipe indices to the pipe sizes available in the market.

3. Update the input file of the simulator (only the diameters are changed).

4. Perform hydraulic analysis of each network. EPANET 2 is used to analyze the network and check the pressure at some nodes which are required to meet certain nodal pressures.

5. Compute violations (∑Δ), if the nodal head at any node is less than the required minimum.

6. Calculate the total cost of the network (F_cost) using the network cost and the penalty found in steps 2 and 5, respectively.

7. The total cost found in step 6 is used as the fitness value for each of the trial networks.

The important parameters of the proposed algorithm are α (scaling factor), p_a (discovering probability of alien eggs/solutions), PAR (pitch adjusting rate), and HMCR which control the local search and global search behavior of the algorithm. Due to the importance of the first two parameters, an extensive sensitive analysis for α and p_a is performed for the first case study (Hanoi problem) in the next section. The other parameters belong to the HS stage of the algorithm and will be determined in this sub-section.

In order to use the HS mechanism effectively, the CSHS algorithm adopts the parameters PAR and HMCR ∈ (0,1). If HMCR is too low, only a few elite solutions are selected, and the algorithm may converge too slowly. If this rate is extremely high, the solutions in the HM are mostly used, and other ones are not explored well, so that good regions may be missed. Furthermore, small PAR values can cause poor performance of the algorithm and a considerable increase in required iterations. However, large PAR values usually cause the improvement of best solutions in final generations when the algorithm converges to optimal solution vector. To solve this problem, inspired by the SGHS (a self-adaptive global best HS) algorithm (Pan et al. 2010), a self-adaptive strategy is presented in this section. Unlike the standard HS algorithm, this technique employs an adaptive parameter tuning method. By using this strategy HMCR and PAR become the configurable components of the CSHS and can be adjusted to alter the performance of the algorithm. Table 1 lists the parameters (that must be set before algorithm execution) required for major types of meta-heuristics including GA, PSO, ACO, HS, SA, and CSHS.

In this method, the main assumption is that the HMCR (or PAR) value is normally distributed in the range of [0.8, 0.99] ([0.01, 0.5]) with mean HMCR_m (PAR_m) and standard deviation of 0.01 (0.05) (Pan et al. 2010). Initially, HMCR_m (PAR_m) is set at 0.85 (0.25) (Poursalehi et al. 2013); and then the algorithm starts with a HMCR (PAR) value generated according to the normal distribution. During the evolution process, the HMCR (PAR) value of a generated solution that has successfully replaced the worst member in the HM is recorded. After a specified number of generations (known as the learning period), HMCR_m (PAR_m) are re-computed by averaging all the recorded good HMCR (PAR) values during this period. With the new mean and the given standard deviation of 0.01 (0.05), a new HMCR (PAR) value is produced and used in the subsequent iterations. The above procedure is repeated, and so an appropriate HMCR (PAR) value can be gradually learned to suit the specific problem within the specific phases of the search process. In this way, the difficulty of finding suitable values for these parameters is solved and a self-adaptation operator adjusts the HMCR and PAR values automatically. According to previous studies, the choice of the start points for HMCR and PAR does not affect the performance of the algorithm (Pan et al. 2010; Kulluk et al. 2011; Poursalehi et al. 2013).

In this section, four well-known benchmark networks are optimized using the present method. The final results are compared to the solutions of other advanced meta-heuristic methods to demonstrate the performance of the CSHS method. The first network (N₁) is the Hanoi network, first proposed by Fujiwara & Khang (1990), which consists of 32 nodes, 34 pipes, and 3 loops. This network has no pumping station as it is fed by gravity from a reservoir with a 100 m fixed head. The minimum head requirement of the other nodes is 30 m and the Hazen–Williams coefficient for each new pipe is 130. The second design example (N₂), first introduced by Cisty (2010), is the double Hanoi network. Since this network is derived from the basic Hanoi network, its optimal cost is known. All the parameters for the reservoir, nodes, and lines in the double Hanoi water distribution network are the same as in the original Hanoi network. The third test problem is the Balerma irrigation network (N₃) which was presented for the first time by Reca & Martinez (2006). This network has a total of 443 demand nodes supplied by four source nodes. There are 454 pipes, arranged in eight loops. The minimum required pressure head is 20 m for each node. The last problem (N₄) was taken from a town in the southeast of China which has 237 pipes, one reservoir, and 192 demand nodes. The head provided by the reservoir is 65 m. The minimum pressure requirement for each demand node is 18 m. This network was originally presented by Zheng et al. (2013b). Here, the optimization process is terminated after a fixed number of function evaluations, which is 60,000 for N₁, 90,000 for N₂, and 5,000,000 for N₃ and N₄ case studies.

As demonstrated in Figure 4, in terms of the best and average cost solutions based on 10 runs with different starting random number seeds, the CSHS with a larger N and HMS performed better for each case study. However, the number of evaluations required to find optimal solutions with a larger N and HMS are increased remarkably. Considering both the solution quality and computational efficiency, the CSHS obtained good results when {N = 30, HMS = 15}, {N = 300, HMS = 150}, and {N = 200, HMS = 100} were selected for the N₁, N₃, and N₄ case studies, respectively. The best parameters obtained for the Hanoi network was used for the double Hanoi problem. By comparing the numerical investigations (Figure 4) and the selected {N, HMS} for each case study, one approximate heuristic guideline for practical implementations can be for the CSHS algorithm applied to a WDS optimization problem.

Table 2 reports the statistical results and the required number of evaluations to find the best cost by the present algorithm and some selected meta-heuristic algorithms from the literature. The proposed CSHS found the best feasible solution of $6.0811 × 10⁶ after 31,800 function evaluations. As shown in Table 2, the best convergence speed belongs to MBA and IMBA algorithms with 22,450 and 16,400 function evaluations, respectively; however, the average and maximum costs obtained by these methods were higher than that produced by the CSHS. Comparison with other advanced and hybrid algorithms such as HD-DDS, PSO-DE, and NLP-DE demonstrates the effectiveness and efficiency of the proposed method in solving this problem.

In Equation (1), the scale of the random search is controlled by multiplying the generated Lévy flight by a scaling factor (α). Thus, a small α may slow down the convergence of CSHS because of the limitation in the exploration of only a small subspace of the whole search space. On the other hand, a large α may cause the solution to scatter around the optima as in a random search. From Figures 7 and 8, it can be concluded that while (in most cases) α= 0.06, it makes the algorithm more effective; otherwise, when α > 0.06, new solutions are created outside of the design domain and thus waste evaluations, and while α < 0.06, considerable increase in iterations are needed to find the optimum solution (Figure 8).

According to Equation (4), the local search is very intensive with about p_a of the search time, while global search takes about 1–p_a of the total search time. When p_a = 0.25, about 75% of the total search time is spent on the global searching which allows the search space to be explored more efficiently on the global scale, and consequently the global optimality can be found with a higher probability (Figure 7). As a result, p_a= 0.25 and α= 0.06 are suitable values for the CSHS algorithm.

For our solution described in the previous example ($6.081 × 10⁶), the reference global optimum solution (C_DH) of the double Hanoi network should be $12.114 × 10⁶. The best results obtained by the CSHS, standard CS, BB–BC (Tahershamsi et al. 2012), GA (Cisty 2002), OptiDesigner (Cisty 2002), and HS (Geem et al. 2002) are summarized in Table 3. Using the optimum cost of the Hanoi network (C_H), the reference global optimum solutions of the double Hanoi network can be calculated as $12.114 × 10⁶ (for CS, GA, and HS), $12.182 × 10⁶ (for OptiDesigner), and $12.400 × 10⁶ (for BB–BC), respectively.

The CSHS found the best feasible solution of $12.347 × 10⁶ after 81,900 evaluations. Therefore, deviation from reference optimal solution for the CSHS algorithm is 1.92%. As can be seen from Table 3, comparison with other algorithms such as CS, GA, and HS reveals that the CSHS algorithm is better in term of closeness to the reference global minimum. In other words, the CSHS algorithm is less likely to be trapped by local optimal solutions.

Table 4 summarizes and compares statistical results and required number of evaluations for convergence of the proposed algorithm and other selected meta-heuristics. As shown in Table 4, the best solution found by the CSHS algorithm for the N₃ case study was €1.988 × 10⁶, which is 3.38% higher than the currently best known solution (€1.923 × 10⁶) reported by Zheng et al. (2011) using an NLP–DE method, but lower than solutions obtained by GA, HS, and GHEST.

As shown in Table 5, the CSHS algorithm is superior to the original GA, SA, HS, and CS algorithms in terms of finding the best solution in the same computational overhead. Considering the solution quality, comparison with other hybrid algorithms such as PSHS and MSATS (mixed SA TS) demonstrates the effectiveness of the proposed method in solving this problem.

The statistics of the results for the last case study are given in Table 6. These include the best and average cost of solutions, and the number of evaluations required to find the best cost solution. For comparison, Table 6 also lists the results of other optimization techniques that have previously been used to optimize this case study including GA, DE, and DBA.

It is seen from Table 6 that the proposed CSHS gives a best design cost of $11.766 × 10⁶ which is 5.86% lower than the best design cost obtained by the standard CS. Furthermore, the number of evaluations performed using the CSHS algorithm is only 3 million which is less than those of the GA (i.e., 4,546,000), SDE (i.e., 4,730,200), DBA (i.e., 3,215,685), and standard CS (i.e., 4,000,000) algorithms. In addition, the average cost solution generated by the proposed method is $11.80 million, which is 2.37 and 3.56% higher than the average cost solution of the DE and DBA while 1.02 and 6.27% lower than the average cost solutions of the GA and standard CS, respectively. Thus, for this problem the best cost optimized by CSHS is consistent with the literature but the average cost is slightly high.

In this paper, the use and efficiency of the CS algorithm are investigated in the context of WDSs optimization. It is found that a standard CS algorithm is sometimes unable to produce acceptable solutions to WDS optimization problems. To improve the performance of the standard CS, a new hybrid algorithm, namely CSHS, based on the combined concepts of the CS and HS, is proposed to solve design problems of WDSs. The detailed implementation procedure of this hybrid meta-heuristic is also presented. By incorporating HS concept into standard CS, it is intended to enhance the efficiency and global convergence behavior of the standard CS.

For the CSHS algorithm, the improvements include utilizing the memory (HM) that contains information extracted online during a search and the addition of the pitch adjusting operation of HS during evolution process. Here, the HS strategy is used for fine tuning of the best solution obtained by a standard CS algorithm. In fact, the HM vectors become the CS population and then the evolving process is performed as the usual CS procedure. Obviously, the HM is a pool of the elite solutions and plays a key role in the algorithm. Another improvement is applying a self-adaptive strategy for the dynamically adaption of the control parameters (i.e., HMCR and PAR) using a learning mechanism. These unique features work in a combined framework and can ensure the efficiency of the proposed algorithm. Also, a sensitivity analysis is performed for the CSHS algorithm parameters where the population size, scaling factor (α), and discovering probability of alien eggs (p_a) are concerned.

The performance of the CSHS is demonstrated using four well-known WDS case studies and the results are compared to that of standard CS and previously applied optimization methods. For almost all case studies, CSHS is shown to outperform standard CS both in terms of the ability to find the minimum solution and computational efficiency.

The current best-known solution for the N₁ case study is $6.081 × 10⁶ which was found using the GENOME by Reca & Martinez (2006). This solution has also been found by the proposed CSHS method. As can be observed from the results of the N₁ case study (Table 2), CS and CSHS exhibit similar performance in terms of finding the best-known solution, while CSHS was found to show better performance in terms of computational efficiency in comparison with other standard and hybrid algorithms (e.g., DE, GHEST, NLP–DE, and PSO–DE). It should be noted that the best convergence speed belongs to MBA and IMBA algorithms (Sadollah et al. 2014); however, the average and worst costs obtained by these algorithms were slightly higher than that obtained by the proposed CSHS. Moreover, in order to determine the ability of finding solutions close to the global minimum, the proposed hybrid method has been tested on the N₂ case study. As can be seen from Table 3, deviation of the CSHS from global optimal solution is 1.92% while it is 6.23% for the original CS which indicates the greater ability of the CSHS to effectively explore the search space by incorporating the HS operators. Finally, the CSHS performance can be valuable in the large-scale optimization problems as evidenced by results on the optimization of the last two real-world and complex WDSs (N₃ and N₄ case studies), where the new algorithm is clearly competitive with the other advanced optimization methods. Our results (Table 5) confirm that in solving the N₃ benchmark problem, with a similar computational budget, CSHS performs better than almost all algorithms and for the N₄ benchmark problem the best cost obtained by CSHS is consistent with the literature but the average cost is slightly high (Table 6). One of the main challenges that exists in the use of meta-heuristics is that these algorithms are not adequately tested and no solid conclusions can be drawn on their searching behavior. Therefore, the proposed hybrid method would need to be further examined in order to determine its behavioral characteristics, including (1) the effectiveness of the search in finding optimal or near-optimal solutions, (2) convergence properties, and (3) the amount of the search effort spent in the feasible and infeasible regions of the solution space.

The first author gratefully acknowledges the facilities provided by the Global Institute for Water Security (GIWS), University of Saskatchewan, to carry out this research project.