AbstractIn recent years, pressure-driven analyses of water distribution systems using pressure-flow relationships (PFRs) were shown to be more computationally efficient than those based on flow-pressure relationships (FPRs), such as the widely used Wagner equation and various spline approximations. Using a PFR enhances the convergence properties of the Newton–Raphson method used by most water distribution network (WDN) solvers. This paper derives a new way of incorporating a PFR into the classical Todini and Pilati global gradient algorithm (GGA). The convergence properties of the resulting solution algorithm are compared with those from two other existing PFR algorithms, EPANET 2.2 and the active-set method (ASM), as well as from a conventional flow-pressure–based algorithm on a number of different networks of varying size.IntroductionSolvers 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: 1.The FPR is generally a concave function of H and also its first-order continuous derivatives approximations (Fujiwara and Ganesharajah 1993; Tanyimboh and Templeman 2010; Piller et al. 2003; Elhay et al. 2016; Deuerlein et al. 2019) unavoidably incorporate a concave portion instead of being a fully convex function as required by the N–R approach to rapidly converge;2.The first-order derivatives are noncontinuous for some suggested FPRs in the technical literature (e.g., Wagner et al. 1988) for H=Hm and H=H*, where H is the actual hydraulic head, H* is the desirable or nominal head, Hm the elevation at and below which no flow can be delivered, while Z is the geodetic elevation of the node.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.GGA “Direct” Flow-Pressure d(H) ApproachThe 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) (1) [A11A12A210][QH]=[−A10H0d*]where: QT=[Q1,Q2,…,Qnp] unknown pipe flows (np is the number of WDN pipes); HT=[H1,H2,…,Hnn] unknown nodal heads (nn is the number of WDN unknown head nodes); H0T=[Hnn+1,Hnn+2,…,Hnt] known nodal heads (nt–nn is the number of WDN nodes with known head); and d*T=[d1*,d2*,…,dnn*] known nodal desired demands.The matrices A12=A21T and A10 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 A11 is a diagonal matrix, where the non-null elements are given by (2) 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 form (3) d(i)={0Hi≤Himdi*(Hi−Him)1/c(Hi*−Him)1/cHimdi*.The idea of using an inverted form of a PFR in the GGA solver originates from the way in which emitters were modeled as unconstrained pressure-dependent demands in version 2.0 of EPANET released back in 2000 (Rossman 2000). In 2007, Rossman also described how the concept could be extended to model-bounded PFRs as well (Rossman 2007).The EPANET 2.2 PDA algorithm can be derived as follows. Starting from the basic GGA system of equations (12) [A11A12A210][QH]=[−A10H0d*]whose definition is the same as per Eq. (1), the introduction of a set of new virtual links requires the addition of a vector d of nn unknowns representing the nodal demand flows. This leads to the revised system of equations (13) [A110A120A33−IA21−I0][QdH]=[−A10H0−Hm0]where (14) A33(i,i)={0di≤0(Hi*−Him)di*c|d|ic−10di*I is an identity matrix of size nn×nn, while the nn vector of nodal desired demands is now set equal to zero being replaced by d in the system of equations. Please note that not all the diagonal elements of matrix A33 will be non-null.For a given value for unknowns Q,d,H at a generic iteration, the system of equations will not necessarily be satisfied, namely (15) [A110A120A33−IA21−I0][QdH]+[A10H0Hm0]≠[000]Linearizing the problem by differentiating Eq. (15), the following system is obtained: (16) [D110A120D33−IA21−I0][dQdddH]=[df1df2df3]where D11 is defined as in Eq. (6) while the diagonal matrix D22 is defined as (17) D33(i,i)={0di≤0c(Hi*−Him)di*c|d|ic−10di*and with (18) [df1df2df3]=[A11Q+A12H+A10H0A33d−H+HmA21Q−d]The system to be solved at each iteration is (19) [dQdqdH]=[D110A120D33−IA21−I0]−1[df1df2df3]Which when solved by partitioning leads to (20) H(+)=(A′)−1F′Q(+)=Q−D11−1(A11Q+A12H(+)+A10H0)d(+)=d−D33−1(A33d−H(+)+Hm)with (21) A′=A21D11−1A12+D33−1(22) F′=F−D33−1[(D33−A33)d−Hm]and F being the classical GGA expression (23) F=−d+A21D11−1[(D11−A11)Q−A10H0]Please note that this formulation does not increase the size of the A′ matrix used to solve for new heads: it remains at nn×nn as in the original GGA. However, at this point, a direct solution of Eq. (20) is not possible because the diagonal matrix D33 may contain several zero values in its principal diagonal as per its definition given in Eq. (17), which makes it not invertible and therefore impossible to evaluate A′ or F′ using Eqs. (21) and (22). To solve this problem, EPANET 2.2 uses a barrier-function technique that adds large weights to the main diagonal of matrix A33 which is now defined as (24) A^33(i,i)={BIGdi≤0(Hi*−Zi)di*c|d|ic−10di*with BIG a large positive number, such as 108. This leads to (25) D^33(i,i)={BIGdi≤0c(Hi*−Him)di*c|d|ic−10di*with the result that D^33 will now be invertible.Having modified A33 into A^33 an additional correction must be introduced in the second of Eq. (18) to match the original equations, by modifying Hm into H^m defined as (26) H^m(i)={Himdi≤0Him0di*The new equations to be solved are finally (27) H(+)=(A^′)−1F^′Q(+)=Q−D11−1(A11Q+A12H(+)+A10H0)d(+)=d−D^33−1(A^33d−H(+)+H^m)with (28) A^′=A21D11−1A12+D^33−1and (29) F^′=F−D^33−1[(D^33−A^33)d−H^m]Empirical testing has shown that the PDA performance of EPANET 2.2 could be improved by reducing the magnitude of the change in di produced by Eq. (27) (with its sign preserved) to di*/2 whenever it exceeds di*. It is this slightly modified EPANET, referred to as EPANET 2.2.1, that was used in the examples subsequently presented in this paper.Active-Set MethodIn 2019, Deuerlein et al. presented an inverse FPR approach based on the minimization of the Collins et al. (1978) “content.” As opposed to the scalar form of EPANET 2.2, this approach introduces, in matrix form, three sets of new constraints describing the inverse FPR. The resulting system of equations, which includes the Lagrange multipliers λi relevant to the nodes for which di≤0 and μi relevant to the nodes for which di≥di*, is then simplified to lead to a set of equations, Eq. (14) in Deuerlein et al. (2019). The resolving algorithm, the ASM, is then described in the “iteration loop” section and in Eqs. (15)–(17) of the aforementioned work.In the ASM, the decision on which nodes to use the PFR is not only based on demand as in EPANET 2.2, but also on the Lagrange multipliers, namely either when 0Him) or, additionally, when di=di* and μi<0 (corresponding to Hidi* are set to d^i=di*.It is not difficult to show that the equations in Eq. (43) are exactly the same equations given as Eqs. (15)–(17) in Deuerlein et al. (2019) and that they lead to the same solution of Eq. (27) of the EPANET 2.2 algorithm. Thus, the basic differences in the three algorithms are essentially limited to the selection criteria used to set the variables at their lower, intermediate, or upper limit, which can be summarized as in Table 1.Table 1. Selection criteria for variables at the lower, intermediate, or upper limit used by the three approachesTable 1. Selection criteria for variables at the lower, intermediate, or upper limit used by the three approachesEPANETActive-set methodGGA-PDALower limitd^i≤0(di*=0)∨(di*>0)∧[(d^i<0)∨(d^i=0∧Hi≤Him)](di*=0)∨(di*>0)∧(d^i≤0∧Hi≤Him)Upper limitd^i>di*(di*>0)∧[(d^i>di*)∨(d^i=di*∧Hi≥Hi*)](di*>0)∧(d^iHi≥di*Hi*)Within limits00)∧{[(d^i>0)∧(d^iiHim)∨(d^i=di*∧Hi0)∧(d^i>0∨Hi>Him)∧(d^iHi

Source link

Leave a Reply

Your email address will not be published. Required fields are marked *