Algorithms for optimization of branching gravity-driven water networks

The design of a water network involves the selection of pipe diameters that satisfy pressure and flow requirements while considering cost. A variety of design approaches can be used to optimize for hydraulic performance or reduce costs. To help designers select an appropriate approach in the context of gravity-driven water networks (GDWNs), this work assesses three cost-minimization algorithms on six moderate-scale 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 set. The backtracking algorithm finds the global optimum for all but the largest of cases tested, for which its long runtime makes it an infeasible option. The calculus-based algorithm's discrete-diameter solution produced slightly higher-cost results but was more scalable to larger network cases. Furthermore, the new calculus-based algorithm's continuous-diameter and mapped solutions provided lower and upper bounds, respectively, on the discrete-diameter global optimum cost, where the mapped solutions were typically within one diameter size of the global optimum. The genetic algorithm produced solutions even closer to the global optimum with consistently short run times, although slightly higher solution costs were seen for the larger network cases tested. 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 suboptimal candidates without missing the global optimum, and algorithm run time. We also extend an existing closed-form model of Jones (2011) to include minor losses and 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.

1 Introduction

A gravity-driven water network (GDWN) is commonly used to deliver potable water from a source at a high elevation, such as a natural spring or reservoir, to households or public tap stands (Fig. 1). When feasible, gravity-driven water networks are attractive in comparison to pumped networks because of their simplicity and lower capital, operational, and maintenance costs. In addition, in many locations where GDWN are considered, there may be little or no access to reliable grid-based electrical power for pumps. To improve reliability, networks may be designed with loops or multiple water sources, although often material cost considerations restrict attention to single-source branched networks.

Water networks are modeled as a collection of nodes, each representing a point of water demand or supply, which are connected with links representing pipes. The geometrical layout of the site is known and fixed, including water source and demand locations and elevations of all 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, 1983) 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 (positive) pressure heads are maintained at each node given the design flow rate at the node. Furthermore, application of the energy equation to each link in the network demonstrates that the design problem is nonunique; i.e., choosing different pressure heads at the nodes will result in a different pipe diameter solution for the network, and thus different network costs. Minimizing network cost will produce a unique solution to the design problem, i.e., unique link diameters and nodal pressure heads.

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 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 cost optimized 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 et al., 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).

Figure 1 Element schematic of a GDWN.

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 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 (Bhave, 1983; Kessler and Shamir, 1989; Swamee and Sharma, 2000; Samani and Mottaghi, 2006). LP can also get stuck in a local optimum (Zhao et al., 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 a 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 (turbulent flow in a smooth pipe) 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 development of this algorithm, as well as solution methodology, differs from that of Bhave (1978), which assumes constancy in several terms and thus requires iteration to solve. The Jones algorithm 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 is obtained which produces minimum network cost, as opposed to the distribution main cost as in Swamee and Sharma (2000). The method of Jones 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 modeling (Taormina and Chau, 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 et al., 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 the cost optimization of water distribution networks, and compare their performance on five cases adapted from real GDWNs. These algorithms include (1) the calculus-based (CB) optimization model of Jones (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 (CB and BT) versus stochastic (GA) search process, 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.

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 suboptimal candidates without missing the global optimum, and also computational time. We present two pre-processors which discrete-diameter search methods can use to reduce the search space without pruning the global optimum. To the authors' knowledge, this is the first implementation of “pre-processor 1” in enumeration methods and the first implementation of “pre-processor 2” in any water network design 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.

2 Problem formulation

Branching networks are considered (Fig. 1), where all branches connect a distribution main node with a delivery node, shown as tap stands or houses. For each link in a network of N_L links, pipe length (L) and the net elevation change (Δz) are considered known and fixed. Steady-state flow rates (Q) 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 K or a dimensionless equivalent pipe length, (L_e∕D, or in symbol form, L_ebyD), where L_e 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 (D) from a set of commercially available diameters such that the network's material cost is minimized. With N_D choices of diameters for N_L links, the problem has $N_{\mathrm{D}}^{N_{\mathrm{L}}}$ candidate solutions.

For all nodes, pressure head, h, is greater than or equal to a chosen minimum, h_min. The value for h_min 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 pressure head, Δh, across each link is calculated with the energy equation for pipe flow as follows:

$$\begin{array}{}\text{(1)}& {\displaystyle }\mathrm{\Delta }h=-\mathrm{\Delta }z+\left(\mathit{\alpha }+K+f\left({\displaystyle \frac{L}{D}}+L_{\mathrm{ebyD}}\right)\right){\displaystyle \frac{\mathrm{8}{Q}^{\mathrm{2}}}{{\mathit{\pi }}^{\mathrm{2}}g{D}^{\mathrm{4}}}},\end{array}$$

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 g 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,

$$\mathit{\alpha }=\left\{\begin{array}{ll}\mathrm{2}& \mathit{\text{Re}}\le \mathrm{2300}\\ \mathrm{1.05}& \mathit{\text{Re}}>\mathrm{2300}\end{array},\right.$$

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.

3 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 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 but the focus for the present work is on branching 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 be addressed below.

Consider the three-pipe network shown in Fig. 2. Pipes 1–2, 2–3, and 2–4 meet where head h₂ is unknown. Each pipe has a prescribed volume flow rate and length and unknown diameter D as shown. The change in elevation between the top and bottom of each pipe is Δz, and Δh is the change in pressure head. There is a prescribed head at each outlet for pipes 2–3 and 2–4.

Figure 2 Three-pipe branch network.

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=\mathrm{0.316}{\mathit{\text{Re}}}^{-\mathrm{1}/\mathrm{4}}$, and the Swamee–Jain correlation (Swamee and Jain, 1976), $f=\mathrm{0.175}{\mathit{\text{Re}}}^{-\mathrm{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 aReⁿ). The Blasius equation has higher accuracy (2 % for low Re and 3 % for high Re) in the range 10⁴ < Re < 10⁵, over which most of the GDWNs in this work operate, compared with the Swamee–Jain correlation of +8 %/−3 %, thus the Blasius equation is chosen for this work. A combination of the Blasius equation with the energy equation gives explicit formulas for D 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 therefore obtain

$$\begin{array}{ll}\text{(2)}& {\displaystyle }& {\displaystyle }{D}_{\mathrm{12}}=\mathrm{0.741}{\left({\displaystyle \frac{\mathrm{\Delta }{z}_{\mathrm{12}}+\mathrm{\Delta }{h}_{\mathrm{12}}}{L_{\mathrm{1}}}}\right)}^{-\mathrm{4}/\mathrm{19}}{\left({\displaystyle \frac{Q_{\mathrm{12}}{\mathit{\nu }}^{\mathrm{1}/\mathrm{7}}}{g^{\mathrm{4}/\mathrm{7}}}}\right)}^{\mathrm{7}/\mathrm{19}}{\displaystyle }& {\displaystyle }{D}_{\mathrm{23}}=D_{\mathrm{24}}=\mathrm{0.741}{\left({\displaystyle \frac{\mathrm{\Delta }{z}_{\mathrm{23}}+\mathrm{\Delta }{h}_{\mathrm{23}}}{L_{\mathrm{2}}}}\right)}^{-\mathrm{4}/\mathrm{19}}{\left({\displaystyle \frac{Q_{\mathrm{23}}{\mathit{\nu }}^{\mathrm{1}/\mathrm{7}}}{g^{\mathrm{4}/\mathrm{7}}}}\right)}^{\mathrm{7}/\mathrm{19}}.\end{array}$$

With our assumptions and inspection of Fig. 2, $\mathrm{\Delta }{h}_{\mathrm{12}}=-h_{\mathrm{2}}$ and $\mathrm{\Delta }{h}_{\mathrm{23}}=\mathrm{\Delta }{h}_{\mathrm{24}}=h_{\mathrm{2}}-h_{\mathrm{3}}=h_{\mathrm{2}}-h_{\mathrm{4}}$, we furthermore obtain

$$\begin{array}{ll}\text{(3)}& {\displaystyle }& {\displaystyle }{D}_{\mathrm{12}}=\mathrm{0.741}{\left({\displaystyle \frac{\mathrm{\Delta }{z}_{\mathrm{12}}-h_{\mathrm{2}}}{L_{\mathrm{1}}}}\right)}^{-\mathrm{4}/\mathrm{19}}{\left({\displaystyle \frac{Q_{\mathrm{12}}{\mathit{\nu }}^{\mathrm{1}/\mathrm{7}}}{g^{\mathrm{4}/\mathrm{7}}}}\right)}^{\mathrm{7}/\mathrm{19}}{\displaystyle }& {\displaystyle }{D}_{\mathrm{23}}=D_{\mathrm{24}}=\mathrm{0.741}{\left({\displaystyle \frac{\mathrm{\Delta }{z}_{\mathrm{23}}-h_{\mathrm{3}}+h_{\mathrm{2}}}{L_{\mathrm{2}}}}\right)}^{-\mathrm{4}/\mathrm{19}}{\left({\displaystyle \frac{Q_{\mathrm{23}}{\mathit{\nu }}^{\mathrm{1}/\mathrm{7}}}{g^{\mathrm{4}/\mathrm{7}}}}\right)}^{\mathrm{7}/\mathrm{19}}.\end{array}$$

The single-part pipe-cost model can be assumed to follow a power-law relationship (Swamee and Sharma, 2008)

$$\begin{array}{}\text{(4)}& C=a{\left({\displaystyle \frac{D}{D_{\mathrm{u}}}}\right)}^b,\end{array}$$

where C is cost per unit length of pipe, a is a constant coefficient, b is a constant exponent, and D_u 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, such that the material cost for the network is of prime importance. 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, C_T, is obtained by summing over all links ij,

$$\begin{array}{}\text{(5)}& C_{\mathrm{T}}=a\sum _{ij}{L}_{ij}{\left({\displaystyle \frac{D_{ij}}{D_{\mathrm{u}}}}\right)}^b,\end{array}$$

which, for the present problem, becomes

$$\begin{array}{ll}{\displaystyle }{C}_{\mathrm{T}}& {\displaystyle }=a\left[L_{\mathrm{12}}{\left({\displaystyle \frac{D_{\mathrm{12}}}{D_{\mathrm{u}}}}\right)}^b+L_{\mathrm{23}}{\left({\displaystyle \frac{D_{\mathrm{23}}}{D_{\mathrm{u}}}}\right)}^b+L_{\mathrm{24}}{\left({\displaystyle \frac{D_{\mathrm{24}}}{D_{\mathrm{u}}}}\right)}^b\right]\\ \text{(6)}& {\displaystyle }& {\displaystyle }=a\left[L_{\mathrm{12}}{\left({\displaystyle \frac{D_{\mathrm{12}}}{D_{\mathrm{u}}}}\right)}^b+\mathrm{2}{L}_{\mathrm{23}}{\left({\displaystyle \frac{D_{\mathrm{23}}}{D_{\mathrm{u}}}}\right)}^b\right].\end{array}$$

The mathematical basis for a unique solution for h₂ with cost minimization is now presented. In addition to the fixed pipe lengths, the total cost depends on the diameters for all pipes in the network. For the case of Fig. 2, where we now allow pipe 2–3 and pipe 2–4 to be different, we get

$$\begin{array}{}\text{(7)}& C_{\mathrm{T}}=C_{\mathrm{T}}\left(D_{\mathrm{12}}\left(h_{\mathrm{2}}\right),D_{\mathrm{23}}\left(h_{\mathrm{2}}\right),D_{\mathrm{24}}\left(h_{\mathrm{2}}\right)\right).\end{array}$$

Using the chain rule from the calculus, the total differential of Eq. (7) is

$$\begin{array}{ll}{\displaystyle }\mathrm{d}{C}_{\mathrm{T}}=& {\displaystyle }{\displaystyle \frac{\partial C_{\mathrm{T}}}{\partial D_{\mathrm{12}}}}{\displaystyle \frac{\partial D_{\mathrm{12}}}{\partial h_{\mathrm{2}}}}\mathrm{d}{h}_{\mathrm{2}}+{\displaystyle \frac{\partial C_{\mathrm{T}}}{\partial D_{\mathrm{23}}}}{\displaystyle \frac{\partial D_{\mathrm{23}}}{\partial h_{\mathrm{2}}}}\mathrm{d}{h}_{\mathrm{2}}\\ \text{(8)}& {\displaystyle }& {\displaystyle }+{\displaystyle \frac{\partial C_{\mathrm{T}}}{\partial D_{\mathrm{24}}}}{\displaystyle \frac{\partial D_{\mathrm{24}}}{\partial h_{\mathrm{2}}}}\mathrm{d}{h}_{\mathrm{2}}.\end{array}$$

The minimum value of C_T is found once dC_T=0 (and once it is verified that the second derivative of C_T is positive thus indicating that C_T is indeed a minimum). Requiring this, we obtain

$$\begin{array}{}\text{(9)}& {\displaystyle }\mathrm{0}={\displaystyle \frac{\partial C_{\mathrm{T}}}{\partial D_{\mathrm{12}}}}{\displaystyle \frac{\partial D_{\mathrm{12}}}{\partial h_{\mathrm{2}}}}+{\displaystyle \frac{\partial C_{\mathrm{T}}}{\partial D_{\mathrm{23}}}}{\displaystyle \frac{\partial D_{\mathrm{23}}}{\partial h_{\mathrm{2}}}}+{\displaystyle \frac{\partial C_{\mathrm{T}}}{\partial D_{\mathrm{24}}}}{\displaystyle \frac{\partial D_{\mathrm{24}}}{\partial h_{\mathrm{2}}}}.\end{array}$$

The cost C_T is from Eq. (5), so the derivatives like $\partial C_{\mathrm{T}}/\partial D_{\mathrm{12}}$ in Eq. (9) are written in general as

$$\begin{array}{}\text{(10)}& {\displaystyle \frac{\partial C_{\mathrm{T}}}{\partial D_{ij}}}=ab{\displaystyle \frac{D_{ij}^{b-\mathrm{1}}}{D_{\mathrm{u}}^b}}{L}_{ij}\end{array}$$

for any link ij.

The derivatives like $\partial D_{\mathrm{12}}/\partial h_{\mathrm{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. For the full energy equation, where D appears in a nonlinear way in more than one location, this would be done using numerical methods. However, if we assume minor-lossless, smooth-turbulent flow as noted above, we can use the energy equations like Eq. (3). We therefore obtain the following for pipe 1–2:

$$\begin{array}{}\text{(11)}& {\displaystyle \frac{\partial D_{\mathrm{12}}}{\partial h_{\mathrm{2}}}}=\mathrm{0.156}{\left({\displaystyle \frac{\mathrm{\Delta }{z}_{\mathrm{12}}-h_{\mathrm{2}}}{L_{\mathrm{12}}}}\right)}^{\frac{-\mathrm{23}}{\mathrm{19}}}{\left({\displaystyle \frac{{\mathit{\nu }}^{\mathrm{1}/\mathrm{7}}{Q}_{\mathrm{12}}}{g^{\mathrm{4}/\mathrm{7}}{L}_{\mathrm{12}}^{\mathrm{19}/\mathrm{7}}}}\right)}^{\frac{\mathrm{7}}{\mathrm{19}}};\end{array}$$

for pipe 2–3, we get

$$\begin{array}{ll}{\displaystyle }& {\displaystyle \frac{\partial D_{\mathrm{23}}}{\partial h_{\mathrm{2}}}}=\\ \text{(12)}& {\displaystyle }& {\displaystyle }-\mathrm{0.156}{\left({\displaystyle \frac{\mathrm{\Delta }{z}_{\mathrm{23}}+h_{\mathrm{2}}-h_{\mathrm{3}}}{L_{\mathrm{23}}}}\right)}^{\frac{-\mathrm{23}}{\mathrm{19}}}{\left({\displaystyle \frac{{\mathit{\nu }}^{\mathrm{1}/\mathrm{7}}{Q}_{\mathrm{23}}}{g^{\mathrm{4}/\mathrm{7}}{L}_{\mathrm{23}}^{\mathrm{19}/\mathrm{7}}}}\right)}^{\frac{\mathrm{7}}{\mathrm{19}}};\end{array}$$

and for pipe 2–4,

$$\begin{array}{}\text{(13)}& {\displaystyle \frac{\partial D_{\mathrm{24}}}{\partial h_{\mathrm{2}}}}=-\mathrm{0.156}{\left({\displaystyle \frac{\mathrm{\Delta }{z}_{\mathrm{24}}+h_{\mathrm{2}}-h_{\mathrm{4}}}{L_{\mathrm{24}}}}\right)}^{\frac{\mathrm{23}}{\mathrm{19}}}{\left({\displaystyle \frac{{\mathit{\nu }}^{\frac{\mathrm{1}}{\mathrm{7}}}{Q}_{\mathrm{24}}}{g^{\frac{\mathrm{4}}{\mathrm{7}}}{L}_{\mathrm{24}}^{\frac{\mathrm{19}}{\mathrm{7}}}}}\right)}^{\frac{\mathrm{7}}{\mathrm{19}}}.\end{array}$$

Equations (10)–(13) are combined with Eq. (9) to produce a single algebraic equation that depends on h₂, as well as D₁₂, D₂₃, and D₂₄. Introducing D₁₂,D₂₃, and D₂₄ from Eq. (3) into this algebraic equation, we get

$$\begin{array}{ll}{\displaystyle }\mathrm{0}=& {\displaystyle }{Q}_{\mathrm{12}}^{\mathrm{7}b/\mathrm{19}}{\left({\displaystyle \frac{\mathrm{\Delta }{z}_{\mathrm{12}}-h_{\mathrm{2}}}{L_{\mathrm{12}}}}\right)}^{-(\mathrm{1}+\mathrm{4}b/\mathrm{19})}\\ {\displaystyle }& {\displaystyle }-Q_{\mathrm{23}}^{\mathrm{7}b/\mathrm{19}}{\left({\displaystyle \frac{\mathrm{\Delta }{z}_{\mathrm{23}}+h_{\mathrm{2}}-h_{\mathrm{3}}}{L_{\mathrm{23}}}}\right)}^{-(\mathrm{1}+\mathrm{4}b/\mathrm{19})}\\ \text{(14)}& {\displaystyle }& {\displaystyle }-Q_{\mathrm{24}}^{\mathrm{7}b/\mathrm{19}}{\left({\displaystyle \frac{\mathrm{\Delta }{z}_{\mathrm{24}}+h_{\mathrm{2}}-h_{\mathrm{4}}}{L_{\mathrm{24}}}}\right)}^{-(\mathrm{1}+\mathrm{4}b/\mathrm{19})}.\end{array}$$

The general form of Eq. (14), written at any internal node is

$$\begin{array}{}\text{(15)}& \mathrm{0}=\sum _{ij,\mathrm{in}}{Q}_{ij}^{\mathrm{7}b/\mathrm{19}}{S}_{ij}^{-(\mathrm{1}+\mathrm{4}b/\mathrm{19})}-\sum _{ij,\mathrm{out}}{Q}_{ij}^{\mathrm{7}b/\mathrm{19}}{S}_{ij}^{-(\mathrm{1}+\mathrm{4}b/\mathrm{19})},\end{array}$$

where the hydraulic gradient, S_ij, is

$$\begin{array}{}\text{(16)}& S_{ij}={\displaystyle \frac{\mathrm{\Delta }{z}_{ij}+\mathrm{\Delta }{h}_{ij}}{L_{ij}}}.\end{array}$$

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 D_ij for all links and h_j 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 h_min prescribed for the network. For nodes that do not satisfy this condition, the pressure head is set equal to h_min, as part of the CB algorithm. Thus, $h_j\ge h_{min}$.

Minor losses using the equivalent-length method can be included in the above developments by artificially extending the length of the link by L_e in which minor loss occurs, thus contributing a non-zero L_ebyD term in Eq. (1). We also extend the cost model of Eq. (5) from Swamee and Sharma (2008) 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 D_co, at and below which the cost model for the small (subscript s) pipe sizes applies, and discrete pipe size D_co+1, at and above which the cost model for the large (subscript l) pipe sizes applies. The cutoff diameter, D_co is chosen by the designer based on inspection of cost vs. diameter data. Thus,

$$\begin{array}{}\text{(17)}& C_{ij}=L_{ij}\left\{\begin{array}{l}{a}_{\mathrm{s}}{\left({\displaystyle \frac{D_{ij}}{D_{\mathrm{u}}}}\right)}^{b_{\mathrm{s}}},\\ D_{ij}\le D_{\mathrm{co}}\\ c_{\mathrm{1}}+c_{\mathrm{2}}{\displaystyle \frac{D_{ij}}{D_{\mathrm{u}}}}+c_{\mathrm{3}}{\left({\displaystyle \frac{D_{ij}}{D_{\mathrm{u}}}}\right)}^{\mathrm{2}}+c_{\mathrm{4}}{\left({\displaystyle \frac{D_{ij}}{D_{\mathrm{u}}}}\right)}^{\mathrm{3}},\\ D_{\mathrm{co}}<D_{ij}<D_{\mathrm{co}+\mathrm{1}}\\ a_{\mathrm{l}}{\left({\displaystyle \frac{D_{ij}}{D_{\mathrm{u}}}}\right)}^{b_{\mathrm{l}}},\\ D_{ij}\ge D_{\mathrm{co}+\mathrm{1}}.\end{array}\right.\end{array}$$

In Eq. (17), a_s and a_l are the coefficients for the small (designated by subscript s) and large (subscript l) pipe size regions, respectively, and b_s and b_l are the exponents for the small and large pipe size regions, respectively. A cubic spline is fit between pipe sizes D_co and D_co+1 to complete the transition between small and large pipe sizes. The coefficients of this polynomial are c₁, c₂, c₃, and c₄ as seen in Eq. (17). These coefficients are evaluated by matching the cubic polynomial and pipe data at D_co and D_co+1 and the first derivative of the polynomial with respect to D_ij∕D_u to $a_{\mathrm{s}}{b}_{\mathrm{s}}{\left(\frac{D_{\mathrm{co}}}{D_{\mathrm{u}}}\right)}^{b_{\mathrm{s}}-\mathrm{1}}$ at D_ij=D_co and to $a_{\mathrm{l}}{b}_{\mathrm{l}}{\left(\frac{D_{\mathrm{co}+\mathrm{1}}}{D_{\mathrm{u}}}\right)}^{b_{\mathrm{l}}-\mathrm{1}}$ at $D_{ij}=D_{\mathrm{co}+\mathrm{1}}$. An example of data for polyvinyl chloride (PVC) pipe and the curve fit is shown in Fig. 3. The results of the curve fit are as follows: D_co=2.067 in., $D_{\mathrm{co}+\mathrm{1}}=\mathrm{2.469}$ in., a_s = USD 1.349 m⁻¹, b_s=1.157, a_l = USD 1.381 m⁻¹, b_l=1.344, c₁ = USD 237.516 m⁻¹, c₂ = USD 316.125 m⁻¹, c₃ = USD 140.450 m⁻¹, c₄ = −USD 20.499 m⁻¹. It is clear from inspection of Fig. 3 that a one-part cost model would not have produced an acceptable curve-fit to pipe-cost data.

With the inclusion of the two-part cost model and minor loss term, Eq. (15) becomes

$$\begin{array}{}\text{(18)}& \mathrm{0}=\begin{array}{l}\sum _{ij,\mathrm{in}}{\displaystyle \frac{C_{ij}^{\prime }{A}_{ij}^{\frac{\mathrm{4}}{\mathrm{19}}}{\left(\mathrm{1}+{\mathit{ϵ}}_{ij}\right)}^{\frac{\mathrm{4}}{\mathrm{19}}}{S}_{ij}^{-\frac{\mathrm{23}}{\mathrm{19}}}{\left(\frac{Q_{ij}^{\mathrm{7}}\mathit{\nu }}{g^{\mathrm{4}}{D}_{\mathrm{u}}^{\mathrm{19}}}\right)}^{\frac{\mathrm{1}}{\mathrm{19}}}}{\mathrm{1}-B{A}_{ij}^{\frac{\mathrm{4}}{\mathrm{19}}}{\mathit{ϵ}}_{ij}^{\prime }{\left(\mathrm{1}+{\mathit{ϵ}}_{ij}\right)}^{-\frac{\mathrm{15}}{\mathrm{19}}}{S}_{ij}^{-\frac{\mathrm{4}}{\mathrm{19}}}{\left(\frac{Q_{ij}^{\mathrm{7}}\mathit{\nu }}{g^{\mathrm{4}}{D}_{\mathrm{u}}^{\mathrm{19}}}\right)}^{\frac{\mathrm{1}}{\mathrm{19}}}}}\\ -\sum _{ij,\mathrm{out}}{\displaystyle \frac{C_{ij}^{\prime }{A}_{ij}^{\frac{\mathrm{4}}{\mathrm{19}}}{\left(\mathrm{1}+{\mathit{ϵ}}_{ij}\right)}^{\frac{\mathrm{4}}{\mathrm{19}}}{S}_{ij}^{-\frac{\mathrm{23}}{\mathrm{19}}}{\left(\frac{Q_{ij}^{\mathrm{7}}\mathit{\nu }}{g^{\mathrm{4}}{D}_{\mathrm{u}}^{\mathrm{19}}}\right)}^{\mathrm{1}/\mathrm{19}}}{\mathrm{1}-B{A}_{ij}^{\frac{\mathrm{4}}{\mathrm{19}}}{\mathit{ϵ}}_{ij}^{\prime }{\left(\mathrm{1}+{\mathit{ϵ}}_{ij}\right)}^{-\frac{\mathrm{15}}{\mathrm{19}}}{S}_{ij}^{-\frac{\mathrm{4}}{\mathrm{19}}}{\left(\frac{Q_{ij}^{\mathrm{7}}\mathit{\nu }}{g^{\mathrm{4}}{D}_{\mathrm{u}}^{\mathrm{19}}}\right)}^{\mathrm{1}/\mathrm{19}}}},\end{array}\end{array}$$

where B=0.1989 and

$$\begin{array}{ll}\text{(19)}& {\displaystyle }& {\displaystyle }{\mathit{ϵ}}_{ij}=\sum _k{\left({\displaystyle \frac{L_{\mathrm{e}}}{D}}\right)}_{k,ij}{\displaystyle \frac{D_{ij}}{L_{ij}}}{\displaystyle }& {\displaystyle }{\mathit{ϵ}}_{ij}^{\prime }=\sum _k{\left({\displaystyle \frac{L_{\mathrm{e}}}{D}}\right)}_{k,ij}{\displaystyle \frac{D_{\mathrm{u}}}{L_{ij}}}\\ \text{(20)}& {\displaystyle }& {\displaystyle }{A}_{ij}=\left\{\begin{array}{ll}\mathrm{0.318},& \mathrm{smooth}\mathrm{pipe}\\ \mathrm{0.420},& \mathrm{steel}\mathrm{pipe}\end{array}\right.\end{array}$$

and A accounts for the effect of pipe roughness (smooth and commercial steel). The term $C_{ij}^{\prime }$ is the derivative of the cost function per unit length with respect to D∕D_u. For the two-part cost model from above, we obtain

$$\begin{array}{}\text{(21)}& C_{ij}^{\prime }=\left\{\begin{array}{l}{a}_{\mathrm{s}}{b}_{\mathrm{s}}{\left({\displaystyle \frac{D_{ij}}{D_{\mathrm{u}}}}\right)}^{b_{\mathrm{s}}-\mathrm{1}},\\ D_{ij}\le D_{\mathrm{co}}\\ c_{\mathrm{2}}+\mathrm{2}{c}_{\mathrm{3}}\left({\displaystyle \frac{D_{ij}}{D_{\mathrm{u}}}}\right)+\mathrm{3}{c}_{\mathrm{4}}{\left({\displaystyle \frac{D_{ij}}{D_{\mathrm{u}}}}\right)}^{\mathrm{2}},\\ D_{\mathrm{co}}<D_{ij}<D_{\mathrm{co}+\mathrm{1}}\\ a_{\mathrm{l}}{b}_{\mathrm{l}}{\left({\displaystyle \frac{D_{ij}}{D_{\mathrm{u}}}}\right)}^{b_{\mathrm{l}}-\mathrm{1}},\\ D_{ij}\ge D_{\mathrm{co}+\mathrm{1}}.\end{array}\right.\end{array}$$

Figure 3 PVC pipe cost from 2011 data.

Equation (18), and its simpler form Eq. (15) for minor-lossless flow and a single-part pipe-cost model (it is easy to show that Eq. (18) regresses to Eq. (15) for these conditions), is the root of the calculus-based optimization in this work and is applied at all internal nodes to uniquely determine h_j. 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 31 January 2018). 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) using slightly different notation than here. For clarity, we re-present Eq. (15) using Bhave's notation as

$$\mathrm{0}=\sum Q_{ij}^{\mathrm{7}b/\mathrm{19}}{S}_{ij}^{-(\mathrm{1}+\mathrm{4}b/\mathrm{19})}-\sum Q_{jk}^{\mathrm{7}b/\mathrm{19}}{S}_{jk}^{-(\mathrm{1}+\mathrm{4}b/\mathrm{19})},$$

where the ij and jk notation are shown in Fig. 4. Index j spans all internal nodes along the distribution main. A quantitative comparison between Eq. (18) and the method of Bhave is presented below.

4 Backtracking algorithm and genetic algorithm

Figure 4 Bhave (1978) index notation at an internal node, j.

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 N_L elements where each element corresponds to a commercially available diameter of a network link. 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 L, Δz, K, L_ebyD, and α, the change in pressure head Δh is evaluated for all N_D×N_L combinations of pipe diameter and link index:

$$\begin{array}{}\text{(22)}& \mathbf{\Delta }\mathbf{h}=\left[\begin{array}{ccc}\mathrm{\Delta }{h}_{\mathrm{11}}& \mathrm{\cdots }& \mathrm{\Delta }{h}_{\mathrm{1}{N}_{\mathrm{L}}}\\ \mathrm{⋮}& \mathrm{\ddots }& \mathrm{⋮}\\ \mathrm{\Delta }{h}_{N_{\mathrm{D}}\mathrm{1}}& \mathrm{\cdots }& \mathrm{\Delta }{h}_{N_{\mathrm{D}}{N}_{\mathrm{L}}}\end{array}\right].\end{array}$$

While an algorithm evaluates a candidate solution, the 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 h_i for a given node by adding the head at the upstream node, h_i−1 to the change in head for that link i_L and the diameter i_D under consideration. Thus,

$$\begin{array}{}\text{(23)}& h_i=h_{i-\mathrm{1}}+\mathbf{\Delta }\mathbf{h}\left(i_{\mathrm{D}},i_{\mathrm{L}}\right).\end{array}$$

Along with the hydraulic evaluation of a candidate solution, the cost of the partial candidate is found through the use of a lookup table C:

$$\begin{array}{}\text{(24)}& \mathbf{C}=\left[\begin{array}{ccc}{C}_{\mathrm{11}}& \mathrm{\cdots }& C_{\mathrm{1}{N}_{\mathrm{L}}}\\ \mathrm{⋮}& \mathrm{\ddots }& \mathrm{⋮}\\ C_{N_{\mathrm{D}}\mathrm{1}}& \mathrm{\cdots }& C_{N_{\mathrm{D}}{N}_{\mathrm{L}}}\end{array}\right],\end{array}$$

where C(i_Di_L) returns the additional cost of assigning a diameter with index i_D to link i_L. 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 pressure head and cost upon consideration of each partial candidate, where GA calculates these values on full candidates as part of the objective function.

4.1 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. For the BT algorithm in particular, 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 pre-processor is used to provide a maximum diameter (D_max) 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 pressure head below h_min, 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 D_max to allow the algorithm to select a larger than necessary diameter if this is able to save cost elsewhere. 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.

4.2 BT and GA pre-processor 2: adjusted minimum pressure head

A second pre-processor adjusts the minimum pressure head requirement for each internal node by considering the total head required at downstream nodes. It can be recognized that, without the use of a pump, the total head cannot increase at nodes downstream of a given node i. Furthermore, the total head must decline at a minimum grade that is determined by the demand volume flow rate and the largest pipe diameter available (D_max) for selection. This energy constraint is utilized to reduce the number of candidates to be considered by increasing the minimum pressure head at nodes where these rules produce a higher minimum head than the original h_min. For example, nodes upstream of a local network high point can have their minimum pressure head increased beyond the normal minimum, since the pressure head must be great enough to ensure adequate flow to the higher elevation downstream node. To begin this process, each node i is initialized with a baseline minimum total head:

$$\begin{array}{}\text{(25)}& t{h}_{min,i}=z_i+h_{min}.\end{array}$$

$t{h}_{min,i}$ is thus initialized by considering only the node's hydraulic requirements in isolation, i.e., without acknowledging the neighboring downstream nodes. The pre-processor then considers updating $t{h}_{min,i}$ by checking the following condition, which is false when the minimum pressure head at downstream nodes produces further constraints on an upstream node i. Thus, for all nodes i which are upstream of some node j, the following inequality can be evaluated:

$$\begin{array}{ll}\text{(26)}& {\displaystyle }& {\displaystyle }t{h}_{min,i}-t{h}_{min,j}\ge {\displaystyle }& {\displaystyle }\left({\mathit{\alpha }}_{i-j}+K_{i-j}+f_{i-j}\left({\displaystyle \frac{L_{i-j}}{D_{i-j}}}+L_{{\mathrm{ebyD}}_{i-j}}\right)\right){\displaystyle \frac{\mathrm{8}{Q}_{i-j}^{\mathrm{2}}}{{\mathit{\pi }}^{\mathrm{2}}g{D}_{max}^{\mathrm{4}}}}.\end{array}$$

Also, consider that when flow rate Q_i−j is small and D_max is large, the right-hand side of Eq. (26) approaches zero, representing the simple statement that upstream total head must always be greater than downstream total head. When the condition in Eq. (26) is false, the minimum total head can be updated in node i such that the maximum diameter size in link i−j is able to meet the downstream node's minimum total head, or

$$\begin{array}{ll}\text{(27)}& {\displaystyle }& {\displaystyle }t{h}_{min,i}=t{h}_{min,j}+{\displaystyle }& {\displaystyle }\left({\mathit{\alpha }}_{i-j}+K_{i-j}+f_{i-j}\left({\displaystyle \frac{L_{i-j}}{D_{i-j}}}+L_{{\mathrm{ebyD}}_{i-j}}\right)\right){\displaystyle \frac{\mathrm{8}{Q}_{i-j}^{\mathrm{2}}}{{\mathit{\pi }}^{\mathrm{2}}g{D}_{max}^{\mathrm{4}}}}.\end{array}$$

In this way, $t{h}_{min,i}$ may be updated for each node until the condition in Eq. (26) is true for all nodes i with a downstream node j connected by a single link.

After the values for $t{h}_{min,i}$ are updated, they are converted back into minimum pressure head values by subtracting the elevation z_i from $t{h}_{min,i}$. This pre-processor serves to narrow the search for viable candidate solutions by potentially increasing the minimum pressure head. Since backtracking and GA consider network links in the downstream direction, these algorithms are otherwise blind to future downstream pressure head requirements. This limitation is alleviated by the pre-processor, which allows these algorithms some implicit information about what local diameter choices will be viable for the full network solution. Note that both pre-processors discussed will not prune the global optimum from the solution.

4.3 Backtracking algorithm (BT)

The backtracking algorithm is a partial enumeration method that 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 in the literature are those by Gessler (1985) and González-Cebollada et al. (2011). The algorithm proposed by Gessler proposes a pipe-grouping strategy which speeds up the algorithm but risks pruning the global optimum. Additionally, pipe grouping represents its own optimization problem (Raad, 2011). The González-Cebollada algorithm does not include such pipe-grouping criteria, although it does halt its search after finding the first feasible solution, thus it does not guarantee a global optimum. The present study's BT algorithm, once run to completion, does guarantee a global optimum. It operates similarly to the method presented by González-Cebollada et al. (2011), with the major differences being that the algorithm continues searching once it has found its first feasible solution and uses pre-processors 1 and 2 to further reduce the search space. This implementation of BT, however, scales poorly with larger network sizes and would not be appropriate for use on large urban networks. Its appropriateness is shown here for many of the GDWNs encountered in practice, as evidenced by its use on real-world GDWN test cases in this paper. Moreover, it serves as a benchmark against which other algorithms can be compared.

BT uses two rejection criteria to discard candidate solutions from further consideration. The first rejection criterion is that when a candidate violates 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 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 second 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.

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 pressure head and the partial candidate cost at the outlet node are calculated with the Δh and C lookup tables. If this partial candidate meets 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 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. In this way, the algorithm continues to extend and reject candidate solutions until a full candidate satisfies the pressure head requirements. 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. While the minimum pressure head criterion tends to prune candidates with diameters that are too small, the cost-based criterion tends to prune candidates of diameters that are too large.

4.4 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, a network designer would not consider such designs, 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 downstream from smaller ones. For this reason, the BT-NoUp algorithm, unlike the BT algorithm, may miss the global optimum at the expense of its greater computational efficiency.

4.5 Genetic algorithm (GA)

Genetic algorithms are stochastic optimization techniques that mimic the process of natural selection, and numerous variations of GAs have demonstrated acceptable performance on WDN design (Nicklow et al., 2010). Given their popularity, the GA included in this study is meant to provide a point of comparison to the BT and CB algorithms when applied to GDWNs.

When implemented in water network design, each candidate solution represents a selection of pipe diameters. The algorithm is initialized with a population of candidates of size N_c that repeatedly undergoes the processes of mutation, crossover, and selection

$$\begin{array}{}\text{(28)}& c_i=\left[\begin{array}{cccc}{D}_{\mathrm{1},i}& D_{\mathrm{2},i}& \mathrm{\dots }& D_{N_{\mathrm{L}},i}\end{array}\right],\end{array}$$

where each candidate in the population c_i contains N_L diameters. In the present work, candidates are represented as a string of natural numbers, which is used over a binary representation to improve the ease of encoding (Vairavamoorthy and Ali, 2000). The mutation operator replaces pipe diameters with a diameter from a uniform random distribution, where each link diameter has a probability of p_mut of mutating on each generation. The crossover operator randomly pairs individuals in the population with probability p_xover and performs a single-point crossover of the two individuals, where the point of crossover is randomly chosen. The fitness, f_i, of each candidate is assessed with penalties associated with the solution's pipe cost, C_pipe,i, and hydraulic cost, C_hyd,i, which is assigned when violations of the pressure head requirements occur:

$$\begin{array}{}\text{(29)}& f_i={\displaystyle \frac{\mathrm{1}}{C_{\mathrm{pipe},i}+C_{\mathrm{hyd},i}}}.\end{array}$$

The hydraulic cost is found for each individual by identifying nodes in which the pressure head is less than h_min and multiplying the total amount of head violation by a hydraulic penalty coefficient, a_hyd:

$$\begin{array}{}\text{(30)}& C_{\mathrm{hyd},i_C}=a_{\mathrm{hyd}}\sum _{\mathrm{1}}^{N_{\mathrm{L}}}\left(h_{min}-h_{i_N}\right)\mathrm{|}{h}_{i_N}<h_{min}.\end{array}$$

To allow for a hydraulic penalty coefficient to produce similar results in both small-scale (inexpensive) network and a large-scale (more expensive) cases, the hydraulic penalty coefficient is made directly proportional to the average solution cost. With each generation, a_hyd is updated by multiplying the normalized penalty coefficient, a_hyd,norm, by the average pipe cost of the population,

$$\begin{array}{}\text{(31)}& a_{\mathrm{hyd}}=a_{\mathrm{hyd},\mathrm{norm}}{\displaystyle \frac{\sum _{\mathrm{1}}^{N_{\mathrm{c}}}{C}_{\mathrm{pipe},i_C}}{N_{\mathrm{c}}}}.\end{array}$$

The algorithm then selects candidates to be carried into the next generation with a tournament selection method, where N_c groups of s individuals are randomly assigned and the fittest candidate among each group is selected, thus replacing the previous population with an equally sized population of N_c individuals.

Table 1 Characteristics of test cases.

In this study, the genetic algorithm parameters used were N_c=200, p_mut=0.05, p_xover=1, N_gen=500, $a_{\mathrm{hyd},\mathrm{norm}}=\mathrm{0.05}$, and s=10. These parameters were chosen by systematically varying parameter values until the optimum cost of a network, case 2, could no longer be significantly improved. The first four of these values are in line with typically used values from Simpson et al. (1994) of N_c (30–200), p_mut (0.01–0.05), p_xover (0.7–1.0), and N_gen (100–1000).

5 Cases studied

Six 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 4a–f. Each network is a branching type without loops. The total lengths of the networks range from less than 1 to over 15 km. Two serial networks are tested to demonstrate the effect of a local high point on the algorithm solutions. Elevation plots for each case are shown in Fig. 5.

The choice of h_min 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 h_min in GDWN design is typically in the range of 5–20 m (Arnalich, 2010; Bouman, 2014; Swamee and Sharma, 2008). In the present study h_min = 7 m, although this requirement was reduced at selected nodes at the beginning of a network where changes in elevation are still small (case 2, where the pressure head at node 2 is relaxed to 2 m). At the source node, the 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 described above. All algorithms were run in a late-version of MATLAB (or Mathcad for CB) on an Intel i5 processor at 2.50 GHz.

Figure 5 Network elevation (z) and hydraulic grade lines (HGLs) of algorithm solution for main distribution links.

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 in 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 in a step-by-step manner. A nearby plot of the pressure heads compared with the theoretical D_ij from the CB approach (e.g., on the same Mathcad page for our solution) will highlight the acceptability or unacceptability of any change. This exercise also gives the designer valuable understanding of the sensitivity of the design to small changes in pipe sizes.

2. Based on the theoretical D_ij from the CB approach, a split-pipe can be created for each link. That is, the lengths for the two discrete pipes sizes that bound the theoretical D_ij 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 match the CB solution in terms of cost.

7 Results

The current study evaluated three types of algorithms that optimize the design of gravity-driven water networks (GDWN). The algorithms include the calculus-based (CB) algorithm (Eq. 18), a backtracking algorithm (BT) and its modified version (BT-NoUp), and a genetic algorithm (GA). The algorithms were applied to six test cases that are based on real GDWNs. Our results show that the CB, GA, and BT-NoUp algorithms could find solutions to the GDWNs within 25 % of the BT global optimum. All cases assume minor-lossless flow and a two-part pipe-cost model. Solution costs from each algorithm are shown in Table 2 and runtime statistics are shown in Table 3. BT could run to completion in < 1 min in all but the largest case (case 6 with 59 links), which did not complete after 7 days. As such, cost comparisons to BT are not made for case 6.

Table 2 Solution costs for each algorithm.

^∗ Note: BT did not complete after 7 days of runtime.

Table 3 Runtime and size of solution space for each algorithm.

^a BT, BT-NoUp, and GA algorithm run times do not include approximately 2 s of pre-processing time. ^b BT did not complete case 6 after 7 days of runtime.

Figure 6 Diameter sizes from calculus-based (CB-Disc) solutions compared with global optimum solutions (from backtracking, BT). A global optimum for case 6, Los Mangos, is not included since BT did not complete after 7 days of runtime.

The CB algorithm based on Eq. (18), 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 BT discrete-diameter global optimum (5.5 to 2.6 % lower cost than BT). In fact, because of the discrete pipe sizes needed for an actual network, the continuous model will always produce the smallest theoretical network cost. 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 known global optimum solutions, as shown in Fig. 6. From all the combined test cases with known global optima, 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 be more likely to produce results identical to the global optimum (see Sect. 6). This was performed in the current study but the results are not presented because of space constraints. The CB-Disc solution costs were, in all cases, larger than the global optimum, with costs ranging from 3.9 to 22.6 % above the global optimum. 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 but case 6, the CB algorithm converged on a solution in 3 min or less.

Table 4 Case network properties, diameter (D) results (inch nominal sizes, with CB-Theor in inches), and nodal h (in meters) for (a) Case 1, Kiangan, (b) Case 2, Los Modulos, (c) Case 3, Cañazas, (d) Case 4, San Miguel, (e) Case 5, El Guabo, and (f) Los Mangos.

BT-NoUp, a modified version of BT which does not consider smaller diameters upstream of large diameters, completed itself within 4 s for all cases, and found solutions which matched or came very close to the BT global optimum. BT-NoUp missed the global optimum in cases 2 and 3, although by a small percentage increase in cost (2.1 and 0.4 % respectively). BT-NoUp, however, finished its search in a shorter amount of time in comparison to BT, a benefit that becomes relevant on problems with larger solution spaces, such as cases 3 (1.0 × 10²³ candidate solutions) and case 6 (7.3 × 10⁴⁹ candidate solutions).

GA was run on each case a total of 100 times, each run itself evolved 200 candidates for 500 generations. The lowest-cost candidate amongst the final population that did not violate the pressure head condition was chosen as the GA solution. 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. Overall, GA costs came close to the global optima (within 3 %) for cases 1–5 where the global optimum was known from BT. GA solution costs increased with larger network sizes, with its solution cost 18 % higher than CB-Theor for case 6, the largest case run. Each GA run finished consistently within 1–5 s, not including about 2 seconds of pre-processor time. We note that variations of GAs have been reported in the literature and several of these may 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). One reported improvement, the scaling of the fitness function to magnify the rewards towards slightly fitter candidates at later generations (Dandy et al., 1996), was attempted for case 2 but did not result in a noticeable effect on performance.

Table 5 Optimization results from Bhave (1978) algorithm. LHS sum and RHS sum are the left and right sides of his Eq. (19), which should be equal.

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 h_min, which was set to 7 m in the present work.

We compared the CB results for the Los Mangos network with those from the Bhave (1978) optimization algorithm (see Table 5). Like Eq. (15) in the present work, Bhave's optimality equation (his Eq. 19) equates the sum of a weighted term for all links entering and leaving each internal node in the distribution main. In the present work the term is proportional to the hydraulic gradient and the weighting factor is proportional to flow rate. In Bhave's case the term is the ratio of pipe cost to head loss, where the weighting factor is pipe-cost exponent b. There are 60 nodes for this network, including 30 nodes in the distribution main. The rest are delivery nodes (note there are 2 branches from node 30 of the distribution main). The terms required for the calculations include b, pipe cost, and head loss in the main and branches. The designation LHS refers to nodes in the distribution main entering, and RHS to those leaving, the node at the far-left side of Table 5. The exponent b comes from curve fitting pipe-cost data to the two-part pipe-cost model. Linear interpolation was used between diameters D_co and D_co+1 to obtain b in this range. Except for a few nodes, agreement between the two CB algorithms is very good. Although Bhave's Eq. (19) and Eq. (18) in the present work, appear quite different due to the different ways each was developed, both produce optimality for the networks considered in this paper. The key distinction between the two developments is the assumption of constancy in terms that comprise the coefficient A in Bhave (his Eq. 13), mainly the exponent b (m in his paper). In general, b depends on pipe diameter, thus making b=b(D) for multi-part cost models. When taking derivatives to obtain the final algorithms in both works, this dependence must be included, which produces additional terms in the optimization equation (see our Eq. 18 above). However, if the system of equations is solved by an iterative method, as Bhave proposed, the dependency may be neglected (though issues with convergence of the numerical solution may arise because of this). It is very important to note that if a non-iterative method is used to solve the system of equations as done in the present work (using a commercial program like Mathcad), all terms in the governing equations must be treated as continuous, not discrete, and the b=b(D) dependence must be explicitly included. It should also be noted that the optimization algorithm of Eq. (18) in this paper includes minor loss, which is not included in the Bhave (1978) work.

8 Conclusions

Algorithms to optimize the cost of branching gravity-driven water networks are evaluated on six 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, when run to completion, 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 rejecting all candidates that included a small diameter directly upstream of a larger diameter but allowed for the possibility of missing the global optimum. The third type of algorithm evaluated was a genetic algorithm (GA) that used single-point crossover and tournament selection.

BT could find the global optimum in most test cases with relatively little computational effort, although its poor scaling to larger networks is evidenced by its inability to find a solution to case 6, a network with 60 nodes and 59 links. The BT-NoUp completed its search in less time than BT and could find a solution to case 6. Based on case 1–5 results, the extra pruning condition adopted in BT-NoUp sacrificed only small cost increases. Both BT and BT-NoUp, however, could become prohibitively time-consuming when dealing with networks with significantly more links, diameter choices, or an unfavorable layout. 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 larger 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, like Kadu et al. (2008), by only including the two diameters above and below the CB-Theor diameters, producing four diameter choices per link. The calculus-based methodology provides an additional benefit to the designer by explicitly revealing the sensitivities to cost for a design. The calculus-based algorithm requires greater computational effort than backtracking for smaller networks, however, this effort scales more linearly with the number of network links, while backtracking scales exponentially. Furthermore, backtracking's computational time is sensitive to the number of available diameters. Still, when applied to the present study's GDWN test cases with a modest number of links (23), backtracking quickly found a global optimum. In addition, because it is guaranteed to find the global optimum, it can be useful for benchmarking the performance of other algorithms which scale better with more network links. While the genetic algorithm produced solutions with good proximity to the global optimum, its solution costs tended to be further from the global optimum in cases with more links.

For all test cases, the calculus-based algorithm's theoretical diameter solutions (CB-Theor) produced a lower cost than the discrete-domain global optimum. This result is made possible because it is not constrained to a discrete set of diameters. As such, the CB-Theor results represent a lower-bound on the optimum solution within the problem formulation, which could be approached with a finer selection of pipe diameters. We also demonstrated good agreement between the CB-based optimization algorithm developed here and that of Bhave (1978). Though Bhave's algorithm and Eq. (18) in the present work appear quite different due to the different ways each was developed, both produce optimality for the networks considered in this paper. The key distinction between the two developments is that Bhave assumed exponent b constant in the pipe-cost model, which was justified based on his iterative method of solution. In the present work, which uses a commercial program to solve the nonlinear governing equations for D and h, b(D) dependence is explicitly included for multi-part cost models. Contrasted with Bhave, minor losses are included in the CB optimization algorithm in the present work.