## 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

*k*th 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:where 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*.

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

*Z*is the elevation of the node. The exponent 1/

_{i}*c*is typically taken as 0.5 to mimic flow through an orifice opening.

### Direct PDA method

*et al.*(1988) equation of Equation (6), an diagonal matrix must then be defined, where the

*i*th diagonal element is given by:and the resolving equations (see Todini

*et al.*2021) become:withwhere the matrices

**A**and

**F**are given by Equations (3) and (4), respectively.

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

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

*BIG*is a large positive number, such as 10

^{9}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

*d*of the

_{i}*i*th 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.

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

*et al.*(1978) ‘content’, while in the GGA-based approach they are directly obtained by expanding the original GGA set of equations with the introduction of an additional set of equations describing the Wagner

*et al.*(1988) inverse FPR, to link the unknown delivered demands 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:with ; and where

*j*indicates the position of the

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

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.

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 .

*δ*belonging to the active set of nodes then becomes:The detailed derivation of these equations can be found in Todini

*et al.*(2021).

*et al.*(2021). A very small positive number is added to in the estimation of , to avoid matrix inversion problems when the estimated demand value approaches zero, leading to:This has the advantage of guaranteeing the continuity all along the Wagner equation while marginally modifying its value.

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.

Each network was run under four different pressure ranges (*P _{req}–P_{min}*) 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.

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 |

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

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

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.

## REFERENCES

*Minima of Functions of Several Variables with Inequalities as Side Constraints*

*MSc Thesis*