- Email: [email protected]
A software tool for assessing the performance of and implementing water distribution system solution methods
Accepted Manuscript A software tool for assessing the performance of and implementing water distribution system solution methods Mengning Qiu, Bradley Alexander Angus R, Simpson Sylvan Elhay PII:
S1364-8152(18)30416-X
DOI:
https://doi.org/10.1016/j.envsoft.2018.10.016
Reference:
ENSO 4335
To appear in:
Environmental Modelling and Software
Received Date: 2 May 2018 Revised Date:
8 October 2018
Accepted Date: 27 October 2018
Please cite this article as: Qiu, M., Angus R, B.A., Elhay, S.S., A software tool for assessing the performance of and implementing water distribution system solution methods, Environmental Modelling and Software (2018), doi: https://doi.org/10.1016/j.envsoft.2018.10.016. This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT
A Software Tool for Assessing the Performance of and Implementing Water Distribution System Solution Methods Mengning Qiu∗
RI PT
Bradley Alexander Angus R. Simpson Sylvan Elhay University of Adelaide, South Australia, 5005
SC
Highlights
M AN U
• A library for the steady-state analysis of a water distribution system (WDS) • An open-source C++ software implementation of a number of WDS solution methods. • A fast simulation platform for both once-off and multi-run simulations • Several improvements to the existing solution methods have been made.
TE D
Abstract
WDSLib is an extensible simulation toolkit for the steady-state analysis of a water distribution system. It includes a range of solution methods: the forest-core partitioning algorithm, the global gradient algorithm, the reformulated co-tree flow method, and also combinations of these methods.
EP
WDSLib has been created using a modularized object-oriented design and implemented in the C++ programming language, and has been validated against a reference MATLAB implementation. WDSLib has been designed: (i) to avoid unnecessary computations by hoisting each of the modules to
AC C
its appropriate level of repetition, (ii) to perform the computations independently of measurement units using scaled variables, (iii) to accurately report the execution time of all the modules in that it is possible to produce a timing model to parameterize multiple simulation times (such as in an optimization using a genetic algorithm) from a series of sampling simulation runs and (iv) to guard against numerical failures. Two example applications, a once-off simulation and a network optimization design application simulation, are presented. This toolkit can be used (i) to implement, test and compare different solution methods, (ii) to focus the research on the most time-consuming parts of a solution method and (iii) to guide the choice of solution method when multiple simulation runs are required. ∗ Corresponding Author Email: [email protected]
1
ACCEPTED MANUSCRIPT Keywords: Water Distribution System; C++ toolkit; Object-Oriented design; Forest-Core Partitioning Algorithm; Reformulated Co-tree Flows Method; Global Gradient Algorithm; open source software
Software availability
2
Name of the Software: WDSLib
3
Version: 1.0
4
Available from: https://github.com/a1184182/WDSLib
5
Language: C++
6
Supported System: Windows, MacOS, Linux, Unix
7
Year first available: 2018
8
1 Introduction
9
Hydraulic simulation has been used to model water distribution systems (WDSs) for several decades
10
and is an essential tool for the design, operation, and management of WDSs in industry and research.
11
Hydraulic simulation allows users (1) to optimize WDS network parameters, such as pipe diameters,
12
in a design setting, (2) to calibrate network parameters, such as demand patterns, in a conventional
13
operational setting, (3) to conduct real-time monitoring and calibration of the network elements in a
14
supervisory control and data acquisition (SCADA) operational setting, and (4) to adjust control devices,
15
such as valves, in a management setting. In the design setting and both the above operational settings,
16
repeated hydraulic assessment is required on a network with fixed topology. In the management setting,
17
repeated hydraulic assessment is required on a network with flexible network parameter settings. With
18
ever-increasing network sizes and the need for real-time management using a SCADA system, it is
19
important to have a robust simulation package which can be configured to be maximally efficient whatever
20
the setting.
AC C
EP
TE D
M AN U
SC
RI PT
1
21
In the field of hydraulic simulation, the system of equations can be formulated as a large and sparse
22
non-linear saddle point problem. Saddle point problems are defined to be a function of two variables (say
23
f(q,h)) and the objective is to find a pair of variables (q*,h*) such that the value f(q*,h*) is minimized
24
with respect to the first variable and maximized with respect to the second variable or equivalently,
25
f (q∗, h) < f (q∗, h∗) < f (q, h∗)∀q, h. There are several well-known iteration methods for solving the non-
26
linear saddle point problem. These include: range space methods, a method that operates in the subspace
27
defined by the rows of the unknown-head node-arc incidences matrix, (Global Gradient Algorithm (Todini
28
and Pilati 1988)), Null space methods, a method that operates in the subspace defined by a null-space that
29
is orthogonal to the column of the unknown-head node-arc incidences matrix (Co-Tree flow formulation 2
ACCEPTED MANUSCRIPT variations (Rahal 1995; Elhay et al. 2014)), loop-based methods, a method is based on the sum of head
31
losses around a loop must be zero (Loop flow correction (Cross 1936)), and pre-and-post-processing
32
methods (forest-core partitioning algorithm (Simpson et al. 2014), domain decomposition (Diao et al.
33
2014), network clustering (Perelman and Ostfeld 2011)). Their relative performance in terms of speed,
34
rate-of-convergence, and accuracy depends among other things on the topology of the target network:
35
size of the forest component, the number of network loops, and the density of these network loops.
36
It is difficult to evaluate the impact of these topology factors by only examining the incidence matrix
37
that describes the pipe network connectivity. As a result, the best method to use for a particular
38
network cannot be easily determined a priori. Moreover, extra complexity is introduced when a multi-
39
run hydraulic assessment is required. During a multi-run hydraulic simulation, the elapsed computation
40
time of each method can be broken down into two parts: the components that are only required to be
41
performed once at the very beginning for the same network, called the overhead, and the components that
42
are required to be carried out repeatedly for each separate run until the required number of iterations has
43
been met, called the hydraulic-phase. It is desirable to have a simulation platform, given the different
44
levels of repetition, to implement these alternative algorithms efficiently. Equipped with such a platform
45
a user would be able to easily benchmark the performance of alternative methods on a small number
46
of evaluations for a given network and use that performance to inform the choice of algorithm to use
47
for either a once-off simulation setting or for a multiple simulation setting (such as for an evolutionary
48
algorithm (EA)).
TE D
M AN U
SC
RI PT
30
This work describes an extensible WDS simulation platform called WDSLib. WDSLib is a numerically
50
robust, efficient and accurate C++ library that implements many WDS simulation methods. WDSLib is
51
written using a modular object-oriented design which allows users to easily mix and interchange solution
52
components, thereby enabling users to avoid redundant computations. It has been optimized to use
53
sparse data structures which are oriented to the pattern of access required for each solution method.
54
WDSLib has been validated for accuracy on a range of realistic benchmark water distribution networks
55
against reference implementations and tested for speed. The program accepts the input file formats of
56
the industry standard EPANET2 (Rossman 2000) toolkit and its performance is faster than EPANET2
57
in all tested settings and benchmarks.
AC C
EP
49
58
The remainder of this paper is structured as follows. The next section describes related methodologies
59
and implementations. A general description of the WDS demand-driven steady-state problem is given
60
in the next section. Section 3 presents a mathematical formulation of the network and the solution
61
methods that are used in WDSLib. The tool-kit structure is then given in section 4. This is followed,
62
in section 5, by the toolkit implementation details. Section 6 provides some examples of how the toolkit
63
can be utilized in a simulation work flow. The results are discussed in Section 7. Finally, section 8
64
summarizes the results of this paper and describes future extensions to the toolkit.
3
ACCEPTED MANUSCRIPT
2 Background
66
This section describes related water distribution system network solution methods and implementations.
67
The first sub-section describes solution methods, including those used by WDSLib. This is followed by
68
a description of currently available implementations and compares these with WDSLib.
69
2.1 Related Methods
70
This research considers a water distribution model made up of energy conservation equations and the
71
demand driven model continuity equations. The Hardy Cross method (Cross 1936), also known as the
72
loop flow corrections method, is one of the oldest methods and uses successive approximations, solving for
73
each loop flow correction independently. It is a method that was widely used for its simplicity at the time
74
when it was introduced. More than three decades later, Epp and Fowler (1970) developed a computer
75
version of Cross’s method and replaced the numerical solver with the Newton method, which solves for
76
all loop flow corrections simultaneously. However, this method has not been widely used because of the
77
need (i) to identify the network loops, (ii) to find initial flows that satisfy continuity and (iii) to use
78
pseudo-loops.
M AN U
SC
RI PT
65
The GGA is a range space method that solves for both flows and heads. It was the first algorithm, in
80
the field of hydraulics, to exploit the block structure of the Jacobian matrix to reduce the size of the key
81
matrix in the linearization of the Newton method. The GGA has gained popularity through its rapid
82
convergence rate for a wide range of starting values. This is the result of using the Newton method on
83
an optimizations problem that has a quadratic surface. However, it was reported by Elhay and Simpson
84
(2011) that the GGA fails catastrophically in the presence of zero flows in a WDS when the head loss
85
is modeled by the Hazen-Williams formula. Regularization methods have been proposed by both Elhay
86
and Simpson (2011) and Gorev et al. (2012) to deal with zero flows when the head loss is modeled by
87
the Hazen-Williams formula.
EP
TE D
79
The GGA as it was first proposed, applied only for the WDSs in which the head loss is modeled by the
89
Hazen-Williams formula, where the resistance factor was independent of flow. Rossman (2000) extended
90
the GGA to allow the use of the Darcy-Weisbach formula. It has been pointed out in Simpson and Elhay
91
(2010), however, that Rossman incorrectly treated the Darcy-Weisbach resistance factor as independent
92
of the flow. They introduced the correct Jacobian matrix to deal with this. It has been demonstrated
93
that once the correct Jacobian matrix is used, the quadratic convergence rate of the Newton method is
94
restored. Furthermore, Elhay and Simpson (2011) reported that the GGA does not fail in the presence
95
of zero flows when the derivatives of the Darcy-Weisbach Jacobian matrix are correctly computed for
96
laminar flows.
97
AC C
88
The co-trees flow method (CTM) (Rahal 1995) is a null space method that solves for the co-tree
4
ACCEPTED MANUSCRIPT flows and spanning tree flows separately. The CTM, unlike the loop flow corrections method, does
99
not require the initial flows to satisfy continuity. However, it does require: (i) the identification of
100
the associated circulating graph; (ii) the determination of the demands that are to be carried by tree
101
branches; (iii) finding the associated chain of branches closing a circuit for each co-tree chord; (iv)
102
computing pseudo-link head losses. The reformulated co-trees flow method (RCTM) (Elhay et al. 2014)
103
is also a null space method that solves for co-tree flows and spanning trees flows separately. It represents
104
a significant improvement on the CTM by removing requirements (i) to (iv) above. It uses the Schilders’
105
factorization (Schilders 2009) to permute the node-arc incidence matrix into an invertible spanning tree
106
block and a co-tree block. This permutation reduces the dimension of the Schur complement from being
107
equal to the number of junctions (as in the GGA) to approximately the number of loops in the network.
108
Abraham and Stoianov (2015) proposed a novel idea to speed-up the solution process when using a
109
null space method to solve a WDS network. Their idea exploits the fact that a significant proportion of
110
run-time is spent computing the head losses. At the same time, flows within some pipes exhibit negligible
111
changes after a few iterations. As a result, there is no point in wasting computer resources to re-compute
112
the pipe head losses for the pipes that have little or no change in flows. This partial update can be used
113
to economize the computational complexity of the GGA, the RCTM and their variations.
M AN U
SC
RI PT
98
The forest-core partitioning algorithm (FCPA) (Simpson et al. 2014) speeds up the solution process
115
in the case where the network has a significant forest component. This algorithm permutes the system
116
equations to partition the linear component of the problem, which is the forest of the WDS, from the
117
non-linear component, which is the core of the WDS. It can be viewed as a method that simplifies the
118
problem by solving for the flows and the heads in the forest just once instead of at every iteration. The
119
FCPA reduces the number of pipes, number of junctions, and the dimension of the Jacobian matrix in
120
the core by the number of forest pipes (or nodes).
EP
TE D
114
The graph matrix partitioning algorithm(GMPA) (Deuerlein et al. 2015) exploited the linear relation-
122
ships between flows of the internal trees within the core and the flows of the corresponding super-links
123
after the forest of the network has been removed. This was a major breakthrough. The GMPA permutes
124
the node-arc incidence matrix in such a way that all of the nodes with degree two in the core can be
125
treated as a group. By partitioning the network this way, the network can be solved by a global step,
126
which solves for the nodes with degree greater than two (super nodes) and the pipes which connect to
127
them (path chords), and a local step, which solves for the nodes with degree two (interior path nodes)
128
and pipes connected to them (path-tree links).
129
2.2 Related Implementations
130
EPANET 2 (Rossman 2000) is a widely used WDS simulation package. EPANET 2 implemented the
131
GGA to provide a demand-driven steady-state solution of a WDS. The code for EPANET 2 is in the
AC C
121
5
ACCEPTED MANUSCRIPT 132
public domain, allowing many extensions to be developed. Currently available extensions include: the
133
implementation of a pressure-dependent model (Cheung et al. 2005; Morley and Tricarico 2008; Siew and
134
Tanyimboh 2012; Jun and Guoping 2012) and a real-time simulation capability (Vassiljev and Koppel
135
2015). The EPANET 2 implementation is not explicitly designed to necessarily be easy to understand or ac-
137
commodate alternative solution methods (Guidolin et al. 2010). The elements that are used in EPANET
138
2 are stored by the variables that describe their graph properties. For example, (1) junctions, reservoirs,
139
and tanks are stored as a C struct called Node and (2) all valves, pipes, and pumps are stored as a C
140
struct called Link. The abundant use of global variables limits the reusability and the possibility of the
141
thread-safe design (Guidolin et al. 2010).
RI PT
136
Consequently, it is difficult to cleanly incorporate new solution methods into EPANET 2 in a manner
143
that allows a fair comparison of performance between these methods. Moreover, because there are no
144
clearly defined interfaces for the incorporation of third-party code components in EPANET 2, there is
145
no guarantee that independently authored extensions will be easy to combine with each other.
M AN U
SC
142
In the absence of a popular easy-to-modify WDS simulation platform there is currently no straight-
147
forward means for comparing different solution methods. To date, when new solution methods have
148
been developed they have been compared using different research systems, on different platforms with
149
different implementation languages. This leads to difficulty in comparing methods, limits the reusability
150
of code, and creates a barrier for researchers to reproduce and replicate results. To address these is-
151
sues, an extensible framework is required that allows implementation of new methodologies to be easily
152
incorporated without an adverse impact on the performance of the rest of the system.
TE D
146
To this end, a number of attempts have been made to implement an object-oriented wrapper to
154
encapsulate the EPANET 2 solver (openNet (Morley et al. 2000) and OOTEN(van Zyl et al. 2003)).
155
However, these two systems were focused on providing more flexibility in the processing of input to the
156
core EPANET solver. They did not address any issues relating to the solution process. CWSnet, a C++
157
implementation in object-oriented style, was produced by Guidolin et al. (2010) as an alternative to
158
EPANET 2. In CWSnet, more attention has been given to the hydraulic elements of the WDS network.
159
In addition, CWSNet provides a pressure driven model, and takes advantage of the computing power of
160
the computer’s Graphics Processing Unit (GPU). However, in CSWnet the data structures representing
161
the network are specialized to the solution methods that it uses. These data structures are not easily
162
adapted to work efficiently with the different traversal orders, and graph algorithms used by newly
163
developed solution methods. However, CWSnet still uses the same hydraulic solver and the same linear
164
solver techniques implemented in EPANET 2 (Guidolin et al. 2010).
AC C
EP
153
165
To accommodate the deficiencies referred to above, this paper presents a new hydraulic simulation
166
toolkit WDSlib. WDSlib is coded in C++, and incorporates a number of recently published techniques.
6
ACCEPTED MANUSCRIPT This toolkit offers users the ability to: (i) choose from, or modify, different approaches and implementa-
168
tions of different WDS model analyses, and (ii) extend the toolkit to include new developments. These
169
features have been implemented using fast and modularized code. A focus of attention in this research
170
has been program correctness, robustness and code efficiency. The correctness of the toolkit has been
171
validated against a reference MATLAB implementation. The differences between all results (intermedi-
172
ate and final) produced by the C++ toolkit and the MATLAB implementation were shown to be smaller
173
than 10−10 . In the interest of toolkit robustness, special attention has been paid to numerical processes
174
to guard against avoidable failures, such as loss of significance through subtractive cancellation, and nu-
175
merical errors, such as division by zero. The data structures and code libraries in WDSLib are shared and
176
all implementations have been carefully designed to ensure fairness of performance comparisons between
177
algorithms. WDSLib uses a pluggable architecture where solution-methods, and their accompanying
178
pre-processing and post-processing code are easily substituted. In addition, different numerical linear
179
algebra techniques can be incorporated using a well-defined interface. This concludes the discussion of
180
related work. The mathematical formulations of the solution methods used in WDSLib are presented in
181
the next section.
182
3 General WDS Demand-Driven Steady-State Problem
183
This section describes the general WDS demand-driven steady-state problem. The following starts with
184
the basic definitions and notation, followed by the system equations. Finally, the relevant equations are
185
shown for each of the different solution methods that are implemented in WDSLib. All variables are
186
described in the nomenclature section in Appendix E.
187
3.1 Definitions and Notation
188
Consider a water distribution system that contains np pipes, nj junctions, nr fixed head nodes and nf
189
forest pipes and nodes. The j − th pipe of the network can be characterized by its diameter Dj , length
190
Lj , resistance factor rj . The i − th node of the network has two properties: its nodal demand di and its
191
elevation zi .
SC
M AN U
TE D
EP
AC C
192
RI PT
167
Let q = (q1 , q2 , ....qnp )T denote the vector of unknown flows, h = (h1 , h2 , ....hnj )T denote the vector
193
of unknown heads, r = (r1 , r2 , ....rnp )T denote the vector of resistance factors, d = (d1 , d2 , .....dnj )T
194
denote the vector of nodal demands, el = (el1 , el2 ....elnr )T denote the vector of fixed head elevations.
195
The head loss exponent n is assumed to be dependent only on the head loss model: n = 2 for
196
the Darcy-Weisbach head loss model and n = 1.852 for Hazen-Williams head loss model. The head loss
197
within the pipe j, which connects the node i and the node k, is modelled by hi −hk = rj qj |qj |n−1 . Denote
198
by G(q) ∈ Rnp ×np , a diagonal square matrix with element [G]jj = rj |qj |n−1 for j = 1, 2, ....np . Denote
7
ACCEPTED MANUSCRIPT 199
by F (q) ∈ Rnp ×np , a diagonal square matrix where the j-th element on its diagonal [F ]jj =
200
A1 is the full rank, unknown head, node-arc incidence matrix, where [A1 ]ji is used to represent the
201
relationship between pipe j and node i; [A1 ]ji = −1 if pipe j enters node i, [A1 ]ji = 1 if pipe j leaves
202
node i, and [A1 ]ji = 0 if pipe j is not connected to node i. A2 is the fixed-head node-arc incidence matrix,
203
where [A2 ]ji is used to represent the relationship between pipe j and fixed head node i, [A2 ]ji = −1 if
204
pipe j enters fixed head node i, [A2 ]ji = 1 if pipe j leaves fixed head node i, and [A2 ]ji = 0 if pipe j is
205
not connected to fixed head node i.
206
3.2 System of Equations
207
The steady-state flows and heads in the WDS system are modeled by the demand-driven model (DDM)
208
continuity equations (1) and the energy conservation equations (2):
RI PT
SC
(1)
M AN U
−A1 T q − d = O 209
G(q)q − A1 h − A2 el = O, 210
which can be expressed as
where its Jacobian matrix is
F (q) J = −A1 T
214
−A1 O
(4)
This non-linear system is normally solved by the Newton method, in which q (m+1) and h(m+1) are
EP
213
(3)
and it is sometimes referred to as a nonlinear saddle point problem (Benzi et al. 2005).
repeatedly computed from q (m) and h(m) by
AC C
212
(2)
−A1 q A2 el = 0, − d h O
TE D
G(q) −A1 T 211
d dqj [G]jj qj .
F
(m) (m)
q
−A1 T
(m+1)
(m)
−A1 q −q G = − O h(m+1) − h(m) ||q (m+1) −q (m) || ||q (m+1) ||
(m) (m)
q
(m)
− A1 h
− A2 e l , −A1 T q (m) − d
(5)
(m+1) (m) −h || are sufficiently small. and ||h (m+1) ||h ||
215
until the relative differences
216
3.3 Global Gradient Algorithm
217
Todini and Pilati (1988) applied block elimination to Eq. (5) to yield a two-step Hazen-William solver:
218
Eq. (6) for the heads and Eq. (7) for the flows.
n o h(m+1) = U −1 −nd + A1 T [(1 − n)q (m) − G−1 A2 el ] 8
(6)
ACCEPTED MANUSCRIPT 219
where U = A1 T G−1 A1
q (m+1) =
h i V h(m+1) = −d + A1 T F −1 (G − F ) q (m) − A2 el
221
222
(7)
Later, Simpson and Elhay (2010) proposed
RI PT
220
o 1n (n − 1)q (m) + G−1 (A2 el + A1 h) n
where V = A1 T F −1 A1
223
(8)
224
h i q (m+1) = q (m) + F −1 A1 h(m+1) − F −1 Gq (m) − A2 el
225
as the generalized equations that can be applied when the head-loss is modeled by the Hazen-Williams
226
equation or the Darcy-Weisbach equation. The correct Jacobian matrix with the formula for F , when
227
head loss is modeled by Darcy-Weisbach equation, can be found in Simpson and Elhay (2010). They
228
showed that the use of the correct Jacobian matrix restores the quadratic rate of convergence.
M AN U
SC
(9)
It is important to note that the GGA, as it was originally proposed, solves the entire network by a
230
non-linear solver, and this can include some unnecessary computations which can be avoided by exploiting
231
the structural properties of the WDS graph composition. The methods described below exploit these
232
structural properties to potentially improve the speed of the solution process.
233
3.4 Forest-Core Partitioning
234
Associated with a WDS is a graph G = (V, E), where the elements of V are the nodes (vertices) of the
235
graph G and elements of E are the pipes (links) of the graph G. The graph G can be partitioned into
236
smaller subgraphs with special properties. The special properties that are exploited in WDSLib and their
237
formulations are described in this subsection. The concept of partitioning the WDS network was proposed
238
by Deuerlein (2008) in order to simplify the WDS solution process. Simpson et al. (2014) extended the
239
idea of the network partitioning of Deuerlein (2008) and introduced the forest-core partitioning algorithm
240
(FCPA), which partitions the network into a treed component and a looped or core component. The
241
FCPA starts with a searching algorithm which identifies the forest subgraph, Gf = (Vf , Ef ), in which
242
S ∈ Nnf ×np is the permutation matrix which identifies the pipes in the forest, Ef , as distinct from the
243
pipes in the core , Ec , and T ∈ Nnf ×nj is the permutation matrix which identifies the nodes in the
244
forest, Vf , as distinct from the nodes in the core, Vc , as distinct from the core subgraph, Gc = (Vc , Ec ),
245
in which P ∈ Nnpc ×np is the permutation matrix for Ec and C ∈ Nnjc ×nj is the permutation matrix for
246
Vc .
AC C
EP
TE D
229
9
ACCEPTED MANUSCRIPT 247
The flows of the pipes in the forest, Sq, can be found directly from
248
−1 Sq = − T A1 T S T T d.
249
Finally, once the iterative solution process for the core has stopped, the forest heads can be found by
250
solving a linear system:
RI PT
(10)
−1 SA2 el − SGS T Sq + SA1 C T Ch . T h = −SA1 T T
251
(11)
3.5 Reformulated Co-Tree flows method
253
A graph, with or without forest, can be partitioned into two sub-graphs: a spanning tree subgraph and
254
a complementary co-tree subgraph. The reformulated co-tree flow method (RCTM) (Elhay et al. 2014)
255
exploited the relationship between the spanning tree pipes and the co-tree pipes. The RCTM starts with
256
a spanning tree search algorithm which identifies a spanning tree subgraph, Gst = (V, Est ), in which
257
K1 ∈ Nnpst ×np is the permutation matrix that identifies the pipes in the spanning tree, Est , as distinct
258
from the pipes in the co-tree, Ect . R is the permutation matrix for the nodes which traverse the same
259
sequence as the corresponding spanning tree pipes, Est , and K2 ∈ Nnpct ×np is the permutation matrix
260
for the pipes in the complementary co-tree edges, Ect . It is important to note that there are many
261
choices of spanning tree for any cyclic graph. The choice of spanning tree and co-tree combination does
262
not affect the correctness of the method.
264
M AN U
TE D
263
SC
252
By exploiting the relationship between the spanning tree and cotree, Elhay et al. (2014) proposed
EP
(m+1)
q1
265
(m+1)
and second for the co-tree flows q2
AC C
266
(m+1)
267
268
269
270
(m+1)
the following equations to solve the WDS for the flows: first for the spanning tree flows q1
W (m+1) q 2
(m)
= L21 T q 2
b − R−T 1 d
272
273
(12)
:
(m) (m) (m+1) (m+1) (m+1) (m) − G1 q1 + F 2 − G2 q 2 + a2 = L21 F 1
(13)
(m) (m) b b where: R1 = K1 Aˆ1 RT ; R2 = K2 Aˆ1 RT ; L21 = −R2 R−T = K1 F (m) K1 T ; F 2 = K2 F (m) K2 T 1 ; F1 (m) (m) (m) (m) b b G1 = K1 G(m) K1 T ; G2 = K2 G(m) K2 T ; W (m) = L21 (F 1 )−1 L21 T + (F 2 )−1 ; a1 = K1 Aˆ2 el ;
a2 = L21 K1 Aˆ2 el + K2 Aˆ2 el . (0)
271
,
Note that in Eq. (12), an initial set of the co-tree flows q2
is needed to commence the solution
process. The heads are found after the iterative process of the RCTM has been completed by using a linear
10
ACCEPTED MANUSCRIPT 274
solution process: (m+1)
R1 h = F 1 q 1
(m)
− (F 1 − G1 ) q 1
− a1
(14)
This partitioning of the network equations reduces the size of the non-linear component of the solver
276
to np − nj (the number of co-tree elements in the network). It has been proven by Elhay et al. (2014)
277
that the RCTM and the GGA have identical iterative results and solutions if the same starting values
278
are used. However, for the RCTM, the user only needs to set the initial flow estimates for the co-tree
279
pipes, q2 , in contrast to GGA where initial flow estimates are required for all pipes. The flows in the
280
complementary spanning pipes are generated by Eq.(12) in the RCTM.
281
4 WDSLib Structure
282
WDSLib is a WDS simulation toolkit consisting of a set of C++ member functions, which henceforth
283
will be referred to just as functions, that can be composed to solve for the steady state solution of a
284
WDS. WDSLib can be used for a once-off simulation or a multi-run simulation. Pre-packaged driver
285
code is provided to perform once-off simulations using a choice of solver methods. For a multi-simulation
286
setting, where the use-cases are very diverse, the user is able to select the desired components of WDSLib
287
to compose and compile their own driver.
M AN U
SC
(0)
RI PT
275
Individual functions in WDSLib are classified according to their role in the simulation workflow. In
289
any simulation workflow, there will be functions that will only have to be executed once. For example,
290
functions to read the input file or partition the network will only have to execute once at the start of the
291
simulation (or of all simulations). Likewise, code to reverse the network partitioning and write simulation
292
results will only have to execute once at the end of the simulation. In this work, these functions that are
293
only required to be run once are called level one (L1 ) functions. L1 functions relate to network topology,
294
which is invariant for the whole simulation. In a multi-simulation setting, certain functions will need to
295
be run once for every hydraulic-phase. An example of such a module is the module making the initial
296
guesses of pipe flow rates for the updated network configuration. In this work, these, once-per-assessment
297
functions, are called level two (L2 ) functions.
AC C
EP
TE D
288
298
Finally, for every hydraulic assessment there is a non-linear iterative phase in the solution process.
299
The functions in this phase run many times for each hydraulic assessment until the stopping test has
300
been satisfied. Examples of these include the functions to calculate the G and F matrices (see Eqs. (3)
301
and (4)) and running the Cholesky solver. These iterative-phase functions are called level three (L3 )
302
functions.
303
Fig. 1 illustrates the global structure of WDSLib under a once-off simulation setting and a multi-run
304
simulation setting. The modular setup of WDSLib allows each module to be run the minimum number
305
of times determined by its simulation setting. Under the module structure described above a once-off
11
ACCEPTED MANUSCRIPT simulation setting can be viewed as a special case where the L1 functions and L2 functions are both run
307
once. Note that after running the initial L1 functions it is possible to run hydraulic assessments of the
308
network in parallel. This mode of execution might be used in a design setting such as using a genetic
309
algorithm (GA) to optimize pipe diameter sizes.
TE D
M AN U
SC
RI PT
306
Figure 1: Global structure of WDSLib for both simulation settings
L1 and L2 functions are classified into parts a and b according to whether they run before or after
311
the lower level processing that they embed. These functions are detailed in Fig. 2. The L1 functions
312
that run at the start of the simulation are called L1a functions. These include the module (1) to read
313
the configuration file; (2) to parse the EPANET .inp file using EPANET programmer’s toolkit; (3) to
314
partition the network; and (4) to solve the linear part(s) of the network. The corresponding L1b functions
315
are run at the end of the simulation. These include tasks such as reversing the network partitioning.
316
Note that certain L1a functions require their corresponding L1b functions to be used. For example the
317
forest search module needs to be paired with the reverse FCPA permutation. There is a similar structure
318
for L2 functions. L2a functions are run at the start of each hydraulic assessment and L2b functions run
319
at the end. The functions that must be included for the FCPA method are denoted with single asterisks.
320
Likewise the functions that must be included with the RCTM method are denoted with double asterisks.
321
For these methods to work correctly all affiliated functions must be included in the simulation workflow.
322
Note that it is also possible to run both the RCTM and FCPA in the same workflow. Also note that
323
the user cannot run both GGA and RCTM in the same workflow – the user must choose between these
AC C
EP
310
12
ACCEPTED MANUSCRIPT Table 1: Key function descriptions, names, their classes, inputs and outputs. The affiliated functions are shown in sub-tables (1a) (1b) (1c) (1d). (a) Shared Modules
Module name readConfig getInputData scale AMD getRf init getGF stopTest rScale
Class Simulation Input Solver Suitesparse Solver Solver Solver Solver Solver
Input config file name EPANET .inp file void void net diameter net, resistance constant result void
(b) Global gradient algorithm (GGA)
Module name runH
Class GGASolver(Solver)
Input void
Output void
SC
Description GGA Solver
Output void EPANET err code void void resistance constant flow rate void norm void
RI PT
Description Read the configuration file Read EPANET input file Variables scaling AMD bandwidth reduction Calculate the resistance constants Generate initial guesses of flows Calculate the head loss coefficients Stopping test Recover scaled variables
(c) Forest-core algorithm (FCPA)
Module name forestSearch forestFlow forestHead rFCPA
Class topology Solver Solver Solver
Input SN, EN demands result void
M AN U
Description Forest search Calculate flows in forest Calculate heads in forest Reverse FCPA permutation
Output void flows in forest pipes heads in forest pipes void
(d) Reformulated cotree flows method (RCTM)
324
Module name STSearch runH RCTMHead rRCTM
Class topology RCTMsolver (Solver) RCTMsolver RCTMsolver
Input SN, EN void flows in ST and CT void
Output void void void void
TE D
Description Spanning tree search RCTM solver Calculate heads in ST and CT Reverse RCTM permutation
solution methods.
Table 1 provides a mapping from the function descriptions in Fig. 2 to the function names in WDSLib.
326
In addition, the dependencies between functions for each solution method are shown in Table 1a, Table
327
1b, Table 1c and Table 1d. The columns in each table list, respectively, the description of the function,
328
its name in WDSLib, the C++ class in which it appears, its input parameters, and its output values.
329
Note, that void is used in these latter two columns to denote that the function interacts with the class
330
variables rather than through its parameters and return value. Examples of how these functions can be
331
coded are presented in section 6. The key data-access functions in WDSLib are described next.
332
Getter and Setter methods
333
parameters and retrieving the results of the WDS network. These methods allow the user to reconfigure
334
the network before and during simulation runs. The names of the setter methods all start with a prefix
335
set and the names of the getter methods all start with a prefix get. For example, a user can set (write-to)
336
the diameter of pipe index to value by calling pipe->setD(index,value) and get (read-from) the head
337
of node index by calling h[index]=result->getHsol(index). A summary of the variables that can be
AC C
EP
325
Each class in WDSLib has various methods available for setting the network
13
AC C
EP
TE D
M AN U
SC
RI PT
ACCEPTED MANUSCRIPT
Figure 2: Function classification in WDSLib. The functions marked with single asterisks must be used for the FCPA method. The functions marked with double asterisks must be used for the RCTM method. Note that it is possible to use both methods at the same time.
14
ACCEPTED MANUSCRIPT Table 2: The getter and setter functions of each class and the variables they access Description Basic network properties, & Pipe and Node Node properties Pipe properties Flag information Parameter information Manage hydraulic simulation Parent class of solution methods GGA solution method RCTM solution method Network topology information Results of the simulations
Read-Access Write-Access Node,Pipe, np , nj , ns d, zs , zu SN , EN , L, D, R, pipe ID getFlag(”flagN”,flagV) setFlag(”flagN”) getPara(”paraN”,paraV) setPara(”paraN”) Result Result getCore, getForest qIter, hIter, GIter , FIter, numIter, CresIter, EresIter, Time
RI PT
Class Name Net Node Pipe Flag Parameter Simulation Solver GGASolver RCTMSolver Topology Result
read-from (read-access through getter methods) and written-to (write-access through setter methods for
339
each key classes is specified in Table 2. This concludes the discussion of the the broad structure of the
340
WDSLib package. The next section describes key aspects of the implementation of the package.
341
5 WDSLIB: Toolkit Implementation
342
This section outlines key implementation details of WDSLib. As previously mentioned, the overall aim of
343
WDSLib is to provide a clearly-structured, flexible and extensible hydraulic simulation toolkit that allows
344
testing, evaluation, and use, in production settings, of both existing and new WDS solution techniques.
345
These aims require WDSLib be implemented so that it is fast to execute, flexible to configure, robust to
346
challenging input data cases, and easy to understand and modify. The following describes aspects of the
347
implementation of WDSLib that enable it to meet these requirements. The next subsection describes the
348
general considerations that informed the design of the whole toolkit. This general discussion is followed
349
by a summary of key improvements to the solution processes encoded in forest searching and spanning
350
tree searching in the WDSLib package.
351
5.1 General capabilities and properties
352
This sub-section describes design aspects underpinning the utility and performance of WDSLib. In-
353
turn, the following outlines measures taken to: (1) maximize code clarity and modularity; (2) increase
354
the efficiency of memory access and storage; (3) maximize numerical robustness; (4) facilitate accurate
355
timing of code execution; and (5) maximize simulation speed for different settings.
356
Design Considerations 1: Modularity
357
The modular design of WDSLib is central to the evaluation and testing of different WDS solution
358
methods. All methods have been defined to perform a single, well–defined, function and each class can
AC C
EP
TE D
M AN U
SC
338
15
ACCEPTED MANUSCRIPT be compiled, used and tested independently. These features allow users to assemble the methods of
360
interest from independently developed components to create a customized WDS solution method in a
361
reliable way. WDSLib’s modular design also allows the users to profile the computation time of each
362
individual component of an algorithm. Functions communicate through well-defined interfaces and the
363
function code has been factored to minimize development and testing cost. This architecture allows
364
customized simulation applications (i) to combine the functions of interest and (ii) to implement new
365
solution algorithms to extend the functionalities of WDSLib.
366
Design Considerations 2: Memory Considerations
367
Care was taken to minimize the memory footprint of executing code (in order to reduce memory re-
368
quirements and prevent memory leaks) in the interest of the toolkit efficiency and toolkit robustness.
369
Reducing memory requirements allows the solution of larger WDS problems for a given memory ca-
370
pacity. In WDSLib, memory reduction was achieved through both, using sparse matrix representations
371
and the systematic allocation and deallocation of working structures in the C++ code. The matrices
372
used in WDS simulation are often sparse, with the density of the full node-arc incidence matrix being
373
only 2/nj . Consequently, it is more efficient to store these matrices using sparse storage schemes which
374
store only the non-zero elements of the matrix and pointers to their locations (Davis et al. 2014). It is
375
important to note that the choice of a sparse matrix representation is made based on (1) the storage
376
requirements of the matrix and (2) common search orders to column elements and row elements. This
377
latter factor means that the best format for sparse matrix representation varies with the preponderant
378
orders of search, (row-wise, column-wise, or both), employed by each method. There is a number of
379
common storage formats for sparse matrices (Compressed column storage (CCS) of Duff et al. (1989)),
380
Compressed row storage (CRS), Block Compressed column storage (BCCS), Block Compressed row stor-
381
age (BCRS), and Adjacency lists). As will be described shortly, WDSLib, uses a modified adjacency-list
382
representation.
AC C
EP
TE D
M AN U
SC
RI PT
359
383
Other implementations use a variety of storage schemes. In EPANET 2, the A1 matrix is stored
384
as two arrays of node indices, which represent start nodes (SN ) and the end nodes (EN ) of each pipe.
385
The i − th entry of the SN and EN arrays represent the start node and end node of i − th pipe of the
386
network. This storage format minimizes the memory required to store the A1 matrix because only the
387
indices are required to be stored because [A1 ](i,SNi ) = −1 and [A1 ](i,ENi ) = 1. As shown in Table 4,
388
searching through rows (pipes) of matrices that are stored in this format is efficient. However, searching
389
though the columns (nodes) is relatively inefficient. This storage format is also used in CWSnet.
390
Both CCS and CRS are used in the FCPA implementation reported in Simpson et al. (2014), and the
391
RCTM implementation reported in Elhay et al. (2014). The partial update null space method (Abraham
392
and Stoianov 2015) used CCS. The memory requirement for storing the A1 matrix in CCS is 2 × nnz +
16
ACCEPTED MANUSCRIPT 393
nj + 1 as shown in Table 4. This storage scheme is fast for searching through columns (nodes) of matrices
394
that are stored in CCS and slow for searching though rows (pipes). Table 3: The adjacency-list matrix presentation
nj
adjacent to {(vi , ej )|vi ∈ N (v1 ) ej connects v1 and vi } {(vi , ej )|vi ∈ N (v2 ) ej connects v2 and vi } {(vi , ej )|vi ∈ N (v3 ) ej connects v3 and vi } .. .
Size Deg(v1 ) Deg(v2 ) Deg(v3 ) .. .
RI PT
Node Index 1 2 3 .. .
{(vi , ej )|vi ∈ N (vnj ) ej connects vnj and vi } Deg(vnj ) Note: (1) ej is negative when vi is the start node and (2) ej is positive when vi is the end node
In WDSLib, a modified adjacency list, described in Table 3, tailored for WDS hydraulic simulation,
396
is used. An adjacency list for an undirected and unweighted graph consists of nj unordered lists for each
397
vertex ni , which contains all the vertices to which vertex ni is adjacent. The network that is shown in
398
Fig. 3 has one source, three nodes, and four pipes. The adjacency list for this network can be described
399
by four lists {{2, 3}, {1, 4}, {1, 4}, {2, 3}}. Each list describes the set of adjacent vertices of a vertex in
400
the graph. For example, the first list, {2, 3}, represents that the vertex 1 is adjacent to the vertex 2 and
401
vertex 3.
AC C
EP
TE D
M AN U
SC
395
Figure 3: A simple sample network. Numbers denote junction and pipe indices in the network.
402
The adjacency list is modified to include a directed and weighted graph for WDSLib. This mod-
403
ified adjacency list for a directed and weighted WDS graph consists of nj unordered lists for each
404
vertex ni . This list contains all the vertex and edge pairs to which vertex ni is adjacent. For ex-
405
ample, the adjacency list for the same network that is shown in Fig. 3 can be described by four lists
406
{{(2, 1), (3, −4)}, {(1, −1), (4, 2)}, {((1, 4), (4, 3)}, {(2, −2), (3, −3)}}. Each list represents the set of adja-
407
cent vertex and edge pair of a vertex in the graph. For example, the first list, {(2, 1), (3, −4)}, describes
408
that the vertex 1 is adjacent to the vertex 2 by edge 1 (from vertex 1 to vertex 2) and the vertex 3 by
409
edge 4 (from vertex 3 to vertex 1). It is fast to search through both the rows and columns of the A1 17
ACCEPTED MANUSCRIPT Table 4: Different sparse representations for A1
CCS CRS EPANET WDSlib
410
size(A1 )
size(A2 )
size([A1 A2 ])
Column Search
Row Search
2 × nnz + nj + 1 2 × nnz + np + 1 -
2 × nnz + nf + 1 2 × nnz + np + 1 -
4 × np + nn + 2 6 × np + 2 2 × np 4 × np
O(n) O((np )n) O(n) O(n)
O((nj )n) O(n) O((nj )n) O(n)
RI PT
Types
matrices that are stored in this format.
411
In addition to these optimized encodings, both G and F are diagonal square matrices, which require
412
less storage when stored as vectors than in sparse matrix form. The storage methods used for the
413
variables in WDSLib and their associated memory usage are given in Table 5.
As a final note, to offer further assurance of the correctness of memory management in WDSLib,
415
Valgrind (Nethercote and Seward 2007), a programming debugging tool, was deployed during testing to
416
detect any memory leaks, memory corruption, and double-freeing.
M AN U
SC
414
Table 5: Vectors and matrices in WDSLib type vector vector matrix matrix matrix
size np × 1 nj × 1 np × np np × nj (np − nj ) × nj
storage method vector vector vector sparse matrix sparse matrix
memory requirements np × double nj × double np × double (2 × np )× integer ≤ (np − nj ) × nj × integer
TE D
variables q, L, D, r h, d G, F A1 , A2 L21
Design Considerations 3: Numerical Considerations
418
The calculations in WDSLib are performed in C++ under IEEE-standard double precision floating point
419
arithmetic with machine epsilon mach = 2.22 × 10−16 . Invariant terms and parameters in every equation
420
were evaluated in advance and replaced by full 20-decimal digit accuracy constants. Intermediate results
421
of calculations, (which are not easily accessible in EPANET), can be output at the user’s request. The
422
stopping tolerance and stopping test can be set by the user either through the configuration file or by
423
the relevant setter method in the Parameter class.
AC C
EP
417
424
In the construction of any numerical solver, there are two primary dangers that are associated with
425
floating point arithmetic that cannot be ignored: (i) subtractive cancellation and (ii) overflow and
426
underflow. To avoid problems associated with these, all input variables are scaled to a similar range
427
to minimize the risk of avoidable computational inaccuracy or failure in floating point arithmetic. It is
428
important to note that unscaled or poorly scaled variables can unnecessarily confound a computation.
429
These scaled input variables are physically dimensionless, which allows computation which is independent
430
of the system of measurement units.
18
ACCEPTED MANUSCRIPT Table 6: WDS variables and units Variables Length Diameter Nodal head Source elevation flow demand G, F
SI unit m m m m m3 /s m3 /s s/m2
US Customary unit ft ft ft ft f t3 /s f t3 /s s/f t2
Scaling factor L0 = max (L) D0 = max (D) h0 = max (el ) el0 = max (el ) q0 = max (d) d0 = max (d) L0 n−1 G0 = D p |q0 |
RI PT
0
The variables, that are provided in EPANET input file for the package, and their corresponding units
432
in US customary and SI units are shown in Table 6. As with the input variables, the system equations
433
were modified to use dimensionless variables. Once the stopping test is satisfied, the original variables
434
are then recovered by reversing the initial scaling. Details of the scaling are shown in Appendix A.
SC
431
To help ensure that WDSLib solution methods are both fast and reliable. The sparse matrix oper-
436
ations are implemented using SuiteSparse (Davis et al. 2014). SuiteSparse is a state-of-the-art sparse
437
arithmetic suite with exceptional performance, from which the approximate minimum degree permuta-
438
tion (AMD) and the sparse Cholesky decomposition routines have been used.
439
Design Considerations 4: Timing Considerations
440
When executing WDSLib, each function reports the time spent in it by sampling wall clock time at
441
the start and end of its execution. Although the overhead for sampling wall clock time is small, there
442
are at least two special considerations involved in the interpretation of these timings: (i) the operating
443
system, at its own discretion, may launch background processes (for example anti-virus software), which
444
will distort the timings and (ii) extrapolating the timing for multiple hydraulic simulations from a single
445
analysis (as may be required, for example, in a genetic algorithm or other evolutionary algorithm run)
446
must be done with care because the relationship between the different settings is not linear.
EP
TE D
M AN U
435
This concludes the discussion of the main considerations concerning the global design of WDSLib. In
448
the following, key details of the implementation of selected parts of the solution processes are described.
449
5.2 Key Improvements to Solution Processes
450
The WDSLib implementation makes several improvements to extant solution processes. This section
451
focuses on the improvement of the network partitioning processes in FCPA and RCTM.
452
Key Optimization 1: Improvements to partitioning in Forest Search
453
The forest-core partitioning algorithm (FCPA) in this paper is a substantial improvement over the
454
algorithm of the original paper (Simpson et al. 2014). Specifically the original FCPA algorithm almost
455
always required many sweeps of the columns (nodes) of the A1 matrix in order to reduce the forest
AC C
447
19
ACCEPTED MANUSCRIPT 456
component down to the core component. The refined algorithm exploits the adjacency list representation
457
of the A1 matrix so that the partitioning process is achieved in a single sweep. This improves the speed
458
of the partitioning process from being O(anp ) to O(np + nf ) where a is the depth of the deepest tree-
459
component in the forest. This can lead to substantial time savings in the case when a is relatively
460
large. The pseudo-code for this refined forest search algorithm is shown in Appendix B This algorithm
462
traverses each tree component in turn from its leaf nodes, which maximizes the locality of operations
463
with respect to the graph representation. In this algorithm, a node is identified as a leaf node when
464
its node degree is one. Every time a leaf node, node k, is identified, the node pointer is moved to its
465
adjacent node, node k, and the node degree of node k is reduced by one. This process repeats if the
466
adjusted node degree of node k is one. Otherwise, node k is the root node for this tree and the algorithm
467
progresses to the next tree in the forest.
468
Key Optimization 2: Improvements to Spanning Tree Search
469
The reformulated co-tree flows method (RCTM) in this paper is also a substantial improvement over
470
the algorithm of the original paper (Elhay et al. 2014). The original spanning tree search algorithm
471
sweeps the rows of the A1 matrix (pipes) in order to identify the singleton rows and their corresponding
472
columns. The spanning tree search in the original RCTM required a sweep of the A1 matrix to identify
473
the next pipe in the spanning tree. This algorithm is O(np nj ), which is relatively inefficient.
TE D
M AN U
SC
RI PT
461
The pseudo-code for the refined spanning tree/forest search algorithm is shown in Appendix C. This
475
improved algorithm takes as input the adjacency list describing the network and the pipe indexes of
476
the core component of the network from the Algorithm 1 (if the FCPA is used). In this algorithm, all
477
water sources are the starting point of the search process, SN , and marked as visited. The nodes in SN
478
are then used as to identify a spanning tree within the WDS. This is achieved by repeatedly finding all
479
adjacent pairs, node t and pipe s, of and removing the first node in SN by using the adjacency list. If
480
the adjacent node t is not visited then node t is inserted into the spanning-tree node vector, ST N , and
481
search node vector, SN , and node t is marked as visited and pipe s to the spanning-tree pipe vector,
482
ST P , and pipe s is marked as visited. If the adjacent node t is visited and the pipe s is not visited
483
then the pipe s is inserted into the co-tree pipe vector, CT P and mark pipe s as visited. This process
484
is repeated until SN is empty. The overall time-complexity of this algorithm is O(np + nj ) (compared
485
to O(np nj ) as mentioned above) is the same as the best asymptotic complexity of breadth-first search
486
on a graph.
AC C
EP
474
20
ACCEPTED MANUSCRIPT
6 Example Applications
488
WDSLib consists of a collection of functions which can be used either as a standalone application for fast
489
one-off simulations or as a library of software components that can be integrated into a user’s own WDS
490
solution processes. This section presents two example applications. The first application is the setup for
491
a basic one-off simulation of a WDS. The second application (described in subsection 6.1) presents an
492
example using WDSLib to implement a simple 1+1 Evolutionary Strategy (Beyer and Schwefel 2002)
493
(1+1-ES or, more commonly, 1+1EA) for sizing pipes in a WDS.
494
Example 1 - Once-off Simulation
495
The setup for WDSLib as a standalone application is straightforward. The user provides a configuration
496
text file that specifies input and output filenames; the name of the solver; the desired output variables;
497
and simulation parameters. These values have sensible defaults so the user can set up the solver by using
498
a minimal configuration such as that shown in Fig. 4. By using this config file, WDSLib is configured to
499
run a single hydraulic analysis of the network that is stored as say ”hanoi.inp”, an EPANET- formatted
500
input file, under ”Network/” sub-directory, using the reformulated co-tree flows method with the forest-
M AN U
SC
RI PT
487
core partitioning algorithm. The full set of configuration parameters for once off simulations is shown in
EP
TE D
[ InputFile ] < directory > % the input file d i r e c t o r y Network / < file > % the input file name hanoi . inp % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -% [ controlFlag ] < SolverFlag > % 1 for GGA 2 for RCTM 2 < FCPAFlag > 1 % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -%
AC C
1 2 3 4 5 6 7 8 9 10 11 12 13 14
Figure 4: A minimal configuration file to run the GGA in WDSLib
501
502
Fig. 11 in Appendix D.
503
6.1 Example 2 - A Simple Network Design Application
504
As a minimalist example of the application of WDSLib to a WDS network design problem, the following
505
example uses 1+1EA for optimally sizing pipe diameters. This algorithm takes an existing network
506
with randomly generated pipe diameters and optimizes the network to minimize cost, subject to given
507
pressure head constraints. A 1+1EA is a very simple evolutionary strategy (Beyer and Schwefel 2002)
508
which starts with a randomly generated individual (in this case a WDS diameter configuration). This
21
ACCEPTED MANUSCRIPT
M AN U
SC
RI PT
# include " Simulation . h " # include " result . h " using namespace std ; # define GTN 100000 // the total number of g e n e r a t i o n s /* a v a i l a b l e pipe d i a m e t e r s ( inches ) */ vector < int > ADiameter ={36 ,48 ,60 ,72 ,84 ,96 ,108 ,120 ,132 ,144 ,156 ,168 ,180 ,192 ,204}; /* dollar per unit length */ vector < int > unitCost ={93.6 ,133.7 ,176.3 ,221 ,267.6 ,315.8 ,365.4 ,416.5 ,468.7 , 522.1 ,576.6 ,632.1 ,688.5 ,745.1 ,804.1}; vector < int > g e n e r a t e i n i t i a l D i a m e t e r s ( int ) ; // see fig . 6 vector < int > mutate ( vector < int >) ; // see fig . 7 double evaluate ( Net * , Result *) ; // see fig . 8 int main ( int argc , char * argv []) { srand (1) ; char * config = argv [1]; Result * result = new Result () ; Simulation * simulation1 = new Simulation () ; // i n i t i a l i z e the s i m u l a t i o n class Net * net = simulation1 - > L1a ( result , config ) ; // perform the L1a f u n c t i o n s vector < int > p1 = g e n e r a t e i n i t i a l D i a m e t e r ( net - > getNp () ) ; // initial guesses of d i a m e t e r vector < int > p2 ; // work space for current best double cost , currentbest ; for ( int i =0; i < GTN ; i ++) { simulation1 - > setD (& p1 ,& ADiameter , net - > getPIPE () -> Diascale () ) ; simulation1 - > solve ( result ) ; // perform the L2 and L3 f u n c t i o n s simulation1 - > L1b ( result ) ; // reverse the p e r m u t a t i o n cost = evaluate ( result , net ,& p1 ) ; // e v a l u a t e the network cost if ( cost < currentbest || i ==0) { // s e l e c t i o n o p e r a t o r currentbest = cost ; p2 = p1 ; cout < dispResult ( result ) ; delete simulation1 ; delete result ; return 0; }
TE D
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38
Figure 5: Example code for 1+1EA for Pipe Sizing 1+1EA then progresses by applying a mutation to a random pipe diameter size, and then evaluating the
510
new individual. If the new individual is better it replaces the old network. This process continues in a
511
loop until a given number of evaluations is reached.
513
514
515
516
517
The C++ code for this example is shown in Figs. 5, 6, 7, and 8. If the name of the file containing this code is: simpEA.cc then the simplest command to compile this code is:
AC C
512
EP
509
g++ simpEA.cc -o simpEA -Llib -lWDSLib To run this code the user would type: ./simpEA config.txt where config.txt contains the same configuration text as for the previous example.
518
Starting with the main function in Fig. 5, line 15 points to the config file specified by the command
519
line. The next two lines initialize the result and the simulation according to the configuration file.
520
This is followed by the L1a module to perform the user selected L1a functions. Line 19 generates the
521
initial pipe diameters of the network and line 20 initializes the workspace for the mutated string. Line 22
ACCEPTED MANUSCRIPT
vector < int > g e n e r a t e i n i t i a l D i a m e t e r s ( int np ) { vector < int > Diameter = vector < int >( np ) ; for ( int i =0; i < np ; i ++) Diameter [ i ]= rand () % ADiameter . size () ; // set the index for pipe d i a m e t e r s return Diameter ; }
RI PT
1 2 3 4 5 6 7
vector < int > mutate ( vector < int >* string ) { vector < int > string1 ; string1 =* string ; int aa = rand () %( string - > size () ) ; // choose which pipe to mutate int a = rand () % ADiameter . size () ; // choose a pipe d i a m e t e r after the m u t a t i o n ( string1 ) [ aa ]= a ; // set the pipe index return string1 ; }
M AN U
1 2 3 4 5 6 7 8
SC
Figure 6: Implementation code for pipe size initialization
double evaluate ( Result * result , Net * net , vector < int > * p1 ) { PIPE * pipe = net - > getPIPE () ; double np = net - > getNp () ; double nj = net - > getNj () ; double P =1 e7 ; double cost1 =0; penaltyCost =0; vector < double > hsol = result - > getHsol () ; // get the vector of nodal heads vector < double >* L = pipe - > getL () ; // get the vector of pipe lengths vector < double >* zu = net - > getNODE () -> getZU () ; // get the vector of nodal e l e v a t i o n s double Lscale = pipe - > Lenscale () ; // get the scaling factor for length for ( int i = 0; i < np ; i ++) cost1 +=(* L ) [ i ]* Lscale * unitCost [(* p1 ) [ i ]]; // c a l c u l a t e the m a t e r i a l cost for ( int i = 0; i < nj ; i ++) if (( hsol ) [ i ] -(* zu ) [ i ] < minP ) penaltyCost +=( minP -( hsol ) [ i ]+(* zu ) [ i ]) * scaleP ; // c a l c u l a t e penalty cost return cost1 + penaltyCost ; }
AC C
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
EP
TE D
Figure 7: Implementation code for the mutation operator
Figure 8: Implementation code for the evaluation function
23
ACCEPTED MANUSCRIPT 23 sets the pipe diameters of the network. Line 24 evaluates the current network configuration. The
523
permutation and scaling for the current individual is reversed by L1b in line 25 of Fig. 5. Line 26
524
calculates the fitness of the current network configuration by using the evaluate function in Fig. 8. This
525
function applies a penalty for pressure head constraint violations and pipe material costs. The body
526
of the 1+1EA is contained in the selection operator and mutation operator that follow. Lines 27 to 31
527
compare the string in the current generation with the current best string if the individual p1, as measured
528
by evaluate is better than the individual p2 then p1 replaces p2. Line 32 mutates the current network,
529
p2, using mutate (see Fig. 7). The mutate function changes the diameter of a randomly selected pipe in
530
the network to a randomly selected diameter, chosen from a set of commercially available pipe diameters.
531
The mutated individual, stored in the workplace p1, is used as the network configuration for the next
532
iteration. Until the total number of generations is reached, the user selected information about the best
533
individual is outputted by dispResult in line 34 of Fig. 5.
SC
RI PT
522
It should be noted that the algorithm described above can be used to design a simple WDS but is
535
not optimal in terms of speed of convergence. Other EA’s such as genetic algorithms (Simpson et al.
536
1994) will perform better. However the above example has the advantage of simplicity and contains all
537
the basic elements that a GA would use when interacting with WDSLib.
538
M AN U
534
This concludes the presentation of examples in this work. The next section presents a case study that illustrates the performance of WDSLib in a multi-simulation setting.
540
7 Case Study
541
The following presents timing results for WDSLib running the 1+1EA described in the previous section.
542
The results below compare the four different solvers plus EPANET2. Note, that detailed timings for
543
once-off simulations comparing the four methods can be found in Qiu et al. (2018). Three networks
544
were benchmarked in these experiments. These were the N1 , N3 , and N4 case-study networks used in
EP
TE D
539
AC C
Simpson et al. (2014). Table 7 summarizes the characteristics of these networks.
Network N1 N3 N4
Full np 934 1975 2465
Table 7: Benchmark networks summary Network nj ns 848 8 1770 4 1890 3
Forest & Core Networks nf (nf /n# njc n pc p ) 361 (38%) 573 487 823 (42%) 1152 947 429 (17%) 2036 1461
Co-tree Network nct 84 205 757
545
546
Table 8 shows the results of the 1+1EA from Fig. 5 for the GGA, GGA with FCPA, RCTM, RCTM
547
with FCPA and the EPANET2 solvers. For each of the four WDSLib solvers above, the timings are
548
given for running the EA with and without the L1 modules hoisted out the main EA loop. Each
24
ACCEPTED MANUSCRIPT 549
experiment evaluates the WDS network 100,000 times. And the best performing method for each network
550
is highlighted in bold. It is important to note that 1+1EA using both the GGA and the WDSLib Table 8: The actual 1+1EA run-time with 100,000 evaluations (min.) for each of the four solution methods applied networks N1 , N3 , N4 GGA with FCPA min. 4.64 9.79 16.29
RCTM min. 4.53 13.75 23.92
RCTM with FCPA min. 4.13 10.30 21.93
EPANET min. 9.81 26.43 67.11
RI PT
N1 N3 N4
GGA min. 6.73 15.21 21.14
The results show that the EA runs using WDSLib are substantially faster than the runs using the
552
EPANET2 solver. This is, in part, due to the fact that the EPANET2 solver is designed as a standalone
553
solver which does not facilitate lifting out of invariant computations from the EA loop.
SC
551
As a demonstration of how the performance of an EA can be traced Fig. 9 shows the evolution of the
555
fitness values of the N1 network. These traces were extracted from a file written to in line 30 in Fig. 5.
556
As can be seen, the cost and the pressure head violation terms drop during the EA run. Note that there
M AN U
554
will be considerable variation between 1+1EA runs due to its highly stochastic nature.
3.4
10 8
TE D
3.2
3
Cost
2.8
2 10 0
AC C
2.2
EP
2.6
2.4
Total Cost Pressure Violation Penalty
10 1
10 2
10 3
10 4
10 5
Number of evaluations
Figure 9: The evaluation of the fitness value of network N1 557
558
8 Conclusions
559
This paper has described WDSLib, a library for steady-state hydraulic simulation of WDS networks.
560
WDSLib is fast, modular, and portable with implementation of several standard and recently published 25
ACCEPTED MANUSCRIPT hydraulic solution methods. We have outlined the supported solution methodologies, the structure of the
562
package and key aspects of WDSLib’s implementation. Two example applications have been presented
563
including a design case study using a simple EA. The EA results were benchmarked for different solvers
564
demonstrating a substantial improvement in speed over the industry standard EPANET2 package. These
565
benchmarks also have illustrated how the modular structure of WDSLib could be exploited to speed up
566
execution time in a design setting.
RI PT
561
As well as providing a fast simulation platform for both once-off and multi-run simulations, WD-
568
SLib also provides a testbed for comparing different solution methods in different settings and network
569
topologies. As such WDSLib is designed with a pluggable architecture which can extended to efficiently
570
incorporate new solution methods as they are created. This will enhance the capability of the research
571
community to demonstrate the efficacy of new methods without having to re-engineer the content of
572
shared WDSLib functions and data representations.
SC
567
WDSLib accepts the input file format used by EPANET2 which allows for simple deployment. The
574
amount of coding required to use the provided solvers in WDSLib is minimized by the use of a simple
575
and extensible configuration file.
M AN U
573
WDSLib is currently limited to finding solutions of demand-driven WDS systems without control
577
devices. The development of software components to simulate control devices and pressure driven models
578
is future work.
579
References
580
Abraham, E. and I. Stoianov (2015). “Sparse null space algorithms for hydraulic analysis of large-scale
581
water supply networks”. Journal of Hydraulic Engineering, p. 04015058. doi: 10.1061/(ASCE)HY.
582
1943-7900.0001089.
584
585
586
EP
Benzi, M., G. H. Golub, and J. Liesen (2005). “Numerical solution of saddle point problems”. Acta Numerica 14, pp. 1–137. doi: 10.1017/S0962492904000212.
AC C
583
TE D
576
Beyer, H. G. and H. P. Schwefel (2002). “Evolution strategies–A comprehensive introduction”. Natural Computing 1.(1), pp. 3–52. doi: 10.1023/A:1015059928466.
587
Cheung, P. B., J. E. Van Zyl, and L. F. R. Reis (2005). “Extension of EPANET for pressure driven
588
demand modeling in water distribution system”. Computing and Control for the Water Industry 1,
589
pp. 311–316.
590
591
592
Cross, H. (1936). “Analysis of flow in networks of conduits or conductors”. University of Illinois. Engineering Experiment Station. Bulletin; no. 286. Davis, Tim, WW Hager, and IS Duff (2014). “SuiteSparse”. http://faculty.cse.tamu.edu/davis/suitesparse.html.
26
ACCEPTED MANUSCRIPT 593
594
Deuerlein, J. W. (2008). “Decomposition model of a general water supply network graph”. Journal of Hydraulic Engineering 134.(6), pp. 822–832. doi: 10.1061/(ASCE)0733-9429(2008)134:6(822).
595
Deuerlein, J. W., S. Elhay, and A. R. Simpson (2015). “Fast Graph Matrix Partitioning Algorithm
596
for Solving the Water Distribution System Equations”. Journal of Water Resources Planning and
597
Management 142.(1), pp. 04015037-1–04015037-11. Diao, Kegong, Zhengji Wang, Gregor Burger, Chien-Hsun Chen, Wolfgang Rauch, and Yuwen Zhou
599
(2014). “Speedup of water distribution simulation by domain decomposition”. Environmental Mod-
600
elling & Software 52, pp. 253–263. issn: 1364-8152. doi: https://doi.org/10.1016/j.envsoft.
601
2013.09.025. url: http://www.sciencedirect.com/science/article/pii/S1364815213002247.
602
Duff, I. S., R. G. Grimes, and J. G. Lewis (1989). “Sparse matrix test problems”. ACM Transactions on Mathematical Software (TOMS) 15.(1), pp. 1–14.
SC
603
RI PT
598
Elhay, S. and A. R. Simpson (2011). “Dealing with zero flows in solving the nonlinear equations for
605
water distribution systems”. Journal of Hydraulic Engineering 137.(10), pp. 1216–1224. doi: 10 .
606
1061/(ASCE)WR.1943-5452.0000561.
M AN U
604
607
Elhay, S., A. R. Simpson, J. Deuerlein, B. Alexander, and W. H. Schilders (2014). “Reformulated co-
608
tree flows method competitive with the global gradient algorithm for solving water distribution sys-
609
tem equations”. Journal of Water Resources Planning and Management 140.(12), pp. 04014040-1–
610
04014040-10. doi: 10.1061/(ASCE)WR.1943-5452.0000431.
612
Epp, R. and A. G. Fowler (1970). “Efficient code for steady-state flows in networks”. Journal of the
TE D
611
Hydraulics Division 96.(1), pp. 43–56.
Estrada, Carlos, C´esar Gonz´ alez, Ricardo Aliod, and Jara Pa˜ no (2009). “Improved pressurized pipe
614
network hydraulic solver for applications in irrigation systems”. Journal of irrigation and drainage
615
engineering 135.(4), pp. 421–430.
EP
613
Gorev, N. B., I. F. Kodzhespirov, Y. Kovalenko, E. Prokhorov, and G. Trapaga (2012). “Method to Cope
617
with Zero Flows in Newton Solvers for Water Distribution Systems”. Journal of hydraulic engineering
618
139.(4), pp. 456–459. doi: 10.1061/(ASCE)HY.1943-7900.0000694.
AC C
616
619
Guidolin, M., P. Burovskiy, Z. Kapelan, and D. A. Savic (2010). “CWSNet: An object-oriented toolkit for
620
water distribution system simulations”. Proceedings of the 12th Annual Water Distribution Systems
621
Analysis Conference, WDSA, pp. 12–15.
622
Jun, L. and Y. Guoping (2012). “Iterative methodology of pressure-dependent demand based on EPANET
623
for pressure-deficient water distribution analysis”. Journal of Water Resources Planning and Man-
624
agement 139.(1), pp. 34–44. doi: 10.1061/(ASCE)WR.1943-5452.0000227.
625
Morley, M. S., R. M. Atkinson, D. A. Savic, and G. A. Walters (2000). “OpenNet: An applicationin-
626
dependent framework for hydraulic network representation, manipulation, and dissemination”. Hy-
627
droinformatics 2000 Conference, University of Iowa, Iowa City, USA, pp. 23–27.
27
ACCEPTED MANUSCRIPT 628
629
630
631
Morley, M. S. and C. Tricarico (2008). Pressure-Driven Demand Extension for EPANET (EPANETpdd). Tech. rep. University of Exeter. Nethercote, N. and J. Seward (2007). “Valgrind: a framework for heavyweight dynamic binary instrumentation”. ACM Sigplan notices. Vol. 42. 6. ACM, pp. 89–100. doi: 10.1145/1273442.1250746. Perelman, Lina and Avi Ostfeld (2011). “Topological clustering for water distribution systems analysis”.
633
Environmental Modelling & Software 26.(7), pp. 969–972. issn: 1364-8152. doi: https://doi.org/
634
10.1016/j.envsoft.2011.01.006. url: http://www.sciencedirect.com/science/article/
635
pii/S1364815211000259.
RI PT
632
Qiu, M., S. Elhay, A. R. Simpson, and B. Alexander (2018). “A Benchmarking Study of Water Distribu-
637
tion System Solution Methods”. Manuscript accepted for publication to Journal of Water Resources
638
Planning and Management.
640
Rahal, H. (1995). “A co-tree flows formulation for steady state in water distribution networks”. Advances in Engineering Software 22.(3), pp. 169–178. doi: 10.1016/0965-9978(95)00020-W.
M AN U
639
SC
636
641
Rossman, L. A. (2000). “Epanet 2 users manual, us environmental protection agency”. Water Supply and
642
Water Resources Division, National Risk Management Research Laboratory, Cincinnati, OH 45268.
643
Schilders, W. H. (2009). “Solution of indefinite linear systems using an LQ decomposition for the linear
644
constraints”. Linear Algebra and its Applications 431.(3), pp. 381–395. doi: 10.1016/j.laa.2009.
645
02.036.
647
Siew, C. and T. T. Tanyimboh (2012). “Pressure-dependent EPANET extension”. Water resources man-
TE D
646
agement 26.(6), pp. 1477–1498. doi: 10.1061/41203(425)9. Simpson, A. R., G. C. Dandy, and L. J. Murphy (1994). “Genetic algorithms compared to other tech-
649
niques for pipe optimization”. Journal of water resources planning and management 120.(4), pp. 423–
650
443. doi: 10.1061/(ASCE)0733-9496(1994)120:4(423).
EP
648
Simpson, A. R. and S. Elhay (2010). “Jacobian matrix for solving water distribution system equations
652
with the Darcy-Weisbach head-loss model”. Journal of hydraulic engineering 137.(6), pp. 696–700.
653
doi: 10.1061/(ASCE)HY.1943-7900.0000341.
AC C
651
654
Simpson, A. R., S. Elhay, and B. Alexander (2014). “Forest-Core Partitioning Algorithm for Speeding
655
Up Analysis of Water Distribution Systems”. Journal of Water Resources Planning and Management
656
140.(4), pp. 435–443. doi: 10.1061/(ASCE)WR.1943-5452.0000336.
657
Todini, E. and S. Pilati (1988). “A gradient algorithm for the analysis of pipe networks”. Computer
658
applications in water supply: vol. 1—systems analysis and simulation. Research Studies Press Ltd.,
659
pp. 1–20.
660
Vassiljev, A. and T. Koppel (2015). “Estimation of real-time demands on the basis of pressure mea-
661
surements by different optimization methods”. Advances in Engineering Software 80, pp. 67–71. doi:
662
10.1016/j.advengsoft.2014.09.023.
28
ACCEPTED MANUSCRIPT
TE D
M AN U
SC
RI PT
EPANET”. Advances in Water Supply Management (CCWI 2003).
EP
664
van Zyl, J. E., J. Borthwick, and A. Hardy (2003). “Ooten: An object-oriented programmers toolkit for
AC C
663
29
ACCEPTED MANUSCRIPT 665
Appendix A Scaling
666
In WDSlib, all input variables are scaled to a similar range to minimize the risk of avoidable computa-
667
tional inaccuracy or failures in floating point arithmetic. ˆ = h/h0 , eˆl = el /el0 , d ˆ = d/d0 where ˆ = G/G0 , q ˆ = q/q0 , h The variables are scaled as following: G
669
ˆ eˆl , ˆ h, G, h, el , and d are the original input vectors, G0 , h0 , el0 , and d0 are the scaling factors and G,
670
and dˆ are the scaled input vectors.
671
RI PT
668
ˆ 0 , el = eˆl el0 , d = dd ˆ 0 , Eq. (1) and Eq. (2) become: ˆ 0, q = q ˆ q0 , h = hh By substituting G = GG
ˆ − el0 A2 eˆl = 0 ˆ q − h0 A1 h G0 q0 Gˆ
673
ˆ=0 ˆ − d0 d q 0 A1 T q
674
Eq. (15) and Eq. (16) can be further simplified by introducing the notation: a = h0 / (G0 ∗ q0 ), b =
675
el0 / (G0 ∗ q0 ), and c = d0 /q0 :
M AN U
676 677
678
679
with the following matrix form:
680
683
A1 T qˆ − cdˆ = 0
(18)
ˆ bA2 eˆl −A1 q − = 0. ˆ ˆ O ah cd
(19)
Finally, the network can be solved by using Eq. (20) and Eq. (21): ˆ − A1 F ˆ (m+1) = 1 U −1 [−cd ˆ −1 [(F ˆ − G)ˆ ˆ q (m) + bA2 eˆl ]] h a
(20)
ˆ m+1 − F ˆ −1 A1 h ˆ −1 (Gˆ ˆ q (m) − bA2 eˆl ) ˆ (m+1) = q ˆ (m) + aF q
(21)
EP
682
(16)
(17)
AC C
681
(15)
ˆ − bA2 eˆl = 0 ˆ q − aA1 h Gˆ
TE D
ˆ G −A1 T
SC
672
684
Choice of scaling factors
685
The choice of the scaling factor, despite much research, is not well understood. In this subsection, a
686
choice for each scaling factor, based on the experience of the authors, is recommended. There are two
687
types of variables and parameters that need to be scaled: invariants and variants.
688
Data sets that have very wide range of values can confound numerical accuracy. As a result, it may
689
be preferable to scale the data to a narrower range. The default scaling factor for each of the input data
690
is chosen to be its maximum absolute value. For example, the scaling factor for demand is max(d), so
691
that its values range from zero to one.
692
In contrast, it is more difficult to choose a scaling factor a priori for values that vary between iterations 30
ACCEPTED MANUSCRIPT 693
(variants). This is because the range of variants can change as the iteration progresses. As a result, the
694
intermediate and the final results might not be within the same range as the initial guesses.
698
connected to a water source and a good choice of the scaling factor for nodal head is max(el ) because
699
the maximum nodal head cannot exceed the maximum elevation head of the fixed nodes.
696
RI PT
697
There are two variants that need to be scaled: q, h. A good choice of the scaling factor for the P d flow rate is because the demand at each node must be satisfied by the water sources in the WDS nf and it is a reasonable assumption that the all demands are equally carried by each pipe that is directly
695
700
During the process of the computation, the matrices G and F are scaled because their input variables
701
are scaled. For the Darcy-Weisbach head loss model, the diagonal elements of the matrix G are modelled
702
by: [G]jj = diag
703
where the friction factor f is modelled by the Swamee-Jain formula:
M AN U
704
) 8 Lj ( 2 ) 5 fj |qj | for j = 1, . . . , np π g Dj
SC
(
64 Re j fj = P3 (α + β /θ )(Re /2000)k k k j j k=0 0.25
705
log2 (θj )
709
710
711
712
713
if Rej ≥ 4000
TE D
708
(22)
In order to make sure the Reynolds number, a dimensionless variable, is not affected by scaling, D0 4ˆ qj q0 ν ˆ= ν is introduced, Reynolds number Rej becomes , which can be further simplified ˆ j D0 q0 πˆ ν (q0 /D0 )D 4ˆ qj to , where the scaling factors are canceled. ˆj πˆ νD In order to make sure f is also not affected by the scaling, ˆ = D0 is introduced, θj becomes
EP
707
if 2000 < Rei < 4000
5.74 4qj j + and Rej = for j = 1, . . . , np . Note that αk and βk are constant, the 3.7Dj πνD Re0.9 j j values of which can be found in Elhay and Simpson (2011) where θj =
AC C
706
if Rej ≤ 2000
θj =
ˆj D0 5.74 + for j = 1, . . . , np ˆ Re0.9 3.7Dj D0 j
which can be further simplified to
714
θj =
ˆj ˆj 3.7D
+
5.74 for j = 1, . . . , np Re0.9 j
715
where the scaling factors are canceled. It is evident that the friction factors remain the same because
716
the values for the only two variables Re and θ are unchanged.
717
718
ˆ where G0 = ( 8 ) L0 |q0 | and [G] ˆ jj = Finally, diagonal elements of G can be expressed as: G0 G, π 2 g D05 ( ) ˆj L diag f |ˆ q | for j = 1, . . . , np . ˆ5 j j D j
31
ACCEPTED MANUSCRIPT
721
For the Hazen-Williams head loss ) model, the diagonal elements of the matrix G are modelled by: ( 10.67Lj |qj |(n−1) for j = 1, . . . , np , where Cj is the Hazen-Williams coefficient for [G]jj = diag Cj1.852 Dj4.871 the j-th pipe. The Hazen-Williams coefficient, unlike the friction factor in the Darcy-Weisbach head
722
loss model, is independent of flow rate, pipe wall condition, and flow regimes, which means it is an
720
723
724
725
726
independent ( variable. As a result, the scaling factor for)Hazen-Williams G can simply be derived by: ˆ j L0 10.67L [G]jj = diag |ˆ q |(n−1) |q0 |(n−1) and the equation for diagonal elements of G 1.852 1.852 ˆ 4.871 D4.871 j ˆ C0 D Cj 0 j n o ˆ , where G0 = 10.67L0 |q0 |(n−1) for Hazen-Williams equation can be expressed as: G = diag G0 G C01.852 D04.871 ) ( ˆ Lj |ˆ q |(n−1) for j = 1, . . . , np . and [G]jj = diag 1.852 ˆ 4.871 j ˆ D C j
RI PT
719
j
Why scaling should be applied
728
Case 1 from Estrada et al. (2009) is used here to demonstrate the efficacy of using scaled variables. The
729
layout of the Case 1 from Estrada et al. (2009) is given in Fig. 10. Details of the pipe and node values
730
are given in Table 9 and Table 10.
TE D
M AN U
SC
727
EP
Figure 10: The layout of the Case 1 from Estrada et al. (2009)
Table 9: The nodal values of the Case 1 from Estrada et al. (2009)
AC C
Node ID Reservoir A Reservoir B Junc C Junc D
Elevation (m) 10 10 0 0
Demand (LPS) 2.7 0
Table 10: The pipe values of the Case 1 from Estrada et al. (2009)
Pipe ID Pipe 1 Pipe 2 Pipe 3 Pipe 4
731
Start Node Reservoir A Reservoir B Junc D Reservoir A
End Node Junc C Junc C Junc C Junc D
Length (m) 0.25 0.5 100 100
Diameter (mm) 1200 1200 1200 1200
Running this example network from EPANET and WDSLib, we get the following results:
32
ACCEPTED MANUSCRIPT Table 11: The nodal heads found by using EPANET and WDSLib
Node ID Junc C Junc D
EPANET WDSLib Head (m) 10 10 10 10
M AN U
735
In contrast, the energy residual and the continuity residual of the WDSLib result are 1.8723 × 10−12 and 1.7849 × 10−9 respectively which are both acceptable.
TE D
734
respectively. These two residuals are too large which indicates an inaccurate solution.
EP
733
The energy residual and the continuity residual of the EPANET result are 2.6973 and 5.3575 × 10−4
AC C
732
EPANET WDSLib Flow (LPS) 1.3446 1.5924 1.3446 1.5924 0.0108 0.0124 0.0108 0.0124
SC
Link ID Pipe 1 Pipe 2 Pipe 3 Pipe 4
RI PT
Table 12: The flows found by using EPANET and WDSLib
33
ACCEPTED MANUSCRIPT
Appendix B Forest Search Algorithm
EP
TE D
SC
M AN U
k ← 1; for i ← 1 to nj do n = i; while Deg(n) == 1 do for (Adjp ,Adjv ) ∈ Adj(n) do if Deg(Adjv ) == 1 then F orest[k] ←Union(F orest[k], (Adjp , n)) else q(Adjp ) ← d(n); d(Adjv ) ← d(Adjv ) + d(n); Deg(Adjv ) ← Deg(Adjv ) − 1; n = Adjv ; if Deg[n] ≥ 2 then RootN ode[k] = n; k ← k + 1; end if end if end for end while end for
RI PT
Algorithm 1: Forest Search Algorithm input : Adjacency list, d and Deg output: Forest, RootNode, q and d
AC C
736
34
ACCEPTED MANUSCRIPT
Appendix C Spanning Tree Search Algorithm
SC
// mark the source i as visited this node into the search node vector
// mark pipe as a spanning tree pipe as a corresponding spanning tree pipe
EP
TE D
M AN U
for i ← nj to nf do V N [i] = true ; SN ← i ; // insert end for while SN is not empty do i ← the first element in SN ; foreach (Adjp ,Adjv ) ∈ Adj[i] do if VN[Adjv ] == false then ST P ← Adjp ; ST N ← Adjv ; // mark node SN ← Adjv ; VN[Adjv ]=true; VP=[Adjp ]=true; end if else if VN[Adjp ]==false then CT P ← Adjp ; VP[Adjp ]=true; end if end foreach remove i from SN ; end while
35
pipes nodes pipes nodes pipes nodes
RI PT
Algorithm 2: Spanning Tree Search Algorithm input : adjList output: Spanning Tree and Co-Tree elements ST P ← {} ; // initialize an empty vector for spanning tree ST N ← {} ; // initialize an empty vector for spanning tree CT P ← {} ; // initialize an empty vector for co-tree V N ← {} ; // initialize a boolean vector for visited V P ← {} ; // initialize a boolean vector for visited SN ← {} ; // initialize an empty set of searching
AC C
737
// removed pipe j as a singleton pipe
ACCEPTED MANUSCRIPT
Appendix D Complete configuration files
M AN U
SC
RI PT
[ InputFile ] < directory > % R E Q U I R E D I N F O R M A T I O N the input file d i r e c t o r y WDSnetwork / < file > % R E Q U I R E D I N F O R M A T I O N the input file name n1 . inp % end of input file % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -% [ Parameter ] < maxIter > % maximum number of i t e r a t i o n 50 < initV > % initial veloc ity 0.3048 % 1 ft / s < StopTol > % s t o p p i n g t o l e r a n c e used 1.000 e -6 < NormTyp > % norm type 1 for 1 - norm , 2 for 2 - norm , 3 for inf - norm 3.0 < StopTest > % s t o p p i n g test used 1 for q - norm , 2 for h - norm 3 for q &h - norm 1 < SerP > % service p r e s s u r e 20.0 < MinP > % Minimum p r e s s u r e 40.0 < Demandfuc > % type of c o n s u m p t i o n fu nc t i o n 1.00 < MaxNpIterResult > % np t r e s h o l d for disping i t e r a t e s result 50.0 < MinImproTol > 1.0000 e -3 % end of p a r a m e t e r % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -% [ dispFlag ] < BasicFlag > false < NetInfoFlag > false < ConvergenceFlag > false < StatFlag > false < ScalingFlag > 0 < NodalResultflag > 0 < LinkResultflag > false < QitersFlag > false < HitersFlag > false < timeFlag > false % end of d i s p F l a g % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -% [ controlFlag ] < SolverFlag > % 1 for GGA 2 for RCTM 3 for SMPA 4 for PDM 1 < FCPAFlag > 0 % end of c o n t r o l F l a g % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -%
TE D
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
EP
Figure 11: A configuration file to run the RCTM in WDSLib
AC C
738
36
ACCEPTED MANUSCRIPT
Appendix E Nomenclature
740
Acronyms
741
CT
Co-tree
742
CTP
Co-tree pipes
743
DW
Darcy-Weisbach head loss formula
744
FCPA Forest-core partitioning algorithm
745
GGA
Global gradient algorithm
746
HW
Hazen-William head loss formula
747
WDS Water distribution system
748
RCTM Reformulated co-tree flows method
749
nnz
Number of non-zeros
750
STP
Spanning tree pipes
751
STN
Spanning tree nodes
752
ST
Spanning tree
753
Constants
754
g
Gravitational acceleration constant
755
ν
Kinematic viscosity of water
756
FCPA variables
757
C
758
Ec
759
Ef
760
Gc
Core subgraph
761
Gf
Forest subgraph
762
P
Permutation matrix for the pipes in the core, Ec , of the network
763
S
Permutation matrix for the pipes in the forest, Ef ,of the network
764
T
Permutation matrix for the nodes in the forest, Ef , of the network
EP
TE D
M AN U
SC
RI PT
739
AC C
Permutation matrix for the nodes in the core, Ec , of the network Set of core pipes (edges) in Gc Set of forest pipes (edges) in Gf
37
ACCEPTED MANUSCRIPT Vc
Set of core nodes (vertices) in Gc
766
Vf
Set of forest nodes (vertices) in Gf
767
Hydraulic variables for GGA
768
A1
Unknown-head node-arc incidence matrix
769
A2
Fixed-head node-arc incidence matrix
770
d
Vector of nodal demands
771
di
Demand of node i
772
D
Vector of pipe diameters
773
Dj
Diameter of pipe j
774
el
Vector of Fixed-head nodes elevation heads
775
elk
Fixed-head nodes elevation heads at node k
776
E
Set of pipes in graph G
777
EN
Vector of end nodes
778
EN j
End nodes of pipe j
779
f
Vector of Darcy-Weisbach friction factors
780
fj
Darcy-Weisbach friction factor of pipe j
781
F
Diagonal Matrix of generalized headloss derivatives when the headloss is modelled by either the
M AN U
TE D
EP
HW and the DW
782
SC
RI PT
765
[F ]jj
Generalized headloss derivatives for pipe j
784
G
785
G
786
[G]jj
rj |qj |n−1
787
h
Vectors of unknown heads
788
hi
Heads at node i
789
J
Jacobian matrix
790
L1a
Functions run once before multiple simulation
AC C
783
Full WDS graph
Diagonal matrix with elements rj |qj |n−1
38
ACCEPTED MANUSCRIPT L2a
Functions run once before hydraulic assessment
792
L3
Functions run every iteration
793
L2b
Functions run once after hydraulic assessment
794
L1b
Functions run once after multiple simulation
795
L
Vector of pipe lengths
796
Lj
Length of pipe j
797
n
Head loss exponent
798
nf
Number of forest pipes and nodes
799
nj
Number of junctions
800
np
Number of pipes
801
nr
Number of fixed-head nodes
802
nst
Number of ST pipes and nodes
803
nct
Number of CT pipes
804
q
Vector of unknown flows
805
qj
Flow in pipe j
806
R
Vector of Reynolds numbers
807
Rj
Reynolds number for pipe j
808
SN
Vector of start nodes
809
SN j
Start nodes of pipe j
810
U
811
V
812
V
Set of node in graph G
813
zi
Elevation at node i
814
αk
Interpolating spline coefficient
815
βk
Interpolating spline coefficient
816
θ
Vector as defined in Eq. (22)
AC C
EP
TE D
M AN U
SC
RI PT
791
Diagonal matrix of Schur Complement when headloss is modelled by HW Generalized Schur Complement when the headloss is modelled by both the HW and the DW
39
ACCEPTED MANUSCRIPT
Vector of pipe roughness heights
818
j
Roughness height for pipe j
819
RCTM variables
820
Est
Set of ST pipes (edges)
821
Ect
Set of complementary CT pipes (edges)
822
K1
Orthogonal permutation matrix for pipes in the ST
823
K2
Orthogonal permutation matrix for pipes in the CT
824
L21
A part of a basis for the null space of the permuted node-arc incidence for the RCTM
825
R
Orthogonal permutation matrix for nodes in the ST
826
Vst
Set of ST nodes (vertices)
827
W
Schur complement for the RCTM
AC C
EP
TE D
M AN U
SC
RI PT
817
40
ACCEPTED MANUSCRIPT
Highlights
EP
TE D
M AN U
SC
RI PT
A library for the steady-state analysis of a water distribution system (WDS) An open-source C++ software implementation of a number of WDS solution methods A fast simulation platform for both once-off and multi-run simulations Several improvements to the existing solution methods have been made
AC C
• • • •
S1364-8152(18)30416-X
DOI:
https://doi.org/10.1016/j.envsoft.2018.10.016
Reference:
ENSO 4335
To appear in:
Environmental Modelling and Software
Received Date: 2 May 2018 Revised Date:
8 October 2018
Accepted Date: 27 October 2018
Please cite this article as: Qiu, M., Angus R, B.A., Elhay, S.S., A software tool for assessing the performance of and implementing water distribution system solution methods, Environmental Modelling and Software (2018), doi: https://doi.org/10.1016/j.envsoft.2018.10.016. This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT
A Software Tool for Assessing the Performance of and Implementing Water Distribution System Solution Methods Mengning Qiu∗
RI PT
Bradley Alexander Angus R. Simpson Sylvan Elhay University of Adelaide, South Australia, 5005
SC
Highlights
M AN U
• A library for the steady-state analysis of a water distribution system (WDS) • An open-source C++ software implementation of a number of WDS solution methods. • A fast simulation platform for both once-off and multi-run simulations • Several improvements to the existing solution methods have been made.
TE D
Abstract
WDSLib is an extensible simulation toolkit for the steady-state analysis of a water distribution system. It includes a range of solution methods: the forest-core partitioning algorithm, the global gradient algorithm, the reformulated co-tree flow method, and also combinations of these methods.
EP
WDSLib has been created using a modularized object-oriented design and implemented in the C++ programming language, and has been validated against a reference MATLAB implementation. WDSLib has been designed: (i) to avoid unnecessary computations by hoisting each of the modules to
AC C
its appropriate level of repetition, (ii) to perform the computations independently of measurement units using scaled variables, (iii) to accurately report the execution time of all the modules in that it is possible to produce a timing model to parameterize multiple simulation times (such as in an optimization using a genetic algorithm) from a series of sampling simulation runs and (iv) to guard against numerical failures. Two example applications, a once-off simulation and a network optimization design application simulation, are presented. This toolkit can be used (i) to implement, test and compare different solution methods, (ii) to focus the research on the most time-consuming parts of a solution method and (iii) to guide the choice of solution method when multiple simulation runs are required. ∗ Corresponding Author Email: [email protected]
1
ACCEPTED MANUSCRIPT Keywords: Water Distribution System; C++ toolkit; Object-Oriented design; Forest-Core Partitioning Algorithm; Reformulated Co-tree Flows Method; Global Gradient Algorithm; open source software
Software availability
2
Name of the Software: WDSLib
3
Version: 1.0
4
Available from: https://github.com/a1184182/WDSLib
5
Language: C++
6
Supported System: Windows, MacOS, Linux, Unix
7
Year first available: 2018
8
1 Introduction
9
Hydraulic simulation has been used to model water distribution systems (WDSs) for several decades
10
and is an essential tool for the design, operation, and management of WDSs in industry and research.
11
Hydraulic simulation allows users (1) to optimize WDS network parameters, such as pipe diameters,
12
in a design setting, (2) to calibrate network parameters, such as demand patterns, in a conventional
13
operational setting, (3) to conduct real-time monitoring and calibration of the network elements in a
14
supervisory control and data acquisition (SCADA) operational setting, and (4) to adjust control devices,
15
such as valves, in a management setting. In the design setting and both the above operational settings,
16
repeated hydraulic assessment is required on a network with fixed topology. In the management setting,
17
repeated hydraulic assessment is required on a network with flexible network parameter settings. With
18
ever-increasing network sizes and the need for real-time management using a SCADA system, it is
19
important to have a robust simulation package which can be configured to be maximally efficient whatever
20
the setting.
AC C
EP
TE D
M AN U
SC
RI PT
1
21
In the field of hydraulic simulation, the system of equations can be formulated as a large and sparse
22
non-linear saddle point problem. Saddle point problems are defined to be a function of two variables (say
23
f(q,h)) and the objective is to find a pair of variables (q*,h*) such that the value f(q*,h*) is minimized
24
with respect to the first variable and maximized with respect to the second variable or equivalently,
25
f (q∗, h) < f (q∗, h∗) < f (q, h∗)∀q, h. There are several well-known iteration methods for solving the non-
26
linear saddle point problem. These include: range space methods, a method that operates in the subspace
27
defined by the rows of the unknown-head node-arc incidences matrix, (Global Gradient Algorithm (Todini
28
and Pilati 1988)), Null space methods, a method that operates in the subspace defined by a null-space that
29
is orthogonal to the column of the unknown-head node-arc incidences matrix (Co-Tree flow formulation 2
ACCEPTED MANUSCRIPT variations (Rahal 1995; Elhay et al. 2014)), loop-based methods, a method is based on the sum of head
31
losses around a loop must be zero (Loop flow correction (Cross 1936)), and pre-and-post-processing
32
methods (forest-core partitioning algorithm (Simpson et al. 2014), domain decomposition (Diao et al.
33
2014), network clustering (Perelman and Ostfeld 2011)). Their relative performance in terms of speed,
34
rate-of-convergence, and accuracy depends among other things on the topology of the target network:
35
size of the forest component, the number of network loops, and the density of these network loops.
36
It is difficult to evaluate the impact of these topology factors by only examining the incidence matrix
37
that describes the pipe network connectivity. As a result, the best method to use for a particular
38
network cannot be easily determined a priori. Moreover, extra complexity is introduced when a multi-
39
run hydraulic assessment is required. During a multi-run hydraulic simulation, the elapsed computation
40
time of each method can be broken down into two parts: the components that are only required to be
41
performed once at the very beginning for the same network, called the overhead, and the components that
42
are required to be carried out repeatedly for each separate run until the required number of iterations has
43
been met, called the hydraulic-phase. It is desirable to have a simulation platform, given the different
44
levels of repetition, to implement these alternative algorithms efficiently. Equipped with such a platform
45
a user would be able to easily benchmark the performance of alternative methods on a small number
46
of evaluations for a given network and use that performance to inform the choice of algorithm to use
47
for either a once-off simulation setting or for a multiple simulation setting (such as for an evolutionary
48
algorithm (EA)).
TE D
M AN U
SC
RI PT
30
This work describes an extensible WDS simulation platform called WDSLib. WDSLib is a numerically
50
robust, efficient and accurate C++ library that implements many WDS simulation methods. WDSLib is
51
written using a modular object-oriented design which allows users to easily mix and interchange solution
52
components, thereby enabling users to avoid redundant computations. It has been optimized to use
53
sparse data structures which are oriented to the pattern of access required for each solution method.
54
WDSLib has been validated for accuracy on a range of realistic benchmark water distribution networks
55
against reference implementations and tested for speed. The program accepts the input file formats of
56
the industry standard EPANET2 (Rossman 2000) toolkit and its performance is faster than EPANET2
57
in all tested settings and benchmarks.
AC C
EP
49
58
The remainder of this paper is structured as follows. The next section describes related methodologies
59
and implementations. A general description of the WDS demand-driven steady-state problem is given
60
in the next section. Section 3 presents a mathematical formulation of the network and the solution
61
methods that are used in WDSLib. The tool-kit structure is then given in section 4. This is followed,
62
in section 5, by the toolkit implementation details. Section 6 provides some examples of how the toolkit
63
can be utilized in a simulation work flow. The results are discussed in Section 7. Finally, section 8
64
summarizes the results of this paper and describes future extensions to the toolkit.
3
ACCEPTED MANUSCRIPT
2 Background
66
This section describes related water distribution system network solution methods and implementations.
67
The first sub-section describes solution methods, including those used by WDSLib. This is followed by
68
a description of currently available implementations and compares these with WDSLib.
69
2.1 Related Methods
70
This research considers a water distribution model made up of energy conservation equations and the
71
demand driven model continuity equations. The Hardy Cross method (Cross 1936), also known as the
72
loop flow corrections method, is one of the oldest methods and uses successive approximations, solving for
73
each loop flow correction independently. It is a method that was widely used for its simplicity at the time
74
when it was introduced. More than three decades later, Epp and Fowler (1970) developed a computer
75
version of Cross’s method and replaced the numerical solver with the Newton method, which solves for
76
all loop flow corrections simultaneously. However, this method has not been widely used because of the
77
need (i) to identify the network loops, (ii) to find initial flows that satisfy continuity and (iii) to use
78
pseudo-loops.
M AN U
SC
RI PT
65
The GGA is a range space method that solves for both flows and heads. It was the first algorithm, in
80
the field of hydraulics, to exploit the block structure of the Jacobian matrix to reduce the size of the key
81
matrix in the linearization of the Newton method. The GGA has gained popularity through its rapid
82
convergence rate for a wide range of starting values. This is the result of using the Newton method on
83
an optimizations problem that has a quadratic surface. However, it was reported by Elhay and Simpson
84
(2011) that the GGA fails catastrophically in the presence of zero flows in a WDS when the head loss
85
is modeled by the Hazen-Williams formula. Regularization methods have been proposed by both Elhay
86
and Simpson (2011) and Gorev et al. (2012) to deal with zero flows when the head loss is modeled by
87
the Hazen-Williams formula.
EP
TE D
79
The GGA as it was first proposed, applied only for the WDSs in which the head loss is modeled by the
89
Hazen-Williams formula, where the resistance factor was independent of flow. Rossman (2000) extended
90
the GGA to allow the use of the Darcy-Weisbach formula. It has been pointed out in Simpson and Elhay
91
(2010), however, that Rossman incorrectly treated the Darcy-Weisbach resistance factor as independent
92
of the flow. They introduced the correct Jacobian matrix to deal with this. It has been demonstrated
93
that once the correct Jacobian matrix is used, the quadratic convergence rate of the Newton method is
94
restored. Furthermore, Elhay and Simpson (2011) reported that the GGA does not fail in the presence
95
of zero flows when the derivatives of the Darcy-Weisbach Jacobian matrix are correctly computed for
96
laminar flows.
97
AC C
88
The co-trees flow method (CTM) (Rahal 1995) is a null space method that solves for the co-tree
4
ACCEPTED MANUSCRIPT flows and spanning tree flows separately. The CTM, unlike the loop flow corrections method, does
99
not require the initial flows to satisfy continuity. However, it does require: (i) the identification of
100
the associated circulating graph; (ii) the determination of the demands that are to be carried by tree
101
branches; (iii) finding the associated chain of branches closing a circuit for each co-tree chord; (iv)
102
computing pseudo-link head losses. The reformulated co-trees flow method (RCTM) (Elhay et al. 2014)
103
is also a null space method that solves for co-tree flows and spanning trees flows separately. It represents
104
a significant improvement on the CTM by removing requirements (i) to (iv) above. It uses the Schilders’
105
factorization (Schilders 2009) to permute the node-arc incidence matrix into an invertible spanning tree
106
block and a co-tree block. This permutation reduces the dimension of the Schur complement from being
107
equal to the number of junctions (as in the GGA) to approximately the number of loops in the network.
108
Abraham and Stoianov (2015) proposed a novel idea to speed-up the solution process when using a
109
null space method to solve a WDS network. Their idea exploits the fact that a significant proportion of
110
run-time is spent computing the head losses. At the same time, flows within some pipes exhibit negligible
111
changes after a few iterations. As a result, there is no point in wasting computer resources to re-compute
112
the pipe head losses for the pipes that have little or no change in flows. This partial update can be used
113
to economize the computational complexity of the GGA, the RCTM and their variations.
M AN U
SC
RI PT
98
The forest-core partitioning algorithm (FCPA) (Simpson et al. 2014) speeds up the solution process
115
in the case where the network has a significant forest component. This algorithm permutes the system
116
equations to partition the linear component of the problem, which is the forest of the WDS, from the
117
non-linear component, which is the core of the WDS. It can be viewed as a method that simplifies the
118
problem by solving for the flows and the heads in the forest just once instead of at every iteration. The
119
FCPA reduces the number of pipes, number of junctions, and the dimension of the Jacobian matrix in
120
the core by the number of forest pipes (or nodes).
EP
TE D
114
The graph matrix partitioning algorithm(GMPA) (Deuerlein et al. 2015) exploited the linear relation-
122
ships between flows of the internal trees within the core and the flows of the corresponding super-links
123
after the forest of the network has been removed. This was a major breakthrough. The GMPA permutes
124
the node-arc incidence matrix in such a way that all of the nodes with degree two in the core can be
125
treated as a group. By partitioning the network this way, the network can be solved by a global step,
126
which solves for the nodes with degree greater than two (super nodes) and the pipes which connect to
127
them (path chords), and a local step, which solves for the nodes with degree two (interior path nodes)
128
and pipes connected to them (path-tree links).
129
2.2 Related Implementations
130
EPANET 2 (Rossman 2000) is a widely used WDS simulation package. EPANET 2 implemented the
131
GGA to provide a demand-driven steady-state solution of a WDS. The code for EPANET 2 is in the
AC C
121
5
ACCEPTED MANUSCRIPT 132
public domain, allowing many extensions to be developed. Currently available extensions include: the
133
implementation of a pressure-dependent model (Cheung et al. 2005; Morley and Tricarico 2008; Siew and
134
Tanyimboh 2012; Jun and Guoping 2012) and a real-time simulation capability (Vassiljev and Koppel
135
2015). The EPANET 2 implementation is not explicitly designed to necessarily be easy to understand or ac-
137
commodate alternative solution methods (Guidolin et al. 2010). The elements that are used in EPANET
138
2 are stored by the variables that describe their graph properties. For example, (1) junctions, reservoirs,
139
and tanks are stored as a C struct called Node and (2) all valves, pipes, and pumps are stored as a C
140
struct called Link. The abundant use of global variables limits the reusability and the possibility of the
141
thread-safe design (Guidolin et al. 2010).
RI PT
136
Consequently, it is difficult to cleanly incorporate new solution methods into EPANET 2 in a manner
143
that allows a fair comparison of performance between these methods. Moreover, because there are no
144
clearly defined interfaces for the incorporation of third-party code components in EPANET 2, there is
145
no guarantee that independently authored extensions will be easy to combine with each other.
M AN U
SC
142
In the absence of a popular easy-to-modify WDS simulation platform there is currently no straight-
147
forward means for comparing different solution methods. To date, when new solution methods have
148
been developed they have been compared using different research systems, on different platforms with
149
different implementation languages. This leads to difficulty in comparing methods, limits the reusability
150
of code, and creates a barrier for researchers to reproduce and replicate results. To address these is-
151
sues, an extensible framework is required that allows implementation of new methodologies to be easily
152
incorporated without an adverse impact on the performance of the rest of the system.
TE D
146
To this end, a number of attempts have been made to implement an object-oriented wrapper to
154
encapsulate the EPANET 2 solver (openNet (Morley et al. 2000) and OOTEN(van Zyl et al. 2003)).
155
However, these two systems were focused on providing more flexibility in the processing of input to the
156
core EPANET solver. They did not address any issues relating to the solution process. CWSnet, a C++
157
implementation in object-oriented style, was produced by Guidolin et al. (2010) as an alternative to
158
EPANET 2. In CWSnet, more attention has been given to the hydraulic elements of the WDS network.
159
In addition, CWSNet provides a pressure driven model, and takes advantage of the computing power of
160
the computer’s Graphics Processing Unit (GPU). However, in CSWnet the data structures representing
161
the network are specialized to the solution methods that it uses. These data structures are not easily
162
adapted to work efficiently with the different traversal orders, and graph algorithms used by newly
163
developed solution methods. However, CWSnet still uses the same hydraulic solver and the same linear
164
solver techniques implemented in EPANET 2 (Guidolin et al. 2010).
AC C
EP
153
165
To accommodate the deficiencies referred to above, this paper presents a new hydraulic simulation
166
toolkit WDSlib. WDSlib is coded in C++, and incorporates a number of recently published techniques.
6
ACCEPTED MANUSCRIPT This toolkit offers users the ability to: (i) choose from, or modify, different approaches and implementa-
168
tions of different WDS model analyses, and (ii) extend the toolkit to include new developments. These
169
features have been implemented using fast and modularized code. A focus of attention in this research
170
has been program correctness, robustness and code efficiency. The correctness of the toolkit has been
171
validated against a reference MATLAB implementation. The differences between all results (intermedi-
172
ate and final) produced by the C++ toolkit and the MATLAB implementation were shown to be smaller
173
than 10−10 . In the interest of toolkit robustness, special attention has been paid to numerical processes
174
to guard against avoidable failures, such as loss of significance through subtractive cancellation, and nu-
175
merical errors, such as division by zero. The data structures and code libraries in WDSLib are shared and
176
all implementations have been carefully designed to ensure fairness of performance comparisons between
177
algorithms. WDSLib uses a pluggable architecture where solution-methods, and their accompanying
178
pre-processing and post-processing code are easily substituted. In addition, different numerical linear
179
algebra techniques can be incorporated using a well-defined interface. This concludes the discussion of
180
related work. The mathematical formulations of the solution methods used in WDSLib are presented in
181
the next section.
182
3 General WDS Demand-Driven Steady-State Problem
183
This section describes the general WDS demand-driven steady-state problem. The following starts with
184
the basic definitions and notation, followed by the system equations. Finally, the relevant equations are
185
shown for each of the different solution methods that are implemented in WDSLib. All variables are
186
described in the nomenclature section in Appendix E.
187
3.1 Definitions and Notation
188
Consider a water distribution system that contains np pipes, nj junctions, nr fixed head nodes and nf
189
forest pipes and nodes. The j − th pipe of the network can be characterized by its diameter Dj , length
190
Lj , resistance factor rj . The i − th node of the network has two properties: its nodal demand di and its
191
elevation zi .
SC
M AN U
TE D
EP
AC C
192
RI PT
167
Let q = (q1 , q2 , ....qnp )T denote the vector of unknown flows, h = (h1 , h2 , ....hnj )T denote the vector
193
of unknown heads, r = (r1 , r2 , ....rnp )T denote the vector of resistance factors, d = (d1 , d2 , .....dnj )T
194
denote the vector of nodal demands, el = (el1 , el2 ....elnr )T denote the vector of fixed head elevations.
195
The head loss exponent n is assumed to be dependent only on the head loss model: n = 2 for
196
the Darcy-Weisbach head loss model and n = 1.852 for Hazen-Williams head loss model. The head loss
197
within the pipe j, which connects the node i and the node k, is modelled by hi −hk = rj qj |qj |n−1 . Denote
198
by G(q) ∈ Rnp ×np , a diagonal square matrix with element [G]jj = rj |qj |n−1 for j = 1, 2, ....np . Denote
7
ACCEPTED MANUSCRIPT 199
by F (q) ∈ Rnp ×np , a diagonal square matrix where the j-th element on its diagonal [F ]jj =
200
A1 is the full rank, unknown head, node-arc incidence matrix, where [A1 ]ji is used to represent the
201
relationship between pipe j and node i; [A1 ]ji = −1 if pipe j enters node i, [A1 ]ji = 1 if pipe j leaves
202
node i, and [A1 ]ji = 0 if pipe j is not connected to node i. A2 is the fixed-head node-arc incidence matrix,
203
where [A2 ]ji is used to represent the relationship between pipe j and fixed head node i, [A2 ]ji = −1 if
204
pipe j enters fixed head node i, [A2 ]ji = 1 if pipe j leaves fixed head node i, and [A2 ]ji = 0 if pipe j is
205
not connected to fixed head node i.
206
3.2 System of Equations
207
The steady-state flows and heads in the WDS system are modeled by the demand-driven model (DDM)
208
continuity equations (1) and the energy conservation equations (2):
RI PT
SC
(1)
M AN U
−A1 T q − d = O 209
G(q)q − A1 h − A2 el = O, 210
which can be expressed as
where its Jacobian matrix is
F (q) J = −A1 T
214
−A1 O
(4)
This non-linear system is normally solved by the Newton method, in which q (m+1) and h(m+1) are
EP
213
(3)
and it is sometimes referred to as a nonlinear saddle point problem (Benzi et al. 2005).
repeatedly computed from q (m) and h(m) by
AC C
212
(2)
−A1 q A2 el = 0, − d h O
TE D
G(q) −A1 T 211
d dqj [G]jj qj .
F
(m) (m)
q
−A1 T
(m+1)
(m)
−A1 q −q G = − O h(m+1) − h(m) ||q (m+1) −q (m) || ||q (m+1) ||
(m) (m)
q
(m)
− A1 h
− A2 e l , −A1 T q (m) − d
(5)
(m+1) (m) −h || are sufficiently small. and ||h (m+1) ||h ||
215
until the relative differences
216
3.3 Global Gradient Algorithm
217
Todini and Pilati (1988) applied block elimination to Eq. (5) to yield a two-step Hazen-William solver:
218
Eq. (6) for the heads and Eq. (7) for the flows.
n o h(m+1) = U −1 −nd + A1 T [(1 − n)q (m) − G−1 A2 el ] 8
(6)
ACCEPTED MANUSCRIPT 219
where U = A1 T G−1 A1
q (m+1) =
h i V h(m+1) = −d + A1 T F −1 (G − F ) q (m) − A2 el
221
222
(7)
Later, Simpson and Elhay (2010) proposed
RI PT
220
o 1n (n − 1)q (m) + G−1 (A2 el + A1 h) n
where V = A1 T F −1 A1
223
(8)
224
h i q (m+1) = q (m) + F −1 A1 h(m+1) − F −1 Gq (m) − A2 el
225
as the generalized equations that can be applied when the head-loss is modeled by the Hazen-Williams
226
equation or the Darcy-Weisbach equation. The correct Jacobian matrix with the formula for F , when
227
head loss is modeled by Darcy-Weisbach equation, can be found in Simpson and Elhay (2010). They
228
showed that the use of the correct Jacobian matrix restores the quadratic rate of convergence.
M AN U
SC
(9)
It is important to note that the GGA, as it was originally proposed, solves the entire network by a
230
non-linear solver, and this can include some unnecessary computations which can be avoided by exploiting
231
the structural properties of the WDS graph composition. The methods described below exploit these
232
structural properties to potentially improve the speed of the solution process.
233
3.4 Forest-Core Partitioning
234
Associated with a WDS is a graph G = (V, E), where the elements of V are the nodes (vertices) of the
235
graph G and elements of E are the pipes (links) of the graph G. The graph G can be partitioned into
236
smaller subgraphs with special properties. The special properties that are exploited in WDSLib and their
237
formulations are described in this subsection. The concept of partitioning the WDS network was proposed
238
by Deuerlein (2008) in order to simplify the WDS solution process. Simpson et al. (2014) extended the
239
idea of the network partitioning of Deuerlein (2008) and introduced the forest-core partitioning algorithm
240
(FCPA), which partitions the network into a treed component and a looped or core component. The
241
FCPA starts with a searching algorithm which identifies the forest subgraph, Gf = (Vf , Ef ), in which
242
S ∈ Nnf ×np is the permutation matrix which identifies the pipes in the forest, Ef , as distinct from the
243
pipes in the core , Ec , and T ∈ Nnf ×nj is the permutation matrix which identifies the nodes in the
244
forest, Vf , as distinct from the nodes in the core, Vc , as distinct from the core subgraph, Gc = (Vc , Ec ),
245
in which P ∈ Nnpc ×np is the permutation matrix for Ec and C ∈ Nnjc ×nj is the permutation matrix for
246
Vc .
AC C
EP
TE D
229
9
ACCEPTED MANUSCRIPT 247
The flows of the pipes in the forest, Sq, can be found directly from
248
−1 Sq = − T A1 T S T T d.
249
Finally, once the iterative solution process for the core has stopped, the forest heads can be found by
250
solving a linear system:
RI PT
(10)
−1 SA2 el − SGS T Sq + SA1 C T Ch . T h = −SA1 T T
251
(11)
3.5 Reformulated Co-Tree flows method
253
A graph, with or without forest, can be partitioned into two sub-graphs: a spanning tree subgraph and
254
a complementary co-tree subgraph. The reformulated co-tree flow method (RCTM) (Elhay et al. 2014)
255
exploited the relationship between the spanning tree pipes and the co-tree pipes. The RCTM starts with
256
a spanning tree search algorithm which identifies a spanning tree subgraph, Gst = (V, Est ), in which
257
K1 ∈ Nnpst ×np is the permutation matrix that identifies the pipes in the spanning tree, Est , as distinct
258
from the pipes in the co-tree, Ect . R is the permutation matrix for the nodes which traverse the same
259
sequence as the corresponding spanning tree pipes, Est , and K2 ∈ Nnpct ×np is the permutation matrix
260
for the pipes in the complementary co-tree edges, Ect . It is important to note that there are many
261
choices of spanning tree for any cyclic graph. The choice of spanning tree and co-tree combination does
262
not affect the correctness of the method.
264
M AN U
TE D
263
SC
252
By exploiting the relationship between the spanning tree and cotree, Elhay et al. (2014) proposed
EP
(m+1)
q1
265
(m+1)
and second for the co-tree flows q2
AC C
266
(m+1)
267
268
269
270
(m+1)
the following equations to solve the WDS for the flows: first for the spanning tree flows q1
W (m+1) q 2
(m)
= L21 T q 2
b − R−T 1 d
272
273
(12)
:
(m) (m) (m+1) (m+1) (m+1) (m) − G1 q1 + F 2 − G2 q 2 + a2 = L21 F 1
(13)
(m) (m) b b where: R1 = K1 Aˆ1 RT ; R2 = K2 Aˆ1 RT ; L21 = −R2 R−T = K1 F (m) K1 T ; F 2 = K2 F (m) K2 T 1 ; F1 (m) (m) (m) (m) b b G1 = K1 G(m) K1 T ; G2 = K2 G(m) K2 T ; W (m) = L21 (F 1 )−1 L21 T + (F 2 )−1 ; a1 = K1 Aˆ2 el ;
a2 = L21 K1 Aˆ2 el + K2 Aˆ2 el . (0)
271
,
Note that in Eq. (12), an initial set of the co-tree flows q2
is needed to commence the solution
process. The heads are found after the iterative process of the RCTM has been completed by using a linear
10
ACCEPTED MANUSCRIPT 274
solution process: (m+1)
R1 h = F 1 q 1
(m)
− (F 1 − G1 ) q 1
− a1
(14)
This partitioning of the network equations reduces the size of the non-linear component of the solver
276
to np − nj (the number of co-tree elements in the network). It has been proven by Elhay et al. (2014)
277
that the RCTM and the GGA have identical iterative results and solutions if the same starting values
278
are used. However, for the RCTM, the user only needs to set the initial flow estimates for the co-tree
279
pipes, q2 , in contrast to GGA where initial flow estimates are required for all pipes. The flows in the
280
complementary spanning pipes are generated by Eq.(12) in the RCTM.
281
4 WDSLib Structure
282
WDSLib is a WDS simulation toolkit consisting of a set of C++ member functions, which henceforth
283
will be referred to just as functions, that can be composed to solve for the steady state solution of a
284
WDS. WDSLib can be used for a once-off simulation or a multi-run simulation. Pre-packaged driver
285
code is provided to perform once-off simulations using a choice of solver methods. For a multi-simulation
286
setting, where the use-cases are very diverse, the user is able to select the desired components of WDSLib
287
to compose and compile their own driver.
M AN U
SC
(0)
RI PT
275
Individual functions in WDSLib are classified according to their role in the simulation workflow. In
289
any simulation workflow, there will be functions that will only have to be executed once. For example,
290
functions to read the input file or partition the network will only have to execute once at the start of the
291
simulation (or of all simulations). Likewise, code to reverse the network partitioning and write simulation
292
results will only have to execute once at the end of the simulation. In this work, these functions that are
293
only required to be run once are called level one (L1 ) functions. L1 functions relate to network topology,
294
which is invariant for the whole simulation. In a multi-simulation setting, certain functions will need to
295
be run once for every hydraulic-phase. An example of such a module is the module making the initial
296
guesses of pipe flow rates for the updated network configuration. In this work, these, once-per-assessment
297
functions, are called level two (L2 ) functions.
AC C
EP
TE D
288
298
Finally, for every hydraulic assessment there is a non-linear iterative phase in the solution process.
299
The functions in this phase run many times for each hydraulic assessment until the stopping test has
300
been satisfied. Examples of these include the functions to calculate the G and F matrices (see Eqs. (3)
301
and (4)) and running the Cholesky solver. These iterative-phase functions are called level three (L3 )
302
functions.
303
Fig. 1 illustrates the global structure of WDSLib under a once-off simulation setting and a multi-run
304
simulation setting. The modular setup of WDSLib allows each module to be run the minimum number
305
of times determined by its simulation setting. Under the module structure described above a once-off
11
ACCEPTED MANUSCRIPT simulation setting can be viewed as a special case where the L1 functions and L2 functions are both run
307
once. Note that after running the initial L1 functions it is possible to run hydraulic assessments of the
308
network in parallel. This mode of execution might be used in a design setting such as using a genetic
309
algorithm (GA) to optimize pipe diameter sizes.
TE D
M AN U
SC
RI PT
306
Figure 1: Global structure of WDSLib for both simulation settings
L1 and L2 functions are classified into parts a and b according to whether they run before or after
311
the lower level processing that they embed. These functions are detailed in Fig. 2. The L1 functions
312
that run at the start of the simulation are called L1a functions. These include the module (1) to read
313
the configuration file; (2) to parse the EPANET .inp file using EPANET programmer’s toolkit; (3) to
314
partition the network; and (4) to solve the linear part(s) of the network. The corresponding L1b functions
315
are run at the end of the simulation. These include tasks such as reversing the network partitioning.
316
Note that certain L1a functions require their corresponding L1b functions to be used. For example the
317
forest search module needs to be paired with the reverse FCPA permutation. There is a similar structure
318
for L2 functions. L2a functions are run at the start of each hydraulic assessment and L2b functions run
319
at the end. The functions that must be included for the FCPA method are denoted with single asterisks.
320
Likewise the functions that must be included with the RCTM method are denoted with double asterisks.
321
For these methods to work correctly all affiliated functions must be included in the simulation workflow.
322
Note that it is also possible to run both the RCTM and FCPA in the same workflow. Also note that
323
the user cannot run both GGA and RCTM in the same workflow – the user must choose between these
AC C
EP
310
12
ACCEPTED MANUSCRIPT Table 1: Key function descriptions, names, their classes, inputs and outputs. The affiliated functions are shown in sub-tables (1a) (1b) (1c) (1d). (a) Shared Modules
Module name readConfig getInputData scale AMD getRf init getGF stopTest rScale
Class Simulation Input Solver Suitesparse Solver Solver Solver Solver Solver
Input config file name EPANET .inp file void void net diameter net, resistance constant result void
(b) Global gradient algorithm (GGA)
Module name runH
Class GGASolver(Solver)
Input void
Output void
SC
Description GGA Solver
Output void EPANET err code void void resistance constant flow rate void norm void
RI PT
Description Read the configuration file Read EPANET input file Variables scaling AMD bandwidth reduction Calculate the resistance constants Generate initial guesses of flows Calculate the head loss coefficients Stopping test Recover scaled variables
(c) Forest-core algorithm (FCPA)
Module name forestSearch forestFlow forestHead rFCPA
Class topology Solver Solver Solver
Input SN, EN demands result void
M AN U
Description Forest search Calculate flows in forest Calculate heads in forest Reverse FCPA permutation
Output void flows in forest pipes heads in forest pipes void
(d) Reformulated cotree flows method (RCTM)
324
Module name STSearch runH RCTMHead rRCTM
Class topology RCTMsolver (Solver) RCTMsolver RCTMsolver
Input SN, EN void flows in ST and CT void
Output void void void void
TE D
Description Spanning tree search RCTM solver Calculate heads in ST and CT Reverse RCTM permutation
solution methods.
Table 1 provides a mapping from the function descriptions in Fig. 2 to the function names in WDSLib.
326
In addition, the dependencies between functions for each solution method are shown in Table 1a, Table
327
1b, Table 1c and Table 1d. The columns in each table list, respectively, the description of the function,
328
its name in WDSLib, the C++ class in which it appears, its input parameters, and its output values.
329
Note, that void is used in these latter two columns to denote that the function interacts with the class
330
variables rather than through its parameters and return value. Examples of how these functions can be
331
coded are presented in section 6. The key data-access functions in WDSLib are described next.
332
Getter and Setter methods
333
parameters and retrieving the results of the WDS network. These methods allow the user to reconfigure
334
the network before and during simulation runs. The names of the setter methods all start with a prefix
335
set and the names of the getter methods all start with a prefix get. For example, a user can set (write-to)
336
the diameter of pipe index to value by calling pipe->setD(index,value) and get (read-from) the head
337
of node index by calling h[index]=result->getHsol(index). A summary of the variables that can be
AC C
EP
325
Each class in WDSLib has various methods available for setting the network
13
AC C
EP
TE D
M AN U
SC
RI PT
ACCEPTED MANUSCRIPT
Figure 2: Function classification in WDSLib. The functions marked with single asterisks must be used for the FCPA method. The functions marked with double asterisks must be used for the RCTM method. Note that it is possible to use both methods at the same time.
14
ACCEPTED MANUSCRIPT Table 2: The getter and setter functions of each class and the variables they access Description Basic network properties, & Pipe and Node Node properties Pipe properties Flag information Parameter information Manage hydraulic simulation Parent class of solution methods GGA solution method RCTM solution method Network topology information Results of the simulations
Read-Access Write-Access Node,Pipe, np , nj , ns d, zs , zu SN , EN , L, D, R, pipe ID getFlag(”flagN”,flagV) setFlag(”flagN”) getPara(”paraN”,paraV) setPara(”paraN”) Result Result getCore, getForest qIter, hIter, GIter , FIter, numIter, CresIter, EresIter, Time
RI PT
Class Name Net Node Pipe Flag Parameter Simulation Solver GGASolver RCTMSolver Topology Result
read-from (read-access through getter methods) and written-to (write-access through setter methods for
339
each key classes is specified in Table 2. This concludes the discussion of the the broad structure of the
340
WDSLib package. The next section describes key aspects of the implementation of the package.
341
5 WDSLIB: Toolkit Implementation
342
This section outlines key implementation details of WDSLib. As previously mentioned, the overall aim of
343
WDSLib is to provide a clearly-structured, flexible and extensible hydraulic simulation toolkit that allows
344
testing, evaluation, and use, in production settings, of both existing and new WDS solution techniques.
345
These aims require WDSLib be implemented so that it is fast to execute, flexible to configure, robust to
346
challenging input data cases, and easy to understand and modify. The following describes aspects of the
347
implementation of WDSLib that enable it to meet these requirements. The next subsection describes the
348
general considerations that informed the design of the whole toolkit. This general discussion is followed
349
by a summary of key improvements to the solution processes encoded in forest searching and spanning
350
tree searching in the WDSLib package.
351
5.1 General capabilities and properties
352
This sub-section describes design aspects underpinning the utility and performance of WDSLib. In-
353
turn, the following outlines measures taken to: (1) maximize code clarity and modularity; (2) increase
354
the efficiency of memory access and storage; (3) maximize numerical robustness; (4) facilitate accurate
355
timing of code execution; and (5) maximize simulation speed for different settings.
356
Design Considerations 1: Modularity
357
The modular design of WDSLib is central to the evaluation and testing of different WDS solution
358
methods. All methods have been defined to perform a single, well–defined, function and each class can
AC C
EP
TE D
M AN U
SC
338
15
ACCEPTED MANUSCRIPT be compiled, used and tested independently. These features allow users to assemble the methods of
360
interest from independently developed components to create a customized WDS solution method in a
361
reliable way. WDSLib’s modular design also allows the users to profile the computation time of each
362
individual component of an algorithm. Functions communicate through well-defined interfaces and the
363
function code has been factored to minimize development and testing cost. This architecture allows
364
customized simulation applications (i) to combine the functions of interest and (ii) to implement new
365
solution algorithms to extend the functionalities of WDSLib.
366
Design Considerations 2: Memory Considerations
367
Care was taken to minimize the memory footprint of executing code (in order to reduce memory re-
368
quirements and prevent memory leaks) in the interest of the toolkit efficiency and toolkit robustness.
369
Reducing memory requirements allows the solution of larger WDS problems for a given memory ca-
370
pacity. In WDSLib, memory reduction was achieved through both, using sparse matrix representations
371
and the systematic allocation and deallocation of working structures in the C++ code. The matrices
372
used in WDS simulation are often sparse, with the density of the full node-arc incidence matrix being
373
only 2/nj . Consequently, it is more efficient to store these matrices using sparse storage schemes which
374
store only the non-zero elements of the matrix and pointers to their locations (Davis et al. 2014). It is
375
important to note that the choice of a sparse matrix representation is made based on (1) the storage
376
requirements of the matrix and (2) common search orders to column elements and row elements. This
377
latter factor means that the best format for sparse matrix representation varies with the preponderant
378
orders of search, (row-wise, column-wise, or both), employed by each method. There is a number of
379
common storage formats for sparse matrices (Compressed column storage (CCS) of Duff et al. (1989)),
380
Compressed row storage (CRS), Block Compressed column storage (BCCS), Block Compressed row stor-
381
age (BCRS), and Adjacency lists). As will be described shortly, WDSLib, uses a modified adjacency-list
382
representation.
AC C
EP
TE D
M AN U
SC
RI PT
359
383
Other implementations use a variety of storage schemes. In EPANET 2, the A1 matrix is stored
384
as two arrays of node indices, which represent start nodes (SN ) and the end nodes (EN ) of each pipe.
385
The i − th entry of the SN and EN arrays represent the start node and end node of i − th pipe of the
386
network. This storage format minimizes the memory required to store the A1 matrix because only the
387
indices are required to be stored because [A1 ](i,SNi ) = −1 and [A1 ](i,ENi ) = 1. As shown in Table 4,
388
searching through rows (pipes) of matrices that are stored in this format is efficient. However, searching
389
though the columns (nodes) is relatively inefficient. This storage format is also used in CWSnet.
390
Both CCS and CRS are used in the FCPA implementation reported in Simpson et al. (2014), and the
391
RCTM implementation reported in Elhay et al. (2014). The partial update null space method (Abraham
392
and Stoianov 2015) used CCS. The memory requirement for storing the A1 matrix in CCS is 2 × nnz +
16
ACCEPTED MANUSCRIPT 393
nj + 1 as shown in Table 4. This storage scheme is fast for searching through columns (nodes) of matrices
394
that are stored in CCS and slow for searching though rows (pipes). Table 3: The adjacency-list matrix presentation
nj
adjacent to {(vi , ej )|vi ∈ N (v1 ) ej connects v1 and vi } {(vi , ej )|vi ∈ N (v2 ) ej connects v2 and vi } {(vi , ej )|vi ∈ N (v3 ) ej connects v3 and vi } .. .
Size Deg(v1 ) Deg(v2 ) Deg(v3 ) .. .
RI PT
Node Index 1 2 3 .. .
{(vi , ej )|vi ∈ N (vnj ) ej connects vnj and vi } Deg(vnj ) Note: (1) ej is negative when vi is the start node and (2) ej is positive when vi is the end node
In WDSLib, a modified adjacency list, described in Table 3, tailored for WDS hydraulic simulation,
396
is used. An adjacency list for an undirected and unweighted graph consists of nj unordered lists for each
397
vertex ni , which contains all the vertices to which vertex ni is adjacent. The network that is shown in
398
Fig. 3 has one source, three nodes, and four pipes. The adjacency list for this network can be described
399
by four lists {{2, 3}, {1, 4}, {1, 4}, {2, 3}}. Each list describes the set of adjacent vertices of a vertex in
400
the graph. For example, the first list, {2, 3}, represents that the vertex 1 is adjacent to the vertex 2 and
401
vertex 3.
AC C
EP
TE D
M AN U
SC
395
Figure 3: A simple sample network. Numbers denote junction and pipe indices in the network.
402
The adjacency list is modified to include a directed and weighted graph for WDSLib. This mod-
403
ified adjacency list for a directed and weighted WDS graph consists of nj unordered lists for each
404
vertex ni . This list contains all the vertex and edge pairs to which vertex ni is adjacent. For ex-
405
ample, the adjacency list for the same network that is shown in Fig. 3 can be described by four lists
406
{{(2, 1), (3, −4)}, {(1, −1), (4, 2)}, {((1, 4), (4, 3)}, {(2, −2), (3, −3)}}. Each list represents the set of adja-
407
cent vertex and edge pair of a vertex in the graph. For example, the first list, {(2, 1), (3, −4)}, describes
408
that the vertex 1 is adjacent to the vertex 2 by edge 1 (from vertex 1 to vertex 2) and the vertex 3 by
409
edge 4 (from vertex 3 to vertex 1). It is fast to search through both the rows and columns of the A1 17
ACCEPTED MANUSCRIPT Table 4: Different sparse representations for A1
CCS CRS EPANET WDSlib
410
size(A1 )
size(A2 )
size([A1 A2 ])
Column Search
Row Search
2 × nnz + nj + 1 2 × nnz + np + 1 -
2 × nnz + nf + 1 2 × nnz + np + 1 -
4 × np + nn + 2 6 × np + 2 2 × np 4 × np
O(n) O((np )n) O(n) O(n)
O((nj )n) O(n) O((nj )n) O(n)
RI PT
Types
matrices that are stored in this format.
411
In addition to these optimized encodings, both G and F are diagonal square matrices, which require
412
less storage when stored as vectors than in sparse matrix form. The storage methods used for the
413
variables in WDSLib and their associated memory usage are given in Table 5.
As a final note, to offer further assurance of the correctness of memory management in WDSLib,
415
Valgrind (Nethercote and Seward 2007), a programming debugging tool, was deployed during testing to
416
detect any memory leaks, memory corruption, and double-freeing.
M AN U
SC
414
Table 5: Vectors and matrices in WDSLib type vector vector matrix matrix matrix
size np × 1 nj × 1 np × np np × nj (np − nj ) × nj
storage method vector vector vector sparse matrix sparse matrix
memory requirements np × double nj × double np × double (2 × np )× integer ≤ (np − nj ) × nj × integer
TE D
variables q, L, D, r h, d G, F A1 , A2 L21
Design Considerations 3: Numerical Considerations
418
The calculations in WDSLib are performed in C++ under IEEE-standard double precision floating point
419
arithmetic with machine epsilon mach = 2.22 × 10−16 . Invariant terms and parameters in every equation
420
were evaluated in advance and replaced by full 20-decimal digit accuracy constants. Intermediate results
421
of calculations, (which are not easily accessible in EPANET), can be output at the user’s request. The
422
stopping tolerance and stopping test can be set by the user either through the configuration file or by
423
the relevant setter method in the Parameter class.
AC C
EP
417
424
In the construction of any numerical solver, there are two primary dangers that are associated with
425
floating point arithmetic that cannot be ignored: (i) subtractive cancellation and (ii) overflow and
426
underflow. To avoid problems associated with these, all input variables are scaled to a similar range
427
to minimize the risk of avoidable computational inaccuracy or failure in floating point arithmetic. It is
428
important to note that unscaled or poorly scaled variables can unnecessarily confound a computation.
429
These scaled input variables are physically dimensionless, which allows computation which is independent
430
of the system of measurement units.
18
ACCEPTED MANUSCRIPT Table 6: WDS variables and units Variables Length Diameter Nodal head Source elevation flow demand G, F
SI unit m m m m m3 /s m3 /s s/m2
US Customary unit ft ft ft ft f t3 /s f t3 /s s/f t2
Scaling factor L0 = max (L) D0 = max (D) h0 = max (el ) el0 = max (el ) q0 = max (d) d0 = max (d) L0 n−1 G0 = D p |q0 |
RI PT
0
The variables, that are provided in EPANET input file for the package, and their corresponding units
432
in US customary and SI units are shown in Table 6. As with the input variables, the system equations
433
were modified to use dimensionless variables. Once the stopping test is satisfied, the original variables
434
are then recovered by reversing the initial scaling. Details of the scaling are shown in Appendix A.
SC
431
To help ensure that WDSLib solution methods are both fast and reliable. The sparse matrix oper-
436
ations are implemented using SuiteSparse (Davis et al. 2014). SuiteSparse is a state-of-the-art sparse
437
arithmetic suite with exceptional performance, from which the approximate minimum degree permuta-
438
tion (AMD) and the sparse Cholesky decomposition routines have been used.
439
Design Considerations 4: Timing Considerations
440
When executing WDSLib, each function reports the time spent in it by sampling wall clock time at
441
the start and end of its execution. Although the overhead for sampling wall clock time is small, there
442
are at least two special considerations involved in the interpretation of these timings: (i) the operating
443
system, at its own discretion, may launch background processes (for example anti-virus software), which
444
will distort the timings and (ii) extrapolating the timing for multiple hydraulic simulations from a single
445
analysis (as may be required, for example, in a genetic algorithm or other evolutionary algorithm run)
446
must be done with care because the relationship between the different settings is not linear.
EP
TE D
M AN U
435
This concludes the discussion of the main considerations concerning the global design of WDSLib. In
448
the following, key details of the implementation of selected parts of the solution processes are described.
449
5.2 Key Improvements to Solution Processes
450
The WDSLib implementation makes several improvements to extant solution processes. This section
451
focuses on the improvement of the network partitioning processes in FCPA and RCTM.
452
Key Optimization 1: Improvements to partitioning in Forest Search
453
The forest-core partitioning algorithm (FCPA) in this paper is a substantial improvement over the
454
algorithm of the original paper (Simpson et al. 2014). Specifically the original FCPA algorithm almost
455
always required many sweeps of the columns (nodes) of the A1 matrix in order to reduce the forest
AC C
447
19
ACCEPTED MANUSCRIPT 456
component down to the core component. The refined algorithm exploits the adjacency list representation
457
of the A1 matrix so that the partitioning process is achieved in a single sweep. This improves the speed
458
of the partitioning process from being O(anp ) to O(np + nf ) where a is the depth of the deepest tree-
459
component in the forest. This can lead to substantial time savings in the case when a is relatively
460
large. The pseudo-code for this refined forest search algorithm is shown in Appendix B This algorithm
462
traverses each tree component in turn from its leaf nodes, which maximizes the locality of operations
463
with respect to the graph representation. In this algorithm, a node is identified as a leaf node when
464
its node degree is one. Every time a leaf node, node k, is identified, the node pointer is moved to its
465
adjacent node, node k, and the node degree of node k is reduced by one. This process repeats if the
466
adjusted node degree of node k is one. Otherwise, node k is the root node for this tree and the algorithm
467
progresses to the next tree in the forest.
468
Key Optimization 2: Improvements to Spanning Tree Search
469
The reformulated co-tree flows method (RCTM) in this paper is also a substantial improvement over
470
the algorithm of the original paper (Elhay et al. 2014). The original spanning tree search algorithm
471
sweeps the rows of the A1 matrix (pipes) in order to identify the singleton rows and their corresponding
472
columns. The spanning tree search in the original RCTM required a sweep of the A1 matrix to identify
473
the next pipe in the spanning tree. This algorithm is O(np nj ), which is relatively inefficient.
TE D
M AN U
SC
RI PT
461
The pseudo-code for the refined spanning tree/forest search algorithm is shown in Appendix C. This
475
improved algorithm takes as input the adjacency list describing the network and the pipe indexes of
476
the core component of the network from the Algorithm 1 (if the FCPA is used). In this algorithm, all
477
water sources are the starting point of the search process, SN , and marked as visited. The nodes in SN
478
are then used as to identify a spanning tree within the WDS. This is achieved by repeatedly finding all
479
adjacent pairs, node t and pipe s, of and removing the first node in SN by using the adjacency list. If
480
the adjacent node t is not visited then node t is inserted into the spanning-tree node vector, ST N , and
481
search node vector, SN , and node t is marked as visited and pipe s to the spanning-tree pipe vector,
482
ST P , and pipe s is marked as visited. If the adjacent node t is visited and the pipe s is not visited
483
then the pipe s is inserted into the co-tree pipe vector, CT P and mark pipe s as visited. This process
484
is repeated until SN is empty. The overall time-complexity of this algorithm is O(np + nj ) (compared
485
to O(np nj ) as mentioned above) is the same as the best asymptotic complexity of breadth-first search
486
on a graph.
AC C
EP
474
20
ACCEPTED MANUSCRIPT
6 Example Applications
488
WDSLib consists of a collection of functions which can be used either as a standalone application for fast
489
one-off simulations or as a library of software components that can be integrated into a user’s own WDS
490
solution processes. This section presents two example applications. The first application is the setup for
491
a basic one-off simulation of a WDS. The second application (described in subsection 6.1) presents an
492
example using WDSLib to implement a simple 1+1 Evolutionary Strategy (Beyer and Schwefel 2002)
493
(1+1-ES or, more commonly, 1+1EA) for sizing pipes in a WDS.
494
Example 1 - Once-off Simulation
495
The setup for WDSLib as a standalone application is straightforward. The user provides a configuration
496
text file that specifies input and output filenames; the name of the solver; the desired output variables;
497
and simulation parameters. These values have sensible defaults so the user can set up the solver by using
498
a minimal configuration such as that shown in Fig. 4. By using this config file, WDSLib is configured to
499
run a single hydraulic analysis of the network that is stored as say ”hanoi.inp”, an EPANET- formatted
500
input file, under ”Network/” sub-directory, using the reformulated co-tree flows method with the forest-
M AN U
SC
RI PT
487
core partitioning algorithm. The full set of configuration parameters for once off simulations is shown in
EP
TE D
[ InputFile ] < directory > % the input file d i r e c t o r y Network / < file > % the input file name hanoi . inp
AC C
1 2 3 4 5 6 7 8 9 10 11 12 13 14
Figure 4: A minimal configuration file to run the GGA in WDSLib
501
502
Fig. 11 in Appendix D.
503
6.1 Example 2 - A Simple Network Design Application
504
As a minimalist example of the application of WDSLib to a WDS network design problem, the following
505
example uses 1+1EA for optimally sizing pipe diameters. This algorithm takes an existing network
506
with randomly generated pipe diameters and optimizes the network to minimize cost, subject to given
507
pressure head constraints. A 1+1EA is a very simple evolutionary strategy (Beyer and Schwefel 2002)
508
which starts with a randomly generated individual (in this case a WDS diameter configuration). This
21
ACCEPTED MANUSCRIPT
M AN U
SC
RI PT
# include " Simulation . h " # include " result . h " using namespace std ; # define GTN 100000 // the total number of g e n e r a t i o n s /* a v a i l a b l e pipe d i a m e t e r s ( inches ) */ vector < int > ADiameter ={36 ,48 ,60 ,72 ,84 ,96 ,108 ,120 ,132 ,144 ,156 ,168 ,180 ,192 ,204}; /* dollar per unit length */ vector < int > unitCost ={93.6 ,133.7 ,176.3 ,221 ,267.6 ,315.8 ,365.4 ,416.5 ,468.7 , 522.1 ,576.6 ,632.1 ,688.5 ,745.1 ,804.1}; vector < int > g e n e r a t e i n i t i a l D i a m e t e r s ( int ) ; // see fig . 6 vector < int > mutate ( vector < int >) ; // see fig . 7 double evaluate ( Net * , Result *) ; // see fig . 8 int main ( int argc , char * argv []) { srand (1) ; char * config = argv [1]; Result * result = new Result () ; Simulation * simulation1 = new Simulation () ; // i n i t i a l i z e the s i m u l a t i o n class Net * net = simulation1 - > L1a ( result , config ) ; // perform the L1a f u n c t i o n s vector < int > p1 = g e n e r a t e i n i t i a l D i a m e t e r ( net - > getNp () ) ; // initial guesses of d i a m e t e r vector < int > p2 ; // work space for current best double cost , currentbest ; for ( int i =0; i < GTN ; i ++) { simulation1 - > setD (& p1 ,& ADiameter , net - > getPIPE () -> Diascale () ) ; simulation1 - > solve ( result ) ; // perform the L2 and L3 f u n c t i o n s simulation1 - > L1b ( result ) ; // reverse the p e r m u t a t i o n cost = evaluate ( result , net ,& p1 ) ; // e v a l u a t e the network cost if ( cost < currentbest || i ==0) { // s e l e c t i o n o p e r a t o r currentbest = cost ; p2 = p1 ; cout < dispResult ( result ) ; delete simulation1 ; delete result ; return 0; }
TE D
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38
Figure 5: Example code for 1+1EA for Pipe Sizing 1+1EA then progresses by applying a mutation to a random pipe diameter size, and then evaluating the
510
new individual. If the new individual is better it replaces the old network. This process continues in a
511
loop until a given number of evaluations is reached.
513
514
515
516
517
The C++ code for this example is shown in Figs. 5, 6, 7, and 8. If the name of the file containing this code is: simpEA.cc then the simplest command to compile this code is:
AC C
512
EP
509
g++ simpEA.cc -o simpEA -Llib -lWDSLib To run this code the user would type: ./simpEA config.txt where config.txt contains the same configuration text as for the previous example.
518
Starting with the main function in Fig. 5, line 15 points to the config file specified by the command
519
line. The next two lines initialize the result and the simulation according to the configuration file.
520
This is followed by the L1a module to perform the user selected L1a functions. Line 19 generates the
521
initial pipe diameters of the network and line 20 initializes the workspace for the mutated string. Line 22
ACCEPTED MANUSCRIPT
vector < int > g e n e r a t e i n i t i a l D i a m e t e r s ( int np ) { vector < int > Diameter = vector < int >( np ) ; for ( int i =0; i < np ; i ++) Diameter [ i ]= rand () % ADiameter . size () ; // set the index for pipe d i a m e t e r s return Diameter ; }
RI PT
1 2 3 4 5 6 7
vector < int > mutate ( vector < int >* string ) { vector < int > string1 ; string1 =* string ; int aa = rand () %( string - > size () ) ; // choose which pipe to mutate int a = rand () % ADiameter . size () ; // choose a pipe d i a m e t e r after the m u t a t i o n ( string1 ) [ aa ]= a ; // set the pipe index return string1 ; }
M AN U
1 2 3 4 5 6 7 8
SC
Figure 6: Implementation code for pipe size initialization
double evaluate ( Result * result , Net * net , vector < int > * p1 ) { PIPE * pipe = net - > getPIPE () ; double np = net - > getNp () ; double nj = net - > getNj () ; double P =1 e7 ; double cost1 =0; penaltyCost =0; vector < double > hsol = result - > getHsol () ; // get the vector of nodal heads vector < double >* L = pipe - > getL () ; // get the vector of pipe lengths vector < double >* zu = net - > getNODE () -> getZU () ; // get the vector of nodal e l e v a t i o n s double Lscale = pipe - > Lenscale () ; // get the scaling factor for length for ( int i = 0; i < np ; i ++) cost1 +=(* L ) [ i ]* Lscale * unitCost [(* p1 ) [ i ]]; // c a l c u l a t e the m a t e r i a l cost for ( int i = 0; i < nj ; i ++) if (( hsol ) [ i ] -(* zu ) [ i ] < minP ) penaltyCost +=( minP -( hsol ) [ i ]+(* zu ) [ i ]) * scaleP ; // c a l c u l a t e penalty cost return cost1 + penaltyCost ; }
AC C
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
EP
TE D
Figure 7: Implementation code for the mutation operator
Figure 8: Implementation code for the evaluation function
23
ACCEPTED MANUSCRIPT 23 sets the pipe diameters of the network. Line 24 evaluates the current network configuration. The
523
permutation and scaling for the current individual is reversed by L1b in line 25 of Fig. 5. Line 26
524
calculates the fitness of the current network configuration by using the evaluate function in Fig. 8. This
525
function applies a penalty for pressure head constraint violations and pipe material costs. The body
526
of the 1+1EA is contained in the selection operator and mutation operator that follow. Lines 27 to 31
527
compare the string in the current generation with the current best string if the individual p1, as measured
528
by evaluate is better than the individual p2 then p1 replaces p2. Line 32 mutates the current network,
529
p2, using mutate (see Fig. 7). The mutate function changes the diameter of a randomly selected pipe in
530
the network to a randomly selected diameter, chosen from a set of commercially available pipe diameters.
531
The mutated individual, stored in the workplace p1, is used as the network configuration for the next
532
iteration. Until the total number of generations is reached, the user selected information about the best
533
individual is outputted by dispResult in line 34 of Fig. 5.
SC
RI PT
522
It should be noted that the algorithm described above can be used to design a simple WDS but is
535
not optimal in terms of speed of convergence. Other EA’s such as genetic algorithms (Simpson et al.
536
1994) will perform better. However the above example has the advantage of simplicity and contains all
537
the basic elements that a GA would use when interacting with WDSLib.
538
M AN U
534
This concludes the presentation of examples in this work. The next section presents a case study that illustrates the performance of WDSLib in a multi-simulation setting.
540
7 Case Study
541
The following presents timing results for WDSLib running the 1+1EA described in the previous section.
542
The results below compare the four different solvers plus EPANET2. Note, that detailed timings for
543
once-off simulations comparing the four methods can be found in Qiu et al. (2018). Three networks
544
were benchmarked in these experiments. These were the N1 , N3 , and N4 case-study networks used in
EP
TE D
539
AC C
Simpson et al. (2014). Table 7 summarizes the characteristics of these networks.
Network N1 N3 N4
Full np 934 1975 2465
Table 7: Benchmark networks summary Network nj ns 848 8 1770 4 1890 3
Forest & Core Networks nf (nf /n# njc n pc p ) 361 (38%) 573 487 823 (42%) 1152 947 429 (17%) 2036 1461
Co-tree Network nct 84 205 757
545
546
Table 8 shows the results of the 1+1EA from Fig. 5 for the GGA, GGA with FCPA, RCTM, RCTM
547
with FCPA and the EPANET2 solvers. For each of the four WDSLib solvers above, the timings are
548
given for running the EA with and without the L1 modules hoisted out the main EA loop. Each
24
ACCEPTED MANUSCRIPT 549
experiment evaluates the WDS network 100,000 times. And the best performing method for each network
550
is highlighted in bold. It is important to note that 1+1EA using both the GGA and the WDSLib Table 8: The actual 1+1EA run-time with 100,000 evaluations (min.) for each of the four solution methods applied networks N1 , N3 , N4 GGA with FCPA min. 4.64 9.79 16.29
RCTM min. 4.53 13.75 23.92
RCTM with FCPA min. 4.13 10.30 21.93
EPANET min. 9.81 26.43 67.11
RI PT
N1 N3 N4
GGA min. 6.73 15.21 21.14
The results show that the EA runs using WDSLib are substantially faster than the runs using the
552
EPANET2 solver. This is, in part, due to the fact that the EPANET2 solver is designed as a standalone
553
solver which does not facilitate lifting out of invariant computations from the EA loop.
SC
551
As a demonstration of how the performance of an EA can be traced Fig. 9 shows the evolution of the
555
fitness values of the N1 network. These traces were extracted from a file written to in line 30 in Fig. 5.
556
As can be seen, the cost and the pressure head violation terms drop during the EA run. Note that there
M AN U
554
will be considerable variation between 1+1EA runs due to its highly stochastic nature.
3.4
10 8
TE D
3.2
3
Cost
2.8
2 10 0
AC C
2.2
EP
2.6
2.4
Total Cost Pressure Violation Penalty
10 1
10 2
10 3
10 4
10 5
Number of evaluations
Figure 9: The evaluation of the fitness value of network N1 557
558
8 Conclusions
559
This paper has described WDSLib, a library for steady-state hydraulic simulation of WDS networks.
560
WDSLib is fast, modular, and portable with implementation of several standard and recently published 25
ACCEPTED MANUSCRIPT hydraulic solution methods. We have outlined the supported solution methodologies, the structure of the
562
package and key aspects of WDSLib’s implementation. Two example applications have been presented
563
including a design case study using a simple EA. The EA results were benchmarked for different solvers
564
demonstrating a substantial improvement in speed over the industry standard EPANET2 package. These
565
benchmarks also have illustrated how the modular structure of WDSLib could be exploited to speed up
566
execution time in a design setting.
RI PT
561
As well as providing a fast simulation platform for both once-off and multi-run simulations, WD-
568
SLib also provides a testbed for comparing different solution methods in different settings and network
569
topologies. As such WDSLib is designed with a pluggable architecture which can extended to efficiently
570
incorporate new solution methods as they are created. This will enhance the capability of the research
571
community to demonstrate the efficacy of new methods without having to re-engineer the content of
572
shared WDSLib functions and data representations.
SC
567
WDSLib accepts the input file format used by EPANET2 which allows for simple deployment. The
574
amount of coding required to use the provided solvers in WDSLib is minimized by the use of a simple
575
and extensible configuration file.
M AN U
573
WDSLib is currently limited to finding solutions of demand-driven WDS systems without control
577
devices. The development of software components to simulate control devices and pressure driven models
578
is future work.
579
References
580
Abraham, E. and I. Stoianov (2015). “Sparse null space algorithms for hydraulic analysis of large-scale
581
water supply networks”. Journal of Hydraulic Engineering, p. 04015058. doi: 10.1061/(ASCE)HY.
582
1943-7900.0001089.
584
585
586
EP
Benzi, M., G. H. Golub, and J. Liesen (2005). “Numerical solution of saddle point problems”. Acta Numerica 14, pp. 1–137. doi: 10.1017/S0962492904000212.
AC C
583
TE D
576
Beyer, H. G. and H. P. Schwefel (2002). “Evolution strategies–A comprehensive introduction”. Natural Computing 1.(1), pp. 3–52. doi: 10.1023/A:1015059928466.
587
Cheung, P. B., J. E. Van Zyl, and L. F. R. Reis (2005). “Extension of EPANET for pressure driven
588
demand modeling in water distribution system”. Computing and Control for the Water Industry 1,
589
pp. 311–316.
590
591
592
Cross, H. (1936). “Analysis of flow in networks of conduits or conductors”. University of Illinois. Engineering Experiment Station. Bulletin; no. 286. Davis, Tim, WW Hager, and IS Duff (2014). “SuiteSparse”. http://faculty.cse.tamu.edu/davis/suitesparse.html.
26
ACCEPTED MANUSCRIPT 593
594
Deuerlein, J. W. (2008). “Decomposition model of a general water supply network graph”. Journal of Hydraulic Engineering 134.(6), pp. 822–832. doi: 10.1061/(ASCE)0733-9429(2008)134:6(822).
595
Deuerlein, J. W., S. Elhay, and A. R. Simpson (2015). “Fast Graph Matrix Partitioning Algorithm
596
for Solving the Water Distribution System Equations”. Journal of Water Resources Planning and
597
Management 142.(1), pp. 04015037-1–04015037-11. Diao, Kegong, Zhengji Wang, Gregor Burger, Chien-Hsun Chen, Wolfgang Rauch, and Yuwen Zhou
599
(2014). “Speedup of water distribution simulation by domain decomposition”. Environmental Mod-
600
elling & Software 52, pp. 253–263. issn: 1364-8152. doi: https://doi.org/10.1016/j.envsoft.
601
2013.09.025. url: http://www.sciencedirect.com/science/article/pii/S1364815213002247.
602
Duff, I. S., R. G. Grimes, and J. G. Lewis (1989). “Sparse matrix test problems”. ACM Transactions on Mathematical Software (TOMS) 15.(1), pp. 1–14.
SC
603
RI PT
598
Elhay, S. and A. R. Simpson (2011). “Dealing with zero flows in solving the nonlinear equations for
605
water distribution systems”. Journal of Hydraulic Engineering 137.(10), pp. 1216–1224. doi: 10 .
606
1061/(ASCE)WR.1943-5452.0000561.
M AN U
604
607
Elhay, S., A. R. Simpson, J. Deuerlein, B. Alexander, and W. H. Schilders (2014). “Reformulated co-
608
tree flows method competitive with the global gradient algorithm for solving water distribution sys-
609
tem equations”. Journal of Water Resources Planning and Management 140.(12), pp. 04014040-1–
610
04014040-10. doi: 10.1061/(ASCE)WR.1943-5452.0000431.
612
Epp, R. and A. G. Fowler (1970). “Efficient code for steady-state flows in networks”. Journal of the
TE D
611
Hydraulics Division 96.(1), pp. 43–56.
Estrada, Carlos, C´esar Gonz´ alez, Ricardo Aliod, and Jara Pa˜ no (2009). “Improved pressurized pipe
614
network hydraulic solver for applications in irrigation systems”. Journal of irrigation and drainage
615
engineering 135.(4), pp. 421–430.
EP
613
Gorev, N. B., I. F. Kodzhespirov, Y. Kovalenko, E. Prokhorov, and G. Trapaga (2012). “Method to Cope
617
with Zero Flows in Newton Solvers for Water Distribution Systems”. Journal of hydraulic engineering
618
139.(4), pp. 456–459. doi: 10.1061/(ASCE)HY.1943-7900.0000694.
AC C
616
619
Guidolin, M., P. Burovskiy, Z. Kapelan, and D. A. Savic (2010). “CWSNet: An object-oriented toolkit for
620
water distribution system simulations”. Proceedings of the 12th Annual Water Distribution Systems
621
Analysis Conference, WDSA, pp. 12–15.
622
Jun, L. and Y. Guoping (2012). “Iterative methodology of pressure-dependent demand based on EPANET
623
for pressure-deficient water distribution analysis”. Journal of Water Resources Planning and Man-
624
agement 139.(1), pp. 34–44. doi: 10.1061/(ASCE)WR.1943-5452.0000227.
625
Morley, M. S., R. M. Atkinson, D. A. Savic, and G. A. Walters (2000). “OpenNet: An applicationin-
626
dependent framework for hydraulic network representation, manipulation, and dissemination”. Hy-
627
droinformatics 2000 Conference, University of Iowa, Iowa City, USA, pp. 23–27.
27
ACCEPTED MANUSCRIPT 628
629
630
631
Morley, M. S. and C. Tricarico (2008). Pressure-Driven Demand Extension for EPANET (EPANETpdd). Tech. rep. University of Exeter. Nethercote, N. and J. Seward (2007). “Valgrind: a framework for heavyweight dynamic binary instrumentation”. ACM Sigplan notices. Vol. 42. 6. ACM, pp. 89–100. doi: 10.1145/1273442.1250746. Perelman, Lina and Avi Ostfeld (2011). “Topological clustering for water distribution systems analysis”.
633
Environmental Modelling & Software 26.(7), pp. 969–972. issn: 1364-8152. doi: https://doi.org/
634
10.1016/j.envsoft.2011.01.006. url: http://www.sciencedirect.com/science/article/
635
pii/S1364815211000259.
RI PT
632
Qiu, M., S. Elhay, A. R. Simpson, and B. Alexander (2018). “A Benchmarking Study of Water Distribu-
637
tion System Solution Methods”. Manuscript accepted for publication to Journal of Water Resources
638
Planning and Management.
640
Rahal, H. (1995). “A co-tree flows formulation for steady state in water distribution networks”. Advances in Engineering Software 22.(3), pp. 169–178. doi: 10.1016/0965-9978(95)00020-W.
M AN U
639
SC
636
641
Rossman, L. A. (2000). “Epanet 2 users manual, us environmental protection agency”. Water Supply and
642
Water Resources Division, National Risk Management Research Laboratory, Cincinnati, OH 45268.
643
Schilders, W. H. (2009). “Solution of indefinite linear systems using an LQ decomposition for the linear
644
constraints”. Linear Algebra and its Applications 431.(3), pp. 381–395. doi: 10.1016/j.laa.2009.
645
02.036.
647
Siew, C. and T. T. Tanyimboh (2012). “Pressure-dependent EPANET extension”. Water resources man-
TE D
646
agement 26.(6), pp. 1477–1498. doi: 10.1061/41203(425)9. Simpson, A. R., G. C. Dandy, and L. J. Murphy (1994). “Genetic algorithms compared to other tech-
649
niques for pipe optimization”. Journal of water resources planning and management 120.(4), pp. 423–
650
443. doi: 10.1061/(ASCE)0733-9496(1994)120:4(423).
EP
648
Simpson, A. R. and S. Elhay (2010). “Jacobian matrix for solving water distribution system equations
652
with the Darcy-Weisbach head-loss model”. Journal of hydraulic engineering 137.(6), pp. 696–700.
653
doi: 10.1061/(ASCE)HY.1943-7900.0000341.
AC C
651
654
Simpson, A. R., S. Elhay, and B. Alexander (2014). “Forest-Core Partitioning Algorithm for Speeding
655
Up Analysis of Water Distribution Systems”. Journal of Water Resources Planning and Management
656
140.(4), pp. 435–443. doi: 10.1061/(ASCE)WR.1943-5452.0000336.
657
Todini, E. and S. Pilati (1988). “A gradient algorithm for the analysis of pipe networks”. Computer
658
applications in water supply: vol. 1—systems analysis and simulation. Research Studies Press Ltd.,
659
pp. 1–20.
660
Vassiljev, A. and T. Koppel (2015). “Estimation of real-time demands on the basis of pressure mea-
661
surements by different optimization methods”. Advances in Engineering Software 80, pp. 67–71. doi:
662
10.1016/j.advengsoft.2014.09.023.
28
ACCEPTED MANUSCRIPT
TE D
M AN U
SC
RI PT
EPANET”. Advances in Water Supply Management (CCWI 2003).
EP
664
van Zyl, J. E., J. Borthwick, and A. Hardy (2003). “Ooten: An object-oriented programmers toolkit for
AC C
663
29
ACCEPTED MANUSCRIPT 665
Appendix A Scaling
666
In WDSlib, all input variables are scaled to a similar range to minimize the risk of avoidable computa-
667
tional inaccuracy or failures in floating point arithmetic. ˆ = h/h0 , eˆl = el /el0 , d ˆ = d/d0 where ˆ = G/G0 , q ˆ = q/q0 , h The variables are scaled as following: G
669
ˆ eˆl , ˆ h, G, h, el , and d are the original input vectors, G0 , h0 , el0 , and d0 are the scaling factors and G,
670
and dˆ are the scaled input vectors.
671
RI PT
668
ˆ 0 , el = eˆl el0 , d = dd ˆ 0 , Eq. (1) and Eq. (2) become: ˆ 0, q = q ˆ q0 , h = hh By substituting G = GG
ˆ − el0 A2 eˆl = 0 ˆ q − h0 A1 h G0 q0 Gˆ
673
ˆ=0 ˆ − d0 d q 0 A1 T q
674
Eq. (15) and Eq. (16) can be further simplified by introducing the notation: a = h0 / (G0 ∗ q0 ), b =
675
el0 / (G0 ∗ q0 ), and c = d0 /q0 :
M AN U
676 677
678
679
with the following matrix form:
680
683
A1 T qˆ − cdˆ = 0
(18)
ˆ bA2 eˆl −A1 q − = 0. ˆ ˆ O ah cd
(19)
Finally, the network can be solved by using Eq. (20) and Eq. (21): ˆ − A1 F ˆ (m+1) = 1 U −1 [−cd ˆ −1 [(F ˆ − G)ˆ ˆ q (m) + bA2 eˆl ]] h a
(20)
ˆ m+1 − F ˆ −1 A1 h ˆ −1 (Gˆ ˆ q (m) − bA2 eˆl ) ˆ (m+1) = q ˆ (m) + aF q
(21)
EP
682
(16)
(17)
AC C
681
(15)
ˆ − bA2 eˆl = 0 ˆ q − aA1 h Gˆ
TE D
ˆ G −A1 T
SC
672
684
Choice of scaling factors
685
The choice of the scaling factor, despite much research, is not well understood. In this subsection, a
686
choice for each scaling factor, based on the experience of the authors, is recommended. There are two
687
types of variables and parameters that need to be scaled: invariants and variants.
688
Data sets that have very wide range of values can confound numerical accuracy. As a result, it may
689
be preferable to scale the data to a narrower range. The default scaling factor for each of the input data
690
is chosen to be its maximum absolute value. For example, the scaling factor for demand is max(d), so
691
that its values range from zero to one.
692
In contrast, it is more difficult to choose a scaling factor a priori for values that vary between iterations 30
ACCEPTED MANUSCRIPT 693
(variants). This is because the range of variants can change as the iteration progresses. As a result, the
694
intermediate and the final results might not be within the same range as the initial guesses.
698
connected to a water source and a good choice of the scaling factor for nodal head is max(el ) because
699
the maximum nodal head cannot exceed the maximum elevation head of the fixed nodes.
696
RI PT
697
There are two variants that need to be scaled: q, h. A good choice of the scaling factor for the P d flow rate is because the demand at each node must be satisfied by the water sources in the WDS nf and it is a reasonable assumption that the all demands are equally carried by each pipe that is directly
695
700
During the process of the computation, the matrices G and F are scaled because their input variables
701
are scaled. For the Darcy-Weisbach head loss model, the diagonal elements of the matrix G are modelled
702
by: [G]jj = diag
703
where the friction factor f is modelled by the Swamee-Jain formula:
M AN U
704
) 8 Lj ( 2 ) 5 fj |qj | for j = 1, . . . , np π g Dj
SC
(
64 Re j fj = P3 (α + β /θ )(Re /2000)k k k j j k=0 0.25
705
log2 (θj )
709
710
711
712
713
if Rej ≥ 4000
TE D
708
(22)
In order to make sure the Reynolds number, a dimensionless variable, is not affected by scaling, D0 4ˆ qj q0 ν ˆ= ν is introduced, Reynolds number Rej becomes , which can be further simplified ˆ j D0 q0 πˆ ν (q0 /D0 )D 4ˆ qj to , where the scaling factors are canceled. ˆj πˆ νD In order to make sure f is also not affected by the scaling, ˆ = D0 is introduced, θj becomes
EP
707
if 2000 < Rei < 4000
5.74 4qj j + and Rej = for j = 1, . . . , np . Note that αk and βk are constant, the 3.7Dj πνD Re0.9 j j values of which can be found in Elhay and Simpson (2011) where θj =
AC C
706
if Rej ≤ 2000
θj =
ˆj D0 5.74 + for j = 1, . . . , np ˆ Re0.9 3.7Dj D0 j
which can be further simplified to
714
θj =
ˆj ˆj 3.7D
+
5.74 for j = 1, . . . , np Re0.9 j
715
where the scaling factors are canceled. It is evident that the friction factors remain the same because
716
the values for the only two variables Re and θ are unchanged.
717
718
ˆ where G0 = ( 8 ) L0 |q0 | and [G] ˆ jj = Finally, diagonal elements of G can be expressed as: G0 G, π 2 g D05 ( ) ˆj L diag f |ˆ q | for j = 1, . . . , np . ˆ5 j j D j
31
ACCEPTED MANUSCRIPT
721
For the Hazen-Williams head loss ) model, the diagonal elements of the matrix G are modelled by: ( 10.67Lj |qj |(n−1) for j = 1, . . . , np , where Cj is the Hazen-Williams coefficient for [G]jj = diag Cj1.852 Dj4.871 the j-th pipe. The Hazen-Williams coefficient, unlike the friction factor in the Darcy-Weisbach head
722
loss model, is independent of flow rate, pipe wall condition, and flow regimes, which means it is an
720
723
724
725
726
independent ( variable. As a result, the scaling factor for)Hazen-Williams G can simply be derived by: ˆ j L0 10.67L [G]jj = diag |ˆ q |(n−1) |q0 |(n−1) and the equation for diagonal elements of G 1.852 1.852 ˆ 4.871 D4.871 j ˆ C0 D Cj 0 j n o ˆ , where G0 = 10.67L0 |q0 |(n−1) for Hazen-Williams equation can be expressed as: G = diag G0 G C01.852 D04.871 ) ( ˆ Lj |ˆ q |(n−1) for j = 1, . . . , np . and [G]jj = diag 1.852 ˆ 4.871 j ˆ D C j
RI PT
719
j
Why scaling should be applied
728
Case 1 from Estrada et al. (2009) is used here to demonstrate the efficacy of using scaled variables. The
729
layout of the Case 1 from Estrada et al. (2009) is given in Fig. 10. Details of the pipe and node values
730
are given in Table 9 and Table 10.
TE D
M AN U
SC
727
EP
Figure 10: The layout of the Case 1 from Estrada et al. (2009)
Table 9: The nodal values of the Case 1 from Estrada et al. (2009)
AC C
Node ID Reservoir A Reservoir B Junc C Junc D
Elevation (m) 10 10 0 0
Demand (LPS) 2.7 0
Table 10: The pipe values of the Case 1 from Estrada et al. (2009)
Pipe ID Pipe 1 Pipe 2 Pipe 3 Pipe 4
731
Start Node Reservoir A Reservoir B Junc D Reservoir A
End Node Junc C Junc C Junc C Junc D
Length (m) 0.25 0.5 100 100
Diameter (mm) 1200 1200 1200 1200
Running this example network from EPANET and WDSLib, we get the following results:
32
ACCEPTED MANUSCRIPT Table 11: The nodal heads found by using EPANET and WDSLib
Node ID Junc C Junc D
EPANET WDSLib Head (m) 10 10 10 10
M AN U
735
In contrast, the energy residual and the continuity residual of the WDSLib result are 1.8723 × 10−12 and 1.7849 × 10−9 respectively which are both acceptable.
TE D
734
respectively. These two residuals are too large which indicates an inaccurate solution.
EP
733
The energy residual and the continuity residual of the EPANET result are 2.6973 and 5.3575 × 10−4
AC C
732
EPANET WDSLib Flow (LPS) 1.3446 1.5924 1.3446 1.5924 0.0108 0.0124 0.0108 0.0124
SC
Link ID Pipe 1 Pipe 2 Pipe 3 Pipe 4
RI PT
Table 12: The flows found by using EPANET and WDSLib
33
ACCEPTED MANUSCRIPT
Appendix B Forest Search Algorithm
EP
TE D
SC
M AN U
k ← 1; for i ← 1 to nj do n = i; while Deg(n) == 1 do for (Adjp ,Adjv ) ∈ Adj(n) do if Deg(Adjv ) == 1 then F orest[k] ←Union(F orest[k], (Adjp , n)) else q(Adjp ) ← d(n); d(Adjv ) ← d(Adjv ) + d(n); Deg(Adjv ) ← Deg(Adjv ) − 1; n = Adjv ; if Deg[n] ≥ 2 then RootN ode[k] = n; k ← k + 1; end if end if end for end while end for
RI PT
Algorithm 1: Forest Search Algorithm input : Adjacency list, d and Deg output: Forest, RootNode, q and d
AC C
736
34
ACCEPTED MANUSCRIPT
Appendix C Spanning Tree Search Algorithm
SC
// mark the source i as visited this node into the search node vector
// mark pipe as a spanning tree pipe as a corresponding spanning tree pipe
EP
TE D
M AN U
for i ← nj to nf do V N [i] = true ; SN ← i ; // insert end for while SN is not empty do i ← the first element in SN ; foreach (Adjp ,Adjv ) ∈ Adj[i] do if VN[Adjv ] == false then ST P ← Adjp ; ST N ← Adjv ; // mark node SN ← Adjv ; VN[Adjv ]=true; VP=[Adjp ]=true; end if else if VN[Adjp ]==false then CT P ← Adjp ; VP[Adjp ]=true; end if end foreach remove i from SN ; end while
35
pipes nodes pipes nodes pipes nodes
RI PT
Algorithm 2: Spanning Tree Search Algorithm input : adjList output: Spanning Tree and Co-Tree elements ST P ← {} ; // initialize an empty vector for spanning tree ST N ← {} ; // initialize an empty vector for spanning tree CT P ← {} ; // initialize an empty vector for co-tree V N ← {} ; // initialize a boolean vector for visited V P ← {} ; // initialize a boolean vector for visited SN ← {} ; // initialize an empty set of searching
AC C
737
// removed pipe j as a singleton pipe
ACCEPTED MANUSCRIPT
Appendix D Complete configuration files
M AN U
SC
RI PT
[ InputFile ] < directory > % R E Q U I R E D I N F O R M A T I O N the input file d i r e c t o r y WDSnetwork / < file > % R E Q U I R E D I N F O R M A T I O N the input file name n1 . inp
TE D
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
EP
Figure 11: A configuration file to run the RCTM in WDSLib
AC C
738
36
ACCEPTED MANUSCRIPT
Appendix E Nomenclature
740
Acronyms
741
CT
Co-tree
742
CTP
Co-tree pipes
743
DW
Darcy-Weisbach head loss formula
744
FCPA Forest-core partitioning algorithm
745
GGA
Global gradient algorithm
746
HW
Hazen-William head loss formula
747
WDS Water distribution system
748
RCTM Reformulated co-tree flows method
749
nnz
Number of non-zeros
750
STP
Spanning tree pipes
751
STN
Spanning tree nodes
752
ST
Spanning tree
753
Constants
754
g
Gravitational acceleration constant
755
ν
Kinematic viscosity of water
756
FCPA variables
757
C
758
Ec
759
Ef
760
Gc
Core subgraph
761
Gf
Forest subgraph
762
P
Permutation matrix for the pipes in the core, Ec , of the network
763
S
Permutation matrix for the pipes in the forest, Ef ,of the network
764
T
Permutation matrix for the nodes in the forest, Ef , of the network
EP
TE D
M AN U
SC
RI PT
739
AC C
Permutation matrix for the nodes in the core, Ec , of the network Set of core pipes (edges) in Gc Set of forest pipes (edges) in Gf
37
ACCEPTED MANUSCRIPT Vc
Set of core nodes (vertices) in Gc
766
Vf
Set of forest nodes (vertices) in Gf
767
Hydraulic variables for GGA
768
A1
Unknown-head node-arc incidence matrix
769
A2
Fixed-head node-arc incidence matrix
770
d
Vector of nodal demands
771
di
Demand of node i
772
D
Vector of pipe diameters
773
Dj
Diameter of pipe j
774
el
Vector of Fixed-head nodes elevation heads
775
elk
Fixed-head nodes elevation heads at node k
776
E
Set of pipes in graph G
777
EN
Vector of end nodes
778
EN j
End nodes of pipe j
779
f
Vector of Darcy-Weisbach friction factors
780
fj
Darcy-Weisbach friction factor of pipe j
781
F
Diagonal Matrix of generalized headloss derivatives when the headloss is modelled by either the
M AN U
TE D
EP
HW and the DW
782
SC
RI PT
765
[F ]jj
Generalized headloss derivatives for pipe j
784
G
785
G
786
[G]jj
rj |qj |n−1
787
h
Vectors of unknown heads
788
hi
Heads at node i
789
J
Jacobian matrix
790
L1a
Functions run once before multiple simulation
AC C
783
Full WDS graph
Diagonal matrix with elements rj |qj |n−1
38
ACCEPTED MANUSCRIPT L2a
Functions run once before hydraulic assessment
792
L3
Functions run every iteration
793
L2b
Functions run once after hydraulic assessment
794
L1b
Functions run once after multiple simulation
795
L
Vector of pipe lengths
796
Lj
Length of pipe j
797
n
Head loss exponent
798
nf
Number of forest pipes and nodes
799
nj
Number of junctions
800
np
Number of pipes
801
nr
Number of fixed-head nodes
802
nst
Number of ST pipes and nodes
803
nct
Number of CT pipes
804
q
Vector of unknown flows
805
qj
Flow in pipe j
806
R
Vector of Reynolds numbers
807
Rj
Reynolds number for pipe j
808
SN
Vector of start nodes
809
SN j
Start nodes of pipe j
810
U
811
V
812
V
Set of node in graph G
813
zi
Elevation at node i
814
αk
Interpolating spline coefficient
815
βk
Interpolating spline coefficient
816
θ
Vector as defined in Eq. (22)
AC C
EP
TE D
M AN U
SC
RI PT
791
Diagonal matrix of Schur Complement when headloss is modelled by HW Generalized Schur Complement when the headloss is modelled by both the HW and the DW
39
ACCEPTED MANUSCRIPT
Vector of pipe roughness heights
818
j
Roughness height for pipe j
819
RCTM variables
820
Est
Set of ST pipes (edges)
821
Ect
Set of complementary CT pipes (edges)
822
K1
Orthogonal permutation matrix for pipes in the ST
823
K2
Orthogonal permutation matrix for pipes in the CT
824
L21
A part of a basis for the null space of the permuted node-arc incidence for the RCTM
825
R
Orthogonal permutation matrix for nodes in the ST
826
Vst
Set of ST nodes (vertices)
827
W
Schur complement for the RCTM
AC C
EP
TE D
M AN U
SC
RI PT
817
40
ACCEPTED MANUSCRIPT
Highlights
EP
TE D
M AN U
SC
RI PT
A library for the steady-state analysis of a water distribution system (WDS) An open-source C++ software implementation of a number of WDS solution methods A fast simulation platform for both once-off and multi-run simulations Several improvements to the existing solution methods have been made
AC C
• • • •