Pressure Flow–Based Algorithms for Pressure-Driven Analysis of Water Distribution Networks

Solvers for water distribution network (WDN) models, from the original Cross (1936) method to the widely used Todini and Pilati global gradient algorithm (GGA) (Todini and Pilati 1988), were essentially developed for design purposes and run by fixing demands while modifying pipe diameters until the desired pressures were reached at all network nodes. This implied using demand-driven approaches (DDAs) regardless of the resulting nodal pressures. With the advent of reliable codes, such as for instance EPANET (Rossman 2000) and several other research and commercial packages, the interest shifted from mere design to resilience and reliability-based design as well as to extended period simulations (EPS). As already mentioned, the design of WDNs was historically approached using the DDA where pipe diameters were to be found in order to meet desired pressures at nodes for a given topology and demands pattern. If pressures at all nodes were not met, diameters were increased until the objective was reached. In recent approaches, design also takes into account resilience and reliability criteria to allow overcoming critical situations and failures; moreover, leakage control can be improved by simulating in time using EPS leak losses that are usually represented as a function of the pipe average pressure (see for instance Giustolisi et al. 2008); finally, the increasingly extended availability of flow and pressure sensors, including smart meters, used in conjunction with data assimilation techniques (Hutton et al. 2014a, b; Bragalli et al. 2016) allows for offline as well as online EPS using WDN models aimed at increasing management capabilities while reducing risks of failures. Specifically, in EPS one does not change the pipe diameters, which means that pressure during the simulation in time, under certain conditions (Wu and Walski 2006) such as pipe bursts (or isolation) or during excessive water use, may fall below the desired pressure at some of the nodes, with the consequence that demand can’t be fully delivered. In these cases, DDAs fail to realistically represent the actual WDN behavior and, when pressure falls below the desired value, one needs to insert a pressure-based condition at the relevant nodes. This alternative pressure-driven approach (PDA) requires defining a flow-pressure relationship (FPR) to describe the reduced nodal flow delivered as a function of available pressure and demand.

Several analytical formulations of FPR were proposed in the literature (Bhave 1981; Germanopoulos 1985; Wagner et al. 1988; Reddy and Elango 1989, 1991; Chandapillai 1991; Fujiwara and Ganesharajah 1993; Tucciarelli et al. 1999; Tanyimboh et al. 2001; Wu et al. 2009; Tanyimboh and Templeman 2010; Jun and Guoping 2013; Morley and Tricarico 2014) while many other authors have introduced artificial components such as emitters, pipes, reservoirs, and so on, to represent the FPR behavior (Ang and Jowitt 2006; Baek et al. 2010; Suribabu and Neelakantan 2011; Jinesh Babu and Mohan 2012; Gorev and Kodzhespirova 2013; Sivakumar and Prasad 2014, 2015; Abdy Sayyed and Gupta 2013; Abdy Sayyed et al. 2014, 2015; Suribabu 2015; Mamizadeh and Shaoonizadeh 2016; Suribabu et al. 2017; Mahmoud et al. 2017). The interested reader can find a clear synthetic description of all the different approaches available in the literature in Suribabu et al. (2019).

Unfortunately, most of the FPRs were written in the form of demand $d$ as a function of pressure head $H$, namely $d=d(H)$, which, as discussed by Ackley et al. (2001) or by Todini (2006) and also recently demonstrated by Deuerlein et al. (2019), tends to increase the number of iterations required by a gradient approach, such as the Newton–Raphson (N–R), to reach convergence when solving the WDN problem.

The lack of convergence is mostly due to two reasons:

To overcome these problems, several authors introduced underrelaxation coefficients or defined alternative first derivative continuous FPR, as subsequently discussed.

In 2020, the US Environmental Protection Agency (USEPA) released version 2.2 of EPANET (USEPA 2020) which uses an inverse formulation of the FPR by adding a virtual one-way link between a node and a virtual reservoir (fixed head node), to represent an inverse flow-pressure–type relation (FPR) in the form of $H=H(d)$, which will be referred to as a pressure-flow relation (PFR).

This approach, already advocated and tested by Rossman several years before for describing emitters, resulted into a noticeable improvement in the convergence properties without increasing the size or the density of the system matrix. More recently, Deuerlein et al. (2019) found improved convergence performance using the inverse relations in their active-set method (ASM), the equations of which are derived by minimizing the Collins et al. (1978) “content” reaching the same system of equations used in EPANET 2.2, although solved differently, as subsequently discussed.

The aim of this work is to introduce a third derivation of the inverse FPR-based PDA algorithm in line with the original GGA derivation to point out the minor changes needed to convert the original demand-driven GGA to implement the PDA. The resulting equations are the same as the previous two derivations, but the final solution algorithm differs slightly from both.

The PDA problem is commonly solved in its “direct” FPR form to impose a reduction of the desired demand in the nodes where pressure is not sufficient to allow delivering it. From the classical GGA equations relevant to the DDA problem (Todini and Pilati 1988; Todini and Rossman 2013)where: ${\mathbf{Q}}^{\mathrm{T}}=[Q_1,Q_2,\dots ,Q_{n_p}]$ unknown pipe flows ($n_p$ is the number of WDN pipes); ${\mathbf{H}}^{\mathrm{T}}=[H_1,H_2,\dots ,H_{n_n}]$ unknown nodal heads ($n_n$ is the number of WDN unknown head nodes); ${\mathbf{H}}_0^{\mathrm{T}}=[H_{n_n+1},H_{n_n+2},\dots ,H_{n_t}]$ known nodal heads ($n_t–n_n$ is the number of WDN nodes with known head); and ${\mathbf{d}}^{*\mathrm{T}}=[d_1^{*},d_2^{*},\dots ,d_{n_n}^{*}]$ known nodal desired demands.

The matrices ${\mathbf{A}}_{12}={\mathbf{A}}_{21}^{\mathrm{T}}$ and ${\mathbf{A}}_{10}$ are incidence matrices defined as in Todini and Pilati (1988) with the flow considered positive if it is directed toward the relevant node and negative if it leaves the node, while ${\mathbf{A}}_{11}$ is a diagonal matrix, where the non-null elements are given by

As proposed by Salgado-Castro (1988), Salgado-Castro et al. (1994), and subsequently by Todini (2006), one can derive the PDA equations by modifying the right-hand side vector to include the FPR equation. Without loss of generality, for instance, in this work we will use the FPR equation of Wagner et al. (1988) in the formwith $c$ an exponent usually taken as $c=2$, to represent the unknown flow actually released at node $i$ instead of the required $d_i^{*}$.

Eq. (3) is graphically represented in Fig. 1, where $Z$ is the geodetic elevation of the node, $H_i^m$ the height at and below which no flow can be delivered, and $H^{*}$ is the required service height at and above which full demand can be satisfied.

For a generic hypothetical solution $\mathbf{Q}$ and $\mathbf{H}$, the system of Eq. (1) will not be satisfied with the equality sign. Therefore, we can writetaking into account that $\mathbf{d}$, given by Eq. (3), may be in part or entirely different from the desired demand ${\mathbf{d}}^{*}$ being either null or a function of the hydraulic head $\mathbf{H}$. Differentiating the system of Eq. (4), the result can be rearranged aswhereandtogether with

Taking into account the sign of the N–R algorithm, we will define all throughout this paper $\mathbf{d}\mathbf{Q}=\mathbf{Q}-{\mathbf{Q}}^{(+)}$ and $\mathbf{d}\mathbf{H}=\mathbf{H}-{\mathbf{H}}^{(+)}$, leading to the following recursive solution provided by the GGA:with ${\mathbf{H}}^{(+)}$ and ${\mathbf{Q}}^{(+)}$ the updated values of nodal hydraulic heads and pipe flows, and

The iterative process of Eq. (9) will proceed starting from an initial guess solution ${\mathbf{H}}^1$, and ${\mathbf{Q}}^1$.

The estimates of the actual delivered demands ${\mathbf{d}}^{(+)}$ are then updated as a function of ${\mathbf{H}}^{(+)}$ by means of Eq. (3). This solution only differs from the classical DDA solution for a different definition of $\mathbf{d}\ne {\mathbf{d}}^{*}$ and the presence of the non-null diagonal matrix ${\mathbf{D}}_{22}$ which is subtracted from the main diagonal of the system and for a correction factor $-{\mathbf{D}}_{22}\mathbf{H}$ applied to the right-hand side of the DDA problem.

Frequently, when using the direct FPR the convergence of the N–R process may be quite slow due to the concavity of the FPR which, similar to what happens with the N–R nodal approach (Shamir and Howard 1968), may lead to a rather slow alternating convergence (Todini and Pilati 1988; Todini and Rossman 2013).

To overcome this lack of convergence, various alternatives have been used, such as introducing underrelaxation coefficients (e.g., Giustolisi et al. 2008; Siew and Tanyimboh 2012; Elhay et al. 2016; Ciaponi and Creaco 2018). Other approaches, which assumed that the slow convergence was due to the lack of continuity of the first derivatives, defined alternative FPRs (Fujiwara and Ganesharajah 1993; Tanyimboh and Templeman 2010) and cubic functions or spline regularization (Piller et al. 2003; Elhay et al. 2016; Deuerlein et al. 2019) to guarantee the continuity of the first derivatives. Nonetheless, after a number of tests using these regularized functions, which still require underrelaxation to improve convergence, we reached the conclusion that more than the continuity of the first derivative, it is the presence of concavity that affects convergence and, in fact, most of the alternative FPRs introduced to guarantee continuous first-order derivatives are not strictly convex: indeed, they present an inflection point and the resulting functions are convex at one end but concave at the other end.

The proposed approach was tested and compared to ASM and EPANET 2.2.1 on three extremely simple examples, already used as testbeds for other approaches and described in the Appendix, as well as on two larger networks, that can be found in the literature, to test the validity of the approach on more complex networks.

The first network example is a simple no-loop 5-node and 4-pipe example introduced by Gupta and Bhave (1996) and already used by Cheung et al. (2005), Ang and Jowitt (2006), Morley and Tricarico (2014), and others, to analyze the behavior of pressure-driven approaches. The second network example is the two-loop network, with 5 nodes and 6 pipes, recently used by Deuerlein et al. (2019) to test their ASM approach. The third example is based on the classical Hanoi network, introduced by Fujiwara and Khang (1990). A detailed description of pipe length and diameters as well as of nodal elevation, desired demand and hydraulic head for the first three examples can be found in the Appendix.

The two additional tested networks were the Modena water distribution network (Artina et al. 2011) and the Balerma Irrigation Network (Reca and Martinez 2006). The description of both networks can be widely found in the literature but also on the University of Exeter Center for Water Systems web-site (http://emps.exeter.ac.uk/engineering/research/cws/resources/benchmarks/design-resiliance-pareto-fronts/large-problems/). The Modena WDN, originally presented by Artina et al. (2011) includes 4 reservoirs, 268 nodes, and 317 pipes. In order to test the PDD approaches, the desired demand was multiplied by 2 in order to generate unmet demand nodes. The Balerma Irrigation Network example originally presented by Reca and Martinez (2006) includes 4 reservoirs, 443 nodes, and 454 pipes; for testing the PDD approaches, the Hazen–Williams equations with $\mathrm{C}=130$ were used instead of the Darcy–Weisbach equations with $\epsilon =0.0025$.

This work proposes an approach using an inverse FPR in conjunction with the GGA, which has the advantage of leaving in demand-driven mode all the nodes for which the power is equal or higher than the power required to deliver at the node the desired demand at the required head, as well as all the nodes for which the estimated operating pressure and demand result at the same time smaller than or equal to zero, thus failing to deliver water.

The derived solution equations, based on an extension of the DDA GGA, are practically identical to the ones used in EPANET 2.2, which were derived in analogy with the use of a virtual pipe and a virtual reservoir connected to each node. They are also identical to the equations used in the ASM approach, which were derived starting from the constrained minimization of the “content” function. The three solutions only differ in terms of the criteria used to decide, at each N–R iteration, which are the nodes to be considered at the lower or upper limit and which ones fall in between.

The paper also shows the modest modifications that need to be applied to a DDA solver that only affect the main diagonal and the right-hand side of the linear system to be solved at each iteration.

Although demonstrated for the Wagner et al. (1988) equation, the algorithm performs well for any inverse strictly concave PFR, representing hydraulic losses, with the curvature typically ranging from linear to quadratic. On the contrary, inverse PFRs derived from first-order derivative continuous FPR functions [such as the ones proposed by Fujiwara and Ganesharajah (1993), Tanyimboh and Templeman (2010), or the cubic or splines functions] do not bring substantial improvements because they remain both concave and convex also after the inversion.

As demonstrated using a number of test examples, the proposed approach practically shows the same results and convergence properties of ASM and EPANET 2.2.1 over a variety of WDN and a range of desired pressures.