## Abstract

This paper describes an extension to the optimal power use surface (OPUS) methodology, which consists of applying a metaheuristic post-optimization process after a network has been designed through the OPUS algorithm. Four different types of metaheuristics were tested in this study. The OPUS methodology focuses on the setting up of efficient ways in which energy is dissipated and flow is distributed, thereby obtaining deterministic design after a few iterations. On the other hand, the stochastic techniques used for the post-optimization step mimic different phenomena, usually requiring several iterations. The proposed methodology is tested on four networks; three of these are benchmark problems. When compared to the results obtained through other methodologies, this algorithm stands out for allowing designs with constructive costs very close to the lowest found in other investigations. However, when compared to the results obtained through the OPUS algorithm alone, this approach reduces the capital cost of the network by a very small percentage while increasing the number of iterations required. This extension to the OPUS methodology proves that following hydraulic principles allows near-optimal results to be obtained, whose improvement demands a considerable number of iterations, providing minimum benefit as the reduction in cost is only 1% at most.

- hot start
- metaheuristics
- power use surface
- water distribution system (WDS) optimal design

## ABBREVIATIONS

- α
- Exponent of the cost function
- DDP
- Deficit pressure penalty of GA
- Cost of the configuration obtained using a sag of value
*x* - D
_{k} - Decision criteria of Greedy A
- Minimum diameter available at commercial list
- F
- Sag for the HGL quadratic equation
- GA
- Genetic algorithms
- H
- Total head loss
- HGL
_{i} - Hydraulic gradient line at node
*i* - Friction head loss
- HM
- Harmony memory
- HMS
- Harmony memory size
- HS
- Harmony search
- ILP
- Integer linear programming
- K
- Coefficient of the cost function
- L
- Pipe length
- m
_{*} - Probability assigned to a Greedy A criteria
- NG
- Number of generations of GA
- NN
- Number of nodes in the network
- OPUS
- Optimal power use surface
*P*_{cross}- Crossover probability of GA
*P*_{mut}- Mutation probability of GA
- PS
- Population size criteria of GA
- Q
- Flow
- RI
- Resilience index
- SA
- Simulated annealing
*S*_{f}- Friction slope
*T*_{cr}- Factor termed cooling rate at SA
*T*_{ini}- Initial temperature at SA
*T*_{stop}- Threshold at SA
- UC
- Unit capacity
- W
_{*} - Difference in terms of that criteria from the actual network designed of Greedy A
- WDS
- Water distribution systems
- X
- Topological distance between the node and the reservoir

## INTRODUCTION

The problem of optimal design of water distribution systems (WDS) has become much more important in recent decades. This is due to the limitation of funds to solve this issue, the fact that water supply is essential for human life, and the rapid increase in demand as the world population is growing all the time at a higher rate. For over 40 years, many researchers have studied the problem, suggesting several methodologies that offer viable solutions to this difficulty.

One of the most common criterion to compare and validate different design methodologies is the construction cost, but this does not mean that other criteria such as reliability, environmental impact, and water quality should not be taken into account. This type of design consists of determining the set of pipe diameter sizes that offer the minimum capital cost, satisfying flow demands with an adequate pressure. Notwithstanding the fact that pipes are usually manufactured in discrete diameter sizes, the amount of possible pipe configurations is immense, which makes the problem highly indeterminate. Yates *et al*. (1984) proved that it is a NP-HARD problem, which implies that only approximate methods could be successful in finding adequate solutions.

Opening approximations involved traditional optimization techniques such as enumeration, linear, and non-linear programming. Later, different metaheuristic algorithms started to gain popularity due to their ease of implementation and additional advantages such as their broader search of the solution space, a relatively small reliance on the system's initial configuration, and their capability to incorporate the discrete-sized diameter restrictions. Successful attempts include genetic algorithms (GAs) by Savic & Walters (1997), harmony search (HS) by Geem (2006), scatter search by Lin *et al*. (2007), cross entropy by Perelman & Ostfeld (2007), simulated annealing (SA) by Reca *et al*. (2007), and particle swarm by Geem (2009).

These metaheuristics consist of phenomenon-mimicking or evolutionary algorithms that randomly generate a large number of possible solutions and test their fitness in terms of quality and capital costs. The results are progressively improved due to the use of generic learning functions. In the WDS design context, each solution represents an alternative design, namely a different set of pipe diameter sizes. Every time an alternative design is generated, a static hydraulic simulation must be performed in order to evaluate its feasibility, thus a large number of iterations is required to reach convergence. This makes metaheuristics very demanding in terms of computational effort, regardless of their flexibility and their capability of accomplishing near-optimal results. Due to this, besides the cost of the final solution, the number of hydraulic simulations (or iterations) must be taken into account to measure and compare the efficiency of the different methodologies. Even though the learning functions used in metaheuristic algorithms involve testing the hydraulic performance of each of the candidate solutions, they do not make use of additional hydraulic criteria.

As a response to these stochastic algorithms, more recently some researchers have suggested new hydraulic approaches based on the use of available energy along the system. Metaheuristics intend to optimize an objective function treating the optimization variables simply as a series of numbers that must follow certain logic, without understanding the machinery behind that logic; however, these new approaches try to characterize the behavior of the different hydraulic variables and understand the underlying dynamics, focusing on the optimal distribution of the power used throughout the network.

In this sense, Wu (1975) analyzed the behavior of the function that was intended to be minimized, which is the sum of each individual pipe cost. This analysis was carried out for simple systems, namely systems composed of series of pipes with a known uniform demand (i.e., demand-driven model). As a result, Wu found that the total head distribution that minimizes constructive costs follows a quadratic curve, which is concave upward and separated from the straight line that connects the hydraulic grade line (HGL) at the start and end nodes, a value of 15% of the total head-loss (H). Thus, optimal designs could be obtained by computing objective head loss values for each pipe derived from the HGL fabricated using Wu's criterion.

The extension of Wu's basic concept was proposed by Featherstone & El-Jumaily (1983), suggesting its applicability to looped networks which are more complex network systems than those studied by Wu. This concept has been used on demand-driven models with open and closed topologies (Ilemobade & Stephenson 2006), and more recently on open systems with pressure-driven demands (Páez *et al*. 2013). The optimal power use surface (OPUS) methodology, introduced by Takahashi *et al*. (2010), proposes a net hydraulic approach that follows the aforementioned criteria and proves that hydraulic criteria could be used as the basis of WDS design in order to replace the iteration-intensive stochastic approach required by metaheuristics. Furthermore, the application of these hydraulic principles along with linear programming formulations integer linear programming (ILP) has been used in other studies, as presented in Saldarriaga *et al*. (2012). The results achieved using these methodologies are outstanding not only in terms of the objective function, proving to reach near-optimal solutions with small differences compared with the global records, but also in terms of the computational effort required, which is always several orders of magnitude less than most of the comparable solutions reached through different approaches. Additionally, this approach offers a clearer insight into the inner mechanics that govern WDS design, which makes it easy to understand and very versatile to implement, allowing its coupling with tools such as ILP in order to accelerate the process.

Taking these results into account and making a detailed analysis of the differences between the models obtained by means of optimal power use criteria and global optimal designs, it was found that these differences are a consequence of the commercial diameter restriction. This restriction usually poses difficulties in WDS design as hydraulically based methodologies such as OPUS have intermediate steps that compute continuous theoretically optimal designs, which are then completely affected by the rounding of diameters. This happens because the hydraulics of the system may change drastically after modifications in the pipe diameters. In this sense, this research implements designs obtained in the final stage of the OPUS methodology as a ‘hot start’ for metaheuristic methodologies typically used in the design of WDS, such as GA, HS and SA, and an optimization technique as a greedy algorithm. This was first proposed by Saldarriaga *et al*. (2014), and now is extended with the objective of improving the results already obtained with hydraulic criteria jointly with stochastic algorithms. The results obtained in this new stage of the development of methodologies based on the energy use for WDS design are presented in comparative tables that portray the obtained costs, the hydraulic executions (measurement of the computational effort) to obtain such costs, and the resilience index of the solution designs. The proposed methodology is tested on different benchmark problems (Hanoi, Balerma, and Taichung) and a fourth network named R28, giving results very close to the global records, but with an insignificant difference regarding the OPUS methodology. As well, the number of iterations increases substantially with respect to the OPUS methodology, which implies that improving hydraulically based designs requires great additional computational effort for a minimum reduction in costs.

## OPUS METHODOLOGY

The OPUS design methodology is composed of six different steps. These are presented in Figure 1 in the order in which they must be executed. Each of the processes is described below, and the detailed algorithm that each of them follows is explained in the work developed by Saldarriaga *et al*. (2012).

### Sump search or tree structure

This step is based on two principles: the first one is that least-cost WDS should supply water to each demand node using a single route from the water sources; the second one states that as a pipe's design flow increases, its marginal cost decreases. The first principle comes from the fact that redundancy, although favoring reliability, is hydraulically inefficient and therefore open WDS are much cheaper than looped networks. The second principle follows from the flow expression derived from the Darcy–Weisbach and Colebrook–White equations. Leaving all the other parameters constant, the flow presents a relation approximately proportional to the diameter to a power of 2.6. Assuming a standard pipe cost equation and replacing the diameter according to this proportion, the cost of a pipe is defined as a function of its design flow (see Equation (1)). The behavior of this new relation is depicted in Figure 2, where it can be seen that as the design flow for a pipe increases, the marginal cost decreases.
1where is the cost of the pipe when it is transporting a flow *Q*, *K* is the coefficient of the cost function of the network, *x* is the exponent of the cost function of the network, and is the friction slope of the network.

Following these principles, an algorithm was designed in order to obtain the tree structure. This was done by aggregating flow values in the least number of main routes possible. The open network is set up starting from the water sources and then adding adjacent pipe–node pairs, one at a time. The group of available pairs in each iteration is termed the ‘search front’ and each of these pairs is assigned a cost–benefit value ), creating a recursive process.

The pair in the search front with the highest cost–benefit value is incorporated into the tree structure. The cost–benefit function of a pair is calculated by computing the quotient between the demand of the new node and the marginal cost of connecting it to the source. This entails the addition of the total cost of the pair's pipe to the cost difference of transporting the additional flow through all of the upstream pipes (see Equation (2)).
2where is the cost of adding the adjacent pipe–node pair *i*. is the flow demand of node *i*, namely the additional flow that will be transported through the tree structure if pair *i* was added to the tree. is the flow that goes through pipe *j* of the tree before adding any adjacent pipe–node pairs. The first term in Equation (2) represents the cost of the adjacent pipe and the summation is the marginal cost upstream. This equation must be solved together with Equation (1). It can be proved that the construction of the tree using this cost–benefit function has an time complexity, where is the number of nodes. The corresponding pseudo-code is presented in Figure 3.

In Figure 3, the is an object composed of nodes each of which has a list of neighbors associated. The function creates a tree object with only the reservoirs added to the structure. The function receives a tree object and returns a list of neighbors of all nodes from the tree. The is a list of possible nodes that could be added to the tree. The procedure calculates the cost–benefit function of adding the analyzed neighbor to the tree. Finally, the function determines the neighbor with the highest value of the cost–benefit function.

### Optimal power use surface

This step consists of assigning an objective hydraulic head to every node in the network and thus predefining the head losses for each pipe. This sub-process gives the name to the entire methodology as it is here where the I-Pai Wu concept of optimal energy gradient line is applied. By assigning all sump nodes the minimum allowable pressure, and knowing the head at the reservoir, the intermediate nodes' head for each path are calculated with a parabolic HGL, as shown in Figure 4. The sag optimal value can be estimated using a methodology suggested by Ochoa (2009), who found that this value is dependent on the demand distribution, the ratio between flow demands and pipe length, and the cost function. As a result of the sub-process, every node in the network has been assigned an objective head and therefore a design flow is needed in order to calculate the diameter of each pipe in the network.

As the branches converge while going over the open network upstream, the sag must be recalculated at each intersection by weighting the flow on each downstream route. It is important to pay special attention to particularly elevated nodes, to avoid assigning a head value that does not meet the pressure requirement. The HGL at every node is computed following Equation (3):
3where *x* is the topological distance of the node, *F* is the sag at the analyzed node, is the HGL at the source, is the minimum HGL admissible at the analyzed node, and is the total distance from the source to the sump of the route where the analyzed node is located.

In order to obtain the general optimum sag of a network, three different continuous designs are executed and their costs are calculated. Each of the three configurations is obtained by considering the exact same criteria except for the value of the general sag. Subsequently, the optimum sag is estimated by adjusting the three points obtained (cost vs. sag) to a second-order polynomial regression, and finding the sag that corresponds to the minimum point (least cost) of the parabola. In order to calculate the sag that corresponds to the least-cost design, the following equation is used:
4where is the cost of the configuration obtained using a sag of value *x*.

As can be seen in Equation (4), the sag values used for the three designs are 0, 0.1, and 0.25. A minimum value of 0 was selected because it represents the case of a straight HGL; a negative sag would give as a result a concave down parabolic HGL, counter to Wu's criterion. As for the maximum value of the sag, it is obtained by differentiating the expression that defines the HGL at a specific node of the system as a function of its topological distance to the source (Ochoa 2009) (see Equation (3)). This differential equation is then set equal to zero and the topological distance is replaced by the total length of the system. This procedure produces as a result an expression for the sag value that corresponds to the case in which there are no head losses in the last pipe of the system. The result is a value of 25% as the maximum sag for a system with a last pipe of infinitesimal length; sags larger than this would imply the generation of energy at some points of the system. Even though real systems do not have pipes of infinitesimal length, the length of the last pipe is negligible compared to the total system's length; due to this, the maximum sag of any real size network will be a number very close to 0.25. Finally, a third sag value of 0.1 is arbitrarily selected in order to obtain a third design.

### Optimal flow distribution

Taking into account that in a looped network the same hydraulic gradient surface can be obtained through many different diameter configurations when the set of allowable pipe sizes is (positive real numbers) (Takahashi *et al*. 2010), it is necessary to predefine also the objective flow for each pipe in order to obtain a diameter configuration that minimizes costs. Therefore, this sub-process finds a unique flow distribution scheme that respects mass conservation and adjusts to the optimal power use surface previously established. In this case, the process is executed using the original graph instead of the spanning tree, as in the previous step.

Starting from the sumps, the flow demand is divided into the upstream pipes according to one of the following criteria: (1) uniform distribution: the total flow demand of each node is divided into the number of upstream pipes connected to it, which means that all pipes are assigned the same flow; (2) proportional distribution: for each pipe, the flow is distributed proportionally to its hydraulic favorability; and (3) least-cost distribution: only one pipe will be assigned the largest flow value, while the others will be assigned the flow that corresponds to the minimum diameter available . Depending on the criterion used, the reliability of the network varies. It was demonstrated that the uniform distribution favors the reliability of the network in comparison to the other two criteria (Takahashi *et al*. 2010). In this study only the third criterion will be used as it has been shown that the most cost-efficient way to transport a fluid is by using as few routes as possible (Takahashi *et al*. 2010).

In order to determine the principal pipe (i.e., the pipe that will have the largest portion of the total flow demanded), several criteria can be used to evaluate their fitness. Similarly to the tree structure step, a function is used to determine the pipe with the highest flow. Namely, the pipe with the maximum value of hydraulic favorability, defined according to Equation (5), is the one selected. For non-sump nodes, the total demand is calculated by adding its own flow demand and the flow demanded downstream. An iterative-recursive algorithm (IRA) can be used to perform all of the calculations with an time complexity
5where is the HGL of the upstream node of pipe *i*, is the HGL of the downstream node of pipe *i*, and *L* is the length of pipe *i*. At the end of this step, every pipe in the system must have been assigned a design flow. The pseudo-code of this procedure is presented in Figure 5.

The function getChildren() for the type Node returns the list of nodes connected to it through its downstream pipes. The function assignFlows() assigns flow values to all of the node's upstream pipes taking care to respect the mass conservation principle.

### Diameter calculation

This sub-process calculates continuous diameter sizes for every pipe using the outcome of the previous steps. As the objective head losses and the design flow rate for every pipe are already known, the continuous diameter needed can be easily obtained from a straightforward calculation. This computation is explicit when the Hazen–Williams equation is used and iterative when the Darcy–Weisbach and Colebrook–White equations are employed. The resulting continuous design is a full-operational WDS, with a cost very close to the minimum. Due to the limited availability of diameter sizes, a next step is required to transform this ‘optimal’ design into a feasible one.

### Diameter round-off

This step consists of approximating each continuous diameter to a discrete value from the list of commercially available diameter sizes. It was found that rounding to the nearest equivalent flow value offers the best results, even though it can be done following several criteria. It was found that the best way to accomplish this is approximating to the nearest flow, by elevating the diameter values to a power of 2.6 and rounding these values, as explained in the Tree Structure step. Unfortunately, this step affects drastically the system's hydraulic behavior, especially if all the diameter sizes are rounded up or down.

### Optimization

This final step has two main goals: to ensure that every node has a pressure higher than or equal to the minimum allowable; and to look for possible cost reductions. The first goal is reached by increasing diameters (if required) starting with the ones with larger unit head-loss difference between real and objective values, until the whole system meets the pressure restriction. The second goal is reached through a double sweep starting from the reservoirs going toward the sumps in the direction of the flow, and then backwards. The reduction of each pipe's diameter is considered twice, and the diameter is permanently reduced if no constraint is violated after the change. To make sure the minimum pressure is not being violated numerous hydraulic simulations are required as it is necessary to run a hydraulic simulation per diameter modification of every single pipe. This sole process could be used alone to obtain sound designs; in spite of this, it is strongly dependent on the initial diameter configuration.

The OPUS methodology was tested by Saldarriaga *et al*. (2012) on three well-known benchmark systems: Hanoi, Balerma, and Taichung; as well as on the R28 network which is introduced in this paper. The results reached near-optimal costs with a number of hydraulic executions (iterations) about three orders of magnitude smaller than other approaches. For the Hanoi network, the OPUS methodology reached a cost of €6.147 million after 83 iterations, which is only 1.5% above the global record of €6.056 million obtained through different metaheuristics by Cunha & Sousa (1999), Geem (2002), Vairavamoorthy & Ali (2005), Kadu (2008), among others. Furthermore, the OPUS solution for Hanoi was obtained with a number of iterations three orders of magnitude below. In the case of the Balerma network, a cost of €2.040 million was obtained after 711 iterations. This result presents a difference of 5.1% in comparison to the lowest cost reported in previous work, which is of €1.940 million with 30,000,000 iterations obtained by Tolson *et al.* (2009). The discrepancy in the number of iterations in this case is of five orders of magnitude, which represents a substantial computational effort. As for the Taichung network, the OPUS methodology gave a design with a cost of €8,939,900 requiring 49 iterations. This cost represents a difference of 1.9% when compared with the result previously reported of €8,774,900 obtained through Tabu search by Sung *et al*. (2007).

These results become more valuable when one notices that the OPUS methodology is deterministic and thus an identical design would be found by any user after the same number of iterations, contrary to other algorithms with a stochastic component. This methodology is significantly more efficient than heuristic techniques and the efforts made to understand the hydraulic principles are well rewarded. The solution obtained through OPUS will offer a very good initial configuration to non-deterministic algorithms, reducing the number of iterations required to reach near-optimal results through metaheuristics and providing OPUS a broader search of the solution space. This is the reason why it is suitable to be used as a hot start for phenomenon-mimicking methodologies, with the objective of improving the results reached through the application of hydraulic principles.

## HOT START METHODOLOGY

The hot start methodology introduced herein combines the OPUS hydraulically based algorithm with a metaheuristic one. The set of diameters obtained through OPUS is used as the initial configuration of the network to be post-optimized using a stochastic technique. The metaheuristics tested in this investigation were adapted to the case of WDS design and were implemented using different software: REDES, MATLAB, and GANetXL.

REDES is a WDS modeling program developed at the CIACUA of the Universidad de los Andes, Colombia. MATLAB is a well-known high-level language and interactive environment for numerical computation, visualization, and programming, developed by MathWorks. Finally, GANetXL is an optimization add-in for Microsoft Excel, which uses GA to solve complex optimization and search problems, developed at the Centre of Water Systems of the University of Exeter, UK. The main characteristics of the algorithms included in this work are described below.

### Genetic algorithm REDES

The GA implemented in REDES uses generational reproduction with standard crossover, and roulette wheel selection as the reproduction method. In addition, REDES takes into account a series of constraints to avoid becoming stuck at a local minimum, as well as an algorithm to reduce the stochastic error. This software ranks each configuration obtained and saves a limited number of the best solutions from all the generations. This list of solutions is stored by the shake algorithm every time a new generation is created, so that every new solution can be compared with the solutions in the list, and thus it can be properly updated. The ranking and list storing done by REDES is an inverse function of the cost and the minimal pressure. This increases the number of hydraulic simulations, as opposed to other GAs, as every time it compares a new solution with the previous ones, it executes the hydraulics again for the latter. It is fundamental to define the population size (PS), the number of generations (NG), the mutation probability (*P*_{mut}), the probability adjustment constant (PAC), which defines a constraint for the number of descendants for the dominant solution, and the deficit pressure penalty (DPP), as initial parameters of the algorithm.

### Genetic algorithm GANetXL

The software GANetXL allows the user to choose between a multi-objective GA and a mono-objective one. For this case, the second option is chosen and the function to be optimized corresponds to the cost equation of the network. The reproduction algorithm can be one of the following: generational, where all the individuals in the population are the product of random mutations produced on the individuals in the previous population; elitist selection, which is similar to the earlier copying without any mutation of the number of individuals indicated by the EL (which is chosen by the user) to the next generation; and steady state, where only a few individuals mutate before being copied to the new generation. Finally, the way in which the crossover will be made must be selected from the following three options: simple one point, which assigns diameter values of one parent to half of the pipes and diameter values of the other parent to the remaining pipes; simple multi-point, which divides the network into different sections and randomly assigns the diameter values of one parent to some of the sections and the diameter values of the other parent to the remaining sections; and random, which randomly decides pipe by pipe the parent that will transfer its diameter size to the new individual. The crossover probability is defined by the parameter (*P*_{cross}). Before being assigned, the parents' diameters mutate according to a *P*_{mut} that is defined by the user.

### Genetic algorithm MATLAB

The GA implemented in MATLAB makes an elitist selection generation by generation. This means that the population in one generation does not interact with the population in the following. The only information that is shared between subsequent generations is that best individuals are copied from one generation to the other without being modified. The number of individuals copied is determined by the elite individuals parameter. The hot start can be arranged in two ways: taking as the initial population the exact diameter configuration obtained through OPUS, or assigning an initial range of possible diameters to each pipe in the network. The range used in this case includes three diameters: the original diameter assigned by OPUS, the discrete diameter immediately below in the list of available diameters of the network, and the discrete diameter immediately above in the same list. This option was selected as it prevents the algorithm becoming stuck in a local optimum. The mutation process is performed using a uniform function.

### HS REDES

HS is an evolutionary algorithm that mimics the improvization process followed by musicians to create a ‘fantastic’ harmony. It consists of three basic steps, which are described below.

#### Prepare a harmony memory

The harmony memory (HM) is a matrix that stores the best harmonies. In the WDS design context, each line in the HM is a diameter configuration of the network and each column indicates the discrete diameter that has been selected for each pipe. The harmony memory size (HMS), which is defined by the user, determines the number of configurations that can be held by the HM. Each of the possible designs is rated using the following objective function:
6where is the value of the objective function for the diameter configuration *j*, is the number of pipes in the system, is the linear cost of pipe *i*, is the length of pipe *i*, is the number of nodes that have a head deficit, is the pressure at node *k*, and are the penalty function's parameters, which can take values of the order of magnitude of 1,000 and 100,000, respectively. For the hot start methodology, the HM is initialized with the diameter configuration obtained through the OPUS methodology. Then a series of stochastic configurations are generated, until the HM is full. This is done by randomly assigning to each pipe a diameter close to the one assigned by OPUS.

#### Improvize a new harmony

The new harmonies are selected based on the set of configurations in the HM. For this reason, the diameter of a pipe is randomly selected from the configurations in the HM, in order to create new designs.

#### Update the HM

If, according to the objective function the last harmony generated is better than the worst design in the HM, the previous design is replaced by the new one.

### SA MATLAB

SA is a metaheuristic algorithm introduced by Kirkpatrick *et al*. (1983), which was designed to avoid stopping at local optimums, like traditional hill-climbing (also called hill-descending) approaches do. It consists of a stochastic relaxation technique inspired by the physical process of annealing a metal. In contrast with GAs, standard SA uses a single solution during the optimization process. In SA, better solutions are accepted in every case, whereas the acceptance of worsening solutions depends on a parameter, called temperature, which states the probability of this happening. The initial temperature (*T*_{ini}) decreases in the following integrations at a rate given by the factor termed cooling rate (*T*_{cr}). The temperature is included in the Metropolis function introduced by Metropolis *et al*. (1953), simultaneously controlling the number of iterations of the algorithm and defining the probability for a certain solution to be accepted. The process concludes when the temperature reaches a given threshold, *T*_{stop}, which is usually close to zero. MATLAB allows two different ways to generate new individuals for the next iteration: fast annealing and Boltzmann annealing. The first one generates a new population taking random steps, with size proportional to the temperature. The latter takes random steps with size proportional to the square root of the temperature. The way in which temperature diminishes is restricted to one of the following: exponentially, it decreases ; logarithmic, it decreases ; and linear, the temperature decreases , where refers to the corresponding iteration number.

### Greedy algorithm REDES

Greedy algorithm is a metaheuristic first presented by Andrade *et al*. (2013), which was designed to improve other metaheuristic methods results at real WDS networks. The algorithm is only intended to redefine diameters as a post-optimization objective presented in this work. The algorithm consists of two primary steps: (1) the first one consists of defining a pipes set (S) from all the pipes that after reducing the diameter by one commercial size generate no violation of network constrains; (2) the second step consists of determining the best candidate from S that will be permanently decreased, and it is done using Equation (7) as an evaluation index for each pipe. At REDES this equation depends on four parameters: , pressure , resilience index , and the unit capacity defined as the flow head-loss per unit of length on a pipe:
7where is the probability assigned to each criterion α and *W*_{α} is the difference of the criterion α between the network with pipe *k* decreased and the original network. After calculating the decision criteria value *D _{k}* the best candidate corresponds to the pipe that generates the highest result.

## RESULTS

All the hot start methodology was validated with the aid of three benchmark networks: Hanoi, Balerma, and Taichung; and a fourth network known as R28, which will be introduced below. Each of these WDS was designed using all of the five metaheuristics explained before, once the OPUS methodology was applied.

### Hanoi

The Hanoi network was first presented by Fujiwara & Khang (1990). The head-loss equation commonly used is Hazen–Williams with a roughness coefficient (*C*) of 130, the minimum pressure for the design scenario is 30 m and the pipes' costs can be calculated using a potential function of the diameter with a unit coefficient of $1.1/m and an exponent of 1.5. The pipes commercially available are: 12, 16, 20, 24, 30, and 40 inches. The network is formed by 34 pipes and 31 nodes configured in three loops. The whole system is supplied by one reservoir with a constant head of 100 m. The topology of the network is presented in Figure 6.

In Table 1, all the parameters tested for each of the metaheuristics are presented, as well as the specific values which gave the solutions of minimum cost. The results obtained for this network are presented in Table 2 where the costs obtained through OPUS alone and with the global record are compared.

### Balerma

Balerma corresponds to a WDS of an irrigation district in Almería, Spain. It was first introduced by Reca & Martínez (2006). The head-loss expression commonly used is the Darcy–Weisbach equation. The pipe diameter sizes commercially available for its design are manufactured exclusively in PVC, with an absolute roughness of 0.0025 mm. The minimum pressure allowable is 20 m and it has 10 commercially available pipes whose unit costs are listed in pairs (diameter in mm, cost in €/m): 113, 7.22; 126.6, 9.1; 144.6, 11.92; 162.8, 14.84; 180.8, 18.38; 226.2, 28.6; 285, 45.39; and 361.8, 76.32. It has 454 pipes and 443 consumption nodes which are supplied by four reservoirs. The topology of the network is presented in Figure 7.

In Table 3, all the parameters tested for each of the metaheuristics are presented, as well as the specific values which gave the solutions of minimum cost for Balerma. The results obtained for this network are presented in Table 4 and they are compared with the costs obtained through OPUS alone and the least cost reported for the Balerma network in previous works.

### Taichung

The Taichung network was first presented by Sung *et al*. (2007) and it corresponds to a WDS located in Taichung, Taiwan. The network's topology consists of 12 loops formed by 31 pipes and 20 nodes, which are supplied by a single reservoir with a head of 113.98 m. For its design, there are 13 pipe diameter sizes commercially available, whose unit costs are listed in pairs (diameter in mm, cost in NT dollars/m): 100, 860; 150, 1,160; 200, 1,470; 250, 1,700; 300, 2,080; 350, 2,640; 400, 3,240; 450, 3,810; 500, 4,400; 600, 5,580; 700, 8,360; 800, 10,400; and 900, 12,800. The head-loss equation used is Hazen–Williams with a roughness coefficient (C) of 100 and the minimum pressure for the design scenario is 15 m. The network's topology is presented in Figure 8.

In Table 5 all the parameters tested for each of the metaheuristics are presented, as well as the specific values which gave the solutions of minimum cost. The results obtained for the Taichung network are presented in Table 6 and they are compared with the costs obtained through OPUS alone and the least cost reported for the Taichung network in previous works.

### R28

R28 is a hypothetical WDS created for research purposes at the Water Distribution and Sewer Systems Research Center (CIACUA) of the University of Los Andes in Bogotá, Colombia. The network's topology consists of 28 loops formed by 67 pipes and 39 nodes that are supplied by a single reservoir, with a head of 40 m. All the nodes are at the same height except for the reservoir which is 15 m above. This network has the characteristic of having high-demand density at the center of the system and low demands in external nodes. The head loss is computed with the Darcy–Weisbach equation and the pipe material for the entire network is PVC, with an absolute roughness of 0.0015 mm. The minimum pressure for the design scenario is 15 m. The system counts 19 available diameters, which are: 50, 75, 100, 150, 200, 250, 300, 350, 400, 450, 500, 600, 750, 800, 1,000, 1,200, 1,400, 1,500, and 1,800 mm. The pipes' cost in millions of dollars per meter ($/m) can be calculated using a potential function of the diameter (mm) with a coefficient of 1.5 and an exponent of 1.45.

The demands for this network are listed in order of node identification number: 3, 7, 8, 7, 4, 12, 21, 24, 22, 8, 10, 28, 31, 27, 15, 24, 29, 33, 30, 10, 11, 31, 34, 30, 23, 27, 31, 28, 10, 9, 25, 27, 24, 13, 6, 11, 12, 11, 4 L/s.

The lengths of the pipes are listed in order of pipe identification number: 80, 150, 100, 120, 100, 120, 100, 150, 80, 200, 120, 100, 150, 80, 180, 120, 100, 150, 80, 220, 120, 100, 150, 80, 200, 120, 100, 150, 80, 180, 120, 100, 150, 80, 150, 120, 100, 150, 80, 150, 200, 220, 180, 100, 200, 180, 100, 150, 200, 180, 180, 200, 220, 100, 150, 200, 180, 180, 200, 220, 100, 150, 200, 180, 180, 200, 220 m. The network's topology is presented in Figure 9.

In Table 7, all the parameters tested for each of the metaheuristics are presented, as well as the specific values which gave the solutions of minimum cost. The results obtained for the R28 network are presented in Table 8 and they are compared with the costs obtained through OPUS methodology alone and with the least cost achieved for the network using GA MATLAB without hot start.

## CONCLUSIONS

The hot start methodology herein introduced makes use of the OPUS algorithm, which is entirely based on hydraulic principles. The designs obtained through this algorithm are used as a hot start for the application of a metaheuristic approach. In other words, the network is re-designed after the OPUS methodology has been applied, but the second time making use of a metaheuristic algorithm that takes the previous solution as the initial configuration of the network. The hot start methodology combines the use of hydraulic principles to efficiently obtain near-optimal designs, with a broader search of the solution space given by different stochastic approaches.

According to the results obtained, the methodology consistently reduces the costs obtained through the OPUS algorithm by a very little percentage, up to 1%. As for the number of iterations, it increases substantially in every case, around three orders of magnitude. As well, the RI of the solutions reached by the hot start methodology is smaller in most cases in comparison with the OPUS ones.

This methodology clearly proves that considering hydraulic bases allows the optimization of WDS design to significantly reduce the number of iterations required in comparison with metaheuristic approaches, giving near-optimal solutions. The results obtained through the hot start methodology are significantly close to the records, but do not present a considerable improvement with respect to the OPUS methodology. In addition, the computational effort required increases substantially, which implies that it is not worth trying to refine a solution that is already very close to the optimum and required minimum computational and human effort to be reached; given that the use of hydraulic knowledge for the design of WDS allows excellent results to be obtained without requiring expert knowledge of optimization techniques.

Ultimately, notwithstanding the fact that the OPUS algorithm was designed to be applied on entirely gravity-driven WDS, the methodology could be extended to the case of networks with pumping stations. The effect of this has not been investigated yet, but it is possible to speculate about some adjustments to the methodology that would allow the inclusion of pumps. In cases where the water is pumped at the reservoir its effect could be easily considered by increasing the initial HGL of the network. On the other hand, if the water is pumped somewhere downstream of the source, the HGL of every node that is fed by the pumping station could be affected in order to take into account the effect of the sudden pressure increase. The HGL of these nodes would have to be decreased in order to simulate a topological descent of the network and thus mimic the action of gravity. Namely, the power surface would be affected as a consequence of the addition of pressure at localized points of the system.

- First received 31 January 2014.
- Accepted in revised form 31 October 2014.

- © IWA Publishing 2015

Sign-up for alerts