Simulation-optimization Model for Design of Water Distribution System using Ant Based Algorithms

: In this paper, the Ant Colony Optimization Algorithm (ACOA) is applied to solve Water Distribution System design optimization problem proposing two different methods. Considering pipe diameters as decision variables of the problem, Ant System and Max-Min Ant System, referred to ACOA1 and ACOA2 respectively, are applied to determine pipe diameters. In proposed methods, the ant-based models are interfaced with EPANET as simulator for the hydraulic analysis. Three benchmark test examples are solved with proposed methods and the results are presented and compared with those obtained with other existing methods. The results show the capability of the proposed methods to optimally solve the design optimization problem in which best results are obtained with ACOA2 in comparison with other available results. Furthermore, the results show the superiority of the proposed ACOA2 over than the ACOA1 in which the trade-off between the two contradictory search characteristic of exploration and exploitation is managed better by using Max-Min Ant System.


Introduction
A Water Distribution System (WDS) is a manmade system that is used to convey water from the source to the user. The main aim of every WDS is to deliver water to the user when it is necessary, in the correct quantity and in accordance with relevant quality standards. WDSs are therefore considered to be essential infrastructure in every modern community. As WDSs continue to get aged and cities continue to grow, the design of new WDSs and the rehabilitation or upgrading of existing WDSs will continue to be a central problem.
Construction of WDS is very expensive due to cost of WDS components. A quite small change in the WDS components, therefore, leads to an enormous saving. The general WDS design optimization problem, therefore, involves minimizing the system cost subject to hydraulic constraints such as meeting minimum allowable pressure and/or maximum allowable velocity constraints under design demand levels. A significant amount of research has focused on the optimal design, upgrade or rehabilitation of WDSs. In general, the available methods can be classified as 1) Linear programming (LP), 2) Nonlinear programming (NLP), 3) Dynamic programming (DP) and 4) Evolutionary Algorithm (EA) such as Genetic Algorithm (GA), Simulated Annealing (SA), Shuffled Frog Leaping Algorithm (SFLA), Tabu Search (TS), Ant Colony Optimization Algorithm (ACOA), Harmony Search (HS), Particle Swarm Optimization (PSO), Cross Entropy (CE), Honey-Bee Mating Optimization (HBMO) and Differential Evolution (DE) in which they were usually used for optimal design of WDS. Haghighi et al. (2011) listed different methods applied to design of WDS such as enumeration (Gessler 1985;LP (Bai et al. 2007; NLP Xu and Goulter 1999;DP Schaake and Lai 1969; GA Jian and Yanbing 2010;Bi et al. 2015;SA Costa et al. 2000;SFLA Baoyu et al. 2011;TS Lippai et al. 1999;ACOA Afshar 2007;HS Yang et al. 2012); PSO Babu and Vijayalakshmi 2013; HBMO Mohan and Babu 2010; CE Shibua and Janga Reddya 2012; DE Dong et al. 2012 andother EA andhybrid methods Zhou et al. 2016;Sheikholeslami et al. 2016) in which some of them are presented in this paper. Each of these optimization methods has its own limitations. For example, calculating derivatives or requirement of an initial policy to start off the solution process are some important limitations of the conventional mathematical optimization methods (Mays and Tung 1992). In addition, DP is known to suffer from the "curse of dimensionally" when large scale problems are attempted to solve.
Nowadays, EAs, therefore, have originally been proposed to overcome the limitations of the conventional optimization methods. However, some limitations such as the element of randomness used in generation of initial solutions and generating alternative solutions inherent in most of EAs may reduce their efficiency. In many large scale cases, generally, EAs require additional amount of computer memory which is based on the fact that they need more objective function evaluations than conventional optimization methods. In addition, some EAs such as ACOA have discrete nature and some others have continuous nature. Therefore, ACOA is suitable algorithm for solving discrete problem such as WDS design optimization problem considered here.
ACOA is one of the EAs that has specific characteristic for solving NP-hard combinatorial discrete optimization problems with reasonable computational time cost. For WDS design optimization problem, the pipe diameters are discrete variables in design process and, therefore continuing optimization methods cannot determine them properly. The specific characteristic of ACOA is that in ACOA each artificial ant incrementally builds a solution by adding opportunely defined solution components to a partial solution under construction. This specific unique feature, namely incremental solution building capability, is very useful for solving optimization problems of sequential nature. This algorithm has been inspired in the real ant colony behaviour for finding the shortest path from food source to the nest. Starting with Ant System (AS) (Colorni et al. 1991), a number of algorithmic approaches based on this characteristic of real ant was developed and applied with considerable success to a variety of combinatorial optimization problems from academic as well as from real world applications. Many other algorithms have been proposed to improve the performance of AS, such as Ant Colony System (ACS), Elitist Ant System (ASelite), Elitist-Rank Ant System (ASrank) and Max-Min Ant System (MMAS). Advancements have been made on the AS to improve the operation of the decision policy and the manner in which the policy incorporates new information to help explore the search space of the problem. These developments have primarily been aimed at managing the trade-off between the two contradictory search characteristic of exploration and exploitation. Ali et al. 2009;Ostfeld 2011;Stutzle et al. 2011;and Chandar Mohan and Baskaran 2012) reviewed an amount of research work in the field of using ACOA to solve different optimization problems in the last 20 years. A review of the literature reveals that studies on the application of ACOA to different optimization problems has not been satisfactorily fruitful and is worth further investigation. Numerous papers on the ACOA application for optimization of water engineering problems such as sewer network design (Moeini and Afshar 2013a), WDS design (Zeccchin et al. 2007) and reservoir operation (Moeini and Afshar, 2013b) have been published in recent years.
In this paper, one of the EAs with discrete nature means ACOA is used to effectively solve one of the most important discrete problem means optimal design of WDS by proposing two different methods. In the proposed methods, the ACOA is used to determine pipe diameters as the discrete decision variables of the problem. AS and MMAS referred to ACOA1 and ACOA2 respectively, are applied to determine pipe diameters. MMAS is applied here to overcome the premature convergence problem to suboptimal solution which is well trained in AS by managing the trade-off between the two inconsistent search characteristic of exploration and exploitation. By determining the pipe diameters of WDS, a hydraulic analysis is required. For the hydraulic analysis, in this paper, the ACOA-based models are interfaced with EPANET as simulator. Three benchmark test examples are solved here with proposed methods and the results are presented and compared with those obtained with other existing methods. It is worth noting that many different EAs have been used to solve this problem.
However, based on the discrete nature of design of WDS, here, ACOA is used to solve this problem in order to overcome the limitations of pervious works. In some researches, the hydraulic analysis is based on the simplified assumptions such as considering inaccurate value for hydraulic value coefficients equations. However, EPANET simulator is used here for the hydraulic analysis. In addition, in some researches the continues nature of EAs are used to solve discrete nature problem of WDS design and in some other researches the equations of simulation model are solved using iterative numerical or traditional methods and therefore some of the WDS components such as pump cannot be easily modelled. Here, this fact is regarded considering third test example in which the WDS is fed by pump. The efficiency of these innovations is fully presented and highlighted in section "Model application and results" by comparison of the results obtained with proposed methods with other available results.

Ant Colony Optimization Algorithm (ACOA).
The ACOA is based on the natural foraging behaviour of a colony of real ants to find the shortest path between food source and their nest. When ants are travelling, they deposit a substance called pheromone which is more attractive for other ants to follow them. Ants can smell pheromone and they incline to choose probability paths marked by strong pheromone concentrations when choosing their path. It has been shown experimentally that this pheromone trail following behaviour can give rise to the emergence of the shortest paths once employed by a colony of ants. Based on the real ant behaviour, the ACOA was developed (Colorni et al., 1991). The original and simplest version of ACOA is AS. Later, many other algorithms such as MMAS have been developed to overcome the limitations and improve the performance of the AS. Here, AS and MMAS are used to solve WDS design optimization problem proposing two different methods. Definition of proper graph is the first step of process for solving an arbitrary optimization problem using ACOA. A graph G=(DP,OP,CO), therefore, can be defined, where DP=   I dp ,...., 2 dp , 1 dp is the set of decision points,

Ant System (AS)
The original and simplest version of ACOA is AS. The basic steps of solving optimization problem using AS can be defined as follows (Colorni et al. 1991): 1. At the start of the process, some proper values are initialized to the amount of pheromone trail on all options ij op .
2. A colony with the size of M is selected and at one time, each of them is placed randomly on the decision points of the problem. 3. By using a transition rule for selection proper option at each decision point i, (Eq. 1), each arbitrary ant, m, starts its movement from one decision point to the next and the solutions are incrementally constructed. This procedure is repeated until all decision points of the problem are covered.
is the probability that the ant m selects option j of the ith decision point, M, at each iteration t, the pheromone is updated using the general form of Eq. 2 for the pheromone updating rule.
represents the change in pheromone concentration associated with option ij op which is defined as Eqn. 3: is the cost of the solution produced by the ant m and R is pheromone reward factor.
6. Steps 2 to 5 are continued until convergence criterion is met.

Max-Min Ant System (MMAS)
The most important issue that can be experienced by all ACOAs is premature convergence to suboptimal solutions. To overcome this problem, the MMAS was developed by Stutzle and Hoos (2000). The basis of MMAS is the provision of dynamically evolving bounds on the pheromone trail best In Eqs. 4 and 5, ) t ( min  and ) t ( max  represent the lower and upper limit on the pheromone trail strength at iteration t, is the cost of the solution produced by the best ant at each iteration, best p is the probability that the best solution is constructed again; Javg represents the average number of options available at decision points of the problem; and I is the number of decision points.
In MMAS, only the cost of iteration-best ant is used for calculating the change of pheromone In Eqn. 6, best ) ( f  is the cost of the solution produced by the best ant. It is worth noting that application of MMAS to some benchmark optimization problems has shown that MMAS overcomes the stagnation problem, hence improves the performance of the ACOAs (Afshar and Moeini 2008).

Water Distribution System Design Optimization Problem
WDSs consist of interconnected elements such as pipes, tanks, pumps and other components. Such systems are used to transfer treated water from one or more sources to users spread over a wide area. Designing an effective WDS with minimum cost is a quite complex task which can be achieved by formulating and solving it as an optimization problem. Optimal WDS design requires optimal determination of WDS components such as pipe diameters, tank and pumping station positions and heights leading to a complex optimization problem with a large number of hydraulic constraints. Having a suitable and cost-effective WDS is normally interpreted as finding the solution for this problem that minimizes total cost without violating the constraints.
The standard description form of a WDS design optimization problem with a prespecified layout can be presented as follows. The objective function of this problem can be formulated as Eqn. 7 in which it is minimizing the total construction cost of WDS.
In Eqn. 7, T C is the total cost of the pipes in the WDS; l L is length of the lth pipe; l c is per unit cost of the lth pipe and N is the total number of existing pipes. The cost function is to be minimized under the following constraints.

Continuity equation:
A continuity constraint should be satisfied for each node as follows: In Eqn. 8, Qk is the required demand at consumption node k; l q is the flow rate in pipe l; ) k ( in is the list of the pipes incoming to node k; ) (k out is the list of the pipes out coming from node k; and K is the total number of exiting nodes in the WDS.
Energy equation: The energy constraint for each loop in the network of WDS can be written as follows: In Eqn. 9, l J is the head loss in the lth pipe and P is the total number of loops in the WDS.
Several equations can be used to evaluate the hydraulic parameters of WDS such as head losses in the pipes. Here, Hazen-Williams equation is used to find hydraulic parameters. Therefore the head losses can be calculated as follows: In Eqn. 13, l d is the diameter of pipe l and D denotes the set of commercially available pipe diameters.
It should be noted that the penalty method is often used to effectively guide solutions of optimization problems from one infeasible solution area to a feasible solution one. Here, the penalty method is applied for formulation of the WDS as an unconstrained optimization problem. In the used method, therefore, nodal pressure head and flow velocity constraints are added in the original objective function. Therefore, a new problem can be defined by the minimization of the penalized objective function of Eq. 14 subject to the constraints defined in Eqns. (8) to (13): In Eqn. 14, p C is the penalized total cost function of WDS; p  is the penalty parameter with a large enough value to ensure that any infeasible solution will have a higher total cost than any feasible solution; v CSV is a measure of the flow velocity constraint violation of the trial solution and H CSV is a measure of the nodal pressure head constraint violation of the trial solution. It should be noted that in calculating the CSV, the summation ranges over those nodes and pipes at which a violation of constraints (11) and (12) occurs, that is the terms in parentheses are positive. Other hydraulic constraints are satisfied among using simulation program (EPANET). In other words, EPANET satisfies the continuity and energy conservation equations while calculating the flow rate in each pipe and the pressure head at each node. Furthermore, here, the flow velocity and the nodal pressure head constraint violations have the same penalty parameter value due to the fact that these constrains are normalized.

Methodology
The WDS design optimization problem is solved here by proposing the simulationoptimization approach. In the proposed approach, the EPANET software is used as a simulator and ACOA is used as an optimization method with the goal of finding pipe diameters.
Here, the optimization model is coded in MATLAB and linked with EPANET.
In the proposed approach, at first, a proper graph should be defined to formulate the WDS design as an optimization problem using ACOA. This graph consists of a set of nodes referred to as decision points and edges referred to as options available at each decision points.
Here the decision variables of the problem are the pipe diameters of the WDS which are determined in ACOA1 and ACOA2 methods using AS and MMAS, respectively.
In the proposed methods, the pipes of WDS are considered as decision points. The options available at each decision point are, therefore, defined by a finite number of commercially available pipe diameters. Figure 1 represents the problem graph defined for ACOA1 or ACOA2, where vertical lines show the decision points, dashed small lines show the components of commercially available pipe diameters (j=1,…..J) at each decision point i (i=1,…..I), and finally the bold small lines show a trial solution on the graph constructed by an arbitrary ant.
By determining the pipe diameters of WDS, the ACOA-based models are interfaced with EPANET for the hydraulic analysis. For both methods, Eqn. 14 is used to calculate the total cost of the trial solutions generated. Finally, optimal pipe diameters are determined using proposed methods. Figures 2 and 3 illustrate the process of solution constructions in ACOA1 and ACOA2 formulations, respectively.
It is worth noting that the different hydraulic softwares have been developed for pipe network analysis such as EPANET, Water CAD, Water GEMS, MIKEURBAN and so on. Each of these software's has its advantages and limitations. Due to more advantages of EPANET software, this software is used in this paper. In other words, EPANET is freeware, user friendly and free source software that everybody can simply use it. Due to these noted facts and its accuracy and availability, EPANET is so popular in the field of water distribution analysis and therefore it is used here. EPANET water distribution system simulator is the most commercial hydraulic model which has been developed by the US Environment Protection Agency. The EPANET performs extendedperiod simulation of hydraulic and water quality behaviour within pressurized pipe networks. This well-documented and tested public domain simulator provides a convenient platform for implementing the simulationoptimization approach (Rossman 2000). Here, the EPANET software is linked to MATLAB with the help of EPANET toolkit. EPANET toolkit is a Dynamic Link Library (DLL) of functions and an extension of EPANET simulation package that allows developers to customize EPANET software accordingly (Rossman 1999). If other pipe network analysis software can be linked to MATLAB, these software can be easily considered as simulators.

Model Applications and Results
In order to show the ability and performance of the proposed methods for the optimal design of WDS, three benchmark examples are applied here with different scales. Therefore, small (test example I), medium (test example II) and large scale (text example III) test examples are presen-ted and solved here using proposed methods. Furthermore, in the third test example, the WDS is fed by a pump; however, in the first and second test examples the WDS is fed by gravity. The first test example (I) as shown in Fig. 4 is the two-loop WDS (Alperovits and Shamir 1977). This WDS network has 7 nodes, 8 pipes and two loops in which it is fed by gravity from a reservoir with a 689 (ft) fixed head. All WDS network pipes lengths are 3281 (ft) and their Hazen-Williams coefficient are 130. The minimum nodal pressure head is 98.4 (ft). Fourteen commercial candidate pipe diameters for each pipe with corresponding cost values are listed in Table 1. Details of the network including nodal ground elevations and water demand are given in Table 2.
A set of preliminary runs is first conducted to find the proper values of AS and MMAS parameters referred to ACOA1 and ACOA2 respectively, as shown in Table 5 for all three test examples. For both methods of all test examples, a colony size of 100 with maximum number of 1000 iterations amounting to maximum number of 100,000 function evaluations is used. Table 6 shows the results of 10 runs carried out using different randomly generated initial guess for all test examples along with the scaled standard deviation, the number of final feasible solutions and CPU time for each run. It is clearly seen from Table 6 that all measures of the quality of the final solutions such as the minimum, maximum and average cost and the scaled standard deviation are improved when using ACOA2 compared to ACOA1. It should be noted that, due to the fact that premature convergence does not occur in ACOA2 by providing dynamically evolving bounds on the pheromone trail intensities, the ACOA2 has been able to outperform the ACOA1. In other words, the results show the superiority of the MMAS, in conjunction with the ACOA2, over than AS, in conjunction with the ACOA1, in which the trade-off between the two contradictory search characteristic of exploration and exploitation is managed better using MMAS. In order to show the efficiency of proposed methods, these results are compared with those obtained with other existing methods. Comparison of the results shows that the best results are obtained using ACOA2 with no extra computational time cost due to the unique feature of ACOA. Therefore, the proposed method is useful to solve NP-hard combinatorial discrete optimization problems such as optimal design of WDS.  The comparison of the results shows that the optimal solution of 419000 units is obtained at the expense of 4700 function evaluations using ACOA2. It is worth noting that the ACOA2 has been able to outperform the methods of Abebe and Solomatine (1999) and Savic and Walters (1997). Furthermore, this compares favourably with about 250000 function evaluations required by the Savic and Walters (1997) using GA1, about 53000 function evaluations required by Cuncha and Sousa (1999) using SA, about 100000 function evaluations required by of Prasad and Park (2004) using GA, about 11155 function evaluations required by of Eusuff and Lansey (2003) using SFLA, about 5100 function evaluations required by of Afshar (2007) using ACO and about 7467 function evaluations required by Wu et al. (2001) Savic and Walters (1997) (1999) and Eusuff and Lansey (2003) are infeasible. In other words, the minimum required pressure head is not satisfied at nodes 17 and 19 when obtained diameters are simulated in EPANET software. In addition, a smaller value for constant coefficient of Hazen-Williams equation ) ( leads to smaller head losses and therefore smaller pipe diameters in which these solutions are infeasible when they are simulated in EPANET software. Comparison of the results shows that ACOA2 has been able to outperform the other used methods with no extra computational time cost due to the unique feature of this method which is presented before. In addition, the near optimal solution of Maier et al. (2003), Zecchin et al. (2006) and Afshar (2007), 38.64 (M$), is obtained here using proposed ACOA2 with smaller computational time.
Finally, the results of test example III obtained with the proposed method are compared with some other available results. The ACOA2 is able to find the optimal solution of 175783163 (Won) in just 8600 function evaluations. This compares favourably with those obtained with other existing methods. Kim et al. (1994) solved the problem using a projected Lagrangian algorithm supported by GAMS/MINOS and then converted the continuous diameters to discrete commercial diameters. They obtained 179142700 (Won) while the original cost was 179428600 (Won). Furthermore, the present result can be compared with about 10000 function evaluations required by Geem (2006) to get the solutions of 177135800 (Won) using HS. It is worth noting that the ACOA2 has been able to outperform all existing methods with no extra computational time cost due to the unique presented feature of this proposed method.
Superior performance of the ACOA2 compared to ACOA1 is already attributed to the provision of dynamically evolving bounds on the pheromone trail intensities and methodology of calculating the change of pheromone ij   . This fact reflected in the typical convergence curves shown in Figs. 7, 8 and 9 for test example I, II and III, respectively. It is seen that the population of ACOA2 has solution costs way below that of ACOA1 due to the fact discussed before.

Conclusion
In this paper, the Ant Colony Optimization Algorithm (ACOA) was applied to solve Water Distribution System (WDS) design optimization problem using two different methods named ACOA1 and ACOA2. Pipe diameters of the WDS network were determined in both ACOA1 and ACOA2 methods using Ant System (AS) and Max-Min Ant System (MMAS), respectively. For the hydraulic analysis, here, the ACOA-based models were interfaced with EPANET as simulator. Three benchmark test examples were solved using proposed methods and the results were presented and compared.
While two proposed methods showed good performance in solving the problems under consideration, the ACOA2 was observed to produce better results for WDS design optimization problem and to be less sensitive to the randomly generated initial guess required to start the solution process represented by the scaled standard deviation of the solutions produced in ten different runs. In other words, the results showed the superiority of the ACOA2 over other the ACOA1 in which the trade-off between the two contradictory search characteristic of exploration and exploitation was managed better using MMAS which is in conjunction with the ACOA2. It was also shown that the ACOA2 consistently gave better quality solutions than existing traditional and heuristic search methods with less computational effort due to the unique presented feature of this proposed formulation.