## Abstract

This paper presents a novel analysis of the accuracy of quadratic approximations for the Hazen–Williams (HW) head loss formula, which enables the control of constraint violations in optimisation problems for water supply networks. The two smooth polynomial approximations considered here minimise the absolute and relative errors, respectively, from the original non-smooth HW head loss function over a range of flows. Since quadratic approximations are used to formulate head loss constraints for different optimisation problems, we are interested in quantifying and controlling their absolute errors, which affect the degree of constraint violations of feasible candidate solutions. We derive new exact analytical formulae for the absolute errors as a function of the approximation domain, pipe roughness and relative error tolerance. We investigate the efficacy of the proposed quadratic approximations in mathematical optimisation problems for advanced pressure control in an operational water supply network. We propose a strategy on how to choose the approximation domain for each pipe such that the optimisation results are sufficiently close to the exact hydraulically feasible solution space. By using simulations with multiple parameters, the approximation errors are shown to be consistent with our analytical predictions.

- Hazen–Williams
- modelling
- optimisation
- pipe friction
- water supply networks

## NOTATION

- ,
- number of links and nodes, respectively
- number of different loading conditions
*L*[*m*],*D*[*m*],*S*[*m*^{2}]- 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]

## INTRODUCTION

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.

## QUADRATIC APPROXIMATIONS FOR FRICTION HEAD LOSSES

Throughout this paper, the friction head loss across a pipe is represented by the HW formula and it is given by
1with *n* = 1.852 and the positive real number *r*, which is the resistance coefficient of the pipe, is defined as
2where *L*, *D*, *C* are length, diameter and HW roughness coefficient of the pipe, respectively. For a given flow velocity *v* in a pipe, the corresponding flow *q* is given by , where *S* is the cross-sectional area of the pipe. In particular, if is the maximum expected flow velocity in a pipe, we have that is the maximum flow.

Formula (1) is non-smooth due to the rational exponent *n* = 1.852. In various optimisation problems for water supply networks, the friction head losses appear as constraints. In order to apply standard nonlinear programming techniques, it is necessary to use a sufficiently smooth approximation of Equation (1). In particular, given an expected maximum flow , we look for a quadratic polynomial function
3which approximates the original HW model over the range of flows – see as example Figure 1.

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.

## A QUADRATIC MODEL THAT MINIMISES ABSOLUTE ERRORS

In the present formulation we consider a smooth quadratic approximation of the friction losses across a pipe, generated by flow *q* ranging between 0 and some fixed maximum flow . We look for a polynomial function which minimises the integral of square errors defined by
4In the following, we refer to this approximation as . An analytical expression for *I*(*a*, *b*) can be derived as
5
where
6Therefore, a couple minimises the integral if it satisfies the following equations:
7The solution of the above linear system yields:
8

## ANALYSIS OF THE APPROXIMATION ERROR

We study the errors introduced by the considered smooth friction loss approximation formula. Given the quadratic polynomial function with coefficients defined by equations in (8), it holds
9where the function depends only on *n*. In particular, we have
10

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 the following, we further analyse the impact of on the approximation accuracy and we propose a selection strategy for its value. Recall that in our study *n* = 1.852. As shown in Figure 2(a), it holds
11In particular
12suggests that to improve accuracy we should avoid large values of .

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.

## A QUADRATIC MODEL THAT MINIMISES RELATIVE ERRORS

We consider the quadratic approximation that is obtained by minimising the integral of relative errors (Eck & Mevissen 2015). In this approximation scheme, a quadratic function is defined to minimise the integral: 13where and specify the approximation range. We refer to this approach as . Furthermore: 14with 15

We look for which solve the following system of equations: 16therefore: 17

## ANALYSIS OF THE APPROXIMATION ERROR

It can be shown from (17) that the accuracy of the quadratic approximation depends on , and the resistance coefficient *r*. In Eck & Mevissen (2015), a method for choosing was suggested, in order to control the relative error of the proposed approximation. In fact, given that is equal to the maximum considered flow value , we choose so that the minimum of the relative error function is within a given tolerance (see Figure 3).

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 .

Given and , let and be the coefficients defined in (17) where is chosen according to the method proposed in Eck & Mevissen (2015). In this case, we have:
18where function depends only on and *n*. In particular, it is possible to compute a real number such that:
19See Lemma 2 (Appendix 1) for a detailed derivation of the above expressions. In addition, there is an implicit nonlinear relation between function and the couple . From Figure 4, we can see that when , the choice of results in . Therefore, if we expect that most flows are significantly smaller than , we can set . The accuracy of the hydraulic solver for larger flows is improved when smaller values of are used. In the remaining part of the paper we consider the case of .

From Figure 5(b), the value is large for . Therefore, in order to have a small approximation error, should be defined so that the majority of expected feasible flows does not exceed this value. On the other hand, Equation (18) suggests that the use of a unnecessarily large value of can generate high errors. Nonetheless, if is such that for most feasible flows, then the corresponding values of can be very small and, consequently, the approximation is sufficiently accurate. As a result, we cannot conclude that a large value of necessarily results in big approximation errors for , which is contrary to the conclusion for (see also the numerical analysis section).

Now let and be the coefficients of the quadratic frictional loss formulae computed according to and , respectively. From the previous sections, we have:
20In Figure 6, we compare the graphs of and . From Figure 6(a), we conclude that when , results in a smaller absolute error than . In the case when the flow *q* is closer to , the best level of accuracy is achieved with . Note that both approximation schemes can result in large errors when .

Finally, recall that . Therefore, the absolute errors for and friction head loss approximations can be written as:
21Even though the quadratic approximation was observed to be more accurate for rough pipes in the case of DW friction models (Eck & Mevissen 2015), the above formulae demonstrate that this does not hold for HW models. In fact, when *L*, *D* and are fixed, both approximation schemes become less accurate for rough pipes. This is shown also in Figure 7, where a pipe with *L* = 100 m, *D* = 0.25 m and is considered.

## NUMERICAL RESULTS AND DISCUSSION

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.

As a case study, we consider the Smart Water Network Demonstrator operated by Bristol Water, InfraSense Labs at Imperial College London and Cla-Val presented in Wright *et al.* (2014). We refer to this case study model as BWFLnet. BWFLnet consists of 2,374 nodes, 2,434 pipes, two inlets (with fixed known hydraulic heads) and it is simulated under 96 different demand conditions for the extended period hydraulic simulation; its network topology and elevation map is presented in Figure 8.

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.

We define the mean absolute error on nodal pressure as
22where represents the vector of nodal pressures at loading condition *t*, computed by the optimisation process using a quadratic approximation for friction losses; on the other hand, represents the vector of nodal pressures computed by hydraulic simulation with optimised valve settings and HW friction loss model. Analogously, we define
23with and vectors of flows computed in valve optimisation and simulation, respectively. Finally, we formulate two empirical cumulative distribution functions as
24A preliminary investigation of the network model based on hydraulic simulations under different demand conditions has highlighted that the maximum velocities in all links are below 6 *m*/*s* – see Figure 9 for a scatter plot of the maximum velocities across all links.

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*.

Note that the value of (and ) used in computing the quadratic approximation does not need to be equal to the optimisation upper bound on flow (and hence velocity) variables. For example, it is possible to compute the coefficients of the quadratic friction head loss formulae with for all pipes and allow the optimisation to consider velocities up to 10 *m*/*s* – we refer to this scenario as T1. In this case, the feasible flow may be larger than for most and . Therefore, according to Figure 6, we expect to be more precise than . As reported in Figure 10(a), this is verified by our experiment, with being more accurate than .

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*.

The difference in accuracy described above is increased if we consider for all links ; this is scenario T3. In this case, the inaccuracies introduced by in the computation of friction losses, prevent convergence of the nonlinear program solver. On the contrary, in the case of , the mean difference between pressures computed during the optimisation process and those obtained from hydraulic simulation is close to 0.15 *m* (Figure 11(a) and 11(b)).

In the case of , the behaviour of the two quadratic approximation schemes can be described as follows. As shown in the previous section, the inaccuracies due to unnecessarily large on many links can result in high errors of approximations. Note from Figure 12(a) that many feasible flows are much smaller than the expected . Specifically, we observe that for most . This implies that for most ; see also Figure 12(b). In comparison, the values are at least an order bigger for the same links. With reference to Equation (20), we conclude that, in the case of , the value of is small enough to compensate the large on most links . This is not valid for .

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).

However, such inaccuracies do not affect the quality of feasible solutions for the considered application of minimising of AZP. The reason is that the AZP depends on nodal pressures only, which are computed with high accuracy. As shown in Figure 13, the value of AZP computed from the optimisation is close to the one obtained from hydraulic simulation with the original HW friction head loss formula for both and .

We observe that the optimisation with friction head loss model underestimates network pressures, while in the case of the simulated pressures are lower than the values computed from optimisation. For example, Figure 14 plots the pressure profiles at a critical point – i.e., the lowest pressure point within the zone. The minimum allowed pressure at this node has been set by the network operator to 18 *m*. From Figure 14(a), the simulated pressure based on optimisation with satisfies the minimum pressure constraints. On the contrary, the hydraulic solution corresponding to results in a constraint violation.

In the considered case study, most of the feasible flows are smaller than . With reference to Figure 15 and Equation (20), we conclude that overestimates friction losses on the majority of network flows, while underestimates these values. For this reason, pressures computed using are lower than corresponding hydraulically feasible pressures. In comparison, the optimisation with computes higher pressures than obtained from hydraulic simulation. Nonetheless, by appropriately choosing the ranges for the two quadratic approximations, a good level of accuracy is achieved. Consequently, smooth quadratic approximations for friction head loss models enable the application of standard nonlinear programming tools for the mathematical optimisation problems arising in the framework of optimal design and operation of water supply networks.

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.

## CONCLUSIONS

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.

## ACKNOWLEDGEMENTS

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.

- First received 19 July 2016.
- Accepted in revised form 2 January 2017.

- © 2017 The Authors

This is an Open Access article distributed under the terms of the Creative Commons Attribution Licence (CC BY 4.0), which permits copying, adaptation and redistribution, provided the original work is properly cited (http://creativecommons.org/licenses/by/4.0/).

Sign-up for alerts