Quadratic head loss approximations for optimisation problems in water supply networks

⁠,

number of links and nodes, respectively

number of different loading conditions

L[m], D[m], S[m²]

pipe length, diameter and cross-sectional area, respectively

C

Hazen–Williams roughness coefficient

q

pipe flow [] or [l/s]

v

fluid velocity [m/s]

maximum flow value [] or [l/s]

maximum velocity corresponding to [m/s]

⁠,

bounds of the approximation range [] or [l/s]

p

nodal pressure head [m]

⁠,

coefficients in quadratic approximation

⁠, ⁠, ⁠, ⁠,

constants used to compute

⁠, ⁠,⁠,

coefficients in

⁠,

coefficients in quadratic approximation

⁠, ⁠, ⁠, ⁠,

constants used to compute

⁠, ⁠, ⁠, ⁠,

coefficients in

accuracy threshold for pressures [m]

accuracy threshold for flows [] or [l/s]

The optimal management of water supply networks requires the satisfaction of multiple objectives, ranging from leakage reduction to improvements in water quality and system resilience. Consequently, various optimisation problems need to be solved including network design problems (Savic & Walters 1997), pump scheduling problems (Jung et al. 2015) and the optimal placement and operation of control valves (Nicolini & Zovatto 2009; Wright et al. 2015). Any optimisation problem, which is formulated for the hydraulic management of water supply networks, is constrained by the mass and energy conservation laws. These hydraulic expressions take into consideration friction head losses that can be represented by either the Hazen–Williams (HW) or Darcy–Weisbach (DW) formulae. The HW formula is semi-empirical (Liou 1998; Christensen et al. 2000) and it involves a non-smooth fractional exponential function, whose Hessian is unbounded around the origin. Consequently, the HW formula is difficult to handle as a constraint for nonlinear programming solvers, for which second order derivatives are often needed. In DW models, the relation between friction head loss and flow is defined by an implicit equation, which involves non-smooth terms, and it can be numerically calculated through an iterative process (Larock et al. 1999, Section 2.2.2). This complicates the use of such models in a smooth mathematical optimisation framework. As a consequence, optimisation problems for water supply networks are often addressed using heuristic approaches that handle the nonlinear nonsmooth hydraulic equations through highly specialised and customised simulation approaches (Maier et al. 2014). For example, genetic algorithms (Galdiero et al. 2016), ant colony optimisation (Mortazavi-Naeini et al. 2015) and simulated annealing (McCormick & Powell 2003) have been applied to optimal design and operation problems in water supply networks. Heuristic approaches, however, are computationally expensive and do not provide guarantees of optimality nor bounds on the global optimality of the solution. In comparison, this important information can be obtained when a mathematical programming method is applied to the same optimisation problems.

The solution of a wide range of problems for the optimal design and operation of water supply networks (Vairavamoorthy & Lumbers 1998; Bragalli et al. 2012; Gleixner et al. 2012; Menke et al. 2015) by using mathematical programming methods requires rigorously investigated and formulated smooth approximations of the friction head loss formulae. For example, Bragalli et al. (2012) approximated the HW head loss formula with a piecewise function, using a quintic polynomial approximation near zero. As noted in Eck & Mevissen (2015), such approach introduces computational complexities due to the high order polynomial function. Furthermore, in order to model these piecewise approximations, it is necessary to use binary variables. These result in mixed integer nonlinear constraints that are difficult to accomodate.

Various explicit approximations of the DW head loss formula have been published so far. A smooth and asymptotically consistent approximation is presented by Burgschweiger et al. (2009). Different analytical approximations for the DW friction factor are studied by Giustolisi et al. (2011). A smooth quadratic friction loss approximation for both HW and DW models was proposed by Eck & Mevissen (2015), together with a simulation-based analysis of the accuracy. The results in Eck & Mevissen (2015) showed that in practical applications, the use of polynomial quadratic friction loss formulae does not affect significantly the distribution of network pressures and flows. However, these approximations were not analytically discussed and accuracy bounds were not provided. In the case of DW models, the quadratic approximation in Eck & Mevissen (2015) is a least-squares fit on discrete head loss values over a range of flows. In comparison, we show that the explicit HW formula allows a more insightful analysis that better informs solution methods for the optimisation problems in water supply networks.

In this paper, we study the quadratic approximations of the commonly used HW head loss formula. We present a novel error analysis of the different quadratic approximations schemes, and investigate their impact on the solutions of optimisation problems for water supply networks. We propose and analyse an approximation that minimises the integral of square absolute errors, and we also formulate an alternative quadratic approximation to the one described by Eck & Mevissen (2015). We derive novel exact formulae to quantify the absolute errors for both approximations and we show how these errors are related to constraint violations in optimisation problems for water supply networks. Furthermore, we assess the efficacy of the smooth approximation framework for water supply network optimisation using an operational network as a case study. Finally, we propose a strategy for the choice of approximation parameters such that the corresponding optimisation results are sufficiently close to hydraulic feasibility.

Analogously to that observed in Abraham & Stoianov (2015, Appendix I), it is possible to show that is a continuous function with Lipschitz continuous first order derivative. Unlike formula (1), the second order derivative of the approximation is continuous almost everywhere and it is bounded at the discontinuous point q = 0. Note that the condition of Lipschitz continuous first order derivatives is sufficient to prove local convergence properties of numerical methods involving functions like (3); see Nocedal & Wright (2006, Theorem 11.2). As a consequence, and without loss of generality, we can restrict our analysis to positive flows only. In particular, we look for a quadratic polynomial function close to the HW friction head loss formula over the range of flows ⁠. Various mathematical notions of closeness can be used, each one resulting in approximation with particular characteristics. In this paper, we focus on two different approaches and analyse their goodness for optimisation problems in water supply networks.

The reader is referred to Lemma 1 in Appendix 1 (available with the online version of this paper) for a technical proof of the above statements.

In addition, from Figure 2(b) we can conclude that the function is monotone increasing for x > 1. Therefore, the greater the flow is than ⁠, the less accurate the quadratic friction loss approximation is. This is not surprising as the approximation is optimised to minimise errors in ⁠.

In conclusion, we consider an approximation interval large enough to include the majority of expected feasible flows. However, this does not mean that the value of should be unnecessarily large.

Since our aim is to apply the considered approximation scheme to formulate constraints of different optimisation problems for water supply networks, we study the effect of ⁠, and on the absolute errors. In fact, these variables affect the degree of constraints violation for a feasible candidate solution. Therefore, we derive new exact analytical formulae for the absolute error. With the application of these formulae, we provide new insights into the approximation defined in (17). In particular, we show that, as for ⁠, the absolute error is proportional to the resistance coefficient and it is a nonlinear function of ⁠.

Smooth head loss approximations are critical in the formulation and solution of mathematical optimisation problems for water supply networks. This is because frictional head loss formulae appear as nonlinear constraints in many optimisation problems (Burgschweiger et al. 2009; Bragalli et al. 2012; Eck & Mevissen 2012; Menke et al. 2015; Pecci et al. 2016). In this paper, we consider optimisation problems that require the computation of optimal control settings for pressure reducing valves (PRVs) in order to minimise average zone pressure (AZP) (and leakage) (Vairavamoorthy & Lumbers 1998; Wright et al. 2014). In Pecci et al. (2016), the co-design problem of optimal placement and operation of PRVs used to approximate frictional head loss formula within the optimisation constraint. We refer the reader to Appendix 2 (available with the online version of this paper) for a description of the problem formulation.

We model a water supply network as a graph with links and nodes. We define an extended period simulation and optimisation formulation, with different demand conditions. Once the locations of the valves are fixed, the same mathematical framework can be used to optimise valve operational settings, resulting in a nonlinear program. The unknown variables include pressure at demand nodes and flows in pipes while the objective to be minimised is AZP. The optimisation problem is solved using the interior point nonlinear programming solver IPOPT (Waechter & Biegler 2006). All computations were performed within MATLAB 2015a-64 bit for Windows 7, installed on a 2.50 GHz Intel^® Xeon(R) CPU E5-2640 0 with 18 cores.

BWFLnet is composed of two interconnected district metered areas (DMAs) and it is currently operated with dynamic network connectivity (Wright et al. 2014). Two kept shut boundary valves between the DMAs were replaced by two dynamic boundary valves that are autonomously closed at low demand periods (night hours) and opened for the remaining 24 hours of a diurnal operational cycle. Three PRVs are optimally controlled using daily derived flow modulation polynomials to minimise AZP. The network model and control options have been expanded from the model presented in Wright et al. (2014). The HW formula is used to model friction losses within the BWFLnet.

We investigate the effect of the quadratic approximations and on the quality of the solutions. When optimal settings for the PRVs are derived, these are implemented within the hydraulic equations to simulate network pressures and flows, while valves are optimally operated. The hydraulic equations are solved using algorithm 1 in Abraham & Stoianov (2015), which is sometimes called a loop method – see also Appendix 2. The optimisation solution is then compared with simulated pressures and flows.

Let be the vector of maximum expected velocities across all links. We can define the corresponding vector of maximum expected flows as with for all ⁠, where is the cross-sectional area of link j.

We have already observed in the previous section that a better level of accuracy is reached when the values in are larger than all possible feasible flows for the considered case study. However, unnecessarily large maximum flow values can cause significant inaccuracies, see Equation (20). In order to investigate this aspect, we set for all and run the optimisation, we refer to this scenario as T2. In the case of ⁠, the optimisation-simulation difference on nodal pressures reaches large values for most nodes as summarised in Figure 10(c) and 10(d). The expected optimised AZP value is 37 m while the actual simulated AZP with optimised valves is found to be 39 m. When is used, the accuracy is improved, with an average error of 0.13 m.

In order to improve the accuracy of both optimisation strategies, we should tailor the value of the maximum expected flow for each particular link j and avoid unnecessarily large ⁠. Recall that for each ⁠, we have ⁠, where is the cross-sectional area of link j. Therefore, a tailored choice of vector for each link results in a tailored for each link. Given ⁠, we set if the maximum simulated velocity across link j is lower than ⁠, for ⁠. In cases where the optimisation processes result in increased flow velocities, for example, in operational optimisation to improve ‘self-cleaning capacity’ of the network (Vreeburg et al. 2009; Abraham et al. 2016), we can set large enough so that all expected feasible flows are within the approximation interval.

The proposed strategy avoids overestimating values on small pipes with low velocities. We consider tailored maximum velocities with ⁠; this is scenario T4. In this case, we obtain the best pressure accuracy for the two approximations (see Figure 11(c)).

We further analyse the quality of the different solutions provided by the optimisation process with and friction head loss models, using a tailored maximum expected flow for each link. In this case, many feasible flows are such that ⁠. In the previous section, we have shown that this implies and have confirmed by the experimental results. In fact, as shown in Figure 11(c), the optimal solution corresponding to is more accurate than the one related to ⁠. As observed in Figure 11(d), both quadratic approximations cause large errors (more than 1 l/s) in the computation of a small fraction of network flows (less than 2% for the described case study).

In the present analysis, we study eight different optimisation problems; in fact, for each approximation scheme, four different scenarios are investigated. All the considered optimisation problems are large-scale nonlinear programs with 695,232 variables, 934,656 nonlinear constraints and 1,618,368 linear constraints. The computational performance of the solver IPOPT is reported in Table 1. We observe that in the case of T3 the solver failed to converge to a solution of the optimisation problem formulated using ⁠. Although the overall computational performance is satisfactory, improvements can be achieved by simplifying the problem formulation. The study of alternative and computationally efficient formulations for optimal valve operation problems in water supply networks is beyond the scope of this work and is a subject for future research.

Table 1

Computational performance of the solver IPOPT

.	.		.
	CPU .	Iter .	CPU .	Iter .
T1	8,332 s	702	1,196 s	138
T2	2,424 s	295	1,284 s	149
T3	–	–	1,930 s	201
T4	7,445 s	642	7,240 s	590

.	.		.
	CPU .	Iter .	CPU .	Iter .
T1	8,332 s	702	1,196 s	138
T2	2,424 s	295	1,284 s	149
T3	–	–	1,930 s	201
T4	7,445 s	642	7,240 s	590

For each quadratic approximation scheme, four optimisation problems were considered corresponding to different scenarios. In scenario T1, the quadratic approximations of the friction head loss formula is performed over an interval defined by ⁠, while the optimisation allows velocities up to 10 m/s. Scenarios T2, T3 and T4 use maximum velocities of 6 m/s, 10 m/s and tailored ⁠, respectively, both for approximations and optimisation frameworks.

View Large

Quadratic approximations have been effectively used to formulate head loss constraints in different optimisation problems for water supply networks. Therefore, quantifying and controlling their absolute errors, which affect the degree of constraint violations of feasible candidate solutions, has a strong impact on the application of mathematical optimisation. Both the HW and DW friction head loss models need smooth approximations to be posed as explicit constraints in mathematical optimisation. However, we show tight analytic error bounds for the HW case because of its explicit formula, unlike the implicit DW model.

In this paper, we have presented two quadratic approximations that minimise the absolute and relative errors, respectively, for the non-smooth HW friction head loss formula over a range of flows for each pipe. We have derived exact analytical formulae for absolute errors for two quadratic approximations that we investigated. In particular, we have shown that the absolute head loss approximation error for each pipe is proportional to the resistance coefficient and it is a nonlinear function of the approximation domain. Based on the derived explicit formulae, we have provided new insights into quadratic approximations of the HW head loss formula for solving optimisation problems for water supply systems. We have also discussed the critical nonlinear relations between errors and the range of flows. Moreover, in the case of the quadratic approximation that minimises relative errors, our analytical framework allows an efficient strategy for the computation of the quadratic approximation coefficients, which is especially well suited when considering large-scale water supply networks with many pipes.

Friction head loss formulae appear as nonlinear constraints in many optimisation problems for water supply networks. These problems include optimal network design (Bragalli et al. 2012), optimal placement and operation of valves (Pecci et al. 2016) and optimal pump scheduling (Gleixner et al. 2012; Bonvin et al. 2015; Menke et al. 2015). The presented quadratic approximation schemes facilitate the solution of such problems in a more deterministic and rigorous mathematical programming framework.

Furthermore, we have experimentally validated the application of these approximations to the mathematical optimisation problem for the pressure control of an operational network. This case study represents a new generation of intelligent water supply networks that dynamically adapt their connectivity and operational objectives. The constraint violations and approximation errors in numerical experiments are found to be consistent with our analytical predictions. Finally, we have proposed a strategy for the choice of the approximation domain such that the derived optimisation results are sufficiently close to hydraulic feasibility.

The authors would like to acknowledge the support of NEC (the NEC-Imperial Smart Water Systems project) and EPSRC (EP/P004229/1). Dr Abraham was a Post-Doctoral Research Associate at Imperial College London (InfraSense Labs) when the work presented in the manuscript was carried out.