https://orcid.org/0000-0002-3521-6632
Abstract
Modeling of pressure-dependent users’ consumption is mandatory to simulate accurately the hydraulics of water distribution networks (WDNs). Several software solutions already exist for this purpose, but none of them actually permits the easy integration and testing of new physical processes. In this paper, we propose a new Python simulator that implements a state-of-the-art pressure-dependent model (PDM) of users’ consumptions based on the Wagner’s pressure–outflow relationship (POR). We tested our simulator on eight large and complex WDNs, for different levels of users’ demands. The results show similar precision and efficiency to the ones obtained by the authors of the original model with their MATLAB implementation. Moreover, in case of fully satisfied users’ demands, our simulator provides the same results as EPANET 2.0 in comparable computational times. Finally, our simulator is integrated into the open-source, collaborative, multi-platform, and Git versioned Python framework OOPNET (Object-Oriented Python framework for water distribution NETworks analyses); thus, it can be easily reused and/or extended by a large community of WDN modelers. All this work represents a preliminary step before the incorporation of new processes such as valves, pumps, and pressure-dependent background leakage outflows.
HIGHLIGHTS
A new simulator of pressure-dependent consumption in water distribution network.
Coded in Python, based on a state-of-the-art MATLAB model, easy to extend.
Numerical experiments on networks composed of up to 19,647 pipes and 17,986 nodes.
Instabilities handled through regularization and line search along the Newton descent.
Preliminary step for integration of new processes and first contribution to OOPNET.
INTRODUCTION
Users’ consumption in WDNs
Households, schools, hospitals, businesses, food-processing and pharmaceutical industries, fire departments, etc., are all significant users of water distribution networks (WDNs). Each of them has its own need of water; moreover, this need can vary over time. Over a specific period, a need corresponds to a “demand”, expressed as an outflow rate. From this demand, and depending on the capacities of the WDN, the amount of water a user actually consumes is called its “consumption”. A user’s consumption is between 0 (i.e., no water is provided to the user) and the demand (i.e., the demand of the user is fully satisfied). For example, in the last three decades, the worldwide amount of drinking water consumed daily by households and public services is estimated to be 
, which represents 510,995 Olympic-size swimming pools of 
, and 11 percent of total freshwater withdrawals (Ritchie & Roser 2017). In France, this amount is of 
per day, which leads to an average domestic consumption of 1481 per day and inhabitant (Lao et al. 2022).
The first goal of WDNs is to satisfy the demands of the users. To fulfill these needs, WDNs must be well sized at their building, and then extended according to new needs (e.g., population growing or installation of new facilities). Pumps are used to maintain enough pressure to carry water from treatment plant or natural source to users. Also, valves are used to control the pressure, and to divide large networks in several district-metered areas (DMAs) that are easier to monitor with the help of sensors. Appropriate design and equipment permit us to limit the construction and functioning costs, by reducing, for example, the initial number of pipes, maintenance operations, energy consumption of pumps and valves, number of installed sensors, and amount of water losses due to leakages (Alperovits & Shamir 1977; Gupta et al. 1993; Walski et al. 2007; Gupta et al. 2016; Creaco Walski 2017). To optimize the choice of decision-makers, the use of pressure-driven models of users’ consumption is mandatory (Karadirek et al. 2012; Berardi et al. 2015; Laucelli et al. 2017a,b; Monsef et al. 2018; Koşucu Demirel 2022).
Pressure-driven modeling of users’ consumption
One approach to model users’ consumption is by assuming that the demand of the users is always satisfied. These models are called “demand-driven models” (DDMs); they neglect unsatisfied users’ demand due to insufficient pressure. Nevertheless, this solution is implemented in several softwares, such as EPANET 2.0 (Rossman 2000), its object-oriented C++ interface OOTEN (van Zyl et al. 2003), and Porteau (Piller et al. 2011).
Conversely, the “pressure-driven models” (PDMs) permit us to simulate users’ consumptions that depend on the pressure. Different pressure–outflow relationships (PORs) can be used in this kind of model (Wagner et al. 1988; Fujiwara & Ganesharajah 1993; Gupta & Bhave 1996; Piller et al. 2003; Tanyimboh Templeman 2004). Among them, the Wagner’s POR (Wagner et al. 1988) was proven to perform best against data (Shirzad et al. 2013). Simulators using a POR to model PDM consumptions include the extension of EPANET 2.0 from Siew & Tanyimboh (2012), EPANET 2.2 (Rossman et al. 2020), WNTR (Water Network Tool for Resilience) (Klisel et al. 2018), and the MATLAB simulator developed by Elhay et al. (2016).
Pressure-driven analyses (PDAs) can also be conducted using DDM simulators by adding artificial equipment to the network (e.g., flow and pressure-regulating valves) (Mahmoud et al. 2017; Vaidya et al. 2023). However, this approach increases the size (i.e., the number of elements) of the system, and can require specific treatments and convergence criteria to ensure the robustness and the stability of the algorithm for every test case (Gorev et al. 2021; Sivakumar et al. 2023).
The engines and the programmer’s toolkit of EPANET 2.2 are multi-platform, but they are non objected-oriented; also, EPANET 2.2’s graphical user interface, which makes its engines easier to use, works only on Microsoft Windows operating systems. The WNTR package (Klisel et al. 2018) was designed specifically to simulate and analyze resilience of water distribution networks under disaster scenarios; moreover, its software architecture does not permit us to extend it easily, in particular to add new physical processes like pressure-dependent background leakage outflows. Finally, the MATLAB simulator from Elhay et al. (2016) is efficient and works fine, but any upgrade in its code will need to be run in MATLAB, a paid-for software package not available to everyone (although interfacing with other languages is possible). In summary, there currently exists no open-source and easy/cheap to extend simulator of PDM consumption.
The Python language
Python is a free, multi-platform, high-level and object-oriented language (Lutz 2013). Its core philosophy emphasizes readability, making the codes written in Python easy to reuse and extend, even by people who are not software developers. Many efficient libraries for scientific computing are available in Python (Hunter 2007; Chen et al. 2008; Hagberg et al. 2008; McKinney 2011; Virtanen et al. 2020), and Python is today one of the most popular programing languages1 .
Big efforts have already been made by the community to simulate WDNs with Python. For example, the software of Xing Sela (2020) allows the running of transient simulations using an adapted method of characteristics. For extended-period simulations (EPS), the EPANET-based package WNTR (Klisel et al. 2018), and the object-oriented framework OOPNET (Steffelbauer & Fuchs-Hanusch 2015) are stable and active solutions.
The framework OOPNET
OOPNET (Object-Oriented Python framework for water distribution NETworks analyzes) is a convenient Python API (Application Programing Interface) of EPANET. It is aimed at students, engineers and researchers who want to run hydraulic simulations in a high-level way, apply pre and post-processing on input or output data of EPANET, and make the simulation results from EPANET more readable (Steffelbauer & Fuchs-Hanusch 2015).
OOPNET is based on several state-of-the-art Python libraries for scientific computing (Hunter 2007; Hagberg et al. 2008; McKinney 2011; Hold-Geoffroy et al. 2014; Hoyer & Hamman 2017). It permits us in particular to parse and convert EPANET’s input file to Python objects, run an EPANET simulation through a command line interface, and parse EPANET’s output files and convert them to Python objects. OOPNET is free, open-source, object-oriented, and designed to enhance collaborative effort from the community of programmers and users of WDN software. Its source code is publicly accessible on GitHub (https://github.com/oopnet/oopnet). Moreover, its authors and contributors currently work on making OOPNET more flexible, user-focused and efficient through a domain-driven design (DDD) approach (Steffelbauer et al. 2022). Thus, it is well suited for prototyping, integrating new physical processes, or testing new hypotheses.
Hypothesis, objectives, and research strategy
The consumption model proposed by Elhay et al. (2016) seems the most effective one; however, this model is not easy to extend and it needs a MATLAB environment. Thus, our objective is to implement a new simulator that will later permit us to test new physical processes in a more convenient way and to share our developments with a wider community of users and programmers.
We think that using the Python programing language and the framework OOPNET is today the best solution for developing such a simulator. Doing so, we expect to reproduce the results obtained by Elhay et al. (2016) with their MATLAB implementation when running PDA, and the ones obtained with EPANET 2.0 for demand-driven analyses (DDAs). It is important to note here that EPANET 2.2 was not available yet at the start of our study; this is the reason why we decided to validate our Python simulator against the MATLAB implementation of Elhay et al. (2016) rather than EPANET 2.2 when running PDA. Finally, developing this new Python simulator would permit us to extend in a non-regressive way the current capabilities of OOPNET, which by now allows running EPANET simulations only.
In this paper, we will first describe the steady-state equations and the method to solve them, presenting several numerical enhancements to deal with potential sources of instability. Then, we will explain our Python implementation, the networks to simulate, and the tests and metrics proposed for validation. Finally, we will present our results and discuss them.
METHODS
Hereafter, since we do not model any other type of demand and consumption, we will call “user’s demand” just “demand”, and “user’s consumption” simply “consumption”.
(called just “head” hereafter) will be calculated as:
the elevation (in m) and 
the pressure-head (Walski et al. 2007, p. 29). 
and 
, as well as friction head-losses, will be expressed in mH
O (metres water column), which are consistent to m. Also, as in Piller et al. (2011), all flow rates, demand and consumption will be expressed in 
rather than in 
, to avoid problems of stability due to machine precision.Finally, and in the absence of other indications, scalar parameters and variables will be denoted in italic (e.g., 
), vectors in italic bold (e.g., 
), matrices in italic bold upper-case (e.g., 
), scalar functions in upright (e.g., 
), vector functions in upright bold (e.g., 
), matrix functions in upright bold upper-case (e.g., 
), and sets in blackboard style (e.g., 
).
Friction head-loss
the inside diameter (in mm) of the pipe and 
its roughness coefficient (unit-less), and supposing that the pipe is cylindrical, we first compute the friction coefficient as:
(in m) of the pipe as:
is the algebraic flow rate traversing the pipe, and 
is the Hazen-Williams exponent (unit-less).Consumption
the demand, 
the pressure-head, and 
and 
the fixed minimum and service pressure-heads (in mH
O) such that 
, we compute the fraction of pressure-head (unit-less) as:
and 
(Piller et al. 2003; Giustolisi et al. 2008, 2014; Elhay et al. 2016; Deuerlein et al. 2019).Equations of equilibrium in a WDN
the number of cylindrical and longitudinal pipes of lengths ℓ
, 
the number of junction nodes at which the heads are unknown, 
and 
respectively the number of tank and reservoir nodes at which the heads are known and fixed, 
the total number of source nodes, and 
the total number of nodes. Also, we denote 
and 
the unknown flow rates in pipes and heads at junctions, 
and 
the known and fixed heads at tanks and at reservoirs, 
and 
the heads at all source nodes and at all nodes, 
and 
the pressure-heads at junctions and the water levels at tanks, and 
, 
and 
the elevations of junctions, tanks’ bottom, and reservoirs’ surface, respectively. Then, we have the relations:
junction-pipe, the 
tank-pipe, and the 
reservoir-pipe incidence matrices as respectively 
, 
and 
, the 
source-pipe incidence matrix as 
, the 
node-pipe incidence matrix as 
, and 
the vector function of 
to 
to compute the friction head-losses in pipes. 
, we compute the friction head-loss 
along the full pipe 
from Equation (3), as:
the fixed demand at junctions, and 
the vector function of 
to 
to compute consumptions. 
, 
is computed from Equation (5) as:
In this study, for the sake of simplicity, we choose the same 
for all junctions. However, it presents no difficulty to set specific values for some junctions, to describe for example pressure-independent fire demands as in Sivakumar et al. (2023).
and heads 
in the WDN at steady-state, we solve the non-linear system of equations:
found satisfies the principle of the conservation of energy, corresponding, in any pipe 
of the network, to the Bernoulli’s equation:
is the friction head-loss along 
, and 
and 
are respectively the elevations and the pressure-heads at the nodes 
and 
located at the extremities of 
; as a reminder, the kinetic energy is not represented in Equation (14) because it is negligible compared to the pressure-head (Walski et al. 2007, p. 29). Also, the flow rates 
satisfy the conservation of the mass at nodes, derived from the Kirchhoff’s first law: “the algebraic sum of the branch currents towards any node of an electric network is zero” (IEC 2023).Solution algorithm
We solve the system (11) with a Newton’s method that includes a system reduction.
Newton’s method
the initial guesses of flow rates in pipes, and 
the ones of heads at junctions,and
the Jacobian matrix of 
at iteration 
,
is invertible, the Newton’s method consists of computing, at each iteration 
, the flow rates 
and heads 
as
System reduction
and 
by inverting numerically the potentially very large matrix 
can take much computational time. Thus, we rather look for the descent directions on 
and 
, defined respectively as
in function of 
:
is invertible. Next, we substitute (21) into the second equation of the system (20) to obtain:
is the Schur complement of the block 
of the matrix 
, and is defined as:
is symmetric and 
is diagonal. Also, supposing that none of the flow rates in vector 
are zero, then both 
and 
are positive-definite. Thus, 
is symmetric and positive-definite, and we can use a Cholesky factorization to compute 
of Equation (22). Once 
is calculated, we can easily compute 
from Equation (21) because 
is diagonal and supposed invertible.
and 
, we can replace
and 
as
Initial guesses
, the initial guess of the flow rate in the pipe 
is chosen consistent with a velocity 
:
, we set the initial guess of the head at 
as:
Convergence criterion
,
. We discuss briefly the choice of the criterion in the Supplementary Material.Other convergence criteria are used in WDN simulators. For example, the one in EPANET 2.0 checks that the 1-norm of flow changes divided by the 1-norm of flows is less than a given tolerance (Rossman 2000, p. 94). EPANET 2.2 improves this criterion by considering also the 
-norm of flow changes and the 
-norm of energy residuals (Rossman et al. 2020, p. 116). Then, Gorev et al. (2022) enhances EPANET 2.2’s criterion by checking for flow residual in each pipe too. Finally, Sivakumar et al. (2023) add one more criterion that verifies the residuals on the inverse of the POR (e.g., the Wagner’s POR) to improve the convergence of EPANET 2.2’s solution algorithm. In our study, we don’t need the criterion proposed by Sivakumar et al. (2023) because our formulation does not consider the POR inverse function. Therefore, we choose to simply reuse the criterion from Elhay et al. (2016), since preliminary tests showed good performances, and to be consistent with the other technical choices made by Elhay et al. (2016).
Sources of instabilities
Different sources of instabilities can hinder the convergence of the solution algorithm described in Section 2.4. We present them below, along with numerical enhancements to deal with them. These numerical enhancements are of primary importance to ensure the convergence of the solution algorithm in many common situations. However, to preserve readability we choose to describe them in the Supplementary Material.
The numerical enhancements presented in Sections 2.5.1 – 2.5.3 are the same as the ones implemented in the Porteau software (Piller et al. 2011).
Pipes with zero flow rate
, the element 
of the Jacobian matrix is equal to the derivative of the friction head-loss, 
, which can be computed as:
for 
. Thus, since we need to compute the inverse of 
at each iteration (see Equations (21), (20) and (23)), the use of Equation (28) to compute 
could lead to a division by zero. To prevent such error, we choose, like Piller (1995), to regularize the friction head-loss and its derivative for 
close to 
, by respectively a cubic and a quadratic polynomial.Junction nodes with pressure-head close to the minimum or the service pressure-head
, 
is equal to 
, which can be computed as:
and 
defined by Equations (4) and (5). However, Equation (31) is undefined at 
and discontinuous at 
. Thus, using Equation (31) to compute 
could lead to convergence failure when 
is close to 
or 
. To make Equation (31) continuous 
, we choose, like Piller et al. (2003), to regularize the consumption and its derivative by respectively a cubic and a quadratic polynomial, when 
is close to 
and 
.Initial guesses far from the solution and/or objective function with (sub-)linear gradient growth
The bottleneck of the solution algorithm is the solving of the system (22). But the algorithm can need numerous iterations (and therefore numerous solving of system (22)) if the initial guesses of the flow rates and heads at junctions are far from the solutions at equilibrium. In this case, the simulation of large WDNs can become unfeasible in reasonable computational time. Moreover, the convergence of the algorithm is no longer guaranteed if the gradient of the objective function (11) to minimize does not have a super-linear growth (Piller et al. 2017). Thus, to reduce the number of iterations needed by the solution algorithm and to guarantee its convergence for any network configuration, we choose to reuse the damping method proposed by Elhay et al. (2016) (see description of the method in the Supplementary Material).
Some formulations of Equation (11) do not need damping corrections because they consider the inverse of the POR function rather than the POR function (Todini et al. 2021). However, to achieve full convergence, these formulations require a supplementary criterion on the inverse POR equation (Sivakumar et al. 2023).
WDN with highly contrasted values of flow rates and/or heads
In general, at each iteration of a solution algorithm based on Newton’s method, the magnitude of the elements in the Jacobian matrix can vary strongly from one part of a network to another, according to demands, pipe resistances, tank levels, use of pumps and/or valves, etc. In the present study, when demands and/or pipe resistances differ by several orders of magnitude across the network, using the Jacobian matrix 
as it can cause instabilities and increase significantly the number of iterations needed to converge. To limit the difference of magnitude order between all elements of 
, we choose to extend the preconditioning method initially proposed by Elhay & Simpson (2011). Indeed, in Elhay & Simpson (2011), the method dealt with consumptions that do not depend on the pressures; thus, we extend it to now deal with pressure-dependent consumptions too. Using this preconditioning method, the solution algorithm becomes a quasi-Newton method.
Gorev et al. (2022) propose an alternative solution to deal with zero flows and low-resistance pipes in EPANET 2.2. To do so, they replace the head-loss by a linear function in the pipes where the flow rate becomes very small, avoiding the appearance of zero head-loss derivatives in the Jacobian matrix. This regularization approach preserves the convergence rate of the Newton method, which is generally better than that of quasi-Newton methods. Nevertheless, it also changes slightly the head-loss behavior (for low flow rates), it uses an absolute threshold flow rate for the cutoff of the derivatives (while the preconditioning method from Elhay & Simpson (2011) is based on the difference of magnitude order between the elements of the Jacobian), and it needs some additional checking computations. Thus, both the Elhay & Simpson (2011) and Gorev et al. (2022) approaches have their own advantages and limitations. We choose to use the one from Elhay & Simpson (2011) because we master it better.
Implementation and framework
To implement our Python simulator, we use the Python programing language2 through its Anaconda distribution3. In addition to the Python’s standard library4, we use several third-party libraries:
NetworkX (Hagberg et al. 2008), for the creation and manipulation of the graphed data structures of the WDNs,
NumPy (Harris et al. 2020), for the creation and manipulation of the vectors and basic linear algebra operations (e.g., computation of vector norms),
SciPy (Virtanen et al. 2020), to find the coefficients of the polynomials defined in the Supplementary Material5, and to handle sparse matrices6,
Scikit-Sparse7, which itself wraps the CHOLMOD routine (Chen et al. 2008), to compute the Cholesky decomposition of the Schur complement matrix,
XArray (Hoyer & Hamman 2017), to store the outputs in multi-dimensional labeled arrays,
Pandas (McKinney 2011), to convert the outputs to tabular data structures, to post-process them, and to write them to comma-separated values (CSV) files,
Matplotlib (Hunter 2007), to plot the outputs,
and Sphinx8, to extract the documentation from the source code, and generate an HTML file from it.
Simulated networks
, which numbers of nodes and pipes are presented in Table 1. EPANET’s input files of networks 
are freely downloadable from the American Society of Civil Engineers (ASCE)’s library9; they are plotted on Figure 1. Other networks 
can not be plotted here, because of ownership and/or security concerns. 
Networks 
used to test our Python simulator. Reservoirs are represented by 
.
Number of pipes 
, junctions 
and reservoirs 
in the networks 
used to test our Python simulator
| Network |  |  |  |
|---|---|---|---|
 | 934 | 848 | 8 |
 | 1,118 | 1,039 | 2 |
 | 1,976 | 1,770 | 4 |
 | 2,465 | 1,890 | 3 |
 | 2,508 | 2,443 | 2 |
 | 8,584 | 8,392 | 2 |
 | 14,830 | 12,523 | 7 |
 | 19,647 | 17,971 | 15 |
| Network | | | |
|---|---|---|---|
| | 934 | 848 | 8 |
| | 1,118 | 1,039 | 2 |
| | 1,976 | 1,770 | 4 |
| | 2,465 | 1,890 | 3 |
| | 2,508 | 2,443 | 2 |
| | 8,584 | 8,392 | 2 |
| | 14,830 | 12,523 | 7 |
| | 19,647 | 17,971 | 15 |
Tests and metrics
Functioning and performances of the Python simulator
presented in Section 2.7, multiplying the peak demand by a multiplier 
. Like so, we test the Python simulator for different levels of demand satisfaction (in %), defined for each multiplier 
and at any junction 
as:
and 
.From the outputs of the simulations, we first check that the 
-norm of the residuals is less than 
(in 
if the greatest absolute value is a mass residual, or in mH
O if it is an energy one). If so, we then assume that the solver reaches the equilibrium state with enough precision.
), defined for each demand multiplier 
as, respectively:
. Also, we check that the cumulated demand increases quicker than the cumulated consumption, and that the global demand satisfaction, defined for each demand multiplier 
as:
.
, the mean central processing unit (CPU) time elapsed during all runs with all 
, first considering the whole simulations (i.e., pre-processing, solver and post-processing):
is the cardinality of the set 
. We expect logically the elapsed CPU times to increase with the size of the networks. Also, since the critical part of the solver (i.e., the Cholesky decomposition of the Schur complement matrix) uses the optimized CHOLMOD compiled library (Chen et al. 2008), we expect that the times needed by the solver increase slower than those needed for the whole simulations when larger networks are simulated. To check this point, we compute, for each network 
, the ratio of time needed by the solver when compared to the times needed by the whole simulations, defined as:
for each demand multiplier 
, first with no numerical enhancement, and next adding successively each numerical enhancement. Hereafter, we denote this set of enabled enhancements as 
(“
” means “in addition to all previously enabled enhancement(s)”). We then compute, for each run, the convergence order of the solution algorithm when it converges. To do so, for each enabled enhancement(s) 
, each network 
and each demand multiplier 
, denoting 
the residuals at the iteration 
of the solution algorithm, we compute the difference between the 
-norm of 
and 
, 
, where 
is the number of iterations needed to reach convergence, as:
then satisfies the relation:
a constant. After rewriting Equation (39) in 
scale as:
by linear regression of Equation (40).
computed by Equation (38) can sometimes be negative, especially at 
or at an iteration that follows a damping correction. Thus, to avoid error in the computation of 
through Equation (40), we do not consider these negative 
, and we use as a previous iteration the last iteration 
such that 
. Equation (38) then becomes:
, we compute, for each enabled enhancement(s) 
, the mean convergence order as:
Comparison of the Python simulator with the MATLAB one
To validate the ability of our Python simulator to run pressure-driven analyses (PDA), we compare the consumptions computed by our simulator with the ones obtained by Elhay et al. (2016) from their MATLAB implementation, for the same sets of networks 
and demand multipliers 
, and the same minimum and service pressure-heads at all junctions: 
and 
. Then, we assume that the consumptions from both Python and MATLAB simulators are the same as soon as they all differ by less than 
. Supposing that the average flow rate in the tested networks is equal to 
, then 
represents a precision of 
, which corresponds to 1/10th of the best precision expected when measuring flow rates physically in a WDN (Itron 2021, p. 3).
The MATLAB simulator of Elhay et al. (2016) implements the numerical enhancements described in Sections 2.5.1 to 2.5.3.
, the mean convergence orders associated with each simulator, as:
is the convergence order obtained when simulating the network 
with the demand multiplier 
.
, the mean numbers of damping corrections associated with each simulator as:
is the number of damping corrections obtained when simulating the network 
with the demand multiplier 
. Since the damping algorithm implemented in both simulators is the same, we expect the simulators to need (more or less) the same number of damping corrections.Finally, it is not possible to compare the process times elapsed using each simulator, because the MATLAB simulator was run by Elhay et al. (2016) on a different machine than the one we use to run our Python simulator.
Comparison of the Python simulator with EPANET 2.0
To check that our Python simulator is also able to simulate pressure-independent consumptions, we simulate the same set of networks 
with EPANET 2.010 , increasing sufficiently the heads at the source nodes 
, and reducing enough the service pressure-head 
, so that the demand becomes satisfied at every junction. Then, we verify, at each junction of each network, that the consumption computed by the Python simulator and the one computed by EPANET 2.0 differ by less than 
.
Finally, to compare the performances of the Python simulator with those of EPANET 2.0, we compute, for both simulators and for each network, the 
-norms of the residuals at convergence, along with the CPU times elapsed running each simulator on the same machine: an Intel Core i9 with 32 GB of memory. For these comparisons to be relevant, we set the “ACCURACY” parameter of EPANET 2.0 to the same value as the one used in the convergence criterion of the Python simulator Equation (27), that is 
(unit-less). The “ACCURACY” parameter is used by EPANET 2.0’s solver to determine when the convergence is reached. Indeed, according to EPANET 2.0’s documentation (Rossman 2000, p. 153): “The trials end when the sum of all flow changes from the previous solution divided by the total flow in all links is less than this number’.
RESULTS AND DISCUSSION
This section presents the results obtained when running the tests explained in Section 2.8.
Functioning and performances of the Python simulator
When simulating the set of networks for each demand multiplier, the maximal 
-norm of all residuals is equal to 
(in 
or mH
O), which is much less than 
. Thus, the solver reaches the equilibrium with enough precision.
, and that the cumulated demand increases quicker than the cumulated consumption (Figure 2(a)). This last observation is confirmed when computing the global demand satisfaction (Figure 2(b)), which decreases when 
increases. This functioning is completely realistic. 
Cumulated demand and consumption (Figure 2(a)), and global demand satisfaction (Figure 2(b)), simulated in networks 
with the Python simulator, for each demand multiplier 
. (a) Cumulated demand and consumption (in 
); (b) Global demand satisfaction (in %).
Cumulated demand and consumption (Figure 2(a)), and global demand satisfaction (Figure 2(b)), simulated in networks with the Python simulator, for each demand multiplier . (a) Cumulated demand and consumption (in ); (b) Global demand satisfaction (in %).
Mean simulator and solver elapsed CPU times (Figure 3(a)), and ratio of mean solver elapsed CPU time when compared to mean simulator elapsed CPU time (Figure 3(b)). (a) Mean simulator and solver CPU times (in s); (b) Ratio of mean CPU time used by the solver (in %).
Mean simulator and solver elapsed CPU times (Figure 3(a)), and ratio of mean solver elapsed CPU time when compared to mean simulator elapsed CPU time (Figure 3(b)). (a) Mean simulator and solver CPU times (in s); (b) Ratio of mean CPU time used by the solver (in %).
. With just friction head-loss regularization (bars “Head-loss”), 96.88% of the simulations converge, with an order of convergence equal to 
. With friction head-loss and consumption regularization (bars “
Consumption”), still 96.88% of the simulations converge, with an order of convergence equal to 
. With friction head-loss regularization, consumption regularization and damping correction (bars “
Damping”), then 100% of the simulations converge, with an order of convergence equal to 
. Finally, with all numerical enhancements enabled (i.e., friction head-loss regularization, consumption regularization, damping correction and preconditioning; bars “
Preconditioning”), still 100% of the simulations converge with an order of convergence of 
. The “not-averaged” orders of convergence per network and demand multiplier are presented in the Supplementary Material. 
Ratio of converging runs (Figure 4(a)) and mean convergence orders (Figure 4(b)) when simulating all networks 
with the Python simulator, enabling successively each numerical enhancement. (a) Ratios of converging runs (in %); (b) Mean convergence orders (unit-less).
Ratio of converging runs (Figure 4(a)) and mean convergence orders (Figure 4(b)) when simulating all networks with the Python simulator, enabling successively each numerical enhancement. (a) Ratios of converging runs (in %); (b) Mean convergence orders (unit-less).
Globally, these metrics show that for most of the configurations (i.e., one network combined with one demand multiplier), we need to regularize the friction head-loss to reach convergence. Also, friction head-loss regularization, consumption regularization and damping correction must be enabled for the simulation of all configurations to reach convergence. If we consider only the networks and demand multipliers for which all simulations converge, then the mean convergence order is better when none of the numerical enhancements are enabled (i.e., 
). Even so, the mean convergence order remains very good (i.e., 
) when a part or all of the numerical enhancements are enabled. Adding preconditioning when all other numerical enhancements are already enabled does not show any improvement in these test cases. However, this does not mean that there would be no difference for other networks and configurations. Thus, since we do not observe computational overhead, we decide, for future simulations, to systematically enable all numerical enhancements, including the preconditioning.
Comparison of the Python simulator with the MATLAB one
When simulating each network for each demand multiplier, the maximal difference between the consumptions computed with the MATLAB and the Python simulators is equal to 
. This difference is not significant. Thus, both simulators give the same result, which means that our simulator is able to run pressure-driven analyses in a correct way.
Mean convergence orders (Figure 5(a)) and numbers of needed damping corrections (Figure 5(b)), simulating each network 
with the Python and MATLAB simulators. (a) Mean convergence orders; (b) Mean numbers of damping corrections.
Mean convergence orders (Figure 5(a)) and numbers of needed damping corrections (Figure 5(b)), simulating each network with the Python and MATLAB simulators. (a) Mean convergence orders; (b) Mean numbers of damping corrections.
For any network, the number of damping corrections is related to the complexity of the network: the higher it is, the more complex is the network. In our study, we could expect that this number increases with the size of the networks, i.e., that 
needs more damping corrections than 
, 
more than 
, etc. However, Figure 5(b) does not show any obvious proportionality relation between the size of the networks and the number of corrections. For example, even if they are smaller, networks 
and 
need many more corrections than 
or 
.
Independently to their size, another source of complexity in WDNs is their degree of meshing. Indeed, Piller (1995) showed, for pressure-independent consumptions, that heavily meshed networks need more operations to be solved. Adapting the work of Piller (1995) to the case of pressure-dependent consumptions, it is possible to demonstrate numerically that 
are more meshed than 
, and therefore more complex to solve.
Finally, strong variability of hydraulic resistances and/or demands across 
could also explain why these networks are more complex to solve than 
even if they are smaller. Indeed, the presence in a network of pipes with very low resistances together with pipes of very high resistances, and/or of junctions with very low demands together with junctions with very high demands, can increase strongly the stiffness of the Schur complement of the Jacobian matrix (defined by Equation (23)).
Comparison of the Python simulator with EPANET 2.0
We find that increasing the heads at sources 
by 
and reducing the service pressure-head 
to 
leads to fully satisfied demand at every junction of every network. With such values of 
and 
, the maximal difference between the consumptions computed by the Python simulator and the ones computed by EPANET 2.0 is equal to 
. Thus, our Python simulator gives the same results as EPANET 2.0 when the demands are fully satisfied, which means that our simulator is able to correctly run demand-driven analyses.
, the 
-norms of the residuals at convergence of EPANET 2.0 are always greater than those of the Python simulator (Figure 6(a)). Thus, our simulator is more accurate than EPANET 2.0 for these test cases. We believe that this difference is due to EPANET 2.0’s convergence criterion, which considers only the differences between iterates of flow rates while the criterion of our Python simulator considers the differences between iterates of both flow rates and heads. 
Infinity norm of residuals (Figure 6(a)) and simulator elapsed CPU time (Figure 6(b)), simulating each network 
with EPANET 2.0 and the Python simulator. (a) Infinity norms of residuals (in 
or mH
O); (b) Simulator CPU times (in s).
Infinity norm of residuals (Figure 6(a)) and simulator elapsed CPU time (Figure 6(b)), simulating each network with EPANET 2.0 and the Python simulator. (a) Infinity norms of residuals (in or mHO); (b) Simulator CPU times (in s).
EPANET 2.0’s convergence criterion has been improved in EPANET 2.2, by considering also the 
-norm of flow changes and the 
-norm of energy residuals (Rossman et al. 2020, p. 116).
As expected, the elapsed CPU times of the Python simulator are slightly greater than the ones of EPANET 2.0 (Figure 6(b)). This is normal because EPANET 2.0’s solution engine is coded entirely with the C programing language, which is compiled and thus more efficient than the Python language. But these differences remain small, maybe because the CHOLMOD (Chen et al. 2008) library used by the Python simulator to compute the Cholesky decomposition of the Schur complement matrix is more optimized than the method used in EPANET 2.0, which is based on an older work from George & Liu (1981); however, a rigorous comparison of these two methods would require some further work.
CONCLUSIONS
In this paper, we designed and implemented a new Python simulator to model pressure-dependent users’ consumptions in WDNs with no need of artificial elements. Our simulator produces the same results as those computed by the state-of-the-art MATLAB simulator developed by Elhay et al. (2016). Extensive numerical experiments on 
complex, large and realistic networks, composed of up to 19,647 pipes and 17,986 nodes, showed similar precision and efficiency for both simulators. Moreover, the performances of our simulator are close to those of EPANET 2.0 in case of fully satisfied demand.
The community of WDN modelers now has at its disposal a new efficient and flexible simulator of pressure-dependent users’ consumption, integrated into the generic and collaborative framework OOPNET (Steffelbauer & Fuchs-Hanusch 2015). This is the first Python simulator that includes all numerical enhancements from Piller (1995), Piller et al. (2003), Elhay & Simpson (2011) and Elhay et al. (2016). These enhancements were mandatory for the simulation of all tested configurations to converge. Among these enhancements, the damping correction will allow the study of more complex models, especially those for which the gradient of the objective function to minimize can have sub-linear growth.
A complete set of tests were driven to assess the good functioning, stability and robustness of our simulator. Also, the native properties of the Python language make code extension easier. Therefore, our Python simulator can now be reused to study trickier physical processes, such as valves, pumps and pressure-dependent background leakage outflows, while limiting the risk of software regression.
An interesting further work could also consist in comparing the results and performances of our Python simulator with those of EPANET 2.2 (which was not available yet at the beginning of this study) when running pressure-driven analyses. In particular, since the decomposition of the Schur complement matrix is a bottleneck of both our simulator and EPANET 2.2’s solution algorithms, we would like to compare the efficiency of the CHOLMOD library (Chen et al. 2008) used in our simulator, with that of the Cholesky decomposition from George & Liu (1981) and implemented in EPANET 2.2. The comparison of these two implementations, both written in C programing language, should not present any difficulty since our simulator has a modular architecture that permits us to change easily the solver and library to use.
AUTHOR CONTRIBUTIONS
C. C. was involved in conceptualization, formal analysis, investigation, methodology, software, validation, visualization, writing the original draft, writing the review, and editing.
O. P. was involved in conceptualization, formal analysis, methodology, supervision, validation, writing the review and editing.
I. M. was involved in conceptualization, supervision, writing , revising and editing.
FUNDING
This work is part of the Oriented Renewal of Pipes (ROC) French research project and was supported by the Nouvelle Aquitaine region, the Loire-Bretagne and Adour-Garonne water agencies, the Aquitaine regional public health authorities (ARS), and the European Regional Development Fund (ERDF).
DATA AVAILABILITY STATEMENT
Please refer to the online version of this article for figures in color. Additional relevant data to the manuscript are included in the Supplementary Materials. Source code of the Python simulator is publicly accessible on https://github.com/oopnet/oopnet/tree/simulators.
CONFLICT OF INTEREST
The authors declare there is no conflict.
Source code from https://www.epa.gov/sites/default/files/2018-10/en2source.zip.
%20cropped.png?versionId=7142)