Speed factors computed for pumping schedules in Water Distribution Networks: DDA versus PDA formulations

. This paper focuses on a methodology allowing to derive the pumping schedule in Water Distribution Networks (WDN), upon a time dependent water demand. The selected test case is a previously studied WDN. Two pumping algorithms give different pumping rules. By solving the nonlinear system of equations, consisting of energy balance equations, mass balance equations and pumping rules, one gets the pumps speed factors. Solutions attached to the Pressure Driven Analysis (PDA) correspond to energy and cost savings, with respect to the solutions given by the Demand Driven Analysis (DDA). The methodology described in this paper is simple and rapid, but the iterative numerical method used to solve the system of equations is highly dependent on the starting guess.


Introduction
Pumping schedules attached to variable speed driven pumps in Water Distribution Networks (WDN) are a major concern for water companies, because appropriate/ personalized schedules allow reducing significantly the overall energy consumption [1,2].
With respect to nodal demands [8], the hydraulic analysis can be mainly performed as Demand Driven Analysis − DDA (where outflows are pressure independent, being set as fixed boundary conditions at the consumers), or as Pressure Driven Analysis − PDA (where outflows are pressure dependent boundary conditions, varying upon the available pressure at each consumer [9,10]). Hybrid simulation methods combining DDA and PDA also can be used [6].
Among the above listed specialised software, only EPANET lacks the PDA formulation, but it incorporates emitters, which allow modelling pressure-dependent outflows at the consumers [11].
All listed software allow modelling variable speed driven pumps, through their speed factor  (or relative speed), defined as ratio between the actual pump speed and its nominal speed. Speed factor values can be inserted into the model as fixed values over time, or through time dependent speed patterns [3] issued from computations performed outside the current software [1,12], or through rule-based controls [3], e.g. where pump's speed factor is adjusted in imposed steps upon the pressure level available at a certain node [2,13]. Within specialised software, speed factor values are not computed as solver outputs, excepting WDNetXL.
In this paper, we present a methodology to derive the pumping schedule in WDNs (to compute pumps speed factors), for both formulations: DDA and PDA. The selected test case is a simple WDN fed by variable speed driven pumps, previously studied [12]. PDA and DDA results are finally compared. DDA results are also discussed with respect to similar results obtained by Georgescu [12], using another DDA-based methodology.

Speed factors computing method
In order to simplify the presentation, the proposed methodology is developed for a test case, consisting of the WDN from figure 1, proposed in an open access paper by Georgescu [12]. The main components are tagged with ID labels, as in figure 1, with ID = 118 for nodes (where ID = 1 denotes the water source − the reservoir R), ID = 120 for pipes and ID = 2123 for three variable speed driven pumps. The pumping station equipped with those pumps is bordered by the main suction node with ID = 2, and the main discharge node with ID = 3. There are 5 consumers (end-users) within the studied looped network, labelled with ID  {11; 12; 14; 16; 18} and marked by red dots in figure 1.
The hydraulic analysis is an Extended Period Simulation (EPS) analysis [14], conducted over a period of 24 hours, with a time step of one hour, starting at midnight; the time t is the starting moment of each time step (thus, t = 1 starts at mid-night and t = 24 at 11 p.m.). During a time step interval, the system behaviour fits a steady-state regime. The EPS model evaluates successively all 24 steadystate intervals, considering the changes from one time step to the next one, through time dependent boundary conditions (here, the requested water demands).
Pipes geometry (diameter D, length L and equivalent roughness), the water demand pattern coefficients cq(t) and the energy price pattern coefficients ce(t), dependent on time t over the above defined 24 hours period, are fully described in Georgescu [12].
Each consumer requests a daily mean water demand Qrm of 0.088 m 3 /s, so the requested hourly demand is computed as: Qr(t) = cq(t)Qrm. The hourly energy price Ch(t) is computed by multiplying the coefficients ce(t) with a daily mean price of 0.065 €/kWh [15].
The hydraulic system is flat, e.g. with nodes elevation z = 0. The reservoir's free surface is open to the atmosphere, so there the head is HR = 0. The reference gauge pressure pref, requested at each consumer to ensure the demand, is pref = 40 meters of water column (i.e. reference pressure head Href = 40 m). The minimal gauge pressure pmin, for which there is outflow at the consumer, is set here as pmin = 0 mWC (i.e. minimal pressure head Hmin = 0 m); usually, pmin = 23 mWC.
The Darcy-Weisbach formula is used to compute head losses h (in meters) on pipes: where Q is the flow rate (in m 3 /s), R is the hydraulic resistance (in s 2 /m 5 ), g is the gravity (g = 9.81 m/s 2 ), and  is the friction factor, defined by the Swamee and Jain formula [3,6]. Inlet and outlet kinetic terms on pipes, and minor losses are neglected. At each moment t (where t = 124, with hourly time step), the centrifugal pumps (denoted by indexes i corresponding to ID-s from figure 1) operate at a certain speed i(t) upon the following characteristic curves [12], namely the head − flow rate curve, with pumping head H in meters, for the flow rate Q in m 3 /s: (2) and the efficiency curve, with efficiency values  in percents, for Q in m 3 /s: Within DDA, the available demand Qa(t) at each consumer equals the requested demand Qr(t) at each time step. Within the PDA formulation, Qa(t) is computed with the following time and pressure dependent demand relationship [13]: where p(t), in mWC, is the available pressure at the consumer, at time t.
For the studied WDN fed by 3 working pumps, the nonlinear system of equations, consisting of energy balance equations on all 23 links and continuity equations in 17 nodes (with ID from 2 to 18), yields a total of 40 equations with 43 unknowns, namely: 23 flow rate values (Q1 to Q23), 17 nodal head values (H2 to H18, more precisely pressure head values for that flat system) and 3 speed factor values (21, 22 and 23). The resulting underdetermined system of equations corresponds to a stiff problem.
Nonlinear system of equations can be solved in GNU Octave [16] or MATLAB, using the fsolve function [17], developed for minimization and zero finding through algorithms based on model-trust region optimization approach [18]; according to fsolve documentation [17], the system of equations cannot be underdetermined. So, the studied hydraulic problem needs 3 additional equations.
One can add 2 equations related to speed factors settings, depending on the pumping station operating algorithm. Two options will be considered further [6,12], namely:  a classical pumping algorithm, where at a time t, a single pump runs at variable speed with (t)  1, while any other opened pump runs at its nominal speed, where the speed factor equals 1.
 an improved pumping algorithm, where at a time t, all opened pumps run with equal speed factor values.
The missing third equation, needed to obtain a system of 43 equations with 43 unknowns, can be the global continuity equation, where the inflow Q1 is the sum of the five available outflows, defined by (4): The global mass balance equation (5) is not linearly independent with respect to the other 17 continuity equations, but it helps to mimic a determined system of equations, which can be solved iteratively using the fsolve function, starting from an initial guess of the solution. As we will show further within Section 3, that pseudo-determined nonlinear system of equations is highly dependent on the starting guess: the algorithm converges only for a certain initial guess, while it fails for a slightly changed initial guess. The same hydraulic problem was solved in Georgescu [12], within a twosteps methodology and a DDA formulation: there, the fsolve function was used to solve a properly-determined nonlinear system (with 6 equations and 6 unknowns), and the same algorithm proved to be globally convergent.
The nonlinear system of equations derived in this paper is too long to be written explicitly here, but it can be presented conveniently as below, mentioning that the first 23 equations (f1f23) are energy balance equations on links (on pipes with ID = 14, on pumps with ID = 2123, then on pipes with ID = 520); the next 17 equations (f24f40) are continuity equations in nodes 218; equation f41 results from (5); the last 2 equations refer to speed factors (which differ for the classic pumping algorithm, with respect to the improved algorithm); in order to simplify notations, the time dependency was not mentioned here: For the considered pumping station, we adopt further the following pumping rules [12]: while pump2 and pump3 are closed ( 0  (6), the list of unknowns contains 3 type of variables: flow rates, nodal heads and speed factors; as already mentioned, HR=0. To obtain a compact form of the system, to be solved numerically, the unknowns become the components wk (with k = 143) of a column vector w, as: 23  22  21  43  42  41   18  3  2  40  25  24   23  2  1  23  2  1   ,  ,  , Thus, the final nonlinear system of equations can be defined as f(w) = 0, where the column vector f contains 43 components, namely the functions fk(w1,w2,...,w43) = 0, with k = 143. Taking into account (4), with reference pressure head Href = 40 m, and minimal pressure head Hmin = 0, the system (6) can be rewritten as: Obviously, the solution of the system (8) corresponds to the pumping rule no. 4. The system (8) simplifies to: • 37 equations with 37 unknowns for the pumping rules no. 2 and 3; • 31 equations with 31 unknowns for pumping rule no. 1.
Finally, based on the parameters attached to the duty points of each pump: flow rate Qi(t), pumping head Hi(t), and pump efficiency i(t) defined by (3) in percents, the power Pi(t) of each pump can be computed (in watts) at every time step, as: where  is the water density ( = 1000 kg/m 3 ).
Starting from (9) with t = 124, one can compute the hourly energy consumption Ei(t) for each pump, as well as the daily energy (denoted E) consumed by the pumping station. Considering the mentioned hourly energy price Ch(t), one can compute, finally, the total cost C of the daily consumed energy.

Results and discussions
For all time steps, based on previously defined pumping rules, the solution w of the corresponding nonlinear system of equations f(w) = 0, written in (8), proved to be highly dependent on the starting guess (the column vector w0).
We present in tables 1 and 2, all speed factor values {21(t), 22(t), 23(t)}, where t = 124, computed for the classical pumping algorithm, as well as for the improved pumping algorithm, for both formulations: DDA and PDA. As one can notice in table 3, the available nodal heads at the consumers with ID = {11; 12; 14; 16} are greater than the reference head (Href = 40 m), meaning that the pressure − demand relationship (4) fits the DDA case, so for those 4 consumers, the available demand Qa(t) equals the requested demand Qr(t). At the last consumer, with ID = 18, there is a pressure deficit, the available nodal head being less than 40 m (here, a pressure drop of 0.75 mWC is obtained for the classical algorithm, and a drop of 0.45 mWC is obtained for the improved algorithm). For that last consumer, the relationship (4) gives Qa18(t) < Qr(t). The mass balance (5) compared to the requested total demand, i.e. 5Qr(t), shows a small water demand shortage (the relative error with respect to the requested total demand being of -0.19% for the classical algorithm, and -0.11% for the improved one).
The above behaviour of the computed results and the above comments apply for all 24 time steps, and for both pumping algorithms, within the PDA formulation: thus, at any time t, the last consumer experiences a pressure deficit and a slight shortage in available demand: the maximum pressure drop (1.29 mWC) corresponds to the classical algorithm, for cq(t) = 1.32, at t = {8; 10; 14; 15; 17; 19}, where the relative error with respect to the requested total demand is of -0.33%; the minimal pressure drop (0.002 mWC) is negligible and corresponds to the improved algorithm, for cq(t) = 1.2, at t = 20. We present in table 4, the computed values of the daily energy E (in kWh), consumed by the pumping station, and its total cost C (in €), for both DDA and PDA formulations, within each pumping algorithm: the classical one and the improved one. For comparison, we also present the results obtained by Georgescu [12], for the same WDN, with another methodology within the DDA formulation. The speed factor values from tables 1 and 2 attached to DDA, can be compared to the corresponding values from Georgescu [12]. Data synthesised in table 4 reveal the following: • solutions attached to PDA within both pumping algorithms ensure energy and cost savings with respect to solutions given by DDA in this paper or previously [12]; • from the energy consumption and cost points of view, the best solution is the one attached to PDA within the improved pumping algorithm − the same conclusion was drawn in previous studies on the same pumping algorithm [6,7,12]; • comments can arise for the fact that PDA formulation gave a pressure deficit (of less than 0.1 bar) and a slight demand shortage (of about 1 litre/second) at the most disadvantaged consumer (out of 5 consumers): it is true, but PDA is the correct formulation, while DDA is not (in real WDNs, consumers do not receive the requested demand at any available pressure). Moreover, the overpressure at the upstream consumers is smaller in PDA (less than 6 mWC for cq(t) = 1.4 in table 3) than in DDA (where it reaches about 8 mWC for cq(t) = 1.4 for both pumping algorithms); • solutions obtained with DDA in this paper, with both pumping algorithms, gave greater energy and cost values than the corresponding solutions obtained by Georgescu [12] (relative errors vary from 1.67% to 3.88%). The reason relies on methodologies: within this paper, the solution is computed quickly, through the nonlinear system of equations (8), but convergence depends strongly on the starting guess (for the same time step, different starting values yield different solutions, so one must choose the best solution; for some time steps, the selected computed results are convenient, while for other time steps, the selected results deviate from the expected values). The DDA methodology proposed previously [12] is quite accurate, but time consuming, since it involves two steps [6]: firstly, a split of the initial WDN into 2 separate hydraulic networks (meaning a change in topology), allowing to compute conveniently the unknown speed factor values and to establish speed patterns for the next step); the second step consists in restoration of the initial WDN configuration and its modelling in EPANET, based on predefined speed patterns (thus, speed factors become known input data in EPANET).

Conclusions
This paper focuses on a methodology allowing to derive the pumping schedule in Water Distribution Networks (WDN), upon a time dependent demand pattern. The selected test case is a simple WDN fed by variable speed driven pumps, previously studied [12]. Two pumping algorithms, namely a classical one and an improved one, allow selecting different pumping rules. Energy balance equations, mass balance equations and additional pumping rules form a nonlinear system of equations, which allows computing pumps speed factors. The hydraulic analysis is conducted over a period of 24 hours, with a hourly time step, within two different approaches, namely the Demand Driven Analysis (DDA) and the Pressure Driven Analysis (PDA).
The results obtained previously [12], using another DDA-based methodology, proved to be better than the results computed here with the DDA formulation. But differences are small (the maximum relative error being less than 3.9%), so the present methodology can prevail, being more rapid. Solutions attached to PDA correspond to greater savings in terms of energy and cost than any other solution given by DDA here, or previously [12].
The methodology described in this paper is simple and rapid, but the iterative numerical method used to solve the system of equations proved to be highly dependent on the starting guess. Further work will focus on other numerical method(s), appropriate to solve the nonlinear system of equations, if possible, in a globally convergent manner.