Pressure-Driven Simulation of Water Distribution Networks: Searching for Numerical Stability

: EPANET uses a demand-driven approach to compute pressures and flows in the water distribution system. The demand-driven approach (DDA) assumes that the required demand is always fully satisfied no matter the existing pressure. In scenarios of pressure-deficient conditions the DDA results are not accurate, and a pressure-driven approach (PDA) is needed. Frequently, the PDA is accomplished by using equations that compute the available demand/leakage as a function of the current pressure. However, embedding such equations into the solver introduces convergence problems. This paper details the actions taken in WaterNetGen—an EPANET extension—to bring numerical stability to the pressure-driven solver, namely, by smoothing the pressure–demand/leakage relationship and the pump curve.


Introduction
EPANET [1] is a public domain Windows software developed by the U.S. Environmental Protection Agency (USEPA). EPANET performs extended period simulations of the hydraulic and water quality behaviour of drinking water distribution systems. EPANET uses a link-node formulation, where consumption along the pipes is assigned to their end-nodes and defined as demand. To compute the nodal heads and the link flows, EPANET iteratively solves a linear system of equations until some convergence criterion is achieved-see Figure 1 and Listing 1a.
A linear system of equations results from the linearization of the mass conservation and energy equations by the Newton-Raphson method [2]. The iterative computation of Figure 1 stops when some level of accuracy is reached. Thus, the heads (and flows) are computed to satisfy the desired demand, which can lead to solutions with negative pressures. In other words, the demand is always satisfied independently of the available nodal pressure, even in the presence of negative pressures-a physical impossibility. Due to this behaviour, it is said that EPANET performs a Demand-Driven Analysis (DDA).
A better approach is to compute the demand as a function of the current pressure-Pressure-Driven Analysis (PDA). In this case, the iterative process must include the "Update demands" block from Figure 1-see Listing 1b. The computation of the demand inside the loop frequently introduces numerical instability, which requires some corrective actions. In [3], the authors used a pressure-demand/leakage relationship to compute the available demand/leakage and employ relaxation coefficients to deal with these convergence difficulties. However, the authors concluded that in some cases the relaxation coefficients are not enough to guarantee convergence. This paper proposes a smoothing technique to build a differentiable pressure-demand/leakage relationship to improve the convergence of the iterative process.
The next section briefly introduces the Global Gradient Algorithm [2], its pressure-driven version [4] and the proposed smoothing technique to be applied to demands, leakages and pump curves computation. Section 3 includes some results. Section 4 concludes the paper and discusses some future developments.

Materials and Methods
At the core of the EPANET solver is the Global Gradient Algorithm (GGA). Section 2.1 introduces the EPANET solver, Section 2.2 describes the pressure-dependent changes to the GGA and Section 2.3 describes the smoothing technique proposed.

EPANET Solver
Considering a hydraulic network composed of np pipes with unknown flow rates, nn nodes with unknown heads and n0 nodes with known heads, the GGA states:  (1), A11 is a np × np diagonal matrix whose elements correspond to the pipe head losses; A12 = A21 T and A10 = A01 T are the topological incidence matrices of size np × nn and np × n0, respectively.
Assuming a head loss formula as follows: the A11 matrix takes the form Considering D11 as the derivative of A11 with respect to pipe flows, the iterative formulation to find a solution of (1) can be stated as follows: The computation of the new head, H k+1 , corresponds to solving a system of linear equations Ax = b, with A = A k , x = H k+1 and b = F k , where the coefficients matrix A is sparse, symmetric and positive definite-see Figure 2. EPANET solves this linear system through a direct method (LL T -Cholesky factorization) customised to sparse matrices. The EPANET implementation of the linear system solver first applies an ordering (Minimum Degree) algorithm to reduce the fill-in problem linked to the factorization of sparse matrices [6]. Figure 2 shows the matrix sparsity pattern that results from the ordering phase (and before numerical factorization) of the benchmark network Wolf-Cordera Ranch.

Pressure-Dependent Demand
There are two types of pressure dependent demands to be considered: the demand itself and pipe leakage. The pressure-dependent demand for a node i can be computed by the following relationship [7]: where Pi ref is the reference (or service) pressure necessary to fully satisfy the required demand qi req , Pi min is the pressure below which no water can be supplied, Pi is the current pressure at node i and DSR (Demand Satisfaction Ratio) is the ratio between the available and the required demands defined by ?
The exponent of the pressure-demand relationship, α, typically assumes the value 0.5, but should be set by calibration. Figure 3 illustrates the DSR for different values of the exponent and considering P min = 3 m and P ref = 14 m. Leakage also depends on the pressure. In fact, leakage increases when pressure rises and decreases when pressure falls. The pipe leakage can be expressed by the sum of two components: background leakage and bursts leakage. The pressure-leakage relationship can be set as follows [8]: where qk leak is the total leakage flow along pipe k; lk is the length of pipe k; αk and βk are parameters of the background leakage model; Ck and δk are parameters of the bursts leakage model (classical orifice flow formulas); and Pk is the average pressure in pipe k computed as the arithmetic mean of the pressure values of its end nodes. According to Lambert [9], the αk parameter can take values between 0.5 and 2.5, depending on the leak type. The βk parameter is related to the pipe material deterioration, and its value must be set by calibration (initial values can be set around 10 −7 ). The iterative formulation for the PDA is as follows: where DL22 and D22 are the derivatives of qleak and qavl functions, Equations (5)- (7), with respect to pipe pressures and nodal pressures, respectively, and k is the current step. WaterNetGen [4,10] uses relaxation coefficients to improve PDA convergence.

Function Smoothing
Despite the incorporation of relaxation coefficients in the pressure-dependent formulation, the PDA still evidences problems in some cases [3]. One source of such problems is the inclusion of the derivative terms of Equations (5)-(7) in the iterative process. For example, considering expression (6), when the current pressure approaches the minimum required pressure (P min ), the derivative tends to infinity. These derivative terms should be smoothed.
One way of dealing with ill-conditioned behaviour derivative terms is to cut part of the pressure-demand relationship extremity. Figure 4 illustrates the idea for the DSR (6): the extremity is replaced by another well-behaved curve, chosen to be differentiable. The logistic sigmoid function (Equation (9)) is a good candidate to this operation, because it is differentiable and can be parametrised. Thus, the logistic function can be implanted at an appropriate location of the original function, for example, at a point X0 corresponding to a DSR of 20%.
The value of K can be computed to guarantee smooth derivatives too: Figure 5 shows the graph of derivatives of the original and the smoothed functions and yellow and red lines, respectively. The above function smoothing technique is a little elaborate to be applied to pressure-leakage functions. Thus, following the same idea of function smoothing, a pragmatic approach results: use the cubic Hermite interpolation [11]. Cubic Hermite curves provide continuity of first-order derivatives and match the set of (chosen) slopes. They are simple to calculate and only need two control points and two control tangents. Accordingly, two points can be considered x1 = P min and x2 = x1 + Δ, with function values f1 = DSR(x1, P min , P max , 0.5) = 0.0 and f2 = DSR(x2, P min , P max , 0.5). The Δ can be selected to produce a good smooth derivative curve ( Figure 5). The corresponding tangent points, df1 and df2, can be set to a small value (for example, 1 × 10 −10 ) and to the derivative of DSR on x2, respectively. The curve coefficients (a, b, c, d) can be obtained by solving the following linear system: which is equivalent to computing In WaterNetGen, the cubic Hermite curve, with the necessary adjustments, is applied to both extremities of the DSR function, to leakage (background and burst) components and to the pump curves.

Results
WaterNetGen is an EPANET extension that allows both demand-driven and pressure-driven simulations. The PDA is achieved with pressure-demand/leakage relationships (Equations (6)-(8)) and with relaxation coefficients. However, in a previous work [3] the authors found that in certain cases the process did not converge. These cases are related to valve interactions and to problems with the derivative terms. Figure 6 shows the C-Town network, used in the Battle of Background Leakage Assessment for Water Networks (BBLAWN) [12], where WaterNetGen was used to define a pump scheduling in order to minimise the operational costs (energy) and background leakage. In some runs WaterNetGen revealed some difficulties and did not converge.
The tool with the proposed smoothing technique is running again and almost does a decent job. However, there remain some cases where the process does not converge. That is, although there are improvements in the number of successful runs, there are still cases that deserve more research. One final note about the research. The EPANET 2.2 [13], released on December 2019, comes with PDA capability based on Equations (5) and (6). This new version was used to simulate the behaviour of the BBLAWN network, but the result is always "WARNING: System unbalanced at…". Therefore, this new version of EPANET still needs some reworking to deal with more complex networks.

Discussion and Conclusions
EPANET is the outstanding tool to study water distribution networks. However, EPANET is built around the demand-driven approach, which is not suitable to simulate pressure-deficient scenarios or pipe leakage. For these cases, the demand cannot be considered as constant known value, but instead must be considered as pressure dependent. A common way of specifying pressuredependent demand is to consider a pressure-demand relationship and use it in the (iterative) simulation process.
The embedding of the pressure-demand relationship into the iterative process introduces instabilities in the solver. In the past, the authors tried to circumvent these instabilities with the inclusion of relaxation coefficients, which revealed to be essential to the convergence of the solver. Nevertheless, there persist some cases in which the solver does not converge.
The smoothing technique proposed in this work allowed one to overcome some more cases of lack of convergence. However, it is not a finished subject; there remain some cases not solved, so more research is needed.
In the future, the authors plan to add a backtracking algorithm, to assure global convergence for the Newton-Raphson method, and exploit it in combination with the techniques here introduced.
Author Contributions: The authors contributed equally to this work. All authors have read and agreed to the published version of the manuscript.
Funding: This research was partially funded by Portuguese Foundation for Science and Technology under project UIDB/00308/2020.

Conflicts of Interest:
The authors declare no conflict of interest.