OPTIMIZING OPERATIONS OF LARGE WATER SUPPLY
NETWORKS
A Case Study
Derek Verleye, El-Houssaine Aghezzaf and Dimi Defillet
Department of Industrial Engineering, Faculty of Engineering and Architecture, Ghent University,
Technologiepark 903, Zwijnaarde, Belgium
Keywords:
Drinking water supply, Mixed integer nonlinear programming, Network optimization.
Abstract:
In this paper we propose a mathematical programming model for a large drinking water supply network and
discuss some possible extensions. The proposed optimization model is of a real water distribution network,
the largest water supply network in Flanders. The problem is nonlinear, nonconvex and involves some binary
variables, making it belong to the class of NP-hard problems. We discuss a way to convexify the nonconvex
term and show some results on two case instances of the actual network.
1 INTRODUCTION
Optimization of drinking water networks has recently
regained attention of many researchers in mathemat-
ical programming. These models are typically non-
convex MINLP’s and are in general hard to solve.
They are often solved by means of some NLP solvers,
after the binary variables are relaxed or the model re-
formulated (Cembrano et al., 2000; Burgschweiger
et al., 2009). Heuristic methods are also frequently
used to solve these models (Savic and Walters, 1997;
L
´
opez-Ib
´
a
˜
nez et al., 2008; Nicklow et al., 2010).
In this paper, we extend the model used to op-
timize an existing network managed by a company
in Flanders, Belgium (Verleye and Aghezzaf, 2011).
This extended model is a nonconvex MINLP for
which we linearize the power term. It is then solved
by means of an open-source solver and tested on two
recent cases (Defillet, 2011).
2 MODEL FORMULATION
Before discussing the solution approach, we first pro-
vide an overview of the model. The components and
restrictions discussed here are based on an existing
supply network of a large drinking water supplier in
Flanders, Belgium. We model the network using a
directed graph G = (N ,A), where N is the node
set representing junctions (J ), delivery points (D),
buffers (B) and raw water sources (S); and A is the
arc set representing pipes (P i), pipes with pure wa-
ter pumps (Pu) and pipes with raw water pumps (R ).
Furthermore, a day is divided in periods t ∈ [1, T ],
where T represents the number of periods. The length
of each period is denoted by τ
t
.
The main variables in a water supply network are
the flow values, Q
i j
t
on arc (i, j), and the piezometric
heads in each node i, H
i
t
. Piezometric head is the sum
of the static level (h
i
) and the manometric pressure.
Since time is discretized, we interpret the variables as
follows. The flows are assumed to be constant during
the interval [t − 1,t]. The head values will be inter-
preted as values at time t. The volume in a tank, V
i
t
,
is the available volume at the end of period t.
In our model, the amount to be produced during
each period needs to be determined. Moreover, at
each time period we need to know which pumps to ac-
tivate for the period to come and at which frequency.
2.1 Restrictions on Nodes (i ∈ N )
On nodes, some bounds on manometric pressure ap-
ply. Conservation of flow is written for each junction.
In what follows, we describe constraints on buffers.
2.1.1 Buffers (i ∈ B)
We consider buffers with free inflow (I
i
t
(+)) in the
tank, whereas outflow is possible through high pres-
sure because of the tank level (I
i
t
(−)) or by pumping
469
Verleye D., Aghezzaf E. and Defillet D..
OPTIMIZING OPERATIONS OF LARGE WATER SUPPLY NETWORKS - A Case Study.
DOI: 10.5220/0003742704690472
In Proceedings of the 1st International Conference on Operations Research and Enterprise Systems (ICORES-2012), pages 469-472
ISBN: 978-989-8425-97-3
Copyright
c
 2012 SCITEPRESS (Science and Technology Publications, Lda.)

water out of the buffer (O
i
t
(−)). For more details, we
refer the reader to (Verleye and Aghezzaf, 2011).
Flow conservation:
∑
k:(k,i)∈A
Q
ki
t
−
∑
j:(i, j)∈A\Pu
Q
i j
t
= I
i
t
(+) − I
i
t
(−) (1)
∑
j:(i, j)∈P u
Q
i j
t
= O
i
t
(−) (2)
Tank volume:
V
i
t
= V
i
t−1
+ (I
i
t
(+) − I
i
t
(−) − O
i
t
− d
i
t
)τ
t
(3)
V
i
0
≥ V
i
T
(4)
d
i
t
denotes demand.
Mean piezometric level (M
i
t
) in tank:
M
i
t
= l f
i
+
L
i
t
+ L
i
t−1
2
(5)
where L
i
t
=
V
i
t
A
i
is the (approximated) water level in
the tank, with A
i
the cross-sectional area of tank.
Tank in- or outflow:
I
i
t
(+) ≤ A
i
l
i
(max)X
i
t
(6)
I
i
t
(−) ≤ A
i
l
i
(max)Y
i
t
(7)
H
i
t
− li
i
≥ −100 (1 − X
i
t
) (8)
M
i
t
− H
i
t
≥ −100 (1 −Y
i
t
) (9)
Here, li
i
is the level at which water flows (freely) in
the tank. The binary variables X
i
t
,Y
i
t
are used to model
in- and outflow. If I
i
t
(+) > 0 pressure H
i
t
at the en-
trance has to be higher than the inflow level li
i
. In the
same way, outflow can only occur if the water level
is sufficiently high. The number 100 is derived from
the fact that the manometric pressure cannot exceed
10 bar (100 m) in any part of the network.
2.2 Restrictions on Arcs ((i, j) ∈ A)
2.2.1 Pipes ((i, j) ∈ P i)
Flow bounds ((Bragalli et al., 2006)):
−3600
π
4
v
i j
max
(d
i j
)
2
≤ Q
i j
t
≤ 3600
π
4
v
i j
max
(d
i j
)
2
(10)
d
i j
is the diameter and v
i j
max
the maximum flow speed.
Pressure losses:
H
i
t
− H
j
t
= κ
i j
Q
i j
t
|Q
i j
t
| (11)
which is the formula of Darcy-Weisbach and is suit-
able for turbulent flow in pipes. The loss coefficient
κ
i j
depends on pipe length, diameter and the pipe fric-
tion coefficient. To determine the last one, we work
with the (simplified) law of Prandtl-K
´
arm
´
an for hy-
draulically rough pipes.
2.2.2 Pure Water Pumps ((i, j) ∈ P u)
The following constraints apply on constant speed
pumps. In addition to flow bounds, we can again write
the pressure loss constraints which are slightly differ-
ent here because water can flow only in one direction
and the pump increases the head, ∆H
i j
t
:
Pressure losses:
H
i
t
− H
j
t
− κ
i j
(Q
i j
t
)
2
+ ∆H
i j
t
= 0 ∀i ∈ N \B (12)
M
i
t
− H
j
t
− κ
i j
(Q
i j
t
)
2
+ ∆H
i j
t
= 0 ∀i ∈ B (13)
The second constraints apply to pumps extracting wa-
ter from a tank with level M
i
t
.
Pump characteristic:
Using data provided by pump suppliers we ap-
proximate the relation between head, ∆H, and flow,
Q as
∆H
i j
t
= h
i j
1
(Q
i j
t
)
2
+ h
i j
2
Q
i j
t
+ h
i j
3
(14)
where h are the function coefficients. In the same way,
the efficiency E is given by
E
i j
t
≤ e
i j
1
(Q
i j
t
)
2
+ e
i j
2
Q
i j
t
+ e
i j
3
(15)
with coefficients e. We can relax the equality con-
straint because efficiency will always be chosen to be
as high as possible. The power term P is given by
P
i j
t
= 2, 73
∆H
i j
t
Q
i j
t
E
i j
t
(16)
To avoid numerical difficulties, we force e
i j
3
> 0
so to prevent the efficiency from becoming 0 when
the pump is not working, that is when Q
i j
t
= 0. This
will prevent the denominator of the power term, P,
from becoming 0.
Substituting the approximations for efficiency and
head in (16), we can also draw a curve for the power
based on pump data (see fig 1). We can then approxi-
mate (16) with:
P
i j
t
= p
i j
1
Q
i j
t
+ p
i j
2
(17)
2.2.3 Raw Water Pumps ((i, j) ∈ R )
In addition to treatment capacity and daily extraction
limits in ground water treatments, we include bounds
on subsequent production flow rates:
−q
i j
( f luct) ≤ Q
i j
t
− Q
i j
t−1
≤ q
i j
( f luct) (18)
−2q
i j
( f luct) ≤ Q
i j
T
− Q
i j
0
≤ 2q
i j
( f luct) (19)
These bounds will prevent quality deterioration
due to excessive changing of production amount from
one period to the next.
ICORES 2012 - 1st International Conference on Operations Research and Enterprise Systems
470

0 200 400 600
0
20
40
60
80
100
Q (m
3
/h)
P (kW)
Pump data
uQ
2
+ vQ + w
Linear approximation
Figure 1: Fitted power function using pump data.
2.3 Models
As two different cases will be tested, the models will
be slightly different. The preceding restrictions are
applicable for both models. In what follows addi-
tional restrictions for each model will be described.
2.3.1 Decision Support Model for the Region
West-Flanders
Due to a recent network expansion, the question on
how to optimally satisfy the demand by using existing
sources arises.
The function that needs to be minimized is the to-
tal cost, consisting of energy delivered by the pumps
and production/electricity at the water production
centers:
Min
T
∑
t=1
∑
(i, j)∈Pu
P
i j
t
c
t
(e)
1000
+
∑
(i, j)∈R
Q
i j
t
c
i j
(p)
!
τ
t
(20)
2.3.2 High-level Transport Optimization Model
Due to the recent collaboration of different drinking
water companies, a significant extra amount of drink-
ing water will be flowing through the main supply
network. Here arises the need of a high level decision
model that investigates the impact of these changes on
the operation of the drinking water supply network.
The objective function (20) is extended with
T
∑
t=1
∑
j:(i, j)∈A
Q
i, j
t
−
∑
k:(k,i)∈A
Q
k,i
t
!
ap
i
τ
t
∀i ∈ D
(21)
where ap
i
denotes the price (in e/m
3
) for delivering
water to node i.
At the ground water sources in this region, raw
water pumps work individually and are set to a con-
stant flowrate. Furthermore, switching these pumps
on/off shouldn’t occur frequently. For each group of
raw water pumps we define the set C which contains
all possible combinations of discrete total flowrates.
We avoid adding binary variables and let α
m(i j)
be 1
when the raw water pumps at i j operate at combina-
tion m. Then (for 0 ≤ α ≤ 1):
Q
i j
t
=
∑
m∈C(i j)
(α
m
bc
i, j
) ∀(i, j) ∈ R (22)
where bc
i, j
is the treatment capacity. As a conse-
quence, the amount of pumped water remains con-
stant over time. Of course, only one possible combi-
nation can be chosen:
∑
m∈C(i j)
α
m
≥ 1 ∀(i, j) ∈ R (23)
α
m
α
n
≤ tol ∀m,n ∈ C(i j) : (i, j) ∈ R (24)
where tol is a very small value.
3 COMPUTATIONAL RESULTS
In this section we discuss the most important results
obtained from the two proposed models. Both mod-
els, were implemented in AMPL (Fourer et al., 2003)
where the open-source solver BONMIN (Bonami and
Lee, 2006) was used to generate solutions. This
solver, suitable only for smooth and convex MINLP,
will only provide local solutions for our nonconvex
models. In both cases we encountered some numeri-
cal difficulties due to buffer in and outflow constraints
(6-9). We relax these constraints and allow only in-
flow in buffers. This is not a major issue, however
the solution value may be slightly inferior because we
don’t admit outflow.
3.1 Decision Support Model for the
Region West-Flanders
This large supply network consists of 160 nodes and
190 arcs. For a division of 1 day in 5 intervals, we
find an optimal solution after 56 seconds. For a more
accurate division in 24 periods of one hour, a solu-
tion is found in 642 seconds. Although this second
approach is more accurate, the goal function changed
only marginally. The most important results for both
simulations are: First, water production centers that
have a high production cost (e.g. due to ground water
taxes) produce considerably less than low-cost pro-
duction centers. Second, the level of water in pure
water tanks remains as high as possible throughout
OPTIMIZING OPERATIONS OF LARGE WATER SUPPLY NETWORKS - A Case Study
471

the optimization period. A possible explanation for
this phenomenon is the fact that the connected pure
water pumps have to deliver less pressure, which re-
sults in a lower energy cost. Third, the level of the
water in the buffers moves in opposite direction to the
energy tariff during each period. This is a direct con-
sequence of the fact that the pumps will mainly work
during periods where the electricity tariff is low.
3.2 High-level Transport Optimization
Model
Setting the Tolerance Value. As stated before, the
value of the actual production flow in a groundwater
production center is allowed to vary slightly accord-
ing to a tolerance value, tol. For smaller values of tol,
the flow will be accurate but computation time will
rise. It is observed that setting this value to 0,0001
in (24) results in a reasonable computation time and a
realistic approximation of the production flow.
Water Hammer Considerations. Initial test runs
showed that the flow in a particular pipe suddenly
changed direction during one period and then changed
again in the subsequent period. Because we want to
minimize chances on water hammer, we add the ad-
ditional constraint that flow can only go in the same
direction for all periods. Ongoing research is con-
ducted to prevent unwanted direction changes in all
pipes during one or several periods.
Before the collaboration took place, network re-
sults are: (1) All pumps have activity points at values
which indicate high levels of efficiency (60-80%). (2)
In one of the pipe sections, water has to travel more
than 15 km before it reaches the other end. Simulation
shows however that this only happens after more than
24 hours. Since no demand takes place at this pipe
section, water will remain inside for quite some time,
which might lead to quality issues. This is certainly
something to take into consideration in future steps.
In the future collaboration, an additional 4000 m
3
of water will be transferred through this part of the
network. This results in: (1) additional costs. This is
of course expected, but it is worth remarking that the
extra cost is only due to increased electricity costs for
the pure water pumps while production cost slightly
drops. In one of the pump station the electricity cost
is increased by as much as 90 %. (2) The equilibrium
in the long pipe has shifted positively because of the
injection of extra water in the network. The quality
issue has thus decreased in size, which is a positive
side effect of the collaboration. Of course we will still
have to take this issue into account in other situations
or for other parts of the network.
4 CONCLUSIONS
In this paper a drinking water network in Flanders in-
volving constant speed pumps and buffers with free
inflow is modeled. Binary variables, used to model
buffer in and outflows, together with nonsmooth and
nonconvex character of some constraints and the
power term in the objective function complicate the
model. A linear approximation of the power term
was proposed, and the resulting model was solved us-
ing the open-source solver BONMIN. The model pro-
vided some initial good results on two real-life cases.
Some future improvements concern the issue of pre-
venting quality deterioration with changing directions
of water flow. The unwanted water hammer effects
need also to be tackled.
REFERENCES
Bonami, P. and Lee, J. (2006). Bonmin users manual. Tech-
nical report, Coin-Or.
Bragalli, C., D’Ambrosio, C., Lee, J., Lodi, A., and Toth,
P. (2006). An minlp solution method for a water net-
work problem. In LNCS, volume 4168, pages 696–
707. Springer-Verlag Berlin Heidelberg.
Burgschweiger, J., Gn
¨
adig, B., and Steinbach, M. C.
(2009). Nonlinear programming techniques for oper-
ative planning in large drinking water networks. Open
Appl Math J, 3:14–28.
Cembrano, G., Wells, G., Quevedo, J., Pehrez, R., and
Argelaguet, R. (2000). Optimal control of a water
distribution network in a supervisory control system.
Control Engineering Practice.
Defillet, D. (2011). Optimalisatie van het transport van
drinkwater op het niveau van het hoofdverdeelnet
van de vlaamse maatschappij voor watervoorziening.
Master thesis, University of Ghent.
Fourer, R., Gay, D., and Kernighan, B. (2003). AMPL: A
Modeling Language for Mathematical Programming.
Duxbury Press/Brooks/Cole Publishing Co.
L
´
opez-Ib
´
a
˜
nez, M., Prasad, T. D., and Paechter, B. (2008).
Ant colony optimization for optimal control of pumps
in water distribution networks. J Water Resour Plan
Manage, 134(4):337346.
Nicklow, J., Reed, P., Savic, D., Dessalegne, T., Harrell,
L., Chan-Hilton, A., Karamouz, M., Minsker, B., Os-
tfeld, A., Singh, A., and Zechman, E. (2010). State
of the art for genetic algorithms and beyond in water
resources planning and management. J Water Resour
Plan Manage, 136(4):412–432.
Savic, D. A. and Walters, G. A. (1997). Genetic algorithms
for least-cost design of water distribution networks. J
Water Resour Plan Manage, 123(2):67–77.
Verleye, D. and Aghezzaf, E.-H. (2011). Modeling and op-
timizing production and distribution of drinking wa-
ter at vmw. In LNCS, volume 6701, pages 315–326.
Springer-Verlag Berlin Heidelberg.
ICORES 2012 - 1st International Conference on Operations Research and Enterprise Systems
472