See discussions, stats, and author profiles for this publication at: [Link]

net/publication/282249133

Mathematical Programming Methods for Pressure Management in Water

Distribution Systems

Conference Paper  in  Procedia Engineering · September 2015

DOI: 10.1016/[Link].2015.08.974

CITATIONS READS

5 152

3 authors:

Filippo Pecci Edo Abraham

Imperial College London Delft University of Technology
14 PUBLICATIONS   59 CITATIONS    46 PUBLICATIONS   441 CITATIONS   

SEE PROFILE SEE PROFILE

Ivan Stoianov
Imperial College London
60 PUBLICATIONS   1,672 CITATIONS   

SEE PROFILE

Some of the authors of this publication are also working on these related projects:

Sustainable Water Resources Management Considering Climate Change Effects View project

Water Nexus Project View project

All content following this page was uploaded by Edo Abraham on 28 September 2015.

The user has requested enhancement of the downloaded file.

 Available online at [Link]

ScienceDirect
Procedia Engineering 119 (2015) 937 – 946

13th Computer Control for Water Industry Conference, CCWI 2015

Mathematical programming methods for pressure management in

water distribution systems
Filippo Peccia,∗, Edo Abrahama , Ivan Stoianova
a InfraSense Labs, Dept. of Civil and Environmental Engineering, Imperial College London, South Kensington Campus, London SW7 2AZ , UK.

Abstract
In this paper, we survey mathematical programming methods for the management of pressure in water distribution systems through
optimal placement and operation of control valves. The optimization framework addresses the minimization of average zone
pressure under multiple demand scenarios, enforcing hydraulic equations as nonlinear constraints. Binary variables are used to
model the placement of valves. The derived nonlinear optimization problem is a non-convex MINLP. We implement and evaluate
a direct MINLP solver and two different reformulation methods for MINLPs that solve sequences of regular NLPs. Moreover, we
investigate the solutions under different design and operation loading conditions.
©c 2015
2015Published
The Authors. Published
by Elsevier byisElsevier
Ltd. This Ltd. article under the CC BY-NC-ND license
an open access
([Link]
Peer-review under responsibility of the Scientific Committee of CCWI 2015.
Peer-review under responsibility of the Scientific Committee of CCWI 2015
Keywords: Water distribution systems, valve placement, mathematical optimization, mixed integer nonlinear programming, pressure
management.

1. Introduction

Pressure is a critical variable for water distribution systems since it is directly related to leakage and burst fre-
quency [1]. Great benefits can be derived by keeping the operational pressure close to the minimum allowed by
standard regulations. The subdivision of water distribution networks into small areas, called District Metered Areas
(DMAs), has been successfully implemented in the last two decades for the management of pressure and leakage.
The sectorization is realized through boundary valves, that are closed in order to form small metered areas, so as to
easily monitor flow between zones and optimize pressure in a simplified way. This practice has improved leakage
reduction, but at the same time has severely reduced the redundancy of the networks, affecting resilience and water
quality [2]. Previous work on water distribution networks with adaptively reconfigurable topology [2] integrated the
achievements in terms of leakage management provided by the DMAs structure with the benefits of high network
connectivity, resulting in increased network redundancy and resilience. The network is segregated into small areas
during the night, to take advantage of the sectorized topology for leakage localization. DMAs are then aggregated in
order to achieve an effective pressure and resilience management in diurnal operation. To benefit from these emerging

∗ Corresponding author.
E-mail address: f.pecci14@[Link]

1877-7058 © 2015 Published by Elsevier Ltd. This is an open access article under the CC BY-NC-ND license
([Link]
Peer-review under responsibility of the Scientific Committee of CCWI 2015
doi:10.1016/[Link].2015.08.974
938 Filippo Pecci et al. / Procedia Engineering 119 (2015) 937 – 946

and advanced control schemes, the retrofit of existing networks requires the formulation and solution of both design
and operational optimization problems.
In the present work we focus on mathematical optimization for network pressure management, minimizing average
zone pressure through the optimal placement and operation of pressure reducing valves. Valve placement is a chal-
lenging problem because it combines binary variables (whether the valve is placed or not) with continuous variables
representing nodal pressures and pipe flows. Valve control is embedded in the optimal placement problem, since the
operational setting of a valve is identified with the pressure at the downstream node. The use of physically feasible
models for real water distribution networks involves the formulation of hydraulic equations that lead to nonlinear
constraints. Together with discrete decision variables, these result in a non-convex optimization problem with a large
number of integer variables, belonging to the class of mixed integer nonlinear programs (MINLP).
Optimal valve placement has been studied using both mathematical programming and heuristic methods. In [3] a
mixed integer linear programming method was applied to a linear approximation of the original optimization problem.
A direct solution with a mixed integer nonlinear solver has been proposed in the report [4], considering just a steady-
state snapshot of daily network operation. A faithful representation of real world water distribution systems requires
the formulation of an extended time optimization problem, which includes multiple demand scenarios. These have
been considered in [5], wherein is proposed a reformulation approach for mixed integer nonlinear programs. By
relaxing the integer constraints, a penalty method is coupled with a heuristic rounding scheme for binary variables in
order to improve the convergence of the method. Genetic algorithms and meta-heuristic approaches have also been
applied to optimal valve placement in [6–8]. However, these methods present some fundamental drawbacks: they do
not guarantee convergence to optimal solutions (not even to local optima) and they require a large number of function
evaluations of objective and constraints in order to achieve good quality solutions. Therefore, with application to
advanced and dynamic control schemes, it is necessary to study reliable and effective mathematical optimization
methods for optimal valve placement and operation.
MINLPs are particularly challenging and several techniques have been applied to these problems, see for refer-
ence [9–11]. Since the integer variables involved are all binary, in the present work we investigate possible solutions
by reformulating the problem as an equivalent mathematical program with complementarity constraints (MPCC). Th
feasible set of an MPCC present a special structure that prevents solution via direct application of standard nonlinear
programming solvers. Various methods have been proposed to handle this pathological nature; we focus on penaliza-
tion and relaxation approaches. Convergence properties under suitable assumptions have been proved for both penalty
methods [12–15] and relaxation methods [16–18]. We propose an algorithmic implementation of these reformulation
approaches and we apply them to the case study.
Since the problem is non-convex, like most nonlinear programming algorithms, under suitable assumptions all the
methods guarantee convergence only to local minimum points and the quality of the solution would depend on the
initial point. We have to take this into account when comparing different approaches. Therefore, we will consider a
solution qualitatively good if it is obtained with various random initial guesses and provides an average zone pressure
close to the best known solution.

2. Problem formulation

We address the minimization of average zone pressure through the placement of nv control valves in a water
distribution network with n0 water sources (eg. reservoirs or tanks), nn nodes and n p pipes, represented as a directed
graph (V, E), with nn + n0 vertices and 2n p links. Note that the jth physical pipe is modelled by two graph links
j and j + n p , for all j = 1, ..., n p . The optimization problem is solved for extended time settings; we discretize
the considered time interval into nl subintervals. For each time step k ∈ {1, ..nl }, we have vectors of known nodal
demands dk ∈ Rnn and known fixed heads at water sources hk0 ∈ Rn0 . Moreover, each node i ∈ {1, ..., nn } has a
known elevation ei . Unknown pressure heads and flows, at each time step k, are represented by the vectors pk =
j
[pk1 , ..., pknn ]T and qk = [qk1 , ..., qk2n p ]T , respectively. The friction headloss across every link i1 →
− i2 is commonly
represented by the Hazen-Williams (HW) or Darcy-Weisbach (DW) formula. In DW models the relation between
friction headloss and flow is expressed by an implicit equation and has to be calculated through an iterative process.
This complicates the implementation of these models in a mathematical optimization framework. On the other hand,
the semi-empirical HW formula is given by HW f (qkj ) = r j (qkj )1.852 , where r j , the resistance coefficient of the pipe, is
 Filippo Pecci et al. / Procedia Engineering 119 (2015) 937 – 946 939

10.670L j
defined by r j = C 1.852 D4.871
with L j , C j , D j length, roughness and diameter of the pipe, respectively. The HW formula
j j
involves a non-smooth exponential function, hence is difficult to handle for most nonlinear programming solvers.
Therefore, a smooth and sufficiently accurate approximation of the headloss formulae is needed. Consider a generic
pipe with length L, diameter D and H-W coefficient C. We model the friction head loss using a quadratic function
q Ä ä2
h f (q) := a q2 +b q, with (a , b ) the minima of the integral of the square error Jq1 ,q2 (a, b) := q12 aq2 +bq−rq1.852 dq
over an approximation range [q1 , q2 ].
Once the quadratic model is calculated, it is possible to formulate the optimization problem. For each link j we
need to introduce a new unknown, which assume value 0 or 1 depending whether a valve is present on j or not
®
1, if a valve is placed on link j
vj =
0, otherwise.

The optimization problem for optimal valve placement and operation can be formulated as follows, for the sake of
brevity we omit some details on the definition of constraints, which can be found in [3–5,19].

nl
1 
nn
min wi pki (1a)
k=1
W i=1
subject to: AT12 qk − dk = 0, ∀k = 1, ..., nl , (1b)
 
S (qk ) −A12 pk − A12 e − A10 hk0 − h f (qk ) ≥ 0, ∀k = 1, ..., nl , (1c)
− A12 p − A12 e −
k
− h f (q ) − M v ≤ 0,
A10 hk0 k k
∀k = 1, ..., nl , (1d)
v j + vn p + j ≤ 1, ∀ j = 1, ..., n p , (1e)

2n p
v j = nv , (1f)
j=1

pkmin ≤ pk ≤ pkmax , ∀k = 1, ..., nl , (1g)

πD2j
0 ≤ qkj ≤ , ∀ j = 1, ..., 2n p , ∀k = 1, ..., nl , (1h)
4
v ∈ {0, 1} . 2n p
(1i)

The objective is the minimization of average zone pressure over anextended time period, expressed by the weighted
L
sum of nodal pressures (1a), whereweights wi are defined by wi := j∈I(i) 2j , with I(i) the set of indices of links that
nn
are connected to node i and W := i=1 wi . The optimization problem is primarily subject to hydraulic constraints, in
order to ensure the physical feasibility of the solutions. The conservation of flow at each junction node is expressed
by (1b), while the couple (1c), (1d) represents the conservation of energy across all links. The matrices AT12 ∈ Rnn ×2n p
and AT10 ∈ Rn0 ×2n p are the node-edge incidence matrices for the nn nodes and the n0 water sources, respectively.
Moreover, S ∈ R2n p ×2n p is the diagonal matrix of flows qk , S (qk ) j, j = qkj , and h f (qk ) = [h f (qk1 ) ... h f (qk2n p )]T . Finally,
M k ∈ R2n p ×2n p is a diagonal matrix of sufficiently large constants. The optimization problem is also subjected to some
specific physical, economic and operational constraints, expressed by (1e), (1f), (1g) and (1h), respectively.
Problem (1) is a non-convex MINLP, with polynomial nonlinear constraints and N = nl (nn + 2n p ) + 2n p variables.
This class of mathematical programs is particularly challenging, especially when the scale of the problem is large. In
fact, it combines the difficulties of dealing with nonlinear non-convex constraints, with the complication of optimizing
in a discrete variable space. When the number of integer variables becomes large, it is in general impossible to find a
solution. We aim to investigate specialized mathematical programming solution methods for (1).

3. Reformulation approaches
 
Let x = [pT , qT , vT ]T ∈ RN be the vector of unknowns. We set B = nl (nn + 2n p ) + 1, ..., nl (nn + 2n p ) + 2n p , the set
of indices corresponding to the binary variables. Similarly, let I and E be the index set of rows in (1) which correspond
940 Filippo Pecci et al. / Procedia Engineering 119 (2015) 937 – 946

to inequality and equality constraints, respectively. Moreover, let g : RN → Rm with m = |I| and h : RN → Rn with
n = |E| be the functions that group inequality and equality constraints, respectively. Since its integer variables are
binaries, Problem (1) can be equivalently reformulated as a mathematical program with complementarity constraints
(MPCC):

min f (x) (2a)

x∈RN
subject to gi (x) ≥ 0, ∀i ∈ I, (2b)
hi (x) = 0, ∀i ∈ E, (2c)
0 ≤ x j ⊥ 1 − x j ≥ 0, ∀ j ∈ B. (2d)

In the framework of nonlinear programming (NLP), linear independence of the gradients of constraints, also referred
to as linear independence constraint qualification (LICQ), is usually required in order to guarantee the convergence of
NLP solvers to stationary points, see [20]. Complementarity constraints result in a special feasible set, which violates
generic constraint qualifications, causing severe convergence issues to standard NLP solvers. Nonetheless, MPCCs
have been solved using tailored algorithms in other engineering frameworks [15,21]. Various methods have been
proposed in recent years, see [12–14,16–18,22] and the references therein.
We are interested in penalty and relaxation approaches, both employ the solution of sequences of regular nonlinear
programs. The penalty method relies on the solution of a sequence of NLPs, each one obtained dropping the special
constraint (2d) and penalizing in the objective deviations from the complementarity condition. It has already been
applied to optimal valve placement in [5], coupled with an heuristic rounding scheme, in order to improve the conver-
gence of their method. On the other hand, the relaxation method focuses on the geometry of the region of feasibility
and solves a sequence of NLPs with a simplified feasible set. Although this relaxation approach has been adopted in
other engineering problems with complementarity constraints, its application to water distribution systems is novel.
In the following we present them with more details. 
Using the same notation of (2), we introduce the following penalty function Ψ(x) = j∈B x j (1 − x j ). Given ρ > 0,
consider the nonlinear program PEN(ρ):

min f (x) + ρΨ(x) s.t.

x∈RN
gi (x) ≥ 0, ∀i ∈ I, (3)
hi (x) = 0, ∀i ∈ E,
0 ≤ x j ≤ 1, ∀ j ∈ B.

Penalty approaches have been studied in [12–15,21], the overall strategy consists in the generation of a sequence of
k→+∞
stationary points (xk )k of PEN(ρk ), with ρk −−−−−→ +∞. Then it is possible to prove (see [12,15]) that the iterates (xk )
converge to a stationary point of Problem (2). We implement the penalty method is Algorithm  guess g0
1. The initial
is selected randomly and the complementarity violation is evaluated as Vio(x) = max j∈B min (x j , 1 − x j ) .

Algorithm 1 Penalty method

1: Initialization:
Select random g0 , β ← 1.1, α ← 0.01,  ← 10−6 ;
2: Solve PEN(0) with initial guess g0 and get the solution x0 and value of objective function z0 ;
3: k ← 0; ρ0 ← α · z0 ;
4: while Vio(xk ) >  do
5: Solve PEN(ρk ) with initial guess xk and get the solution xk+1 ;
6: ρk+1 ← β · ρk ;
7: k ← k + 1;
8: end while

Alternative approaches for the solution of Problem (2) include relaxation methods, see [14,17,18]. In the present
work we focus on the technique proposed by [17], where theoretical proof for the convergence is provided. For t > 0
 Filippo Pecci et al. / Procedia Engineering 119 (2015) 937 – 946 941

consider the nonlinear program REL(t)

min f (x) s.t.

x∈RN
gi (x) ≥ 0, ∀i ∈ I,
hi (x) = 0, ∀i ∈ E, (4)
0 ≤ x j ≤ 1, ∀ j ∈ B,

x j (1 − x j ) ≤ t.
j∈B

The method generates a sequence of stationary points of (4), for decreasing values of the positive parameter t, solving
the nonlinear problem REL(t). The sequence will convergence to a stationary point of Problem (2), for a rigorous
and general proof see the appendix in [17]. We implement the relaxation method in Algorithm 2, where Vio(x) is the
same as above.

Algorithm 2 Relaxation method

1: Initialization:
Select
 random g0 , x0 ← g0 , t0 ← 1, tmin ← 10−15 ,  ← 10−6 , β ← 10−4 , k ← 0;
  
2: while Vio(x ) >  & tk > tmin OR k = 0 do
k

3: Solve REL(tk ) with initial guess xk and get the solution xk+1 ;
4: tk+1 ← β · tk ;
5: k ← k + 1;
6: end while

All the NLP subproblems involved in the presented methods have sparse nonlinear structures and so can be solved
using tailored techniques for large sparse nonlinear programs, offering a scalable approach for large scale water
distribution systems. However, note that both penalty and relaxation approaches are local methods and so converge to
stationary points of Problem (2). If the nonlinear constraints are non-convex, which is the case in our problem, then
the convergence is highly influenced by initial guess and algorithmic parameters.

4. Case study

We have applied the different outlined mathematical programming approaches to Problem (1) and we compare their
solutions with those provided by the MINLP solver BONMIN (v.1.8.1). The NLP subproblems within penalty and
relaxation methods were solved using the interior point solver for large scale nonlinear optimization IPOPT (v3.11.8).
Moreover, we directly provided to BONMIN and IPOPT sparse gradients and Jacobians, in order to exploit the very
sparse structure of the problem. All computations were performed within MATLAB 2013b-64 bit for Windows 7,
installed on a 2.50GHz Intel Xeon(R) CPU E5-2640 0 with 18 Cores.
The selected benchmarking network has been already considered in [4–8,23,24], see Figure 1a. The network has
22 nodes, 37 pipes and 3 reservoirs, see [5,24] for details on pipes’ characteristics, nodal demands and reservoir levels.
In contrast to the report in [4], our formulation is extended time, we divide network’s daily operation into 24 time
steps, assigning a specific demand scenario to each hour. A similar model was considered in [5]. Daily demand profile
is reported in Figure 1b.
In the code implementation, in order toaccommodate the scale of the problem, we considered the minimization
nl  n n k
of the non-normalized version of (1a) : k=1 i=1 wi pi . We follow the recommendation in [9,25] and we apply
BONMIN with branch and bound algorithm to our non-convex MINLP (1), all other input options for BONMIN
were set to their default values. In [5] BONMIN was reported to fail when the placement of more than 2 valves was
considered. Also in our initial implementation the solver encountered
n p convergence issues which we have overcome
with the substitution of the equality linear constraint (1f) with j=1 v j ≤ nv . Note that this inequality constraint results
in an equivalent optimization problem, since we cannot reduce pressure further by having fewer valves. Moreover,
the direct implementation of sparse gradients and Jacobians has improved the performance of BONMIN, allowing
942 Filippo Pecci et al. / Procedia Engineering 119 (2015) 937 – 946

1.4

1.2

Relative demand
1

0.8

0.6

0.4

0.2

0
00.00 06.00 12.00 18.00 24.00

(a) (b)

Figure 1. (a) network model from [23]; (b) demand profile

convergence to the best known optimal solutions for this network. In Table 1 we summarize the results, while in
Figure 2 we present the average zone pressure profiles, corresponding to different number of optimized valves. Note
that a relatively small decrease in pressure (p) can result in significant reduction in the amount of leakage (L), since
L ∝ pαL with 0.5 ≤ αL ≤ 2.5 [1]. Moreover, the placement of an optimized control valve can provide other benefits,
not directly related to the minimization of AZP. For example, when we move from the configuration with 4 valves to
the optimized set of 5 valves we achieve a substantial decrease in night pressure, reducing pressure variability in the
network (see Figure 2).
In Figure 3 we show the optimal daily operation for 4 optimized valves. Valve on link (23, 1) is always open and
acts as a pressure reducing valve (PRV), maintaining outlet pressure close to the minimum. On the other hand, the
control valve on (22, 15) operates as a fixed closed boundary valve (BV). The remaining two valves have a similar
dynamic behaviour; they are closed BVs during the night (minimum demand period) and active PRVs in diurnal
operation. This benchmarking example demonstrates that our optimization constraints are very flexible and allow the
modelling of control valves that operate as PRVs or BVs depending on hydraulic conditions. Therefore, the presented
mathematical formulation is suitable for use in the design and implementation of the adaptive reconfigurable topology
for water distribution networks proposed in [2].
Generally, the non-convexity of nonlinear constraint affects the convergence properties of tools like BONMIN for
the solution of MINLPs, as noted in [25]. However, the optimal placements in Table 1 are obtained with different
random initial guesses and are the best-known solutions for the case study. Therefore, the robustness of our imple-
mentation of BONMIN is satisfactory. Nonetheless, computational time suggests that more issues can arise when
dealing with large scale networks, where the number of integer variables grows. This is the main motivation for the
study of possible reformulation approaches.

Table 1. Optimal placement of control valves with Bonmin

No. of valves Link AZP/h CPU time
1 (13, 12) 33.90 m 76 s
2 (23, 1), (13, 12) 32.94 m 962 s
3 (23, 1), (13, 12), (16, 17) 32.61 m 1345 s
4 (23, 1), (22, 15), (12, 15), (13, 12) 31.90 m 2124 s
5 (23, 1), (22, 15), (12, 15), (25, 16), (13, 12) 31.48 m 1892 s
 Filippo Pecci et al. / Procedia Engineering 119 (2015) 937 – 946 943

50 No valves
1 valve
2 valves
45
3 valves

AZP (m)
4 valves
40 5 valves

35

30
00.00 06.00 12.00 18.00 24.00

Figure 2. Average zone pressure after the optimal placement of a different number of control valves

56 20 48 20
Hydraulic head (m)

Hydraulic head (m)

46
54
15 15
44
Flow (l/s)

Flow (l/s)
Inlet
52
Inlet Outlet
10 42 10
50 Outlet Flow
Flow 40
5 5
48
38

46 0 36 0
00.00 06.00 12.00 18.00 24.00 00.00 06.00 12.00 18.00 24.00

Valve on link (23, 1) Valve on link (22, 15), always closed

56 20 56
Inlet 100
Hydraulic head (m)

Hydraulic head (m)

52 Outlet 54 Inlet
Flow 15 Outlet 80
52
Flow (l/s)

Flow (l/s)

48 Flow
60
10 50
44
48 40
5
40 20
46

36 0 44 0
00.00 06.00 12.00 18.00 24.00 00.00 06.00 12.00 18.00 24.00

Valve on link (12, 15) Valve on link (13, 12)

Figure 3. Operation of 4 optimized valves

944 Filippo Pecci et al. / Procedia Engineering 119 (2015) 937 – 946

4.1. Reformulation approaches

We applied to the case study the penalty and relaxation methods outlined in Section 3. We minimize average zone
pressure, with 24 time steps, the results are reported in Table 2. A  highlights when the algorithm converged to the
BONMIN solution. The penalty approach results in good quality solutions, although these do not always coincide
with the BONMIN optima; however, our implementation of the penalty method results in improved performance with
respect to a similar approach used in [5], both in terms of quality of solutions and computational time. On the other
hand, the relaxation method was observed to achieve the same solution as BONMIN in most cases, with only slightly
sub-optimal solutions where they are different.

Table 2. Optimal placement of control valves with reformulation approaches

No. of valves Link IDs AZP/h CPU time
Solution with Penalty method
1 (13, 12) 33.90 m 126 s 
2 (24, 10), (13, 12) 33.54 m 178 s
3 (23, 1), (24, 10), (13, 12) 32.65 m 225 s
4 (23, 1), (12, 15), (13, 12), (13, 12) 31.90 m 235 s 
5 (23, 1), (22, 15), (12, 15), (24, 10), (13, 12) 31.62 m 101 s

Solution with Relaxation method

1 (13, 12) 33.90 m 66 s 
2 (23, 1), (13, 12) 32.94 m 29 s 
3 (23, 1), (12, 15), (13, 12) 32.93 m 77 s
4 (23, 1), (22, 15), (12, 15), (13, 12) 31.90 m 15 s 
5 (23, 1), (22, 15), (12, 15), (24, 10), (13, 12) 31.62 m 23 s

Enabled by its superior computational performances, we applied the relaxation method to the problem in (1), with
the number of valves ranging from 0 to 37. We report the average AZP for each number of optimized valves in Figure
4. We can consider Figure 4 as an approximation of the pareto front for a multiobjective optimization problem with
objectives the minimization of AZP and the installation cost, expressed here by the number of valves. Note that, as
reported also in [6], the region corresponding to large numbers of valves exhibits high instability and some of those
points can not be optima; in fact, it is impossible to increase optimal AZP by adding more control valves. However,
using Figure 4, we can realize a cost-benefit analysis: it appears that any solution with more than 7 valves should not
be considered for this network, since the corresponding reductions in AZP are negligible.

4.2. Different loading conditions

Finally, we aim to investigate different demand conditions for the placement of valves and then for the determi-
nation of optimal control settings. Consider the static problem with just one fixed demand condition for the optimal
placement of 3 valves. We solve with BONMIN three static optimization problems with demand fixed to the min-
imum, medium and maximum demand conditions from the diurnal demand profile. Then, enforcing the resulting
optimal valve placement, we use IPOPT to solve the resulting nonlinear programs for the optimal operation of each
valve. In Figure 5 compares the average zone pressure corresponding to the different configurations. Since the optimal
placement found in Table 1 corresponds to the solution with medium demand, we conclude that an optimal solution for
difficult extended time problem (1) can be obtained through the solution of a smaller MINLP using only the medium
demand snapshot. Subsequently, we solve for an extended time NLP with the valve locations fixed.

5. Conclusions

We have presented a rigorous mathematical framework for the minimization of average zone pressure in water
distribution networks, through optimal placement and operation of control valves. The mathematical formulation is
 Filippo Pecci et al. / Procedia Engineering 119 (2015) 937 – 946 945

35
Minimum demand
Medium demand
38 Peak demand
34
Average AZP (m)

AZP (m)
36

34 33

32
32
00.00 06.00 12.00 18.00 24.00
30
0 7 14 21 28 37
Number of valves Figure 5. AZP corresponding to optimal placements ob-
tained with different loading conditions: low demand:
Figure 4. Average daily AZP corresponding to a different number (23, 1), (24, 10), (13, 12); medium demand: (23, 1), (13, 12),
of optimized valves (16, 17); peak demand: (23, 1), (22, 15), (13, 12)

sufficiently flexible to incorporate different physical actuators, such as PRVs and BVs, within the same model. The
optimization problem belongs to the class of non-convex mixed integer nonlinear programs. We have implemented
and evaluated a direct MINLP solver BONMIN, using a benchmarking network model. Although our BONMIN
implementation performed better than previous applications to the same problem, the CPU time required by the
algorithm in this relative small network model suggests that infeasible computational time would arise when dealing
with large scale networks, which would involve a large number of integer variables. We have therefore presented and
implemented two reformulation approaches that reduce the MINLP to an NLP, the penalty and relaxation methods.
The two methods rely on the solution of a sequence of sparse nonlinear problems which can be solved using tailored
techniques for large sparse nonlinear programs. They have been both applied to the case study. The penalty method
results in improved performances with respect previous work [5], where a similar scheme was applied to the problem
of optimal valve placement for the same benchmarking network. However, the relaxation approach is shown to
have superior performance both in quality of the solutions and CPU time. This relaxation method has been adopted
in other engineering frameworks when dealing with complementarity constraints, however, its application to water
distribution networks is novel. The numerical results are promising and suggests that our relaxation approach can be
successfully applied to large-scale non-convex optimization problems with binary variables related to optimal design
and operation of water distribution systems. Finally, we have investigated the influence of different loading conditions
for the placement of valves. Results reveal that optimal pressure management is achieved when the placement of the
valves is optimized under medium demand condition. In conclusion, the presented benchmarking study demonstrates
that our mathematical optimization framework can provide effective tools to support design and operation of smart
water distribution systems.

Acknowledgements

This work was supported by the NEC-Imperial Smart Water Systems project. Project Reference: NEC-Imperial
Smart Water Systems.

References

[1] A. Lambert, What do we know about pressure: leakage relationships in distribution systems?, in: Iwa Conference ’System Approach To
Leakage Control and Water Distribution Systems Management’, 2001.
[2] R. Wright, I. Stoianov, P. Parpas, K. Henderson, J. King, Adaptive water distribution networks with dynamically reconfigurable topology,
Journal of Hydroinformatics 16 (2014) 1280.
946 Filippo Pecci et al. / Procedia Engineering 119 (2015) 937 – 946

[3] K. S. Hindi, Y. M. Hamam, Locating Pressure Control Elements for Leakage Minimization in Water Supply Networks: an Optimization Model,
Engineering Optimization 17 (1991) 281–291.
[4] B. J. Eck, M. Mevissen, Non-Linear Optimization with Quadratic Pipe Friction, Technical Report, 2012.
[5] P. D. Dai, P. Li, Optimal Localization of Pressure Reducing Valves in Water Distribution Systems by a Reformulation Approach, Water
Resources Management 28 (2014) 3057–3074.
[6] M. Nicolini, L. Zovatto, Optimal location and control of pressure reducing valves in water networks, Journal of Water Resources Planning and
Management 135 (2009) 178–187.
[7] S. Liberatore, G. M. Sechi, Location and calibration of valves in water distribution networks using a scatter-search meta-heuristic approach,
Water Resources Management 23 (2009) 1479–1495.
[8] M. E. Ali, Knowledge-Based Optimization Model for Control Valve Locations in Water Distribution Networks, Journal of Water Resources
Planning and Management 141 (2015) 1–7.
[9] C. Bragalli, C. D’Ambrosio, J. Lee, A. Lodi, P. Toth, On the optimal design of water distribution networks: a practical MINLP approach,
Optimization and Engineering 13 (2012) 219–246.
[10] P. Belotti, C. Kirches, S. Leyffer, J. Linderoth, J. Luedtke, P. Belotti, Mixed-integer nonlinear optimization, Acta Numerica 22 (2013) 1–131.
[11] J. Lee, S. Leyffer, Mixed Integer Nonlinear Programming, Springer New York, 2012.
[12] X. M. Hu, D. Ralph, Convergence of a Penalty Method for Mathematical Programming with Complementarity Constraints, Journal of
Optimization Theory and Applications 123 (2004) 365–390.
[13] S. Leyffer, G. López-Calva, J. Nocedal, Interior methods for mathematical programs with complementarity constraints, SIAM Journal on
Optimization 17 (2006) 52–77.
[14] D. Ralph, S. J. Wright, Some properties of regularization and penalization schemes for MPECs, Optimization Methods and Software 19 (2004)
527–556.
[15] M. Herty, S. Steffensen, MPCC Solution Approaches for a Class of MINLPs with Applications in Chemical Engineering, Technical Report,
Aachen Institute for Advanced Study in Computational Engineering Science, 2012.
[16] A. U. Raghunathan, L. T. Biegler, An Interior Point Method for Mathematical Programs with Complementarity Constraints (MPCCs), SIAM
Journal on Optimization 15 (2005) 720–750.
[17] S. Scholtes, Convergence properties of a regularization scheme for mathematical programs with complementarity constraints, SIAM Journal
on Optimization 11 (2001) 918–936.
[18] T. Hoheisel, C. Kanzow, A. Schwartz, Theoretical and numerical comparison of relaxation methods for mathematical programs with comple-
mentarity constraints, Mathematical Programming 137 (2013) 257–288.
[19] H. Sherali, E. Smith, A global optimization approach to a water distribution network design problem, Journal of Global Optimization 11 (1997)
107–132.
[20] J. Nocedal, S. J. Wright, Numerical Optimization, Springer Series in Operations Research and Financial Engineering, 2nd ed., Springer-Verlag,
New York, 2006.
[21] A. U. Raghunathan, L. T. Biegler, Mathematical programs with equilibrium constraints (MPECs) in process engineering, Computers &
Chemical Engineering 27 (2003) 1381–1392.
[22] H. Scheel, S. Scholtes, Mathematical Programs with Complementarity Constraints: Stationarity, Optimality, and Sensitivity, Mathematics of
Operations Research 25 (2000) 1–22.
[23] M. Sterling, A. Bargiela, Leakage reduction by optimised control of valves in water networks, Transactions of the Institute of Measurement
and Control 6 (1984) 293–298.
[24] P. Jowitt, C. Xu, Optimal valve control in water-distribution networks, Journal of Water Resources Planning and Management 116 (1990)
455–472.
[25] C. D’Ambrosio, A. Lodi, Mixed integer nonlinear programming tools: An updated practical overview, Annals of Operations Research 204
(2013) 301–320.

View publication stats