Interactive comment on “ Algorithms for Optimization of Branching Gravity-Driven Water Networks ” by Ian Dardani and Gerard F

Author’s Response: Quantitative differences between the cost-minimization function in our paper and the one from Bhave are not the issue. It was never an intent of our paper to compare the two or encourage the use of one method over the other. Ours is simply an alternate equation, based on Taylor series fundamentals applied to a simple branch network, and embedded in an alternate method. We have proved that it works because the fundamentals are correct. The real issue is the basis for the Bhave equation vs. that for ours. Each method needs to stand on its own merit, which comes from how each was developed. We have made the development of our equation as clear as possible by explaining all assumptions and showing nearly every step in its development. In addition, we have explained in detail the differences that we see between the two, including the iterative method of Bhave vs. the use of a mathpackage-based nonlinear equation solver in ours. If you require further evidence of the correctness of our method, please advise us and we will be happy to include this in the paper.

5. It is mentioned in p.1 that three cost-minimization algorithms are adopted for the design of moderate-scale branching Gravity-Driven Water Networks.What are the other feasible alternatives?What are the advantages of adopting these particular algorithms over others in this case?How will this affect the results?More details should be furnished.
Author's Response: 5. The problem considered in this paper is that of selecting a single pipe diameter for each link in a water network to optimize the material cost.This problem has three major categories of methods that are applicable to it: enumeration methods (including both complete enumeration and partial enumeration), nonlinear programming methods, and metaheuristic methods.For each of these categories we have proposed and tested one representative algorithm: backtracking (partial enumeration), the Jones calculus-based algorithm (nonlinear programming), and a genetic algorithm (metaheuristic).While there are many other types of metaheuristic algorithms (simulated annealing, Tabu search, cellular automata, ant colony optimization, and particle swarm optimization), the genetic algorithm is the most representative of these and is also the most commonly used (Zhao et al. 2016).
Note that we have categorized the backtracking algorithm as a partial enumeration method and not a heuristic algorithm.The backtracking algorithm, like heuristic methods, follows a set of deterministic rules to find better solutions, however, those rules are strictly formulated to find cost-optimal solutions and does so without missing the global optimum.In contrast, heuristic algorithms follow rules which achieve some proxy of an optimum solution but do not guarantee the global optimum.One example of a heuristic algorithm comes from Suribabu (2012), whose algorithm uses the uniformity of a solution's flow velocity as a proxy for its cost-optimality.As such, the Suribabu algorithm increments a pipe diameter when its flow velocity is high and decrements the diameter when its flow velocity is low, in an attempt to approach more optimal solutions.It was not necessary to include such a heuristic algorithm for comparison in this paper, since the Backtracking algorithm presented follows a more strictly formulated set of rules that do guarantee a cost optimum.
We note that other methods that can be applied to water network design do not address the problem of focus in our paper.For example, linear programming methods provide split-pipe solutions for each link, while the problem addressed here calls for each link to have a single-diameter solution.Multi-objective optimization methods, which typically involve new implementations of metaheuristic methods, are also outside the scope of (single-objective) cost optimization, although it should be noted that costoptimization algorithms, such as the ones used in this paper, can be used within multi-objective implementations.In addition, decomposition methods, where networks are broken down into smaller sub-networks, can use any of the methods listed above, and are therefore not an exclusive method category.Given the appropriateness of the presented algorithms to the network sizes of gravity water networks, decomposition was not necessary.
We have added a more thorough description to our introduction section that gives clearer context to the selection of these algorithms, highlighting the category of method to which they belong and the key features that make them distinct from one another.6.It is mentioned in p.4 that five actual GDWNs installed in Panama, Nicaragua, and the Philippines are adopted as test cases.What are the other feasible alternatives?What are the advantages of adopting these particular test cases over others in this case?How will this affect the results?More details should be furnished.
Author's Response: 6.This work was prompted by our GDWN design work in Panama, Nicaragua, and the Philippines.Thus, these networks are most pertinent to our interests and are networks for which we have detailed data.Despite an extensive review of the literature, we did not uncover any common test cases for GDWNs, nor many papers focusing on GDWN design, which is likely due to their use in mostly developing countries.While each network is a branching type without loops, the total lengths range to over 15 km.This is considered a very large network in the scale of the present work.In fact, this is the longest GDWN in Central and South America.Two serial networks were tested to demonstrate the effect of a local high point on the algorithm solutions.
7. It is mentioned in p.5 that ". ..The pressure upper bound is not incorporated into the optimization process. . .."Why this is not incorporated?More justifications should be furnished on this.
7. It is correct that the pressure upper bound is not incorporated into the optimization process.Worstcase pressure conditions occur under hydrostatic conditions, which are directly related to the maximum elevation change in the network and where no flow occurs.Therefore, before the optimization process is undertaken, the selections of appropriate pressure ratings for the pipe and, if needed, break-pressure tanks are left to the correct judgment of the designer under no-flow conditions.In addition, precautions against water hammer are left to the designer since the design process presented in the paper specifically addresses network cost minimization and not this aspect of the hydraulic design.
8. It is mentioned in p.5 that ". ..precautions against water hammer are left to the designer. .."Why this is left to the designer?More justifications should be furnished on this.
Author's Response: 8. See response to 7. 9.It is mentioned in p.7 that a two-part model is adopted in this study.What are the other feasible alternatives?What are the advantages of adopting this particular model over others in this case?How will this affect the results?More details should be furnished.
Author's Response: 9. Both single-part and two-part pipe cost models are considered in the paper, which result in two different equations for network cost minimization, Equations ( 15) and ( 18) respectively.The single-part cost model is simpler than the two-part model and results in the simpler Equation ( 15) compared with Equation (18).However, as Figure 3 in the paper shows, actual pipe-cost data are more accurately fit to the two-part cost model.There appears to not be the need for more elaborate pipe cost models, nor has one been proposed in the literature as far as the authors are aware.
10.It is mentioned in p.9 that a cubic spline is adopted to complete the transition between small and large pipe sizes.What are the other feasible alternatives?What are the advantages of adopting this particular function over others in this case?How will this affect the results?More details should be furnished.
Author's Response: 10.As seen in Equation ( 21) in the paper, a first derivative of the pipe cost w.r.t.pipe diameter is needed in the network cost minimization Equation ( 18), which is based on the two-part pipe cost model.Thus, a smooth, continuous function is fit between the small-D and large-D parts of the cost model to ensure the first derivative is well behaved.The cubic polynomial is the simplest one that ensures continuous function and first-derivative behavior between the small-D and large-D parts of the cost model.
11.It is mentioned in p.13 that two rejection criteria are adopted to reduce the number of solutions to be evaluated.What are the other feasible alternatives?What are the advantages of adopting these particular criteria over others in this case?How will this affect the results?More details should be furnished.
11.In the BT algorithm, the selection of a diameter size results in a corresponding pressure head at a downstream node and a corresponding cost, both of which serve as the basis of the two rejection criteria used.Because the direct relationships between diameter and cost, and between diameter and pressure head, span the entire range of diameters, criteria to prune certain diameter choices can be safely proposed that never skip a global optimum.It can be noticed that the two criteria prune diameter solutions at the competing ends of the range of diameters, with the cost criteria pruning out larger diameters, which add more cost, and the head criteria pruning out smaller diameters, which contribute to greater head loss.These two criteria for enumeration methods have been the only ones in the literature that do not risk pruning a global optimum (Raad 2011, González-Cebollada et al. 2011).Other criteria, such as rejecting partial candidates that include a diameter that is larger than an upstream diameter, risk skipping the global optimum and were therefore not included in the BT algorithm.Instead, these were instead included in the modified algorithm, BT-NoUp, to compare the resulting effect on computation time and the cost of the solution.
These rejection criteria are further strengthened with the use of Pre-processor 1 (section 4.1) to adjust the maximum available diameter and Pre-processor 2 (section 4.2) to adjust the minimum static pressure head at each node (section 4.2), again without risk of pruning the global optimum and without the need to update any such calculations during the search process itself.To our knowledge, is the first implementation of Pre-processor 1 in enumeration methods and the first implementation of Preprocessor 2 in any method.
12. It is mentioned in p.14 that the method presented by González-Cebollada et al. ( 2011) is adopted as the BT algorithm.What are the other feasible alternatives?What are the advantages of adopting this particular method over others in this case?How will this affect the results?More details should be furnished.
Author's Response: 12. The use of backtracking for water network design has not been widely used, and only two other backtracking methods can be found in the literature, those by Gessler (1985) and González-Cebollada et al. (2011).However, the algorithm proposed by Gessler, unlike the present work's BT algorithm, includes a pipe-grouping criteria set that risks pruning the global optimum, and is therefore a less appropriate comparator to the present study's BT algorithm.The González-Cebollada algorithm, on the other hand, does not include such a criterion to potentially prune a global optimum from consideration, however, it halts its search after finding the first feasible solution.Thus, out of the two reported backtracking algorithms in the literature, both do not guarantee a global optimum, while the BT algorithm presented in this work does.The trade-off of this expanded search is an increase in computation time, however, the BT algorithm was able to execute in a satisfactory timeframe for all cases (with a maximum runtime of 79 seconds).More clarifications have been added to the paper to further clarify BT's comparison to the González-Cebollada backtracking algorithm..1943-4774.0000308, 2011.13.It is mentioned in p.15 that a single-point crossover is adopted in this study.What are the other feasible alternatives?What are the advantages of adopting this particular crossover method over others in this case?How will this affect the results?More details should be furnished.
Author's Response: 13.The crossover operation in a genetic algorithm has a number of variations including: single-point, two-point, uniform, multi-parent single point, multi-parent two-point, multi-parent uniform, universal single-point, universal two-point, and universal parent uniform.Unlike a single-point crossover, in which two individuals exchange data to the right (or left) of a randomly selected point in their chromosome, a two-point crossover includes two such random points that distinguish the points at which data is exchanged.The uniform crossover is where each link's diameter has an equally likely probability of being exchanged with the other parent.The references to multi-parent techniques involve schemes where more than two parents are used to exchange information.Finally, references to universal parents means that, at low probabilities, one of the parent chromosomes is the set of all commercial diameters.This serves to add in more variation to the population, but effectively accomplishes the same thing as the mutation operator.There is no standard technique used in the literature, though single-point crossover has a long history of use in genetic algorithms for water distribution network design (Simpson et al. 1994, Krapivka and Ostfeld 2009, Dandy et al. 1996).While we did experiment with a two-point crossover design, we did not find the results to have any benefit over a single-point design, and decided on the single-point technique for its greater simplicity.9496(1994)120:4(423), 1994. 14.It is mentioned in p.15 that a pinwheel lottery (with replacement) is adopted to select candidates to be carried into the next generation.What are the other feasible alternatives?What are the advantages of adopting this particular approach over others in this case?How will this affect the results?More details should be furnished.
Author's Response: 14.The technique used to select candidate solutions is more commonly known as a proportionate selection method, and the paper has been updated to include this term.With this technique, the fitness of each solution directly relates to its probability of being selected for the next generation.The proportionate selection method has a long history of use in genetic algorithms for water distribution network design (Goldberg 1989, Simpson et al. 1994).One feature of this selection process is that less- Author's Response: 15.The purpose of the present study is to offer a comparison of the major methods available for gravitydriven water network optimization in order to help inform designers of GDWNs on what to expect from each algorithm, where one of those is a genetic algorithm (GA).The computational research being done with genetic algorithms, and more broadly metaheuristic methods, is vast and presently developing.The number of variations available to implement in a genetic algorithm are too large to properly list in the present work.A good review of these developments in the context of water resources is given by Nicklow et al. (2010).Notably, while a theoretical basis for why genetic algorithms work has been proposed, the performance of a GA variant cannot be predicted on a given problem set, thus requiring experimental characterization (i.e., test runs on a range of cases).Even so, despite certain techniques showing benefits on water network test cases, it is not possible to generalize these findings to all water network cases.
The GA included in this paper is meant to present a benchmark for the results that could be expected when designing a GA with a reasonable amount of effort put into finding optimal parameters and the selection of proper operator schemes (selection, mutation, and crossover).It was not our intention with this paper to find the optimal GA scheme to implement for GDWNs.In fact, one improvement mentioned in the paper (scaling of the fitness function (Dandy et al. 1996)) was attempted in the present work.However, upon the trial of a few scaling schemes (Eg. a fitness exponent transitioning from 1 on the first generation to 2 on the last generation), no benefit was observed.Nonetheless, it is possible that a more thorough search of different scaling exponents may have improved the results, thus, this technique was not given a full consideration.Another improvement mentioned, the systematic optimization of parameters (Reed et al. 2000), requires quite an extensive process to find optimal search parameters.A less intensive and systematic version of this was performed, as mentioned in the parameter tuning process, although it is possible that a more systematic and thorough search would have improved the GA parameters selected.
Seen in a practical context, such advanced GA techniques are not in line with GDWN design practice.A GDWN designer typically would rarely would have access to commercial water network design software, and would typically not have much time to dedicate to optimization.Instead, designers will often select pipe diameters manually by marching along a distribution main and selecting successive pipe diameters, each time checking for hydraulic acceptability.With this as common practice, the use of any feasible optimization technique represents a substantial improvement.In our experience, and as noted above, a designer rarely has access to commercial water network design software, nor would they have the capability to implement more advanced techniques in the literature.As such, the GA represented in the paper represents a feasible performance benchmark with which to compare with the other algorithms, which are already easier to implement even without more advanced GA techniques.16.Some key parameters are not mentioned.The rationale on the choice of the particular set of parameters should be explained with more details.Have the authors experimented with other sets of values?What are the sensitivities of these parameters on the results?
Author's Response: 16.All parameters chosen for the GA algorithm are listed in Section 4.5.These include   = 0.02,   = 50,   = 0.5,   = 100,  ℎ, = 0.1.The first four of these parameters were originally chosen based on typical values presented in the literature and then tuned with a parameter sensitivity study for the first test case.Simspon et al. (1994) present typical values for   (30 -200),   (0.7 -1.0),   (0.01 -0.05), and   (100 -1000).During the parameter tuning process, each parameter was sequentially varied until a decrease or increase in no parameter resulted in significantly lower-cost results, as determined by an average of 100 runs of each.The last parameter ( ℎ, ) was found through a parameter sensitivity study and then checked to be reasonably in line with literature-based hydraulic penalty coefficients, which normalized based on their respective network's optimum price.A value for the normalized hydraulic penalty coefficient,  ℎ, , was chosen such that the GA converged on solutions that tended to satisfy the minimum static pressure constraint, but still allowed the population to gravitate towards smaller diameters with static pressures close to ℎ  .The least sensitive parameters were   ,   , and   = 100, such that even a 50% change in their values resulted in less than 1-5% change in average solution cost, while  ℎ, was moderately sensitive to changes and   was the most sensitive, thus requiring finer-tuning in increments of 0.005, where the best performance was achieved at   = 0.020.As further validation of the value chosen for   , a generally acceptable standard for GA performance occurs when   ≈ 1/ (Reed 2000).
Reed, P., Minsker, B., and Goldberg, D. E.: Designing a competent simple genetic algorithm for search and optimization, Water Resour.Res., 36(12), 3757-3761, doi:10.1029/2000WR900231, 2000.17.Some assumptions are stated in various sections.Justifications should be provided on these assumptions.Evaluation on how they will affect the results should be made.
Author's Response: 17.The main assumptions used in the problem formulation will be discussed below: turbulent flow, smooth pipe, and minimum static pressure head.The assumption of turbulent flow is reasonable for all practical potable water distribution networks, and is a traditionally-held assumption in the literature for this scale of water networks.The validity of this assumption is also verifiable post-analysis, as is mentioned in page 6 of our paper.It would not be appropriate to speculate on how the results would be different if the flow were laminar or transitional, since such cases do not normally occur in practice.The assumption of smooth pipe is consistent with conditions that exist in the field for nearly all the GDWNs we encounter in our work, which are nearly always constructed out of PVC pipe.A deeper discussion of this can be found in Jones ( 2011) pp.39-44.Finally, there is no current standard for assumptions of minimum static pressure head in gravity water network design.Some examples from gravity water network design references include 5 m (Bouman, 2014), 10 m (Arnalich, 2010), 5-15 m (Action Contre La FAIM, 2008), and 8-20 m (Swamee and Sharma 2008).Thus, 7 m is a reasonable assumption and in line with our own experience of common practice.Including even larger values for hmin will unnecessarily increase pipe sizes and thus cost.This justification has been included in the paper.18.The discussion section in the present form is relatively weak and should be strengthened with more details and justifications.
Author's Response: 18. We have added a number of additional elements to the discussion in order to strengthen it, per your suggestion.A discussion of the CB results in the context of the mapping schemes now appears in Section 6.Even though our tabular results are based on up-sizing, each of the methods described in Section 6 produce a solution closer to the global minimum of the more-computationally intensive and not-soscalable BT method.These methods engage the designer to make decisions on pipe-ordering (e.g., large diameters are chosen upstream of small diameters for the compound-pipe approach) and the acceptability or unacceptability of a diameter change relative to the HGL corresponding to the CB continuous diameter results.In addition, discussion of the rationale behind GA techniques used and parameter tuning has been added, including a clearer explanation of the relevance of these criteria may have on the results of the study.Where appropriate, we have included more details in sections which were a better fit for their inclusion, for example, the comparison of our BT algorithm to the González-Cebollada algorithm (section 4.3).
19.Moreover, the manuscript could be substantially improved by relying and citing more on recent literatures about real-life case studies of contemporary soft computing techniques in water resources engineering such as the followings: 19.We have looked at the literature examples you have provided and reviewed their scope in comparison with that of the current paper.Our paper focuses on water distribution network design, a significant body of work in its own right without including other aspects of the much broader field of water resources.We felt that the inclusion of references on aquifer levels and flows prediction for watersheds and runoff did not fit within the scope of our paper.We have included a reference to Taorimina (2015), however, in an attempt to highlight the many uses of metaheuristic methods in water resources engineering.
20.In the conclusion section, the limitations of this study, suggested improvements of this work and future directions should be highlighted.
Author's Response: 20.The conclusions section has been has been updated to highlight the limitations of the study (e.g., the scalability of the BT algorithm, the lack of inclusion of some GA improvements reported in the literature) and improvements and future directions for the work (e.g., integrating the strengths of the BT to map the CB algorithm continuous diameter solution to a discrete diameter solution and verifying the suitability of the BT and BT-NoUp algorithms on other large GDWNs).
Thank you again for your time spent reviewing our paper. 1

Algorithms for Optimization of Branching Gravity-Driven Water Networks
Ian Dardani 1 , Gerard F. Jones 1 1 Villanova University, College of Engineering, Villanova, 19085, United States Correspondence to: Ian Dardani (ian.dardani@villanova.edu) Abstract.The design of a water network requires the selection of pipe diameters that satisfy pressure and flow requirements while optimizing for cost.This work focuses on the design of moderate-scale branching Gravity-Driven Water Networks (GDWNs), in contrast to large urban-scale looping networks, where budgets are highly constrained and where Polyvinyl chloride (PVC) pipe is typically used.To help designers of GDWNs select an appropriate design approach for a given network problem, three cost-minimization algorithms are presented and results compared with five GDWN test cases.Two algorithms, a backtracking algorithm and a genetic algorithm, use a set of discrete pipe diameters, while a new calculus-based algorithm produces a continuous-diameter solution, which is mapped onto a discrete-diameter solution.The backtracking algorithm produced the overall lowest-cost solutions with relative efficiency for the test cases, while the calculus-based algorithm produced slightly higher-cost results but with greater scalability to networks with more links.Furthermore, the new calculusbased algorithm's continuous-diameter and mapped solutions provided lower and upper bounds, respectively, on the discretediameter global optimum cost, where the mapped solutions were typically within one diameter size of the global optimum.
Overall, the genetic algorithm as implemented did not produce results which deemed it compelling over deterministic methods for GDWN design.However, for more complex networks and problem formulations, a genetic algorithm may be more advantageous, particularly if it incorporates improvements reported in the literature.The results of this study highlight the advantages and weaknesses of each GDWN design method including closeness to the global optimum, the ability to prune the solution space of infeasible and sub-optimal candidates without missing the global optimum, and calculation time.We also extend an existing closed-form model of Jones (210) to include minor losses, a more-comprehensive two-part cost model, which realistically applies to pipe sizes that span a broad range typical of GDWNs of interest in this work, and for smooth and commercial steel roughness values.

Introduction
A gravity-driven water network (GDWN) is commonly constructed to deliver potable water to a community in a developing region.These systems draw water from a source at a high elevation, such as a natural spring or a stream, and deliver it through a branching pipe network to household taps or public tapstands (Fig. 1).In principal, loops and 2 loop/branching constructs may be added to networks for greater reliability, but material cost considerations often restrict attention to just branch networks in GDWNs.The methodologies presented in this paper, however, may be extended to all networks, including those with loops.When feasible, gravity water networks are very attractive compared with pumped networks because of their simplicity and lower capital, operational, and maintenance costs.In addition, in most locations where GDWN are considered, there may be little or no access to reliable grid-based electrical power for pumps.
Water networks are modelled as a collection of nodes, each representing a point of water demand or supply, which are connected with links representing pipes.Typically, the layout of the site is known, including water source and demand locations and elevations of all other nodes.For the present work, design flow rates are determined from community survey data, which are extrapolated for future population growth.Networks in this category are referred to as "demand-driven" designs.Bhave (1978) refers to these as "Q-specified" designs.Thus, to design a network of this type, pipe diameters for each link must be chosen such that acceptable but arbitrary minimum pressure heads are maintained at each node given a design flow rate at the node.Furthermore, application of the energy equation to this network demonstrates that the design problem is non-unique; i.e., choosing different pressure heads at the nodes will result in a different pipe diameter solution for the network, and thus different networks costs.
In practice, gravity-driven water networks are commonly designed by a marching method, where diameters for each link of the network are chosen sequentially.After selecting a reasonable diameter for each link, the designer calculates the static pressure head at the link outlet, and proceeds to the next link if this result is acceptable.In this way, the designer marches through the network until all pipe diameters have been selected.This method produces a feasible solution, but not a costoptimized one.As noted by Bhave (2003), cost savings of 20-30% can result from the use of optimization techniques.In developing regions, the cost of a water network can be prohibitive, adding to the importance of optimizing network design.
Within the provided framework, the global optimum can be found through an exhaustive search of the solution space, known as complete enumeration, although this is infeasible when considering networks with many links and diameter choices (Kadu et al. 2008;González-Cebollada 2011).To reduce the computational time required by enumeration, authors have proposed various partial enumeration methods which prune the search space (Kadu et al. 2008), although some of these techniques may remove the global optimum (Simpson et al. 1994).The most common types of algorithms that have been applied to optimize water network design include traditional deterministic methods, heuristic methods, metaheuristic methods, multi-objective methods, and decomposition methods (Zhao et al. 2016).
Deterministic methods include linear programming (LP), dynamic programming, and nonlinear programming (NLP), and typically involve rigorous mathematical approaches (Zhao et al. 2016).A brief overview and comparison of these algorithms is given in Kansal, et al. (1996), who use a single-part cost correlation for metric pipe diameters between 100 mm and 350 mm.Linear programming techniques have relatively low computational complexity and allow each link to be composed of two diameters, called a split-pipe solution, although these may not always be practical to implement (Kessler and Shamir 1989, Swamee and Sharma 2000, Samani and Mottaghi 2006).LP can also get stuck in a local optimum (Zhao et al. 3 2016), although combining LP with metaheuristic techniques can help with the problem's non-smoothness properties (Krapivka and Ostfeld 2009).Dynamic programming has been used by Yang et al. (1975) and Martin (1980) to optimize networks in stages.This approach begins at the discharge nodes, proceeding to select feasible diameters and joints for upstream stages and storing these partial candidates in memory until the source node is reached.At this point, the algorithm reviews the feasible segment design options and selects a combination of stage solutions producing the lowest cost overall solution.This method, however, requires the designer to allow a relatively narrow range for the design pressure of each node, or otherwise store a large set of feasible candidate solutions in memory and also allow adjoining branches to arrive at different heads at the same node.
Nonlinear programming, a calculus-based method, deals with each link's diameter as a continuous variable.Using Lagrange multipliers and a one-part, pipe-cost model with minor-lossless flow, Swamee and Sharma (2000) developed systems of equations for both continuous and discrete pipe diameters for branch networks, assuming constant friction factor.When solved, the solution gives diameter values that minimize distribution main cost, not network cost.In carrying out the solution, iteration is required to update the value of the friction factor.For the discrete diameter case, large computational times were noted by Swamee and Sharma because of the stiffness of the mathematical system.Cases where one or more nodal pressure heads are not acceptable need to be treated manually by the designer in various ways as discussed by the authors.
For branching networks, Jones (2011) showed that by restricting the focus to smooth-turbulent, minor-lossless flow, and the use of a one-part, pipe-cost model, a simple nonlinear algebraic equation for each internal node in the distribution main could be developed.The equation has been extended in the present work to include minor losses and rough pipe.When solved simultaneously with the energy equation for each link, a unique solution for all link diameters and nodal pressure head values are obtained that produces minimum network cost, as opposed to the distribution main cost as in Swamee and Sharma (2000).
The method of Jones (2011) also applies to serial and loop networks because of its generality.
Heuristic methods follow specific rules to incrementally build better solutions, although the rules are not strictly formulated to trend towards local or global optima.An approach by Monbaliu et al. (1990) sets all network pipes to their minimum size, where the pipe that has a maximum head loss gradient is incremented to its next-highest size until all nodal head requirements are satisfied.Similarly, an algorithm by Keedwell and Khu (2006) selects an initial solution and iteratively responds to nodal head deficits and surpluses by incrementing or decrementing pipe sizes accordingly until a feasible solution is found.Suribabu (2012) proposed a heuristic that identifies pipes to increment or decrement in size based on flow velocity and alternative metrics such as proximity to the source node, achieving acceptable cost results with computational efficiency.
While these algorithms are typically computationally efficient, they do not guarantee a global optimum.
Metaheuristic optimization methods allow for a set of solutions to evolve through random processes that are guided with an objective function which rewards low network costs and penalizes hydraulic insufficiencies.Examples include evolutionary algorithms, which are most commonly genetic algorithms (Krapivka and Ostfeld 2009, Simpson et al. 1994, Kadu et al. 2008, Prasad and Park 2004), simulated annealing (Vasan and Simonovic 2010;Tospornsampan et al. 2007), ant colony optimization (Maier et al. 2003), and differential evolution (Vasan and Simonovic 2010).As reviewed by Nicklow et al. (2010), evolutionary algorithms are an emerging popular alternative to the deterministic methods, and they offer the opportunity to accommodate unique constraints and multiple design objectives.The main challenges for evolutionary algorithms are the difficulty of incorporating constraints into objective functions, the optimum selection of parameters, and a relatively large amount of computational effort.In addition to optimizing for cost, multi-objective methods, often based on evolutionary algorithms, allow the designer to choose from a Pareto-optimal front of objectives, such as cost and reliability (Prasad and Park 2004).In addition to water network design, metaheuristic algorithms have been used for a range of problems in water resources engineering, such as rainfall and runoff modelling (Taormina et al. 2015).
Decomposition methods involve the partitioning of networks into smaller sub-networks which are each optimized using one of many types of techniques and then combined into an overall solution.In some cases, the loops in the sub-networks are removed, producing branching trees which are then optimized individually.Techniques used to optimize the sub-networks can involve multiple methods, including linear programming (Saldarriaga 2013) and differential evolution (Zheng et al. 2013), with a later stage optimizing the network as a whole using the sub-network solutions as inputs.Note that another distinct use of the term 'decomposition' refers to the approach of iteratively solving "inner" and "outer" mathematical problem formulations, and has been used in the literature by Krapivka and Ostfeld (2009) who traces its use in this context back to Alperovits and Shamir (1977).
In the present study, we present three algorithms, each from one of three major categories of methods applied to cost optimization of water distribution networks, and compare their performance on five cases adapted from real GDWNs.These Within the broader context of water network problem formulations, this paper is concerned with finding cost-optimal single-diameter solutions to branching water distribution networks with steady-state demand flows and pre-specified pipe locations..By implication of being gravity-driven, the problem does not involve the use of pumping stations.This problem formulation is directly applicable to typical gravity-driven water networks, and is also useful for multi-objective algorithms, the consideration of sub-networks in a decomposition technique, pumped networks, and looped system optimization, which can involve reformulating the problem into a branching configuration.
The results of this study highlight the advantages and weaknesses of each GDWN design method including closeness to the global optimum, the ability to prune the solution space of infeasible and sub-optimal candidates without missing the global optimum, and also computational time.We present two pre-processors with which discrete-diameter search methods can use to reduce the search space without pruning the global optimum.To the authors' knowledge, is the first implementation of Pre-Processor 1 in enumeration methods and the first implementation of Pre-Processor 2 in any method.We also extend the Jones closed-form model to include minor losses, a more-comprehensive two-part cost model, which realistically applies to pipe sizes that span a broad range typical of GDWNs of interest in this work, and for smooth and commercial steel roughness values.

Problem Formulation
Branching networks are considered (Fig. 1), where all branches connect a distribution main node with a delivery For all nodes, static pressure, ℎ, is greater than or equal to a chosen minimum, ℎ  .The value for ℎ  is selected to eliminate possible leakage of contaminated ground water into the network should the operating conditions change in an unanticipated way.The change in static pressure head, Δℎ, across each link is calculated with the energy equation for pipe flow, (1) where for each link, is the kinetic energy correction factor and f is the Darcy friction factor, calculated with the Colebrook-White equation (Colebrook and White 1937) or Churchill correlation (Churchill 1977), and  is acceleration of gravity.The kinetic energy correction factor, , is considered only in the first link, where acceleration from a zero-velocity source is sometimes non-negligible for the smallest of GDWNs that have been encountered.Thus, where Re is the Reynolds number for pipe flow, 4Q/   D, and  is the kinematic viscosity of water.The possibility of laminar flow (Re ≤ 2300) is permitted since branches from the smallest GDWN observed in practice have been in this regime.
The pressure upper bound is not incorporated into the optimization process.Worst-case pressure conditions occur under hydrostatic conditions, which are directly related to the maximum elevation change in the network and where no flow occurs.Therefore, before the optimization process is undertaken, the selections of appropriate pressure ratings for the pipe and, if needed, break-pressure tanks are left to the correct judgment of the designer under no-flow conditions.In addition, precautions against water hammer are left to the designer.

New Calculus-Based Algorithm
In this section we develop a new calculus-based algorithm for pipe diameters that minimize overall pipe cost for the network.First appearing in the text by Jones (2011), this algorithm is solved simultaneously with the energy equation for each link to produce unique solutions for D and nodal pressure head values that minimize network pipe cost, as opposed to only the distribution main cost as in Swamee and Sharma (2000).The method also applies to serial and loop networks.
First consider the physical basis for the existence of a unique set of pipe diameters and static pressures for the demanddriven design problem with cost minimization included.Several works reviewed in the previous section have considered optimization of GDWN and combined pumped and gravity-driven networks.We assume continuous pipe diameters in this section; values that result from the solution of the energy equation.Mapping between continuous diameters and the discrete nominal sizes, required to complete the design, will not be fully addressed in the present work.However, we will discuss two methods we have used for mapping the continuous D solutions onto the discrete pipe diameter set.
Consider the three-pipe network shown in Fig. 2. Pipes 1-2, 2-3, and 2-4 meet where head ℎ 2 is unknown.Each pipe has prescribed volume flow rate and length and unknown diameter  as shown.The change in elevation between the top and bottom of each pipe is Δ and Δℎ is the change in static pressure head.There is a prescribed head at each outlet for pipes 2-3 and 2-4.
To facilitate insight, we at first assume turbulent flow, which can be verified post-calculation if necessary, in smooth pipe and that minor losses are negligible.Two sources for the friction factor for smooth-turbulent flow are considered, namely the classical Blasius equation (reported in Streeter et al. 1998), f = 0.316 Re -1/4 , and the Swamee-Jain correlation (Swamee and Jain 1976), f = 0.175 Re -0.1923 (though not explicitly appearing in this reference, f from the Swamee-Jain correlation is obtained by writing it for smooth pipe and comparing this with the energy equation, where f is assumed to be in the form a Re n ).The Blasius equation has higher accuracy (2% for low Re and 3% for high Re) in the range 10 4 < Re < 10 5 , over which most of the GDWNs in this work operate, compared with the Swamee-Jain correlation of +8% / -3%, thus prompting the Blasius equation to be chosen for this work.A combination of the Blasius equation with the energy equation gives explicit formulas for  for the three links in Fig. 2. For simplicity, and to reduce the number of free parameters, the conditions for pipes 2-3 and 2-4 are assumed to be identical without loss of generality.We obtain (3) The pipe cost model can be assumed to follow a power-law relationship (Swamee and Sharma 2008) where  is a constant coefficient,  is a constant exponent, and   an assumed unit diameter.A more robust, two-part model, valid for a greater range of pipe sizes than that of Swamee and Sharma (2008), will be used below.The use of pipe material cost as the objective function was assumed because of relevance.In most GDWNs of interest in this work, installation labor comes from the local community and has no well-defined associated cost.The material cost for the network is of prime importance since it normally comes from funds raised by nongovernmental organizations or grants, where there is seldom a required repayment but are always in short supply.For a more-general case, the economics of a GDWN may be more encompassing and include materials, labor, operation and maintenance, depreciation, taxes, and salvage, among others.The time value of money may also need to be considered, which includes interest rates and estimation of the network lifetime.
With Eq. ( 4) the general expression for the total cost for the pipe material,   , is obtained by summing over all links A close inspection of Eq. (3) in combination with Eq. ( 6) will reveal the origin of the existence of an optimal ℎ 2 for the design of the network in Fig. 2. Because of its arbitrariness, we are free to vary the value of ℎ 2 .As ℎ 2 increases, say from a small value like 1 m, the pressure difference between the junction and the bottom of pipe 2-3 (and 2-4) increases.Since the volume flow rates in each pipe are fixed, an increase in pressure drop across pipe 2-3 (and 2-4) requires a reduction in  23 (and 𝐷 24 ).This is evident from our inspection of the second of Eq. ( 3), where we see that  23 and  24 are both proportional to (Δ 12 − ℎ 3 + ℎ 2 ) −4/19 ; that is (Δ 12 − ℎ 3 + ℎ 2 ) −4/19 decreases as ℎ 2 increases.
Because the top of pipe 1-2 is at atmospheric pressure, an increase in ℎ 2 will decrease the pressure drop between the top of pipe 1-2 and the junction.Thus, compared with pipes 2-3 and 2-4, the opposite effect occurs in pipe 1-2;  12 increases with increasing ℎ 2 .For insight on how the energy equation supports this explanation, note that the first of Eqs.
From this discussion it is clear that for an increasing ℎ 2 there is a competition between the decrease of  23 (and  24 ) and an increase in  12 .Once the effect of  on pipe cost is included through Eq. ( 6), as ℎ 2 increases we see that the cost for pipes 2-3 and 2-4 decrease, and the cost for pipe 1-2 increases.A consequence of this competition is the existence of an optimum, in this case an optimal ℎ 2 , which produces the smallest possible cost.
The mathematical basis for a unique solution for ℎ 2 with cost minimization is now presented.In addition to the fixed pipe lengths, the total cost depends on the diameters for all of the pipes in the network.For the case of Fig. 2, where we now allow pipe 2-3 and pipe 2-4 to be different, get Using the chain rule from the calculus, the total differential of Eq. ( 7) is The minimum value of   is found once   = 0 (and once it is verified that the second derivative of   is positive thus indicating that   is indeed a minimum).Requiring this, obtain The cost   is from Eq. ( 5), so the derivatives like   / 12 in Eq. ( 9) are written in general as for any link .
The derivatives like  12 /ℎ 2 in Eq. ( 9) are obtained by taking the partial derivative of the pipe diameter with respect to the relevant pressure head in the appropriate energy equation.(13) Equations ( 10)-( 13) are combined with Eq. ( 9) to produce a single algebraic equation that depends on ℎ 2 , as well as  12 The general form of Eq. ( 14), written at any internal node is where the hydraulic gradient,   , is In Eq. ( 15) the indices ij,in and ij,out on the summations refer to inflows and outflows at the node (e.g., in Fig. 2, ij,in =12 and ij,out = 23 and 24).Equation ( 15), the new CB algorithm proposed in this work, is written for each internal node in the network and solved simultaneously with the energy equation for each link to obtain unique and optimal values of Dij for all links and hj for all internal nodes.It is understood that the nodal pressure heads determined from the solution of this system must be greater than or equal to the hmin prescribed for the network.For nodes that do not satisfy this condition, the pressure head is set equal to hmin, as part of the CB algorithm.Thus hj > hmin.
Minor losses using the equivalent-length method can be included in the above developments by artificially extending the length of the link by   in which minor loss occurs, thus contributing a non-zero   term in Eq. ( 1).We also extend the cost model of Eq. ( 5) from Swamee and Sharma (2008) In Eq. ( 17),   and  ℓ are the coefficients for the small and large pipe size regions, respectively, and   and  ℓ are the exponents for the small and large pipe size regions, respectively.A cubic spline is fit between pipe sizes   and  +1 to complete the transition between small and large pipe sizes.The coefficients of this polynomial are  1 ,  2 ,  3 , and  4 as seen in Eq. ( 17).
Equation ( 18), and its simpler form Eq. ( 15), forms the basis for calculus-based optimization in this work and is applied at all internal nodes to uniquely determine hj.Equation ( 18) is valid over the range of ~4000 < Re < ~300,000.
Algorithms to solve a general set of independent, nonlinear algebraic equations using, for example, the Levenberg-Marquardt, Quasi-Newton, Newton-Raphson, or Conjugate Gradient methods are available in most commercial math packages including Matlab (1 Apple Hill Drive, Natick, MA USA 01760) and Mathcad (http://www.ptc.com).We used the package Mathcad in the present work.Thus, compared with an iterative solution procedure, a solution flowchart is not relevant here.Bhave (1978) first proposed an algorithm like Eq. ( 15).However, Bhave used an iterative method to solve the design problem.As such, there are several qualifications leading up to the cost minimization equation in Bhave.These include the assumption of smallness in variation of the static pressure head between two iterations.This allowed the terms in the cost function to be approximated as constants.In the present work, the cost-function coefficient and exponent are not assumed constant at any node joining two sets of links; see Equation ( 18).Nor do we make any assumptions on the orders of magnitude of the terms in our equations to simplify them.For clarity, we re-present Eq. ( 15) using Bhave's (1978) notation as where the ij and jk notation are shown in Fig. 4. Index j spans all internal nodes along the distribution main.

Backtracking Algorithm and Genetic Algorithm
Backtracking (BT) and genetic algorithm (GA) assess candidate solutions composed of discrete diameters from a commercially available set.These candidates are represented by a vector of size [1,   ] where each element corresponds to a network link.
The values of the vector specify a diameter from the commercially available set that are indexed from smallest (  = 1) to largest (  =   ).To reduce the computational time associated with these evaluations, the constraints imposed by the energy equation and cost minimization may be more efficiently evaluated through lookup tables.With fixed , Δ, ,   , and , the change in static pressure head Δℎ is evaluated for all     combinations of pipe diameter and link index: While an algorithm evaluates a candidate solution, the static pressure head at each node is sequentially calculated by "marching" through the network.Starting with the fixed source pressure head, the algorithm finds the pressure head ℎ  for a given node by adding the head at the upstream node, ℎ −1 to the change in head for that link   and the diameter   under consideration.Thus, Along with the hydraulic evaluation of a candidate solution, the cost of the partial candidate is found through the use of a where (  ,   ) returns the additional cost of assigning a diameter with index   to link   .In this way, the candidate solution's hydraulic performance and cost are incorporated into the genetic algorithm and backtracking approaches.In contrast to GA, the backtracking algorithm evaluates static pressure head and cost upon consideration of each partial candidate, where GA calculates these values on full candidates as part of the objective function.

BT and GA Pre-Processor 1: Maximum Available Diameter
To increase the efficiency of BT and GA, it is advantageous to limit the number of pipe diameters in the available set, especially those outside of the range of the optimal solution.In particular for the BT algorithm, larger diameters can require considerable computational effort, since they tend not to violate static head requirements and require multiple-link partial candidates for the algorithm to reject them once their cost exceeds that of an already-found viable candidate.Therefore, a preprocessor is used to provide a maximum diameter (  ) that should be considered during the optimization process.This procedure, which produces a conservative estimate, finds the smallest diameter at which a network with a single pipe diameter choice produces no nodes with a static pressure head below ℎ  , similar to the technique used by Mohan and Jinesh Babu (2009).After this diameter is found, the next-larger diameter in the set is selected as   to allow the algorithm to select a larger-than-necessary diameter if this is able to save cost elsewhere.If   appears in the optimum solution, the designer can elect to further increase this maximum diameter.It worth noting that Kadu et al. (2008) presents another method to further prune the search space with the critical path concept, where Dongre and Gupta (2011) noted the computational advantages of having just four diameter choices per link.This method, however, may prune the global optimum and may not produce feasible head values at intermediate nodes, as in the case of networks with a local high point.
will be viable for the full network solution.Note that both pre-processors discussed will not prune the global optimum from the solution.

Backtracking Algorithm (BT)
The backtracking algorithm employs a systematic search of candidate solutions to find a global optimum.The algorithm works recursively to incrementally build candidate solutions while checking the candidates for hydraulic and cost acceptability.The strength of the BT is that, upon discovery of an infeasible partial candidate, all extensions of that candidate can be eliminated from consideration.In this way, many solutions can be pruned from the solution tree to achieve greater computational efficiency.
Two backtracking methods can be found in the literature, namely those by Gessler (1985) and González-Cebollada et al. (2011).The algorithm proposed by Gessler, however, also proposes a pipe-grouping criteria that risks pruning the global optimum and represents a tricky optimization problem in and of itself (Raad 2011).The González-Cebollada algorithm, on the other hand, does not include such criteria to potentially prune a global optimum from consideration, but it halts its search after finding the first feasible solution.In contrast, the BT algorithm in the present study guarantees a global optimum by continuing its search of the solution tree even after the first solution has been found.In addition, the present BT algorithm utilizes Pre-Processors 1 and 2 to further reduce the search space, without risk of pruning the global optimum, in advance of its search routine.Thus, out of the two reported backtracking algorithms in the literature, both do not guarantee a global optimum, while the BT algorithm presented in this work does.It should be noted that BT is known to scale poorly with large network sizes and would not be appropriate for use on large urban networks, though its appropriateness is demonstrated here for GDWNs, given that the test cases used in this paper representative of the sizes of GDWNs that would be expected in practice.
Backtracking is a type of partial enumeration method, which Raad (2011) notes can drastically reduce the number of solutions to be evaluated based on two rejection criteria.The first rejection criterion is that when a candidate violates static pressure head constraints, all candidates with equal or lesser diameter sizes can be discarded.This condition is leveraged even more effectively with Pre-Processor 2 above, which can increase static pressure heads at individual nodes by anticipating the head requirements at surrounding nodes.The second rejection criterion is that once a feasible candidate has been found, all other partial candidates with a higher cost can also be discarded.The BT algorithm further extends this criterion by considering that the links yet to be considered in a partial candidate, an "extension" to the partial candidate, will cost at a minimum that of the entire extension being composed of the smallest available diameter.Thus, when considering whether the partial candidate will necessarily be more expensive than the running optimum, this minimum extension cost can be added to the partial candidate cost.
The backtracking algorithm begins its search of the solution tree by considering the partial candidate with the smallest diameter size assigned to the first network link.The static pressure head and the partial candidate cost at the outlet node are calculated with the  and  lookup tables.If this partial candidate meets static pressure head and cost requirements, the algorithm extends this partial candidate by assigning the smallest diameter to the downstream link.If a partial candidate produces a node that is rejected on the basis of static pressure head, the next largest larger diameter is chosen for the link upstream of the node.If no diameter satisfies the pressure head condition, the algorithm backtracks to the upstream link and assigns a larger diameter to the link.If a node is connected to a delivery node by a single link, the smallest feasible diameter for that link is found, and if no such diameter exists, the partial candidate is rejected.In this way, the algorithm continues to extend and reject candidate solutions until a full candidate satisfies the static pressure head requirements.Once this has been achieved, the diameter choices and cost of the network are stored as a running optimum.
Once a working solution has been found, candidate solutions may be rejected based on cost.For each new candidate, cost is calculated by adding the cost of diameters that have already been assigned to the cost of assigning all downstream links with the smallest diameter available.If this cost exceeds the cost of the running optimum, the partial candidate is rejected.
Unlike a candidate rejection based on static pressure head, a rejection based on cost does not consider siblings with larger diameters, since these would only further add to cost, rather, the algorithm backtracks immediately to re-assigning the upstream link.
In this way, the rejection criteria based on a minimum static pressure head and cost are used to prune the solution tree and largely reduce the number of non-optimal candidates that need to be considered.The minimum static pressure head criterion tends to prune candidates with diameters that are too small, while the running optimal cost criterion tends to prune candidates of diameters that are too large.
Another pruning technique noted by Raad (2011) is to group together adjacent links that are sized identically.This technique, in contrast to the former two mentioned, cannot guarantee an optimal result, and is therefore not included in the present study's BT algorithm.The present study's BT algorithm operates similarly to the method presented by González-Cebollada et al. (2011), with the major difference being that the BT algorithm continues searching once it has found its first feasible solution and its use of Pre-Processors 1 and 2. The BT algorithm could also be used in this way to find an initial solution very quickly, and then continue as normal to find progressively better solutions until the end of the search space, or a predefined condition such as calculation time, are met.

Modified Backtracking Algorithm (BT-NoUp)
A modification to the BT algorithm was made to further improve its computational speed, although at the risk of pruning the global optimum from the search.This modified algorithm (BT-NoUp) rejects all candidates that feature a smaller diameter that is upstream of a larger diameter when an equal or smaller flow rate is present in the downstream link.Typically, optimal networks would not exhibit this feature, and in cases where a single source feeds into a network with constant-length links, it is advantageous (or equivalent) to place larger diameters upstream of smaller diameters.However, due to the discrete nature of diameter choices and link lengths, an optimization problem may, in fact, have an optimal candidate with a larger diameter cost.With each generation,  ℎ is updated by multiplying the normalized penalty coefficient,  ℎ, , by the average pipe cost of the population, The algorithm then selects candidates to be carried into the next generation through a proportionate selection method, where each candidate has a probability of being selected,  , in direct proportion to its fitness relative to the sum of all fitness values in the population The algorithm replaces the parent generation with a generation of equal size and tends to select more fit individuals in successive generations.In this study, the genetic algorithm parameters used were   = 0.02,   = 50,   = 0.5,   = 100,  ℎ, = 0.1.The first four of these parameters were chosen based on typical values presented in the literature and then tuned with a sensitivity analysis for the first test case.Simspon et al. (1994) present typical values for   (30 -200), (0.7 -1.0),   (0.01 -0.05), and   (100 -1000).The normalized hydraulic penalty coefficient,  ℎ, , was chosen such that the GA converged on solutions which tended to satisfy the minimum static pressure constraint, but still allowed the population to gravitate towards smaller diameters with static pressures close to ℎ  .

Cases Studied
Five cases were studied based on actual GDWN in Panama, Nicaragua, and the Philippines.Global characteristics of each network are presented in Table 1 and the details of each network are presented in Table 3( The choice of ℎ  is not standardized, and should appropriately balance the risk of negative pressure in pipes and the increase in network cost due to the requirement of using larger diameters.The choice of ℎ  in GDWN design is typically in the range of 5 m -20 m (Arnalich 2010;Bouman 2014;Swamee and Sharma 2008).In the present study hmin = 7 m, although this requirement was reduced at selected nodes at the beginning of networks where changes in elevation are still small.At the source node, the static pressure head is fixed at atmospheric pressure.All cases assumed minor-lossless flow, although all algorithms (e.g., Eq. ( 18) for CB-Theor) are capable of handling minor loss coefficients through the equivalent length method as presented above.
6 Mapping the Theoretical D to Discrete Pipe Sizes The mapping between continuous diameters and the discrete nominal pipe sizes was accomplished in our solution by one of the following ways: 1.For small and moderate size networks, the designer may manually adjust the pipe sizes (downward, normally one pipe size) starting from the first link downstream from the source and continuing along the rest of the distribution main to the end.A nearby plot of the static pressure heads compared with the theoretical Dij from our CB approach (on the same Mathcad page) will highlight the acceptability or unacceptability of any change.This exercise also gives the designer an understanding of the sensitivity of the design to small changes in pipe sizes.
2. Based on the theoretical Dij from the CB approach, a composite pipeline can be created for each link.That is, the lengths for the two discrete pipes sizes that bound the theoretical Dij from above and below are calculated such that the pressure drop between two consecutive nodes in the distribution main matches between the composite pipeline and the CB approach.This also provides discrete pipe sizes that nearly matches the CB solution in terms of cost.

Results
The current study evaluated three types of algorithms that optimize the design of gravity-driven water networks (GDWN).The algorithms include a calculus-based (CB) algorithm, a backtracking algorithm (BT), a modified backtracking version (BT-NoUp), and a genetic algorithm (GA).The algorithms were applied to five test cases that are based on real GDWNs.
The global optimum network cost, found with BT, is shown in Table 2.The costs of solutions from all other algorithms are expressed as a percentage difference in cost from the global optimum cost.To visually compare the algorithm solutions, the hydraulic grade lines from BT, BT-NoUp, CB-Theor, and CB-Disc are presented in Fig. 5 along with the network elevation for each test case.For clarity, the hydraulic grade lines of branch links are omitted from the figure.In addition, the GA solutions are omitted since 100 solutions were obtained for each test case.Collectively, the hydraulic grade lines reveal a close alignment of the BT solution (the global optimum) with the CB-Theor solution which utilizes a continuous diameter set.Furthermore, the mapping scheme used to generate a CB-Disc solution is shown to increase pipe sizes in some cases far beyond the limit imposed by ℎ  , which was set to 7 m in the present work.
In practice, a GDWN must be designed with pipe diameters that are selected from a discrete, commercially available set.With a given number of network links,   , and a number of available diameters,   , a total of     candidate solutions exist, yet with only one global optimum (except in the case of no viable solutions or unique solutions with identical costs).For example, a GDWN of 20 links and 13 commercially available pipe sizes will, in principal, produce approximately 1.910 22 candidate solutions.
prune infeasible and sub-optimal candidates.In this study, BT evaluated only a fraction of the candidate solutions, where the fraction ranged from 410 −18 to 710 −4 .To further reduce the number of evaluations required to arrive at a solution, the BT algorithm was modified (BT-NoUp) to prune all solutions that include a smaller diameter that is upstream of a larger diameter.
This criterion, which seems intuitive to the designer, may actually miss the global optimum due to trade-offs associated with discrete solutions.In fact, BT-NoUp missed the global optimum in cases 2 and 3, although by a small percentage increase in cost (2.1% and 0.35% respectively).BT-NoUp, however, finished its search in a shorter amount of time in comparison to BT.
Using a Dell Latitude (i5 CPU at 2.50 GHz), the time to evaluate one candidate for both BT and BT-NoUp was around 0.2 ms.
In addition to approximately 2 s of pre-processing time, the computation times for BT ranged from 0.08 s (case 1) to 79 s (case 3), while BT-NoUp ranged from 0.06 s (case 4) to 0.3 s (case 3).The number of available diameters used in the BT, BT-NoUp, and GA runs are listed in Table 1.
The CB algorithm, unlike the other algorithms in this work, finds a solution with theoretical diameters that are drawn from a continuous domain (CB-Theor).For all test cases, the costs of the CB-Theor solutions was less when compared with the discrete-diameter global optimum (-5.46% to -2.60%).In fact, because of the discrete pipe sizes needed for an actual network, the continuous model will always produce the smallest theoretical cost for the network.The CB algorithm then maps this solution to a commercially-available discrete set (CB-Disc).The mapping process used in this study simply mapped each theoretical diameter to the nearest available diameter of a larger size, thus producing a solution which still satisfies static head requirements but with a higher associated material cost.This tended to oversize the diameters, although the CB-Disc solutions were always within two diameters of the BT global optimum solutions, as shown in Fig. 6.From all of the test cases combined, all but one (71 out of 72) of the diameter selections were within one diameter of the global optimum.More sophisticated mapping schemes, like independently adjusting D for each link in the distribution main in a step-by-step manner starting with the source while ensuring all pressure head constraints are satisfied, would more likely produce results identical to the global optimum.This was not performed in the current study.The CB-Disc solution costs were, in all cases, larger than the global optimum, with a percentage difference ranging from 3.86% to 22.6%.Thus, for all cases, the calculus-based algorithm bounded the cost of the global optima with a lower-cost CB-Theor solution and a higher-cost CB-Disc solution.This trend is a result of the additional constraints imposed by the finite set of diameter choices.If the algorithm is allowed a greater number of discrete diameter choices, i.e., through adding a less-common nominal diameter size to the available set, the cost of the CB-Disc solution would approach the CB-Theor solution.For all cases, the CB algorithm converged on a solution in 5 minutes or less.
GA was run on each case a total of 100 times, each run itself produced 100 generations of 50 candidates.The leastcost candidate in that did not violate the static pressure head condition was chosen as the optimum.Because GA is a stochastic search algorithm producing different results from run-to-run, the costs of the optima from all 100 runs were averaged, with this averaged value presented in Table 2 as a percentage increase from the global optimum.Out of the 100 GA runs for each test case, nearly all runs failed to achieve the global optimum, with the exception of 8 runs of the Kiangan network.On a Dell Latitude (i5 CPU at 2.50 GHz), GA runs took between 0.9 s (case 1) and 1.2 s (case 3), not including about 2 s of pre-processor time.We note that many variations of GAs have been reported in the literature and several of these would likely improve upon the GA results obtained in this study.Potential improvements to the GA a self-adapting penalty function (Wu and Walski 2005), the use of elitism to preserve the best solutions (Kadu et al. 2008), and a reduction in the search space (Kadu et al. 2008).Other reported improvements, including the systematic optimization of operator parameters (Reed et al. 2000) and scaling of the fitness function to magnify the rewards towards slightly fitter candidates (Dandy et al. 1996) were attempted in a less systematic way (Ie.a parameter tuning study and three attempts at scaling the fitness function with an increasing exponent), but these did not result in a noticeable on performance.However, it is possible that if these techniques were followed systematically in full, the GA performance may have been improved.Still, the GA algorithm presented has undergone reasonable attempts to adapt its design and parameters for real-world GDWN cases, and therefore presents a useful point of comparison to the BT and CB algorithms.

Conclusions
Algorithms to optimize the cost of branching gravity-driven water networks are evaluated on five test cases from real networks in the Philippines, Nicaragua, and Panama.A calculus-based algorithm produced a solution composed of theoretical diameters from a continuous set (CB-Theor), which are then mapped onto discrete commercially available diameters (CB-Disc).Backtracking (BT), a recursive algorithm, systematically searches discrete candidate solutions and is guaranteed to find the global optimum by following rules that prune only higher-cost or hydraulically infeasible candidates.The BT algorithm was modified (BT-NoUp) to improve computational speed by also rejecting all candidates that included a small diameter directly upstream of a larger diameter.This criterion allowed BT-NoUp to prune more candidate solutions but allowed for the possibility of missing the global optimum.The third type of algorithm evaluated was a genetic algorithm (GA) that used singlepoint crossover and proportionate selection BT was able to find the global optimum in all test cases with relatively little computational effort, and could be applied to other GDWNs composed of a similar number of links.In addition, while BT-NoUp completed its search in less time than BT, the time required to complete BT would not be burdensome on a designer and therefore BT-NoUp did not produce a compelling relative advantage over BT.BT, however, could become prohibitively time-consuming when dealing with networks with significantly more links, as would be the case with large urban networks.While the test cases represent the range of GDWN sizes encountered in the authors' experience, future work would be needed to verify the suitability of the BT and BT-NoUp algorithms on other large GDWNs.The calculus-based algorithm produced consistently good results for the networks tested, although a more robust mapping scheme from theoretical diameters to discrete diameters would further improve on these results as discussed above.In potential future work, the CB-Theor solutions could be used to prune the BT search space, similar to Kadu et al. (2008), by only including the two diameters above and below the CB-Theor diameters, Δz (m) 25.0 38.9 11.9 42.1 -22.9 32.3 -29.9 40.8 -3.0 -14.7 34.1 -7.6 -5.0 20.0 -15.0 2.0 -12.0 14.0 -6.0 5.0 -1.0 -13.0 9.0 D solutions Suribabu, C. R.: Heuristic-based pipe dimensioning model for water distribution networks, J. Pipeline Syst.Eng.Pract., 3(4), 115-124, doi:10.1061/(ASCE)PS.1949-1204.0000104,2012.Zhao, W., Beach, T., and Rezgui, Y.: Optimization of Potable Water Distribution and Wastewater Collection Networks: A Systematic Review and Future Research Directions, IEEE Transactions on Systems, Man, and Cybernetics: Systems, 46 (5), 659-681, doi:10.1109/TSMC.2015.2461188,2016.
algorithms include (1) the calculus-based (CB) optimization model ofJones (2011), an NLP method, (2) backtracking (BT), a partial enumeration method, and (3) a genetic algorithm (GA), a metaheuristic method.Major distinguishing features of these algorithms include their working use of continuous diameters (CB) versus discrete diameters (BT and GA), their deterministic nature (CB and BT) versus a stochastic nature (GA), and their relative scalability as better (CB, GA) and worse (BT) for larger networks.In terms of their ability to find a global optimum solution for the problem formulation, CB finds a global optimum for continuous diameters but cannot guarantee a discrete diameter global optimum in its mapped solution, BT can guarantee a discrete global optimum, and GA cannot guarantee an optimum.For a direct comparison of techniques, the pipe costs used for all algorithms are found by interpolating a two-part cost formula based on a curve-fit of real cost data for available diameter values.The three algorithms are tested against networks adapted from field data on five actual GDWNs installed in Panama, Nicaragua, and the Philippines.
node, shown as tapstands or houses.For each link in a network of NL links, pipe length () and the net elevation change (Δ) are considered fixed.Steady-state flow rates () are prescribed for each link based on the demand flow data at delivery nodes.As noted above, demand flows are determined by community surveys and extrapolated in time to quantitatively account for population growth.Minor losses are accounted for through a minor loss coefficient  or a dimensionless equivalent pipe length, (  /, or in symbol form,   ), where Le is the pipe length of diameter D whose frictional loss results in the corresponding minor loss.An optimal solution is obtained by selecting pipe diameters () from a set of commercially available diameters such that the network's material cost is minimized.With   choices of diameters for   links, the problem has     candidate solutions.

Figure 6 :
Figure 6: Diameter sizes of calculus-based (CB-Disc) solutions above the global optimum solutions.
fit candidates have ability to be selected into successive rounds, albeit at a lower probability than fitter candidates.This ability to preserve less-fit candidates in the population increases genetic diversity, permitting genetic algorithms to avoid convergence into a local optimum.Alternates to proportionate selection include tournament selection, elitism, and stochastic universal sampling, which have been shown improved performance in genetic algorithms in certain cases of water network problems.
For the full energy equation, where  appears in a nonlinear than one location, this would be done using numerical methods.However, if we restrict our interest to minorlossless, smooth-turbulent flow as noted above, we can use the energy equations like Eq. (3).Obtain for pipe 1-2 ,  23 , and  24 .Introducing  12 ,  23 , and  24 from Eqs (3) into this single algebraic equation, we get to encompass two different ranges of pipe diameters having two different coefficients a and exponents b.The link between the two ranges starts at discrete pipe size Dco, at and below which the cost model for the small (subscript s) pipe sizes applies, and discrete pipe size Dco+1, at and above which the cost model for the large (subscript l ) pipe sizes applies.The cutoff diameter, Dco is chosen by the designer based on inspection of cost vs.
accounts for the effect of pipe roughness (smooth and commercial steel).The term  ′  is the derivative of the cost function per unit length with respect to /  .For the two-part cost model from above, obtain