https://orcid.org/0000-0003-0847-1379
Abstract
This work compares several alternative methods of pressure-driven analysis (PDA) implemented within the code base of the EPANET water distribution system modeling software. The resulting code includes a direct method, the original and a newly modified version of the inverse method currently used by EPANET 2.2, and two recently published methods (the global gradient approach-based PDA and the active set method). The latter four methods solve PDA problems using an inverse flow–pressure relationship, whereas the direct method uses the original relationship. The alternative methods were extensively tested, and their performance was compared over several case study water distribution networks of vastly different sizes and complexities. The results showed that all of the new inverse methods are equally efficient and reliable, whereas the direct method is less reliable by having a higher frequency of failing to converge.
HIGHLIGHTS
A demand-driven approach and five alternative PDAs were incorporated into EPANET.
The EPANET 2.2 barrier function algorithm was modified to improve convergence, a direct PDA and the inverse global gradient approach-based PDA and ASM approaches were also included.
The code, thoroughly tested on the available water distribution network, cases up to over 150,000 nodes, will be made available for further comparisons.
Results confirm the lack of convergence of the direct approach.
INTRODUCTION
Criteria for evaluating the design and operation of water distribution systems typically include cost, resilience and reliability evaluated using extended period simulations (EPSs) of system performance (e.g. Gargano & Pianese 2000; Todini 2006). Moreover, EPSs have become a common practice for real-time control of water distribution network (WDN) leakage control by simulating real-time leak losses (e.g. Almandoz et al. 2005; Giustolisi et al. 2008); finally, using EPS, the increasingly extended availability of flow and pressure sensors, including smart meters, allows for off-line as well as on-line data assimilation (Hutton et al. 2014a, 2014b; Bragalli et al. 2016) and real-time control to increase management capabilities while reducing risks of failures.
In EPS, under certain conditions (Wu & Walski 2006) such as pipe bursts (or isolation) or during excessive water demand and use, pressure at some of the nodes of a WDN may fall below the desired level, with the consequence that full demand cannot be delivered. In these cases, demand-driven approaches (DDAs) not only fail to realistically represent the actual WDN performance but also do not allow estimating such deficiencies (Laucelli et al. 2012). For this reason, 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 mathematically represent the reduced nodal flow delivered as a function of available pressure and demand.
EPANET (US EPA 2020) is a widely used, open-source software package used to model the hydraulic and water quality behavior of water distribution systems. In 2020, pressure-dependent demands were introduced in EPANET 2.2 (Rossman et al. 2020) by using a barrier function approach applied to the inverted form of the FPR, namely a pressure–flow relationship (PFR) that had already proven successful in modeling emitters. The idea 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) and the description on how the concept could be extended to model bounded PFRs as well (Rossman 2007).
It is well known that the direct use of the FPR often requires smoothing the curve and/or dampening of the N–R estimated changes. Elhay et al. (2016) developed in fact the energy and mass residuals (EMRs)-damped Newton method, which was verified guaranteeing convergence when using the direct FPR.
Creaco et al. (2022) developed a third-order Newton–Raphson (N–R) approach combined with a variable coefficient under-relaxation technique, concluding that when using the direct FPR relationship, the implementation of a rigorous dampening and a more sophisticated convergence scheme always guarantee convergence.
However, in the meantime, the same authors who developed the EMR introduced an approach based on the inverse FPR, the active set method (ASM) (Deuerlein et al. 2019) based on the minimization of Collins’ ‘content’, pointing out that the new approach converges much faster than their previous EMR without requiring under-relaxation of the N–R corrections. Same conclusions were drawn for the global gradient algorithm-based PDA (GGA-PDA) (Todini et al. 2021) which was derived as a direct extension of the GGA (Todini & Pilati 1988), showing that the convergence of the N–R approach is better preserved by an inverse form of the FPR, which will be here indicated as PFR. In general, a PFR is a continuous convex function, instead of a concave function as the FPR, although preserving two discontinuities of the first-order derivative at its edges. In the same work, a discussion was presented on the limited benefits for using an alternative first-order derivative continuous approximations aimed at guaranteeing the continuity of the first derivatives, such as the ones proposed by Fujiwara & Ganesharajah (1993) and Tanyimboh & Templeman (2010) or the cubic functions or spline regularization by Piller et al. (2003), Elhay et al. (2016) and Deuerlein et al. (2019).
The inverse form of the Wagner et al. (1988) function, a continuous convex function over the entire domain of existence with the two mentioned discontinuities of the first-order derivative at its edges, was then used to test the three algorithms: the original one used in EPANET 2.2 based on linear barrier functions applied at the edges of the PFR (Rossman et al. 2020), the ASM approach (Deuerlein et al. 2019) derived by minimizing Collins et al. (1978) ‘content’ and the GGA-PDA (Todini et al. 2021) derived by expanding the GGA with a set of constraints describing the PFRs for nodes with pressure laying between a minimum and a desired one. The three PFR-based PDA approaches, the one used in EPANET 2.2, the ASM and the GGA-PDA were implemented in APL (Iverson 1962), a matrix-oriented computer language extremely useful for testing algorithms such as the solvers of WDN problems, but hardly usable on the very large WDN encountered in practice.
This motivated the need for evaluating the behavior of the new algorithms in an already well-developed package such as EPANET 2.2, the main objective being to show that under the same conditions of no under-relaxation, the inverse approaches better agree with the convergence properties of the N–R approach. This does not mean that using higher order N–R schemes and specifically designed under-relaxation approaches, such as the ones proposed by Elhay et al. (2016) or by Creaco et al. (2022), the convergence properties of the direct FPR would not improve. What we want to stress is that the inverse approaches do not require higher order or damped N–R approaches combined to sophisticated under-relaxation schemes to converge.
Moreover, there is no need for inverting smooth direct FPR expressions such as the one by Fujiwara & Li (1998) and Ciaponi et al. (2015) to solve PDA problems involving multi-storey buildings. The three proposed approaches can easily extend the use of the inverse FPR to solve multi-storey buildings’ PDA, by adding to the node of interest several virtual pipes each of which with the appropriate minimum and desired head requirements.
This article also describes the changes from the original code that were deemed necessary, together with the minimal changes needed in the input file, to implement alternative PDA methods into the new EPANET_WB package. Moreover, several test cases, ranging from few to over 150,000 nodes, were used to compare the performance of the alternative PDA approaches.
ALTERNATIVE PDA METHODS
pipes connected at their endpoints to a set of 
nodes at unknown head and 
nodes at fixed head, where the head loss across the kth pipe as a function of its flow rate 
is 
for known constants 
and 
. The conservation of energy and flow that govern this system can be expressed as:
is an 
vector of unknown pipe flows, 
is an [
vector of unknown nodal heads, 
is an 
vector of known nodal heads and 
is an 
vector of known nodal demands. The incidence matrices 
and 
indicate which pipes are connected to which unknown head nodes and to which fixed head nodes, respectively, while 
is a diagonal matrix with entries 
for each pipe k.
is found by solving:
is a diagonal matrix containing the gradient of each pipe's head loss equation with respect to flow: 
. After solving for new heads 
they are used in the following equation to obtain a new set of pipe flows 
:
are no longer fixed to their desired values 
but are instead functions of pressure. As with EPANET 2.2, this paper uses the Wagner et al. (1988) bounded power law expression to represent pressure-dependent demands; however, any other invertible function could be used as well. Assuming there is some minimum pressure 
below which no demand can be delivered and a higher pressure 
above which full demand is satisfied, then the demand realized at some node i is
is the full demand expected, 
is the available head, 
is the head below which no demand can be delivered, 
is the required head above which full demand is satisfied and Zi is the elevation of the node. The exponent 1/c is typically taken as 0.5 to mimic flow through an orifice opening.Direct PDA method
diagonal matrix 
must then be defined, where the ith diagonal element is given by:
It is worthwhile noticing that Equation (9) will remain formally the same for the other PDA solvers while the actual values for matrix 
and vector 
will differ from Equations (10) and (11), as will be shown in the following sections. After solving a N–R iteration for updated head values 
they are used in Equation (5) to obtain a new set of pipe flows 
.
In the direct approach, there is no need for an updating expression for the demands based on the N–R algorithm in analogy to Equation (5). The updated estimates of the actual delivered demands 
are directly obtained via Equation (6) as a function of 
. These iterations continue until some stopping criteria are satisfied.
Constrained emitter method
The PDA method contained in the current release of EPANET 2.2 uses the inverted form of the FPR demand function to treat pressure-dependent demands in an analogous fashion to how it models flow through emitter elements. While emitter flows are unconstrained, demand flows are forced to be between 0 and full expected demand. As with emitters, the demand flow out of some node i at head 
is considered as flow through a virtual pipe connected to a virtual reservoir node with fixed grade 
(see Figure 1).
Pressure-dependent nodal demand schematization using a virtual pipe and a reservoir (from Todini et al. 2021).
Pressure-dependent nodal demand schematization using a virtual pipe and a reservoir (from Todini et al. 2021).
, is computed using the inverse of the Wagner FPR:
.
This scheme, referred to as constrained emitter method (CEM)-1, requires an initial estimate of demands but not of heads and does not affect the size or structure of the solution matrix 
. To avoid matrix inversion problems when the estimated demand value approaches zero the diagonal elements of 
are not allowed to drop below some small value such as 10−7.
Modified CEM-2
, 
, 
, 
and BIG is a large positive number, such as 109 while 
is a small positive number, such as 10−6. This function is smooth with continuous first derivatives that are always non-zero. It is plotted in Figure 2 for the case where 
, 
and 
. Over the range of 0 to full demand, the function faithfully follows the original inverted demand function. Outside this range, the barrier terms quickly increase the head loss across the virtual pipe forcing demands to stay in the range of 0 to 
. For example, although not visible on the figure, for relative demands 
of 0.999, 1.0, 1.00001 and 1.0001 the head loss relative to 
increases smoothly from 0.998 to 1.0005, 10,000 and 100,000, respectively. The same GGA updating formulas apply where now the diagonal elements of 
are the derivatives of Equation (18) with respect to demand di of the ith node. Owing to the introduction of the small positive number 
, there is no longer need to adjust these derivatives to avoid matrix inversion problems. In addition, empirical testing has shown that improved convergence can be obtained when the change in demand at each N–R iteration is limited to be no greater than 40% of full demand.
Dimensionless plot of the inverted FPR used with CEM-2.
The ASM and GGA-based PDA methods
Two additional approaches to solve the PDA problem using the inverse PFR concept were recently introduced in the literature, the ASM due to Deuerlein et al. (2019) and the GGA-based PDA described in Todini et al. (2021). The two approaches are similar, they both solve an identical set of equations, only differing in the controls used to decide which are the nodes to be set at the upper or lower limit and those requiring being treated using the FPR curve.
to the relevant hydraulic head:
,
and 
are the same as in Equation (1); 
is a new diagonal 
matrix, with 
the number of nodes for which the actual demand is unknown, defined as:
; 
and where j indicates the position of the jth node in the vector of 
and 
indicates the corresponding position in the vector of unknown pressure nodes. As opposed to the EPANET 2.2 approach to the PDA, in both ASM and GGA, the size of matrix 
is not necessarily 
but varies with 
which may dynamically differ from one iteration to the following one.Matrix 
is an incidence topological 
matrix filled with zeroes except having a 
at nodes where actual demand is unknown. If 
, then 
, with 
being the identity matrix.
, defined as 
, contains the demands of vector 
for the nodes set at their upper or lower limits, whereas all the other elements of vector 
corresponding to nodes with unknown demands are set to zero, allowing one to define the actual estimated demand vector 
as:
, Equation (22) gives 
.The 
nodes, whose demand 
must be estimated using the inverse FPR in the expanded system of equations, are selected after each N–R iteration as those whose demands are neither at the lower nor at the upper limit.
The selection criteria used by the ASM directly descend from the minimization of Collins et al. (1978) ‘content’ function, whereas the GGA uses the concept of power (energy per time unit) available at a node, to decide where demand must be set at the lower or upper limiting value.
In the ASM, nodes are always set at the lower bound if the desired demand results 
. Alternatively, when 
, nodes are set at the lower bound if the estimated demand is 
or if 
and 
; nodes are instead set at the upper bound if the estimated demand 
or if 
and 
.
Similarly, in GGA-PDA, nodes are always set at the lower bound if the desired demand results 
. Alternatively, when 
, nodes are set at the lower bound if the estimated demand is 
and 
; nodes are instead set at the upper bound if the available power is larger than the desired one, namely 
.
Table 1 provides a summary of the selection criteria for the two approaches. One can clearly see the expressions derived from the Karush–Kuhn–Tucker (Karush 1939; Kuhn & Tucker 1951) conditions in the ASM as opposed to the selection criteria based on 
(a quantity proportional to the power) used by the GGA-PDA.
Selection criteria for variables at the lower, intermediate, or upper limit used in the ASM and the GGA-based PDA approaches
| ASM | |
| Upper limit |  |
| Within limits |  |
| Lower limit |  |
| GGA | |
| Upper limit |  |
| Within limits |  |
| Lower limit |  |
| ASM | |
| Upper limit | |
| Within limits | |
| Lower limit | |
| GGA | |
| Upper limit | |
| Within limits | |
| Lower limit | |
Note that it has been extended to the case 
from its previous version 
.
within the GGA are:
is the 
diagonal matrix of head loss derivatives with respect to demand 
. In addition, according to Equation (20), at each N–R iteration the vector of demands appearing in 
is modified from its original expression 
of Equation (4) to become 
, which contains 
for nodes whose demand has been set to the upper or lower bounds, and to zero for all the nodes in the active set to be updated. The updating formula for new demands δ belonging to the active set of nodes then becomes:
is added to 
in the estimation of 
, to avoid matrix inversion problems when the estimated demand value approaches zero, leading to:
In addition, similar to what was done for CEM-2, in the GGA and ASM approaches the change in demand at each N–R iteration is limited to be no greater than 40% of full demand.
Although not tested in this work, it is worthwhile noting that all the PDA approaches here presented can also be used to impose non-negativity to emitters in a similar manner by dropping the upper bound used for demands and imposing on the emitter flows only the lower non-negativity constraint.
IMPLEMENTATION OF METHODS
The direct, CEM-2, ASM, and GGA-PDA methods were each added to the current EPANET 2.2 engine source code (CEM-1 is already implemented in version 2.2). The choice of which method to use is made through the value assigned to the ‘Demand Model’ option in EPANET's input file: DDA (for demand-driven analysis), DIR (for the direct PDA method), PDA (for CEM-1), CEM-2, GGA or ASM. Initial values assigned to the program's unknowns were chosen to conform to the EPANET 2.2 standard setting with initial pipe flows set to produce a velocity of 
and initial demands equal to the full expected demand. For the direct, ASM, and GGA methods, heads are additionally initialized to the upper pressure limit at which full demand is realized (
).
TEST CASES
Fourteen different WDN models were used to test the different PDA methods. They cover a wide spectrum of well-documented models from classical examples with a small number of nodes and pipes to several real or synthetic cases with over 10,000 including one with over 150,000 nodes and pipes. Table 2 lists the size of each network as well as the URL where full information about each can be found.
Case studies with network sizes and URL
Each network was run under four different pressure ranges (Preq–Pmin) in meters of water column (or the identical quantities expressed in pounds per square inch). They are summarized in Table 3. The first three pressure ranges of 10, 20 and 30 m (14.22, 28.44 and 42.67 psi, respectively) correspond to typical distribution system operating conditions, while the fourth range of 0.1 m (0.142 psi) was included to stress test the algorithms. For several networks it was necessary to increase the full expected demands by a uniform multiplier in order to create pressure-deficient conditions. For each run, the number of trials needed for convergence was recorded as was the number of nodes that could not supply full demand and the percentage shortfall in supplied demand from them.
Testing pressure conditions
| Minimum pressure (m) | Required pressure (m) | Pressure range (m) |
|---|---|---|
| 5 | 15 | 10 |
| 5 | 25 | 20 |
| 5 | 35 | 30 |
| 0 | 0.1 | 0.1 |
| Minimum pressure (m) | Required pressure (m) | Pressure range (m) |
|---|---|---|
| 5 | 15 | 10 |
| 5 | 25 | 20 |
| 5 | 35 | 30 |
| 0 | 0.1 | 0.1 |
is the result of evaluating its FPR form in Equation (6) at the new head value 
.TEST RESULTS
The results of the test runs for the three normal pressure ranges are summarized in Table 4. The first column lists the name of the test network and the multiplier applied to its base demands. Column 2 lists the pressure ranges used in each test. Column 3 contains the number of nodes where pressure was insufficient to deliver the full expected demand and column 4 lists the percentage of reduction in demand at the affected nodes. Column 5 shows the number of GGA trials needed to converge for a DDA, where demands are fixed at their full expected levels. The remaining columns show the number of GGA trials required by each of five different PDA methods: the direct method, the original EPANET 2.2 CEM-1, the modified CEM-2, the GGA-PDA method and the ASM. The entries ‘*’ indicate that convergence was reached with under-relaxation using DAMPLIMIT = 0.1, while the entries ‘**’ indicate that convergence was reached with under-relaxation using DAMPLIMIT = 10 option to reach convergence. The table entries ‘***’ indicate that converge was not reached after 100 trials even after using different under-relaxation weights by setting DAMPLIMIT = 0.1 or DAMPLIMIT = 10. The DAMPLIMIT
option limits flow changes to 60% of the change computed by the GGA at each N–R iteration when Equation (29) convergence criterion drops below the DAMPLIMIT.
Convergence results for tests run at normal pressure ranges
| Network | Pressure range (m) | # Junctions affected | Demand % reduction | DDA | PDA ‘DIRECT’ | PDA EPANET 2.2 | BBA EPANET WB | GGA EPANET WB | ASM EPANET WB |
|---|---|---|---|---|---|---|---|---|---|
| NET1 (DEMAND×5) | 10 | 5 | 27.41 | 4 | *** | 5 | 5 | 5 | 5 |
| 20 | 8 | 26.05 | 10 | 5 | 5 | 5 | 5 | ||
| 30 | 8 | 35.35 | 8 | 5 | 5 | 5 | 5 | ||
| NET4 (DEMAND×1) | 10 | 9 | 40.10 | 5 | 5 | 5 | 5 | 5 | 5 |
| 20 | 9 | 49.21 | 5 | 5 | 5 | 5 | 5 | ||
| 30 | 9 | 55.18 | 5 | 5 | 6 | 6 | 6 | ||
| ANYTOWN (DEMAND×1) | 10 | 9 | 66.47 | 6 | 13 | 9 | 10 | 10 | 10 |
| 20 | 12 | 48.42 | 23 | 9 | 9 | 9 | 9 | ||
| 30 | 16 | 40.35 | 8 | 6 | 7 | 7 | 7 | ||
| NET2 (DEMAND×1) | 10 | 0 | 0 | 5 | 5 | 5 | 5 | 5 | 5 |
| 20 | 2 | 16.68 | 5 | 5 | 5 | 5 | 5 | ||
| 30 | 12 | 8.98 | 5 | 5 | 5 | 5 | 5 | ||
| BWN_2 (DEMAND×1) | 10 | 2 | 67.21 | 10 | 20 | 14 | 9 | 9 | 9 |
| 20 | 8 | 20.07 | 20 | 14 | 9 | 9 | 9 | ||
| 30 | 36 | 13.04 | 29 | 16 | 10 | 10 | 10 | ||
| BALERMA (DEMAND×1) | 10 | 326 | 45.49 | 3 | *** | 11 | 11 | 11 | 11 |
| 20 | 390 | 41.49 | *** | 10 | 9 | 9 | 9 | ||
| 30 | 426 | 41.01 | *** | 7 | 7 | 7 | 7 | ||
| WOLF (DEMAND×5) | 10 | 406 | 63.19 | 9 | 11 | 14 | 14 | 13 | 13 |
| 20 | 566 | 49.85 | 9 | 13 | 13 | 13 | 13 | ||
| 30 | 795 | 33.85 | 9 | 36 | 12 | 13 | 13 | ||
| EXNET (DEMAND×1) | 10 | 346 | 27.05 | 6 | 7 | 6 | 6 | 6 | 6 |
| 20 | 840 | 20.01 | 6 | 7 | 7 | 7 | 7 | ||
| 30 | 1,414 | 19.15 | 5 | 7 | 7 | 7 | 7 | ||
| KY17 (DEMAND×5) | 10 | 184 | 45.85 | 5 | 9 | 8 | 9 | 9 | 9 |
| 20 | 292 | 21.78 | 9 | 8 | 8 | 8 | 8 | ||
| 30 | 495 | 22.66 | 9 | 7 | 7 | 7 | 7 | ||
| NW-MODEL (DEMAND×5) | 10 | 0 | 0 | 8 | 8 | 8 | 8 | 8 | 8 |
| 20 | 910 | 8.56 | 8 | 8 | 8 | 8 | 8 | ||
| 30 | 5,833 | 9.06 | 8 | 8 | 8 | 8 | 8 | ||
| DMAS (DEMAND×1) | 10 | 2 | 93.97 | 6 | 6 | 6 | 6 | 6 | 6 |
| 20 | 14 | 10.02 | 6 | 6 | 6 | 6 | 6 | ||
| 30 | 223 | 6.75 | 6 | 6 | 6 | 6 | 6 | ||
| BWNS_2 (DEMAND×5) | 10 | 234 | 12.93 | 6 | 8 | 6 | 6 | 6 | 6 |
| 20 | 2,072 | 12.40 | 6 | 6 | 8 | 8 | 8 | ||
| 30 | 6,209 | 11.19 | 6 | 6 | 6 | 6 | 6 | ||
| NJ1 (DEMAND×5) | 10 | 11 | 38.43 | 6 | 6 | 6 | 6 | 6 | 6 |
| 20 | 673 | 7.39 | 6 | 6 | 6 | 6 | 6 | ||
| 30 | 2,544 | 11.01 | 6 | 6 | 6 | 6 | 6 | ||
| VIRT_ROME (DEMAND×1) | 10 | 49 | 15.36 | 5 | 6 | 12 | 5 | 5 | 5 |
| 20 | 413 | 17.06 | 6 | 12 | 5 | 5 | 5 | ||
| 30 | 1,330 | 18.83 | 6 | 12 | 5 | 5 | 5 |
| Network | Pressure range (m) | # Junctions affected | Demand % reduction | DDA | PDA ‘DIRECT’ | PDA EPANET 2.2 | BBA EPANET WB | GGA EPANET WB | ASM EPANET WB |
|---|---|---|---|---|---|---|---|---|---|
| NET1 (DEMAND×5) | 10 | 5 | 27.41 | 4 | *** | 5 | 5 | 5 | 5 |
| 20 | 8 | 26.05 | 10 | 5 | 5 | 5 | 5 | ||
| 30 | 8 | 35.35 | 8 | 5 | 5 | 5 | 5 | ||
| NET4 (DEMAND×1) | 10 | 9 | 40.10 | 5 | 5 | 5 | 5 | 5 | 5 |
| 20 | 9 | 49.21 | 5 | 5 | 5 | 5 | 5 | ||
| 30 | 9 | 55.18 | 5 | 5 | 6 | 6 | 6 | ||
| ANYTOWN (DEMAND×1) | 10 | 9 | 66.47 | 6 | 13 | 9 | 10 | 10 | 10 |
| 20 | 12 | 48.42 | 23 | 9 | 9 | 9 | 9 | ||
| 30 | 16 | 40.35 | 8 | 6 | 7 | 7 | 7 | ||
| NET2 (DEMAND×1) | 10 | 0 | 0 | 5 | 5 | 5 | 5 | 5 | 5 |
| 20 | 2 | 16.68 | 5 | 5 | 5 | 5 | 5 | ||
| 30 | 12 | 8.98 | 5 | 5 | 5 | 5 | 5 | ||
| BWN_2 (DEMAND×1) | 10 | 2 | 67.21 | 10 | 20 | 14 | 9 | 9 | 9 |
| 20 | 8 | 20.07 | 20 | 14 | 9 | 9 | 9 | ||
| 30 | 36 | 13.04 | 29 | 16 | 10 | 10 | 10 | ||
| BALERMA (DEMAND×1) | 10 | 326 | 45.49 | 3 | *** | 11 | 11 | 11 | 11 |
| 20 | 390 | 41.49 | *** | 10 | 9 | 9 | 9 | ||
| 30 | 426 | 41.01 | *** | 7 | 7 | 7 | 7 | ||
| WOLF (DEMAND×5) | 10 | 406 | 63.19 | 9 | 11 | 14 | 14 | 13 | 13 |
| 20 | 566 | 49.85 | 9 | 13 | 13 | 13 | 13 | ||
| 30 | 795 | 33.85 | 9 | 36 | 12 | 13 | 13 | ||
| EXNET (DEMAND×1) | 10 | 346 | 27.05 | 6 | 7 | 6 | 6 | 6 | 6 |
| 20 | 840 | 20.01 | 6 | 7 | 7 | 7 | 7 | ||
| 30 | 1,414 | 19.15 | 5 | 7 | 7 | 7 | 7 | ||
| KY17 (DEMAND×5) | 10 | 184 | 45.85 | 5 | 9 | 8 | 9 | 9 | 9 |
| 20 | 292 | 21.78 | 9 | 8 | 8 | 8 | 8 | ||
| 30 | 495 | 22.66 | 9 | 7 | 7 | 7 | 7 | ||
| NW-MODEL (DEMAND×5) | 10 | 0 | 0 | 8 | 8 | 8 | 8 | 8 | 8 |
| 20 | 910 | 8.56 | 8 | 8 | 8 | 8 | 8 | ||
| 30 | 5,833 | 9.06 | 8 | 8 | 8 | 8 | 8 | ||
| DMAS (DEMAND×1) | 10 | 2 | 93.97 | 6 | 6 | 6 | 6 | 6 | 6 |
| 20 | 14 | 10.02 | 6 | 6 | 6 | 6 | 6 | ||
| 30 | 223 | 6.75 | 6 | 6 | 6 | 6 | 6 | ||
| BWNS_2 (DEMAND×5) | 10 | 234 | 12.93 | 6 | 8 | 6 | 6 | 6 | 6 |
| 20 | 2,072 | 12.40 | 6 | 6 | 8 | 8 | 8 | ||
| 30 | 6,209 | 11.19 | 6 | 6 | 6 | 6 | 6 | ||
| NJ1 (DEMAND×5) | 10 | 11 | 38.43 | 6 | 6 | 6 | 6 | 6 | 6 |
| 20 | 673 | 7.39 | 6 | 6 | 6 | 6 | 6 | ||
| 30 | 2,544 | 11.01 | 6 | 6 | 6 | 6 | 6 | ||
| VIRT_ROME (DEMAND×1) | 10 | 49 | 15.36 | 5 | 6 | 12 | 5 | 5 | 5 |
| 20 | 413 | 17.06 | 6 | 12 | 5 | 5 | 5 | ||
| 30 | 1,330 | 18.83 | 6 | 12 | 5 | 5 | 5 |
***No convergence reached in 100 iterations.
Some significant findings revealed by Table 4 are:
All of the PDA methods produced the same demand solution for each pressure range run, which is why there is only a single set of affected node count and percent demand reduction results as shown in the table.
The four inverse PFR methods are consistently able to converge in a number of trials that is not much greater than that required for the DDA.
There is essentially no difference in the performance of the three new inverse methods (CEM-2, GGA-PDA and ASM) and they are slightly more efficient than the CEM-1 method currently used in EPANET 2.2.
For most cases, the direct FPR method requires the same or slightly more trials than the inverse methods but there are five cases where the method has convergence problems making it less reliable than the inverse methods.
For some networks, there appears to be a weak association between the required pressure and the number of trials required to converge.
Table 5 displays the results found for the stress test runs where the pressure range is only 0.1 m (0.142 psi). Over this narrow range, the demand versus pressure curve is almost vertical, approximating the case where nodes have no demand at pressure below some minimum value and full demand at a slightly higher pressure. The contents of the table's columns are the same quantities as those displayed in Table 4. Like the runs made for the larger pressure ranges, all the PDA methods produce the same demand solution (when they converge).
Convergence results for tests run at minimal pressure range
| Network | Pressure range (m) | # Junctions affected | Demand % reduction | DDA | PDA ‘DIRECT’ | PDA EPANET 2.2 | BBA EPANET WB | GGA EPANET WB | ASM EPANET WB |
|---|---|---|---|---|---|---|---|---|---|
| NET1 (DEMAND×5) | 0.1 | 2 | 56.37 | 4 | *** | 10 | 9 | 9 | 9 |
| NET4 (DEMAND×1) | 0.1 | 1 | 17.83 | 5 | *** | 5 | 6 | 5 | 5 |
| ANYTOWN (DEMAND×1) | 0.1 | 7 | 72.10 | 6 | *** | 10 | 13 | 13 | 13 |
| NET2 (DEMAND×1) | 0.1 | 0 | 0 | 5 | 5 | 5 | 5 | 5 | 5 |
| BWN_2 (DEMAND×1) | 0.1 | 0 | 0 | 10 | 20 | 32 | 23 | 22 | 22 |
| BALERMA (DEMAND×1) | 0.1 | 154 | 81.53 | 3 | *** | 63 | 17 | 17 | 17 |
| WOLF (DEMAND×5) | 0.1 | 221 | 97.13 | 9 | *** | 90** | 43 | 41 | 41 |
| EXNET (DEMAND×1) | 0.1 | 30 | 64.93 | 6 | *** | 18 | 12 | 12 | 12 |
| KY17 (DEMAND×5) | 0.1 | 95 | 94.65 | 5 | *** | 13 | 34 | 42 | 42 |
| NW-MODEL (DEMAND×5) | 0.1 | 0 | 0 | 8 | 8 | 8 | 8 | 8 | 8 |
| DMAS (DEMAND×1) | 0.1 | 1 | 100 | 6 | 6 | 6 | 6 | 6 | 6 |
| BWNS_2 (DEMAND×5) | 0.1 | 0 | 0 | 6 | 6 | 6 | 6 | 6 | 6 |
| NJ1 (DEMAND×5) | 0.1 | 0 | 0 | 6 | 6 | 6 | 6 | 6 | 6 |
| VIRT_ROME (DEMAND×1) | 0.1 | 0 | 0 | 5 | 6 | 18 | 8 | 12 | 12 |
| Network | Pressure range (m) | # Junctions affected | Demand % reduction | DDA | PDA ‘DIRECT’ | PDA EPANET 2.2 | BBA EPANET WB | GGA EPANET WB | ASM EPANET WB |
|---|---|---|---|---|---|---|---|---|---|
| NET1 (DEMAND×5) | 0.1 | 2 | 56.37 | 4 | *** | 10 | 9 | 9 | 9 |
| NET4 (DEMAND×1) | 0.1 | 1 | 17.83 | 5 | *** | 5 | 6 | 5 | 5 |
| ANYTOWN (DEMAND×1) | 0.1 | 7 | 72.10 | 6 | *** | 10 | 13 | 13 | 13 |
| NET2 (DEMAND×1) | 0.1 | 0 | 0 | 5 | 5 | 5 | 5 | 5 | 5 |
| BWN_2 (DEMAND×1) | 0.1 | 0 | 0 | 10 | 20 | 32 | 23 | 22 | 22 |
| BALERMA (DEMAND×1) | 0.1 | 154 | 81.53 | 3 | *** | 63 | 17 | 17 | 17 |
| WOLF (DEMAND×5) | 0.1 | 221 | 97.13 | 9 | *** | 90** | 43 | 41 | 41 |
| EXNET (DEMAND×1) | 0.1 | 30 | 64.93 | 6 | *** | 18 | 12 | 12 | 12 |
| KY17 (DEMAND×5) | 0.1 | 95 | 94.65 | 5 | *** | 13 | 34 | 42 | 42 |
| NW-MODEL (DEMAND×5) | 0.1 | 0 | 0 | 8 | 8 | 8 | 8 | 8 | 8 |
| DMAS (DEMAND×1) | 0.1 | 1 | 100 | 6 | 6 | 6 | 6 | 6 | 6 |
| BWNS_2 (DEMAND×5) | 0.1 | 0 | 0 | 6 | 6 | 6 | 6 | 6 | 6 |
| NJ1 (DEMAND×5) | 0.1 | 0 | 0 | 6 | 6 | 6 | 6 | 6 | 6 |
| VIRT_ROME (DEMAND×1) | 0.1 | 0 | 0 | 5 | 6 | 18 | 8 | 12 | 12 |
***No convergence reached in 100 iterations.
*Convergence reached with DAMPLIMIT = 0.1.
**Convergence reached with DAMPLIMIT = 10.
The following points can be drawn from the results shown in Table 5:
Unlike for the tests with larger pressure ranges, some networks require considerably more trials to converge when compared to a DDA. However, the inverse PFR methods also demonstrated good robustness in the stress tests (they always converged).
The three new inverse PFR methods (CEM-2, GGA-PDA and ASM) are about equally efficient and reliable.
The EPANET 2.2 PDA method (CEM-1) can be less reliable than the other inverse PFR methods – for three networks it required the use of under-relaxation to converge in significantly more trials.
The direct method is even less reliable than the others at this small pressure range as 6 of 14 networks failed to converge.
CONCLUSIONS
The objective of this study was to compare the performance of several different PDA methods on a significant number of both small and large pipe networks, using a modified version of the EPANET water distribution system software. Five different PDA methods that can all be implemented within the framework of the N–R-based GGA network solver were evaluated. One of these uses the direct form of the bounded Wagner power function of demand versus pressure. The other four all use the inverted form of this function. Each method was applied to 14 different networks, ranging in size from 9 to over 150,000 nodes, over four different pressure ranges resulting in almost 300 runs of the modified EPANET program.
The results of this analysis show that:
The three new inverse PFR methods (CEM-2, GGA-PDA and ASM) were all equally efficient and reliable (they were all able to converge for all test cases).
The inverse PFR method currently used in EPANET 2.2 (CEM-1) is slightly less efficient than the other inverse methods for typical pressure range requirements and can be less reliable for an extremely narrow pressure range.
The direct PFR method can be just as efficient as the inverse PFR methods when it converges but having failed to converge for a number of test cases makes it less reliable.
The reliability of the direct PFR method can certainly be improved by higher order N–R approaches or by more sophisticated under-relaxation schemes, but in this work we show that the indirect approaches do not need them to converge.
In addition, although the indirect PDA increases the mathematical complexity of the DDA, the new implemented codes in EPANET do not seem to remarkably increase the number of required iterations to converge, particularly for larger WDNs.
Any of the three new PDA methods tested here would be a good replacement for the PDA method currently used by EPANET 2.2.
DATA AVAILABILITY STATEMENT
All relevant data are available from an online repository. The modified EPANET code together with the set of 14 input files describing the WDN used as test examples is made available. Interested users will be able of downloading it from the following website: https://www.unicas.it/siti/laboratori/lia-laboratorio-di-ingegneria-delle-acque/link/epanet-workbench.aspx.
%20cropped.png?versionId=7142)