A Direct Approach for the Near-Optimal Design of Water Distribution Networks Based on Power Use.

The optimal design of Water Distribution Systems (WDS) has been addressed internationally during the last decades

/var/www/pavco.com.co/public/site/pdftohtml/5e81077962f627b67949947c14bfa4ef/index-html.html
background image

water

Article

A Direct Approach for the Near-Optimal Design of
Water Distribution Networks Based on Power Use

Juan Saldarriaga

1,

* , Diego Páez

2

, Camilo Salcedo

1

Paula Cuero

3

, Laura Lunita López

3

,

Natalia León

3

and David Celeita

3

1

Department of Civil and Environmental Engineering, Universidad de los Andes, Bogotá 111711, Colombia;
ca.salcedo959@uniandes.edu.co

2

Department of Civil Engineering, Queen’s University, Kingston, ON K7L 3N6, Canada;
da.paez270@gmail.com

3

Water Supply and Sewer Systems Research Center (CIACUA), Universidad de los Andes, Bogotá 111711,
Colombia; pa.cuero38@uniandes.edu.co (P.C.); ll.lopez28@uniandes.edu.co (L.L.L.);
n.leon40@uniandes.edu.co (N.L.); da.celeita10@uniandes.edu.co (D.C.)

*

Correspondence: jsaldarr@uniandes.edu.co; Tel.:

+57-1-339-49-49 (ext. 2805)

Received: 7 February 2020; Accepted: 24 March 2020; Published: 6 April 2020

Abstract:

In recent years, iterative computational techniques have been considered as the most

e

ffective methods to tackle the problem of Water Distribution System (WDS) minimum-cost design.

Given their stochastic nature, these approaches involve a large number of hydraulic simulations in
order to obtain suitable results. Herein, a WDS design methodology based entirely on hydraulic

principles is presented. This methodology, named Optimal Power Use Surface (OPUS), focuses on

both reaching low-cost designs and diminishing the number of hydraulic executions (iterations),
by establishing e

fficient ways in which energy is dissipated and flow is distributed throughout the

system. The algorithm was tested in four well known benchmark networks, previously reported
in the literature. OPUS proved that following hydraulic principles is a fair choice to design WDS,
showing plenty of potential in other water distribution mathematical modeling applications and
o

ffering an alternative for the extensive search process undertaken by metaheuristics.

Keywords:

minimum cost design; energy-based methods; optimization; WDS; hydraulic behavior

1. Introduction

The optimal design of Water Distribution Systems (WDS) has been addressed internationally

during the last decades. This is as a result of the need to make an e

fficient distribution of scarce economic

resources without compromising the welfare of society. This problem is substantially magnified in
the context of developing countries because the availability of resources is very limited, and there is
a significant deficiency of WDSs that are capable of supplying water in appropriate conditions. To
overcome these di

fficulties, it is important to use methodologies that are both efficient and effective in

addressing the minimum-cost WDS design problem.

For instance, in developing countries it is common to find some cities that have WDS that cannot

supply water continuously for 24 h a day. Thus, water utilities must find alternative ways to solve the

water availability problem with intermittent water supply, which is not well-solved internationally [

1

].

Another problem found in developing countries is related to aging water networks, but few economic

resources are available to perform proper maintenance and rehabilitation programs. In addition, there
is increasing pressure on water supply because of climate change and population growth [

2

]: cities

are requiring more water, but water is becoming scarcer. Hence, water utilities must find strategies
to solve the problem. For instance, it is well-known that pressure control strategies over all of the

Water 2020, 12, 1037; doi:10.3390

/w12041037

www.mdpi.com

/journal/water

/var/www/pavco.com.co/public/site/pdftohtml/5e81077962f627b67949947c14bfa4ef/index-html.html
background image

Water 2020, 12, 1037

2 of 26

network allows for reducing energy losses in water networks, as a workaround for making better use
of energy resources. This reduces both energy consumption at pumping stations and water losses in

pipes [

2

], which will allow water utilities to extend the lifespan of aging assets through the reduction

of head-losses. Furthermore, the optimal selection of diameters for the equipment, pipes and pumps in

WDSs have been found to achieve considerable energy reductions [

3

]. This increases the importance of

making good decisions regarding energy losses in the network, as the objective of water distribution is
to deliver water in the quantity required at adequate pressures with the maximum coverage possible
at a

ffordable costs [

3

].

For all these reasons, we find that it is important for water utilities to have several alternatives

for the design of WDSs. As there is very little research on intermittent water supply, there is also a
general lack in design methods that are straightforward and simple enough to be used by water utilities
that lack the interest, the technology or the understanding required to apply state-of-the art design
methods based on evolutionary algorithms. Thus, we propose this methodology as a workaround to
solve design problems using an energy method, which will allow the designer to fully understand
and control what is happening with the energy losses in the network, which we believe may be very
useful to network designers in cities with few computation resources available, or cities where energy
reduction is the main goal.

It has been shown that finding the optimal configuration of diameters is a NP-Hard problem

(non-deterministic polynomial-time) [

4

], even for the simplest of branched networks, it is a highly

constrained nonlinear optimization problem because of the need to satisfy node pressure and pipe
flow restrictions. This means that finding an optimal solution using exhaustive and deterministic
methods is not viable, hence, approximate methods are required to reach near optimal solutions. Most
approximate methodologies are stochastic in nature, and currently most researchers use evolutionary
algorithms. However, those approximate methodologies will try to search the whole solution space to
find optimal and near-optimal solutions, making it necessary to reduce the search space and guide the
heuristic methods to find the solution. For instance, Wu (2005) proposed self-adaptive penalties with
Genetic Algorithms [

5

] and Montalvo (2010) proposed self-adaptive penalties with Particle Swarm

Optimization [

6

].

So far, several authors have proposed approaches to solve the aforementioned problem, which have

led to less expensive designs. In the 1970s and 1980s, engineers used traditional optimization methods
such as enumeration [

7

], linear and nonlinear programming [

8

]. However, in the mid 1990s and

beginnings of the 20th century, these approaches were slowly replaced by metaheuristic algorithms
due to the greater ease of their implementation, more extensive exploration of the solution space,
less dependence on the initial condition of the system, and the ability to integrate the discrete-sized
diameter restrictions. Among the successful attempts in the literature, there are: Di

fferential Evolution

Algorithms [

9

], Genetic Algorithms [

10

], Simulated Annealing ([

11

,

12

]), Harmony Search [

13

], Scatter

Search [

14

], Cross Entropy [

15

], and Particle Swarm [

16

]. Moreover, other authors have proposed

combined approaches such as binary linear programming with di

fferential evolution algorithm [

17

], as

well as proposed approaches to increase the e

fficiency of the algorithm by reducing the search space,

using afterwards a Genetic Algorithm [

18

].

Despite the advantages of these stochastic algorithms, they have an intrinsic limitation that exist

in the requirement of an appropriate determination of several parameters inherent to each algorithm,

which is a time-consuming and laborious task [

19

]. To accomplish this, metaheuristics use a large

number of randomly-generated alternatives and use generic learning functions that progressively

produce better results. However, these algorithms require the assessment of thousands of di

fferent

design alternatives before they converge on a solution. For WDSs, testing each of these alternatives
involves the use of a calculation engine, such as EPANET [

20

]. Thus, the computational e

ffort required

is considerable in most cases. For this reason, the number of iterations executed to obtain the final
design of the WDS may be considered as an indicator of the e

fficiency of the tested methodologies.

/var/www/pavco.com.co/public/site/pdftohtml/5e81077962f627b67949947c14bfa4ef/index-html.html
background image

Water 2020, 12, 1037

3 of 26

One approach overshadowed by the rapid growth of metaheuristic techniques was based on the

analysis of optimized designs for simple networks, to produce criteria that, combined with hydraulic

principles, could produce optimized results in more complex systems. One of the first examples of this

approach was Wu [

21

], who analyzed the behavior of the cost function in a series of pipes with uniform

demands, establishing a criterion that indicates the target head-loss value for each pipe in the system
to produce a minimum cost design. Featherstone and El-Jumaily [

22

] extended Wu’s criterion by

suggesting its application to looped networks using the Euclidean distance as the main factor to assign
the target Hydraulic Grade Line (HGL) at each node and, therefore, the target head-loss value per pipe.
Based on the contributions of the aforementioned authors, the Optimal Power Use Surface (OPUS)
methodology was first presented by Takahashi et al. [

23

], as a design methodology based on the target

HGL criterion and the use of other hydraulic criteria to extend the concept to looped networks.

This paper describes the OPUS methodology for the optimal design of WDSs and tests its

performance in di

fferent benchmark networks (Hanoi, Balerma, Taichung and Pescara). Several

operational criteria, such as water quality, reliability and environmental impact, can be considered
in an optimization problem related to WDS design and management. However, as a first step in
the development on this alternative type of design method, the objective function considered in this
manuscript is reduced to the construction costs of the system. A comparison against well tested
metaheuristic algorithms shows promising results, both in the performance of the algorithm, as well
as in the understanding of the physical principles of WDS design. The latter makes OPUS a suitable
method for use in decision making at an urban scale, and can be easily understood by engineers.
Finally, based on the obtained results, conclusions and future work guidelines are stated.

OPUS methodology, described in this paper, has become the basis for several applications to

demand-driven models considering di

fferent topologies [

24

], and in pressure-driven models [

25

].

2. Optimal Hydraulic Grade Line (HGL) Criteria

Wu [

21

] analyzed the HGL generated by the global optimum design of pipes in a series system

(drip irrigation main lines). He found that for flat topography, the optimal HGL was always a convex

function. To characterize the function, Wu [

21

] used the straight HGL between the source and the last

demand node, and defined the “sag” as the distance between that straight line and the function at
the middle section. For his case study with uniform demands, the optimal sag was 15% of the total
available head.

Using the linear programming method described by Alperovits and Shamir [

8

], it is possible to

extend Wu’s [

21

] results for di

fferent demand distributions along the series of pipes, as well as different

topography. As can be seen in Figure

1

, the optimal HGL is always a convex function or a combination

of adjacent convex functions with limits at nodes with high elevation (i.e., higher than the straight

HGL). In these figures, the minimum HGL is composed of the elevation along the pipe series and the

minimum pressure allowable. Hence, these figures show di

fferent combinations between topographic

profiles and demand patterns.

/var/www/pavco.com.co/public/site/pdftohtml/5e81077962f627b67949947c14bfa4ef/index-html.html
background image

Water 2020, 12, 1037

4 of 26

Water 202012, x FOR PEER REVIEW 

4 of 26 

 

 

Figure 1. Comparison of Hydraulic Grade Lines (HGL) resulting from Different Demand 
Configurations and Minimum HGL. 

A quadratic equation can be adjusted to the convex function found by Wu [21], as shown in 

Equation (1): 

𝐻𝐺𝐿(𝑥) =  4𝐹 ×

(𝐻𝐺𝐿

− 𝐻𝐺𝐿

)

𝐿

× (𝑥 ) − (1 + 4𝐹) ×

(𝐻𝐺𝐿

− 𝐻𝐺𝐿

)

𝐿

× 𝑥 + 𝐻𝐺𝐿

 

(1) 

where, HGL(x) is the hydraulic grade line at a distance x from the source, HGL

Max

 is the hydraulic 

grade line at the source, HGL

Min

 is the HGL at the final node and represents the minimum value 

allowable, F is the predefined sag as a percentage of (HGL

Max

 − HGL

Min

), and L is the total length of 

the series [26]. 

Featherstone and El-Jumaily [22] extended the concept of optimal hydraulic grade line to lopped 

networks by setting x as the Euclidean distance from each node to the source, and L as the maximum 
value of x in the network (i.e., the furthest node). Using the Euclidean distance is a quick approach to 
extend Wu’s concept, but clearly the water needs to flow through longer distances according to the 
layout of the network. The OPUS methodology introduces a new method to define x and L using the 
topological distance in looped networks according to a cost-benefit optimal spanning tree.  

3. Methodology  

The Optimal Power Use Surface (OPUS) design methodology is composed of six steps, which 

are focused on understanding the flow of energy along the network. These steps are summarized in 
the flowchart shown in Figure 2, and described in detail in the subsequent sections. 

Figure 1.

Comparison of Hydraulic Grade Lines (HGL) resulting from Di

fferent Demand Configurations

and Minimum HGL.

A quadratic equation can be adjusted to the convex function found by Wu [

21

], as shown in

Equation (1):

HGL

(

x

) =

4F ×

(

HGL

max

− HGL

min

)

L

2

×

x

2

(

1

+

4F

)

×

(

HGL

max

− HGL

min

)

L

× x

+

HGL

max

(1)

where, HGL(x) is the hydraulic grade line at a distance x from the source, HGL

Max

is the hydraulic grade

line at the source, HGL

Min

is the HGL at the final node and represents the minimum value allowable, F

is the predefined sag as a percentage of (HGL

Max

− HGL

Min

), and L is the total length of the series [

26

].

Featherstone and El-Jumaily [

22

] extended the concept of optimal hydraulic grade line to lopped

networks by setting x as the Euclidean distance from each node to the source, and L as the maximum

value of x in the network (i.e., the furthest node). Using the Euclidean distance is a quick approach to

extend Wu’s concept, but clearly the water needs to flow through longer distances according to the
layout of the network. The OPUS methodology introduces a new method to define x and L using the
topological distance in looped networks according to a cost-benefit optimal spanning tree.

3. Methodology

The Optimal Power Use Surface (OPUS) design methodology is composed of six steps, which are

focused on understanding the flow of energy along the network. These steps are summarized in the
flowchart shown in Figure

2

, and described in detail in the subsequent sections.

/var/www/pavco.com.co/public/site/pdftohtml/5e81077962f627b67949947c14bfa4ef/index-html.html
background image

Water 2020, 12, 1037

5 of 26

Water 202012, x FOR PEER REVIEW 

5 of 26 

 

 

Figure 2. Optimal Power Use Surface (OPUS) methodology flow chart. 

3.1. Sump Search or Spannin Tree Structure 

The sump search, or the identification of the Spanning Tree structure, is based on two hydraulic 

principles: (1) In a least-cost WDS, the demand at the nodes should be supplied by transporting the 
water from the sources using a single route, and (2) As the pipes flow rates for design purposes 
increase, the marginal costs diminish.  

For the first principle, it can be stated that redundancy is hydraulically inefficient, despite its 

benefits in terms of reliability. As a result, branched WDSs can manage to be less expensive compared 
to looped networks; thus, this step aims to decompose the looped system into an open tree-like 
structure. This leads to the identification of nodes in the original network that can be considered as 
sumps of a fictional branched network, meaning the nodes with the lowest hydraulic head among all 
its neighbors. 

The second principle is based in the analysis of the relationship between the pipe diameter and 

its cost per length unit. A typical equation to represent the constructive cost of a pipe is given by 

Equation (2)

𝐶𝑜𝑠𝑡 = 𝐾 × 𝐿 × 𝑑

 

(2)

 

where, 𝐶𝑜𝑠𝑡  is the cost of pipe 𝑖; 𝐾 is the coefficient of the cost function; 𝐿  is the length of pipe 𝑖; 𝑑  
is the diameter assigned to pipe 𝑖; and 𝑥 is the exponent of the cost function [21]. On the other hand, 
Hazen-Williams friction head losses equation relates to pipe diameters, friction head losses and 
flowrates as shown in Equation (3): 

ℎ =

10.67 × 𝐿 × 𝑄

.

𝐶

.

× 𝑑

.

 

(3)

 

Figure 2.

Optimal Power Use Surface (OPUS) methodology flow chart.

3.1. Sump Search or Spannin Tree Structure

The sump search, or the identification of the Spanning Tree structure, is based on two hydraulic

principles: (1) In a least-cost WDS, the demand at the nodes should be supplied by transporting the
water from the sources using a single route, and (2) As the pipes flow rates for design purposes increase,

the marginal costs diminish.

For the first principle, it can be stated that redundancy is hydraulically ine

fficient, despite its

benefits in terms of reliability. As a result, branched WDSs can manage to be less expensive compared to
looped networks; thus, this step aims to decompose the looped system into an open tree-like structure.

This leads to the identification of nodes in the original network that can be considered as sumps of a

fictional branched network, meaning the nodes with the lowest hydraulic head among all its neighbors.

The second principle is based in the analysis of the relationship between the pipe diameter and

its cost per length unit. A typical equation to represent the constructive cost of a pipe is given by
Equation (2):

Cost

i

=

K × L

i

× d

x
i

(2)

where, Cost

i

is the cost of pipe i; K is the coe

fficient of the cost function; L

i

is the length of pipe i;

d

i

is the diameter assigned to pipe i; and x is the exponent of the cost function [

21

]. On the other

hand, Hazen-Williams friction head losses equation relates to pipe diameters, friction head losses and
flowrates as shown in Equation (3):

h

f i

=

10.67 × L

i

× Q

1.85
i

C

1.85

× d

4.87
i

(3)

/var/www/pavco.com.co/public/site/pdftohtml/5e81077962f627b67949947c14bfa4ef/index-html.html
background image

Water 2020, 12, 1037

6 of 26

where, h

f i

is the head-loss along pipe i (m), Q

i

represents the volumetric flowrate on pipe i (m

3

/s

)

and

C is the pipe roughness coe

fficient (dimensionless). Therefore, by solving Equation (3) for d

i

and then

replacing it in Equation (2), an expression for the cost of the pipe as a function of Q

i

is obtained, as

shown in Equation (4):

CostQ

(

Q

i

) =

Cost

i

=

K

0

× L

i

×

Q

i

S

0.54

f i

x/2.63

(4)

where, K

0

=

K ×

10.67

C

1.85

x/4.87

; and S

f i

is the friction slope for pipe i defined by h

f i

/L

i

. Figure

3

shows

the trend of this new expression, where it is shown that for typical values of x (lower than 2.63),
the marginal cost decreases as the design flowrate of a pipe increases.

Water 202012, x FOR PEER REVIEW 

6 of 26 

 

where, ℎ  is the head-loss along pipe 𝑖 (m), 𝑄  represents the volumetric flowrate on pipe 𝑖 (𝑚 /𝑠) 
and 𝐶 is the pipe roughness coefficient (dimensionless). Therefore, by solving Equation (3) for 𝑑  and 
then replacing it in Equation (2), an expression for the cost of the pipe as a function of 𝑄  is obtained, 
as shown in Equation (4): 

𝐶𝑜𝑠𝑡𝑄(𝑄 ) = 𝐶𝑜𝑠𝑡 = 𝐾

×

𝐿

×

𝑄

𝑆

.

/ .

 

(4)

 

where, 𝐾 = 𝐾 ×

.

.

/ .

; and 𝑆  is the friction slope for pipe 𝑖 defined by ℎ /𝐿 . Figure 3 shows 

the trend of this new expression, where it is shown that for typical values of 𝑥 (lower than 2.63), the 
marginal cost decreases as the design flowrate of a pipe increases. 

 

Figure 3. Schematic relation between pipe cost and flow. 

Considering that WDS are represented as a set of junctions and links, the problem of obtaining 

a tree-like structure for the WDS is similar to the problem of finding a Minimum Spanning Tree (MST) 
for a weighted undirected graph. In this problem, the vertices of the graph are the nodes of the WDS, 
while the edges of the graph are the pipes (links) in the WDS. Prim [27] proposed an algorithm to 
find the MST that minimizes the sum of link lengths. For the OPUS methodology, Prim’s Algorithm 
can be modified to maximize a recursively updated Benefit-Cost (𝐵/𝐶) value for every link in the 
minimum spanning tree. This is equivalent to finding the MST for a graph weighted with the 𝐵/𝐶 
values, but with the 𝐵/𝐶 values updated at each iteration of Prim’s Algorithm. 

The spanning tree is initialized at the water source and grows by incorporating nodes from the 

selection front (set of adjacent nodes that are not yet connected to the spanning tree), ensuring the 
connection performed reaches the highest 𝐵/𝐶  value possible among its set. This procedure is 
repeated until every node (demand and not demand nodes) in the network is connected to the tree. 

If there is more than one water source in the network, there will be the same number of trees as 

sources in the system, since there should exist only one route to supply each node in the spanning 
tree structure. In that case, a tree will be initialized per water source, and the selection front will 
include any node not yet connected to a tree but adjacent to one. 

The 𝐵/𝐶 function of adding a pipe-node pair 𝑖 in the selection front, is calculated as the ratio 

between the demand of the node added to the tree and the marginal cost associated with its 
connection to the source [26], as shown in Equation (5). This marginal cost is computed as the sum of 
the total cost of the pipe 𝑖 and the cost difference related to the transportation of the additional flow 
through all of the upstream pipes [28].  

(𝐵/𝐶) =

𝑂𝑢𝑡𝑓𝑙𝑜𝑤

𝐶𝑜𝑠𝑡𝑄(𝑂𝑢𝑡𝑓𝑙𝑜𝑤 ) + ∑

𝐶𝑜𝑠𝑡𝑄 𝐹𝑙𝑜𝑤𝑅𝑎𝑡𝑒 + 𝑂𝑢𝑡𝑓𝑙𝑜𝑤 − 𝐶𝑜𝑠𝑡𝑄(𝐹𝑙𝑜𝑤𝑅𝑎𝑡𝑒 )

 ∈ 

 

(5)

 

where, (𝐵/𝐶)  is the Benefit–Cost value of adding the attached pipe-node pair 𝑖. 𝑂𝑢𝑡𝑓𝑙𝑜𝑤  is the flow 
demand of node 𝑖, which is the additional flow conveyed through the network, if pair 𝑖 was added 

Figure 3.

Schematic relation between pipe cost and flow.

Considering that WDS are represented as a set of junctions and links, the problem of obtaining a

tree-like structure for the WDS is similar to the problem of finding a Minimum Spanning Tree (MST)
for a weighted undirected graph. In this problem, the vertices of the graph are the nodes of the WDS,

while the edges of the graph are the pipes (links) in the WDS. Prim [

27

] proposed an algorithm to

find the MST that minimizes the sum of link lengths. For the OPUS methodology, Prim’s Algorithm
can be modified to maximize a recursively updated Benefit-Cost (B/C) value for every link in the
minimum spanning tree. This is equivalent to finding the MST for a graph weighted with the B/C

values, but with the B/C values updated at each iteration of Prim’s Algorithm.

The spanning tree is initialized at the water source and grows by incorporating nodes from the

selection front (set of adjacent nodes that are not yet connected to the spanning tree), ensuring the
connection performed reaches the highest B/C value possible among its set. This procedure is repeated
until every node (demand and not demand nodes) in the network is connected to the tree.

If there is more than one water source in the network, there will be the same number of trees as

sources in the system, since there should exist only one route to supply each node in the spanning tree
structure. In that case, a tree will be initialized per water source, and the selection front will include
any node not yet connected to a tree but adjacent to one.

The B/C function of adding a pipe-node pair i in the selection front, is calculated as the ratio

between the demand of the node added to the tree and the marginal cost associated with its connection
to the source [

26

], as shown in Equation (5). This marginal cost is computed as the sum of the total cost

of the pipe i and the cost di

fference related to the transportation of the additional flow through all of

the upstream pipes [

28

].

(

B/C

)

i

=

Out f low

i

CostQ

(

Out f low

i

) +

P

j ∈ Upstream

pipes

h

CostQ

FlowRate

j

+

Out f low

i

− CostQ

FlowRate

j

i

(5)

where,

(

B/C

)

i

is the Benefit–Cost value of adding the attached pipe-node pair i. Out f low

i

is the flow

demand of node i, which is the additional flow conveyed through the network, if pair i was added

/var/www/pavco.com.co/public/site/pdftohtml/5e81077962f627b67949947c14bfa4ef/index-html.html
background image

Water 2020, 12, 1037

7 of 26

to the tree. FlowRate

j

is the flow that goes through pipe j of the tree before inserting any adjacent

pipe-node pairs. The first term in the denominator of Equation (5) stands for the cost of the new

included pipe, while the sum corresponds to the marginal cost upstream [

26

]. This expression is

computed using Equation (4) as the function CostQ.

Given Pescara [

29

], the Italian WDS benchmark shown in Figure

4

a, the sump search procedure

will be developed as follows: Starting from the three sources (identified as 15, 43 and 65), each selection

front is analyzed. In other words, the possible link-node pairs that can be added to Reservoir 65
are

<90,76> (i.e., junction 90 and link 76) and <6,89>. Respectively, the link-node pairs <42,54> and

<14,11> could be added to sources 43 and 15. To decide which source will be the root of the first tree,

and consequently, which link–node pair will be added initially to the tree, the B/C value is computed
for each link–node pair, selecting the one that maximizes this criterion. This procedure is repeated
until every node is connected to a tree. In Figure

4

b, the minimum cost-spanning tree is shown for

Pescara WDS, where each link maximized the B/C value at each selection front. Dashed links are those

who were not included in the spanning tree, but are present in the original WDS.

Water 202012, x FOR PEER REVIEW 

7 of 26 

 

to the tree. 𝐹𝑙𝑜𝑤𝑅𝑎𝑡𝑒  is the flow that goes through pipe 𝑗 of the tree before inserting any adjacent 
pipe-node pairs. The first term in the denominator of Equation (5) stands for the cost of the new 
included pipe, while the sum corresponds to the marginal cost upstream [26]. This expression is 
computed using Equation (4) as the function 𝐶𝑜𝑠𝑡𝑄. 

Given Pescara [29], the Italian WDS benchmark shown in Figure 4a, the sump search procedure 

will be developed as follows: Starting from the three sources (identified as 15, 43 and 65), each 
selection front is analyzed. In other words, the possible link-node pairs that can be added to Reservoir 
65 are <90,76> (i.e., junction 90 and link 76) and <6,89>. Respectively, the link-node pairs <42,54> and 
<14,11> could be added to sources 43 and 15. To decide which source will be the root of the first tree, 
and consequently, which link–node pair will be added initially to the tree, the 𝐵/𝐶 value is computed 
for each link–node pair, selecting the one that maximizes this criterion. This procedure is repeated 
until every node is connected to a tree. In Figure 4b, the minimum cost-spanning tree is shown for 
Pescara WDS, where each link maximized the 𝐵/𝐶 value at each selection front. Dashed links are 
those who were not included in the spanning tree, but are present in the original WDS. 

 

(a) 

 

(b) 

Figure 4. (a) Layout of Pescara Water Distribution System (WDS). Pipe ID shown in the labels above 
the pipes (b) Minimum-Cost Spanning Tree. Node ID shown in the labels. Dashed links correspond 
to the pipes that do not form part of the spanning tree. 

Figure 4.

(a) Layout of Pescara Water Distribution System (WDS). Pipe ID shown in the labels above

the pipes (b) Minimum-Cost Spanning Tree. Node ID shown in the labels. Dashed links correspond to
the pipes that do not form part of the spanning tree.

/var/www/pavco.com.co/public/site/pdftohtml/5e81077962f627b67949947c14bfa4ef/index-html.html
background image

Water 2020, 12, 1037

8 of 26

In particular, it must be noted that as the Pescara network is supplied by three di

fferent sources, it

is expected that this tree-structure identification step will reach three di

fferent branched sub-systems,

each one fed by a unique source. In general terms, for systems with n sources, the nodes in the original
network will be grouped into a maximum number of n di

fferent tree structures. This is accomplished

by seeking to maximize the B/C relationship and ensuring that every node in the original network is
assigned to a tree, except for the reservoirs that can be isolated given their B/C relationship. Hence,
for further steps in this methodology, the obtained trees are considered simultaneously.

Moreover, the values involved in this step are not actual costs, but values obtained proportionally

from the relation presented in Figure

3

. Configuring the tree using this cost-benefit function has an

O(NN

2

) time complexity, where NN is the number of nodes. At the end, the status of the last nodes in

each branch is set to ‘sumps’, which will be useful for the rest of the methodology [

28

].

3.2. Optimal Power Use Surface (OPUS)

The entire methodology is named after this specific step due to the close relation with the work

developed by Wu [

21

]. This step consists of assigning a target HGL to each node in the network,

which determines the head-losses for each pipe in advance [

26

].

In the proposed algorithm, the OPUS is computed using the tree structure, described in the

previous section. Here, the Wu [

21

] concept of target HGL is used by setting the pressure of all of

the sump nodes to the minimum admissible value [

26

], and calculating the intermediate nodes’ head

for each path knowing the head in each reservoir. For this purpose, each node’s topological distance,
measured along the path of pipes that connect the node with the source, must be computed previously.

Then Equation (1) is used to compute the target HGL at every node.

Since Equation (1) depends on the parameter F (sag), a sensitivity analysis to evaluate its e

ffect in

the final design cost was performed. The evaluated range was F ∈

[

0, 0.25

]

. In order to represent a

straight HGL, the minimum value was set to 0, since a negative sag results in a concave down parabolic

HGL, opposite to Wu’s criterion [

26

]. The maximum value of the sag was obtained by finding the

maximum F in Equation (1) that produces a monotonically decreasing HGL for x ∈

[

0, L

]

. It can be

shown that a value of 0.25 is the maximum sag under this consideration (no increasing HGL).

The results of the sensitivity analysis showed that a second-order polynomial could be adjusted to

the cost vs. F relationship. Therefore, the design must be done three times with continuous diameters
and the costs must be calculated to find the optimal target sag of a network. The criteria considered for
the three designs are identical, apart from the general sag used. Subsequently, taking into account that
cost vs. sag is a parabola, the optimal sag is determined by finding the minimum cost of a second-order

polynomial regression obtained by the three cost vs. sag points previously calculated [

26

].

If the three designs are made using F

=

{0, 0.1, 0.25} then the optimal target sag can be computed

using Equation (6):

Optimal sag

=

21C

0

− 25C

0.1

+

4C

0.25

40

(

3C

0

− 5C

0.1

+

2C

0.25

)

(6)

where, C

x

is the cost of the architecture obtained using a sag F

=

x [

26

]. The procedure corresponding

to this step is described in the pseudo-code shown in Figure

5

.

/var/www/pavco.com.co/public/site/pdftohtml/5e81077962f627b67949947c14bfa4ef/index-html.html
background image

Water 2020, 12, 1037

9 of 26

Water 202012, x FOR PEER REVIEW 

9 of 26 

 

 

Figure 5. Algorithm for OPUS subprocess. 

The target sag (𝐹) must be recalculated at each intersection in an upstream analysis through the 

network, where the branches converge. In order to do so, the flow is weighted on each downstream 
course. It is essential to analyze the nodes with a high elevation, in order to avoid the assignment of 
a head value that does not fulfill the pressure demand [26]. 

The optimal power use surface for the Pescara network is shown in Figure 6 as an example. The 

terrain of the network is shown in Figure 6a, while subsequent figures show the level curves for the 
HGL (Figure 6b), and its corresponding surface (Figure 6c). Note that, as mentioned before, for the 
water source, the HGL corresponds to the total head available in the reservoir, while in the sump 
nodes its value was assigned equal to the minimum pressure.  

 

 

(a) (b

 

(c

Figure 6. (a) Topographic profile for Pescara. (b) Level curve for the HGL. (c) Assigned surface for the 
Pescara network. 

Figure 5.

Algorithm for OPUS subprocess.

The target sag

(

F

)

must be recalculated at each intersection in an upstream analysis through the

network, where the branches converge. In order to do so, the flow is weighted on each downstream
course. It is essential to analyze the nodes with a high elevation, in order to avoid the assignment of a
head value that does not fulfill the pressure demand [

26

].

The optimal power use surface for the Pescara network is shown in Figure

6

as an example. The

terrain of the network is shown in Figure

6

a, while subsequent figures show the level curves for the

HGL (Figure

6

b), and its corresponding surface (Figure

6

c). Note that, as mentioned before, for the

water source, the HGL corresponds to the total head available in the reservoir, while in the sump nodes

its value was assigned equal to the minimum pressure.

Water 202012, x FOR PEER REVIEW 

9 of 26 

 

 

Figure 5. Algorithm for OPUS subprocess. 

The target sag (𝐹) must be recalculated at each intersection in an upstream analysis through the 

network, where the branches converge. In order to do so, the flow is weighted on each downstream 
course. It is essential to analyze the nodes with a high elevation, in order to avoid the assignment of 
a head value that does not fulfill the pressure demand [26]. 

The optimal power use surface for the Pescara network is shown in Figure 6 as an example. The 

terrain of the network is shown in Figure 6a, while subsequent figures show the level curves for the 
HGL (Figure 6b), and its corresponding surface (Figure 6c). Note that, as mentioned before, for the 
water source, the HGL corresponds to the total head available in the reservoir, while in the sump 
nodes its value was assigned equal to the minimum pressure.  

 

 

(a) (b

 

(c

Figure 6. (a) Topographic profile for Pescara. (b) Level curve for the HGL. (c) Assigned surface for the 
Pescara network. 

Figure 6.

(a) Topographic profile for Pescara. (b) Level curve for the HGL. (c) Assigned surface for the

Pescara network.

/var/www/pavco.com.co/public/site/pdftohtml/5e81077962f627b67949947c14bfa4ef/index-html.html
background image

Water 2020, 12, 1037

10 of 26

Once this step is concluded, all nodes must have an objective head value assigned. Thus, a flow

value is required in order to calculate the diameter of each pipe in the network, using the process

explained in the next section.

3.3. Optimal Flow Distribution

In a looped network, the same hydraulic gradient surface can be obtained by combining di

fferent

configurations of diameters when the set of allowable pipe sizes is R

+

(positive real numbers) [

24

].

Thus, this step is intended to find a unique flow distribution in the network that fits to the Optimal

Power Use Surface (OPUS) previously defined, ensuring it satisfies the mass conservation at each node.

Since at this point the objective HGL has been assigned using Wu’s [

16

] criteria with the topological

distance, the spanning tree structure is no longer required. Therefore, for further steps, the described

procedures are done using the original graph (network).

The determination of the optimal flow distribution consists of the division of the flow demand into

the upstream pipes, starting from the sumps. Three di

fferent criteria were explored for this procedure

by Saldarriaga et al. [

19

]:

(1)

Uniform distribution: It assumes that all pipes have the same flow, which is calculated by
dividing the total flow demand of each node into the number of upstream pipes connected to
it [

26

].

(2)

Proportional distribution: For each pipe, the flow distributes proportionally to H

/L

2

, where H

stands for the head losses in the pipe and L for its length.

(3)

All-in-one distribution: The conveyance capacity for all the pipes upstream of the analyzed node
is computed assuming they have the minimum diameter available (d

min

), and they are assigned

as the design flows. If this is insu

fficient to transport the total flow demand, the residual flow is

assigned to the one pipe with highest hydraulic favorability.

Depending on the used criterion, the reliability of the network varies. In this context, the reliability

can be defined as the ability of the system to handle failures, and it can be measured through several
approaches such as the Network Resilience Index (NRI), proposed by Prasad and Park [

30

]. Saldarriaga

et al. [

24

] proved that the reliability of the network is magnified when the uniform distribution is used,

contrasted to the other two criteria. However, in this paper only the All-in-one distribution is considered
as it was found to be the one that produced least-cost networks [

21

].

There are multiple criteria to analyze the appropriateness of di

fferent pipes in order to define the

main pipe, which is the one that transports the largest portion of the total flow demanded. Similar to
the tree structure step, the pipe with the highest flow is determined through a function. This expression
is Equation (7), which refers to the Hydraulic Favorability, where the pipe with the maximum value of
this parameter is chosen. For nonsump nodes, the total demand is computed by adding its flow to the
one required downstream. An Iterative-Recursive Algorithm (IRA) can be applied to perform all of the
calculations with an O(NN) time complexity [

26

].

Hydraulic Favorability

i

=

HGL

up

− HGL

down

L

2

(7)

where, HGL

up

is the hydraulic grade line of the upstream node of pipe i, HGL

down

is the HGL of the

downstream node of pipe i, and L is the length of pipe i [

26

]. At the end of the process, all of the pipes

in the system must have been assigned an objective flow value. The subprocess is explained in detail
in Figure

7

.

/var/www/pavco.com.co/public/site/pdftohtml/5e81077962f627b67949947c14bfa4ef/index-html.html
background image

Water 2020, 12, 1037

11 of 26

Water 202012, x FOR PEER REVIEW 

11 of 26 

 

 

Figure 7. Algorithm for the Optimal Flow Distribution. 

3.4. Continuous Diameter Calculation 

Once the assigned head-losses and the objective flow rate for each pipe are defined, the diameter 

required in a continuous range as the hydraulic design can be easily calculated in an explicit manner 
using the Hazen-Williams equation (solving for 𝑄 on Equation (3)) and iteratively with the Darcy-
Weisbach and Colebrook-White equations (Equations (8) and (9), respectively) [26]. Then, the 
obtained continuous diameter must be rounded to an available size to satisfy the constraints of the 
problem and turn an ideal optimal design into a low-cost feasible one. 

In this equation, 𝑄 represents the volumetric flowrate (m /s), ℎ  is the head loss over the length 

of pipe (m), 𝐶 is the pipe roughness coefficient (dimensionless), 𝑑 is the diameter of the pipe (m), and 
𝐿 is the length of the pipe (m): 

𝑄 =

ℎ × 𝑔 × 𝜋 × 𝑑

8𝑓 × 𝐿

/

 

(8)

 

where 𝑔 is the acceleration due to gravity (m/s ), and 𝑓 is the Darcy friction factor (dimensionless). 

1

𝑓

= −2 log

𝑘

3.7𝑑

+

2.51

𝑅𝑒 𝑓

  

(9)

 

where, 𝑘  is the roughness height of the pipe (m), and 𝑅𝑒 is the Reynolds Number (dimensionless). 

3.5. Diameter Round-off 

As described before, the continuous diameters obtained in the previous step must be rounded 

to satisfy the commercial constraints of the problem. Therefore, in this step the continuous diameter 
of each pipe is rounded to a discrete value from the list of commercially available diameter sizes, 
which is represented by the set 

𝑫

=

{

𝐷

1

, … , 𝐷

𝑁𝐷

}

 [26]. 

During this research, it has been found that using approaches for approximating the diameter 

based in either the nearest equivalent flow, as well as the nearest equivalent head-loss led to 
promising results. The latter is accomplished by raising the continuous diameter values to a power 

Figure 7.

Algorithm for the Optimal Flow Distribution.

3.4. Continuous Diameter Calculation

Once the assigned head-losses and the objective flow rate for each pipe are defined, the diameter

required in a continuous range as the hydraulic design can be easily calculated in an explicit
manner using the Hazen-Williams equation (solving for Q on Equation (3)) and iteratively with the
Darcy-Weisbach and Colebrook-White equations (Equations (8) and (9), respectively) [

26

]. Then,

the obtained continuous diameter must be rounded to an available size to satisfy the constraints of the

problem and turn an ideal optimal design into a low-cost feasible one.

In this equation, Q represents the volumetric flowrate (m

3

/s

)

, h

f

is the head loss over the length

of pipe (m), C is the pipe roughness coe

fficient (dimensionless), d is the diameter of the pipe (m), and L

is the length of the pipe (m):

Q

=

h

f

× g ×

π

2

× d

5

8 f × L

1/2

(8)

where g is the acceleration due to gravity (m/s

2

)

, and f is the Darcy friction factor (dimensionless).

1

p f

=

−2 log

10

k

s

3.7d

+

2.51

Re

p f

(9)

where, k

s

is the roughness height of the pipe (m), and Re is the Reynolds Number (dimensionless).

3.5. Diameter Round-o

ff

As described before, the continuous diameters obtained in the previous step must be rounded to

satisfy the commercial constraints of the problem. Therefore, in this step the continuous diameter of
each pipe is rounded to a discrete value from the list of commercially available diameter sizes, which is
represented by the set D

=

{

D

1

,

. . . , D

ND

} [

26

].

During this research, it has been found that using approaches for approximating the diameter

based in either the nearest equivalent flow, as well as the nearest equivalent head-loss led to promising
results. The latter is accomplished by raising the continuous diameter values to a power of 2.6 for

/var/www/pavco.com.co/public/site/pdftohtml/5e81077962f627b67949947c14bfa4ef/index-html.html
background image

Water 2020, 12, 1037

12 of 26

the equivalent flow criterion (as explained in the Tree Structure step), or using a power of −4.87 in
the case of the equivalent head-loss criterion. The aforementioned powers describe the relationship
between the head-loss and the diameter, according to the Hazen-Williams equation calculated using the
International System of Units. However, this action has a negative e

ffect on the hydraulic performance

of the network, requiring an additional step to post-optimize the current network.

3.6. Optimization

Finally, the optimization step has two central purposes: (1) To guarantee that all nodes have a

pressure above the minimum required value; and (2) To find possible reductions in the capital cost of

the network. The first objective is achieved by increasing diameters (if needed), setting as a starting

point the ones with larger unit head-loss di

fferences between real and target values [

26

]. The process is

repeated until all the nodes in the system meet the minimum pressure requirement.

The second goal is reached with a double review over the pipes, decreasing the diameters if no

constraint is violated. This two-way review starts at the reservoirs, heading towards the sumps in the
direction of the flow, and then backwards. This means that the need of reducing the pipe’s diameter
is verified twice. If any variation caused a pressure deficit in the system, it is reversed immediately,
otherwise it remains [

28

].

Given that the diameter of one pipe is increased iteratively while some nodes are below the

pressure demand, this step requires the highest percentage of iterations of the whole method because a

hydraulic simulation must be run per pipe for each diameter variation. This single heuristic can be
used alone to accomplish feasible designs, even though, it is strongly dependent on the initial pipe
configuration [

28

].

4. Results and Discussion

The OPUS methodology was tested and validated using four well-known benchmark networks:

Hanoi [

31

], Balerma [

32

], Taichung [

33

] and Pescara [

29

]. Di

fferent configurations of the parameters

defined on each step of the methodology were tested, looking for optimal designs with discrete
diameters (assigning only diameters commercially available). In some cases, there were assigned
continuous values (D

=

R

+

), where is the diameter of each pipe, to be considered as a theoretical

solution giving the best costs that may be found with the methodology without the restriction of
continuous diameter values. The input files corresponding to the optimal designs reported below are
available in the supplementary files of this research.

4.1. Hanoi

The Hanoi network was first introduced by Fujiwara and Khang [

31

], and it has become a

well-known benchmark WDS for the comparison of design methodologies since its appearance. For

this network the Hazen-Williams equation is used to calculate the friction losses; a Hazen-Williams
coe

fficient of C

=

130 is the value for this benchmark network. For the design process a 30 m minimum

pressure was established and in order to calculate the network total cost, it is necessary to add the costs

of all pipes using a potential function of each diameter. That function has a unit coe

fficient of $1.1/m

and 1.5 as the exponent. The set of available design diameters are 12, 16, 20, 24, 30 and 40 inches [

26

],

and their corresponding unit costs are shown in Table

1

. The network is composed of 34 pipes and

31 nodes configured in 3 loops. The whole system is supplied by 1 reservoir with a constant head of

100 m. A sag value of 0.25 was used in OPUS methodology, while both the layout of the network and

the obtained spanning tree are shown in Figure

8

a.

Table 1.

Unit Costs of Hanoi Network.

Diameter (in.)

12.0

16.0

20.0

24.0

30.0

40.0

Unit Cost ($

/m)

45.73

70.40

98.39

129.33

180.75

278.28

/var/www/pavco.com.co/public/site/pdftohtml/5e81077962f627b67949947c14bfa4ef/index-html.html
background image

Water 2020, 12, 1037

13 of 26

Water 202012, x FOR PEER REVIEW 

13 of 26 

 

 

(a) 

 

(b) 

Figure 8. Cont.

/var/www/pavco.com.co/public/site/pdftohtml/5e81077962f627b67949947c14bfa4ef/index-html.html
background image

Water 2020, 12, 1037

14 of 26

Water 202012, x FOR PEER REVIEW 

14 of 26 

 

 

(c

Figure 8. (a) General layout and Spanning Tree for Hanoi Network. Dashed lines are pipes from the 
original networks that are not included in the Spanning Tree (b) OPUS design considering continuous 
diameters (in inches) (c) OPUS design considering discrete diameters (in inches). 

In OPUS, the Hanoi network was tested using continuous diameters, reaching a total cost of 

$5,535,888 and the pipe configuration shown in Figure 8b. As can be seen in the mentioned figure, 
higher values than the maximum allowable diameter (40 inches according to Fujiwara and Khang 
(1990)), are assigned to some pipes of the network. In order to consider the influence of the diameter 
restrictions on the final cost of the network, two different discrete designs were performed: the first 
one based on the original diameter list, and the second one considering the availability of a 50 inches 
diameter. 

In the first scenario considering the original diameter list, the OPUS algorithm generated a total 

cost of $6,374,525 after 106 iterations, with the resulting pipe diameters (in inches) shown in Figure 
8c. These iterations correspond to the total number of hydraulic simulations executed by the 
algorithm. 

Although this is not the least cost reported, the number of hydraulic simulations needed to reach 

this result is three orders of magnitude smaller than that of other approaches, as is shown in Table 2. 
In regards to Hanoi, the least cost design used as benchmark was $6.081 million, given that lower 
cost results do not satisfy pressure constraints in EPANET2 as shown by [34]. To this end, OPUS 
reached a solution that has a difference of 4.81% with respect to the lowest value reported in a very 
low number of iterations. 

Table 2. Reported costs and number of iterations for the Hanoi WDS. 

Algorithm 

Cost (millions) 

Number of Iterations 

Genetic Algorithm [10] 

$6.073 

1,000,000 

Simulated annealing [11] 

$6.056 

53,000 

Harmony search [35] 

$6.056 

200,000 

Shuffled frog leaping [36] 

$6.073 

26,987 

Shuffled complex evolution [37] 

$6.220 

25,402 

Genetic Algorithm [38] 

$6.056 

18,300 

Figure 8.

(a) General layout and Spanning Tree for Hanoi Network. Dashed lines are pipes from the

original networks that are not included in the Spanning Tree (b) OPUS design considering continuous
diameters (in inches) (c) OPUS design considering discrete diameters (in inches).

In OPUS, the Hanoi network was tested using continuous diameters, reaching a total cost of

$5,535,888 and the pipe configuration shown in Figure

8

b. As can be seen in the mentioned figure,

higher values than the maximum allowable diameter (40 inches according to Fujiwara and Khang

(1990)), are assigned to some pipes of the network. In order to consider the influence of the diameter

restrictions on the final cost of the network, two di

fferent discrete designs were performed: the

first one based on the original diameter list, and the second one considering the availability of a 50
inches diameter.

In the first scenario considering the original diameter list, the OPUS algorithm generated a total

cost of $6,374,525 after 106 iterations, with the resulting pipe diameters (in inches) shown in Figure

8

c.

These iterations correspond to the total number of hydraulic simulations executed by the algorithm.

Although this is not the least cost reported, the number of hydraulic simulations needed to reach

this result is three orders of magnitude smaller than that of other approaches, as is shown in Table

2

.

In regards to Hanoi, the least cost design used as benchmark was $6.081 million, given that lower cost
results do not satisfy pressure constraints in EPANET2 as shown by [

34

]. To this end, OPUS reached a

solution that has a di

fference of 4.81% with respect to the lowest value reported in a very low number

of iterations.

In the second scenario considering the 50” diameter, the total cost of the design obtained following

the OPUS algorithm was of only $5,533,393, due to the closer resemblance between the final hydraulic
gradient surface and the target one, as shown in Figure

9

. In this approach, the unit cost of the 50” pipe

was extrapolated, resulting in a value of $388.91

/m. Finally, the resulting diameters listed by their ID

shown in Figure

8

a are 50, 50, 40, 40, 40, 40, 40, 24, 24, 40, 20, 20, 12, 12, 12, 16, 20, 20, 20, 40, 16, 12, 40,

30, 24, 20, 12, 12, 12, 12, 12, 12, 12 and 20. These diameters are measured in inches.

/var/www/pavco.com.co/public/site/pdftohtml/5e81077962f627b67949947c14bfa4ef/index-html.html
background image

Water 2020, 12, 1037

15 of 26

Table 2.

Reported costs and number of iterations for the Hanoi WDS.

Algorithm

Cost (millions)

Number of Iterations

Genetic Algorithm [

10

]

$6.073

1,000,000

Simulated annealing [

11

]

$6.056

53,000

Harmony search [

35

]

$6.056

200,000

Shu

ffled frog leaping [

36

]

$6.073

26,987

Shu

ffled complex evolution [

37

]

$6.220

25,402

Genetic Algorithm [

38

]

$6.056

18,300

Ant colony optimization [

39

]

$6.134

35,433

Genetic Algorithms [

32

]

$6.081

50,000

Particle Swarm Optimization [

40

]

$6.093

6600

Genetic Algorithms [

12

]

$6.173

26,457

Simulated annealing [

12

]

$6.333

26,457

Simulated annealing with tabu search [

12

]

$6.353

26,457

Local search with simulated annealing [

12

]

$6.308

26,457

Harmony search [

13

]

$6.081

27,721

Cross entropy [

15

]

$6.081

97,000

Scatter search [

33

]

$6.081

43,149

Modified GA 1 [

41

]

$6.056

18,000

Modified GA 2 [

41

]

$6.190

18,000

Particle swarm harmony search [

16

]

$6.081

17,980

Heuristic based approach [

42

]

$6.701

70

Di

fferential evolution [

43

]

$6.081

48,724

Honey-bee mating optimization [

44

]

$6.117

15,955

Heuristic based approach [

45

]

$6.232

259

Self-adaptive di

fferential evolution [

46

]

$6.081

60,582

Pseudo-Genetic Algorithm [

34

]

$6.081

25,000

Linear Programming – Di

fferential Evolution [

17

]

$6.081

33,148

B-GENOME (B-GA) [

18

]

$6.182

26,000

SOGH [

47

]

$6.337

94

OPUS (This research)

$6.374

106

Water 202012, x FOR PEER REVIEW 

15 of 26 

 

Ant colony optimization [39] 

$6.134 

35,433 

Genetic Algorithms [32] 

$6.081 

50,000 

Particle Swarm Optimization [40] 

$6.093 

6600 

Genetic Algorithms [12] 

$6.173 

26,457 

Simulated annealing [12] 

$6.333 

26,457 

Simulated annealing with tabu search [12] 

$6.353 

26,457 

Local search with simulated annealing [12] 

$6.308 

26,457 

Harmony search [13] 

$6.081 

27,721 

Cross entropy [15] 

$6.081 

97,000 

Scatter search [33] 

$6.081 

43,149 

Modified GA 1 [41] 

$6.056 

18,000 

Modified GA 2 [41] 

$6.190 

18,000 

Particle swarm harmony search [16] 

$6.081 

17,980 

Heuristic based approach [42] 

$6.701 

70 

Differential evolution [43] 

$6.081 

48,724 

Honey-bee mating optimization [44] 

$6.117 

15,955 

Heuristic based approach [45] 

$6.232 

259 

Self-adaptive differential evolution [46] 

$6.081 

60,582 

Pseudo-Genetic Algorithm [34] 

$6.081 

25,000 

Linear Programming – Differential Evolution [17] 

$6.081 

33,148 

B-GENOME (B-GA) [18]  

$6.182 

26,000 

SOGH [47] 

$6.337 

94 

OPUS (This research) 

$6.374 

106 

In the second scenario considering the 50” diameter, the total cost of the design obtained 

following the OPUS algorithm was of only $5,533,393, due to the closer resemblance between the final 
hydraulic gradient surface and the target one, as shown in Figure 9. In this approach, the unit cost of 
the 50” pipe was extrapolated, resulting in a value of $388.91/m. Finally, the resulting diameters listed 
by their ID shown in Figure 8a are 50, 50, 40, 40, 40, 40, 40, 24, 24, 40, 20, 20, 12, 12, 12, 16, 20, 20, 20, 
40, 16, 12, 40, 30, 24, 20, 12, 12, 12, 12, 12, 12, 12 and 20. These diameters are measured in inches. 

 

Figure 9. Comparison between the ideal HGL and that obtained by OPUS considering the additional 
50 in. diameter. 

Figure 9.

Comparison between the ideal HGL and that obtained by OPUS considering the additional

50 in. diameter.

/var/www/pavco.com.co/public/site/pdftohtml/5e81077962f627b67949947c14bfa4ef/index-html.html
background image

Water 2020, 12, 1037

16 of 26

4.2. Balerma

The Balerma network corresponds to an irrigation district in Almería, Spain [

32

], and is composed

of a total of 454 pipes and 443 demand nodes which are supplied by 4 reservoirs. The topology of the
network is presented in Figure

10

.

Water 202012, x FOR PEER REVIEW 

16 of 26 

 

4.2. Balerma 

The Balerma network corresponds to an irrigation district in Almería, Spain [32], and is 

composed of a total of 454 pipes and 443 demand nodes which are supplied by 4 reservoirs. The 
topology of the network is presented in Figure 10. 

 

(a) 

 

(b) 

Figure 10. (a) General layout for Balerma network and the spanning tree obtained for each reservoir, 
including the removed pipes. Dashed lines are pipes from the original networks that are not included 
in the Spanning Tree (b) OPUS design for Balerma Network. The diameters are shown in millimeters. 

Figure 10.

(a) General layout for Balerma network and the spanning tree obtained for each reservoir,

including the removed pipes. Dashed lines are pipes from the original networks that are not included
in the Spanning Tree (b) OPUS design for Balerma Network. The diameters are shown in millimeters.

The head-loss expression commonly used is the Darcy-Weisbach equation. The pipe diameters

commercially available for its design are manufactured exclusively in PVC, with an absolute roughness

/var/www/pavco.com.co/public/site/pdftohtml/5e81077962f627b67949947c14bfa4ef/index-html.html
background image

Water 2020, 12, 1037

17 of 26

coe

fficient of 0.0025 mm. The minimum pressure allowable is of 20 m, and it has 10 commercially

available pipes, the costs of which are presented in Table

3

. A sag value of 0.25 was used in the OPUS

methodology, while both the general layout and the obtained spanning tree are shown in Figure

10

a.

In this figure, it is shown that if the network is fed by multiple sources, it will have a maximum of one
spanning tree per reservoir. In this step, all nodes must be assigned to a spanning tree. However, a
reservoir can be isolated given that it is not e

fficient in terms of B/C to feed a part of the network.

Table 3.

Unit costs of Balerma network.

Diameter (mm)

113

126.6

144.6

162.8

180.8

226.2

285

361.8

452.2

581.8

Unit Cost (

€/m)

7.22

9.10

11.92

14.84

18.38

28.60

45.39

76.32

124.64 215.85

Considering continuous diameters for this network, OPUS algorithm generated a design with a

total cost of

€1.755 million. After executing the Round-off and Optimization processes, the optimal

discrete design reached a total cost of

€2.015 million, requiring 1165 hydraulic simulations. The

configuration of diameters in the OPUS design is shown in Figure

10

b. Table

4

shows reported costs

and number of iterations reached by several authors.

Table 4.

Reported costs and number of iterations for the Balerma WDS.

Algorithm

Cost (

€ millions)

Number of Iterations

Genetic algorithm [

32

]

2.302

10,000,000

Harmony search [

13

]

2.601

45,400

Harmony search [

13

]

2.018

10,000,000

Genetic algorithm [

12

]

3.738

45,400

Simulated annealing [

12

]

3.476

45,400

Simulated annealing with tabu search [

12

]

3.298

45,400

Local search with simulated annealing [

12

]

4.310

45,400

Hybrid discrete dynamically dimensioned search [

48

]

1.940

30,000,000

Harmony search with particle swarm [

16

]

2.633

45,400

SOGH [

47

]

2.100

1779

Memetic algorithm [

49

]

3.120

45,400

Genetic heritage evolution by stochastic transmission [

50

]

2.002

250,000

Di

fferential evolution [

46

]

1.998

2,400,000

Self-adaptive di

fferential evolution [

46

]

1.983

1,300,000

OPUS (this study)

2.015

1165

4.3. Taichung

Sung et al. [

28

] introduced the Taichung network, which is a WDS located in a city in Taiwan,

Taichung. This WDS has 31 pipes and 20 nodes, included in 12 di

fferent loops. Water is supplied by

one reservoir with 113.98 m as initial head. Table

5

presents the 13 commercial diameters and their

costs [

26

]. For the design process a minimum head of 15 m was imposed and the friction loss equation

is Hazen-Williams with 100 as the C coe

fficient [

28

]. The network’s topology is shown in Figure

11

a.

The OPUS methodology reached a cost of $8,914,400 after 47 iterations, considering discrete

diameters as shown in Figure

11

b. Reported costs in literature, and its corresponding number of

iterations for Taichung network are shown in Table

6

.

/var/www/pavco.com.co/public/site/pdftohtml/5e81077962f627b67949947c14bfa4ef/index-html.html
background image

Water 2020, 12, 1037

18 of 26

Table 5.

Unit costs of Taichung network.

Diameter (mm)

Unit Cost (NT Dollar

/m)

100

860

150

1160

200

1470

250

1700

300

2080

350

2640

400

3240

450

3810

500

4400

600

5580

700

8360

800

10,400

900

12,800

Water 202012, x FOR PEER REVIEW 

18 of 26 

 

Table 5. Unit costs of Taichung network. 

Diameter (mm) 

Unit Cost (NT Dollar/m) 

100 860 
150 1160 
200 1470 
250 1700 
300 2080 
350 2640 
400 3240 
450 3810 
500 4400 
600 5580 
700 8360 
800 10,400 
900 12,800 

 

 

 

(a) (b

Figure 11. (a) Spanning Tree for Taichung Network. Dashed lines are pipes from the original networks 
that are not included in the Spanning Tree. (b) OPUS design for Taichung Network. The diameters 
are shown in millimeters. 

The OPUS methodology reached a cost of $8,914,400 after 47 iterations, considering discrete 

diameters as shown in Figure 11b. Reported costs in literature, and its corresponding number of 
iterations for Taichung network are shown in Table 6. 

Figure 11.

(a) Spanning Tree for Taichung Network. Dashed lines are pipes from the original networks

that are not included in the Spanning Tree. (b) OPUS design for Taichung Network. The diameters are
shown in millimeters.

Table 6.

Reported costs and number of iterations for the Taichung WDS.

Algorithm

Cost (NT Dollar)

Number of iterations

Tabu search [

33

]

8,774,900

Not Available

OPUS (this study)

8,914,400

47

/var/www/pavco.com.co/public/site/pdftohtml/5e81077962f627b67949947c14bfa4ef/index-html.html
background image

Water 2020, 12, 1037

19 of 26

4.4. Pescara Network

The Pescara network corresponds to the reduced version of the water distribution system of a

medium-size city at Italy [

29

]. This WDS is composed of three reservoirs with a fixed head ranging

between 53.08 m to 57.00 m, 99 pipes and 68 demand nodes.

The available pipes are made of cast iron represented by a uniform Hazen-Williams roughness

coe

fficient of 130. The unitary costs of the available pipes are shown in Table

7

. A sag value of 0.10

was used in the OPUS methodology. However, for this network the general layout and the obtained

spanning tree for each reservoir are shown in Figure

4

a.

Table 7.

Unit costs of Pescara network.

Diameter (mm)

Unit Cost (

€/m)

100

27.7

125

38

150

40.5

200

55.4

250

75

300

92.4

350

123.1

400

141.9

450

169.3

500

191.5

600

246

700

319.6

800

391.1

The minimum pressure head of all demand nodes must be maintained at 20 m. However, this

system has two additional constraints in its design: A set of maximum pressures at each demand node
ranging between 28.5 m and 55.9 m, and a maximum velocity at each pipe of 2 m

/s.

The OPUS methodology reached an optimal discrete design with a total cost of

€2.161 million,

requiring 206 hydraulic simulations. The configuration of the diameters obtained by OPUS methodology
are shown in Figure

12

a. Table

8

presents minimum costs and the number of iterations reached by

several authors for Pescara Network.

Table 8.

Reported costs and number of iterations for the Pescara WDS.

Algorithm

Cost (

€ millions)

Number of Iterations

Mixed Integer Linear Programming—MILP [

29

]

1.820

Not available—Time limit: 7200 s

OPUS (This research)

2.161

206

The cost obtained using OPUS design is 18.7% higher than the reported by Bragalli et al. [

29

];

however, according to the authors, this result was the best solution they could reach within a time
limit of 7200 s. In terms of computational e

fficiency, OPUS methodology requires 210 iterations,

which represents less than one second in time, using an average desktop computer (we used a

fourth-generation Intel Core i7 processor with 16 GB of RAM, however design software requirements
are lower).

/var/www/pavco.com.co/public/site/pdftohtml/5e81077962f627b67949947c14bfa4ef/index-html.html
background image

Water 2020, 12, 1037

20 of 26

Water 202012, x FOR PEER REVIEW 

20 of 26 

 

 

(a

 

 

(b) (c

Figure 12. (a)  OPUS design for Pescara Network. The diameters are shown in millimeters. The 
resulting spanning tree is shown in Figure 4b. (b) Zoom in to the location of Pipe 84. (c) Zoom in to 
the location of Pipe 103. 

The cost obtained using OPUS design is 18.7% higher than the reported by Bragalli et al. [29]; 

however, according to the authors, this result was the best solution they could reach within a time 
limit of 7200 seconds. In terms of computational efficiency, OPUS methodology requires 210 
iterations, which represents less than one second in time, using an average desktop computer (we 
used a fourth-generation Intel Core i7 processor with 16 GB of RAM, however design software 
requirements are lower).  

Finally, OPUS design fulfilled pressure constraints of the network, considering both the 

minimum pressure service as well as the maximum pressures allowable at each node. However, these 
results showed that 2 pipes (ID 84 and ID 103, the latter located downstream the reservoir 43 as shown 
in Figure 12b,c correspondingly) in the system did not meet the maximum velocity constraint. 

Table 8. Reported costs and number of iterations for the Pescara WDS. 

Algorithm 

Cost (€ 

millions) 

Number of Iterations 

Mixed Integer Linear Programming—MILP 

[29] 

1.820 

Not available—Time limit: 7200 

seconds 

OPUS (This research) 

2.161 

206 

Figure 12.

(a) OPUS design for Pescara Network. The diameters are shown in millimeters. The

resulting spanning tree is shown in Figure

4

b. (b) Zoom in to the location of Pipe 84. (c) Zoom in to the

location of Pipe 103.

Finally, OPUS design fulfilled pressure constraints of the network, considering both the minimum

pressure service as well as the maximum pressures allowable at each node. However, these results

showed that 2 pipes (ID 84 and ID 103, the latter located downstream the reservoir 43 as shown in
Figure

12

b,c correspondingly) in the system did not meet the maximum velocity constraint.

4.5. Sensitivity Analysis OPUS Methodology—Sag

As described in the Optimal Power Use Surface (OPUS) step (Section

3.2

) of this methodology,

the optimal HGL is predefined using the sag, extending Wu (1975) criterion. In this context, the sag is a

positive value between 0.0 and 0.50 that represents the di

fference, in percentage, between a straight

HGL and the parabolic HGL that represents the minimum cost design.

To analyze the sensitivity of the resulting OPUS designs regarding the selected sag, several tests

were developed considering values between 0.0 and 0.5, with increments of 0.05. In order to assess the
performance of the methodology, two indexes were considered in addition to the cost (in millions) of

the resulting designs for measuring the reliability of the system: the Resilience Index (RI) proposed
by Todini [

51

], and the Network Resilience Index (NRI) proposed by Prasad and Park [

30

]. In first

place, the RI measures the capacity to be prepared against the occurrence of sudden failures in the

supply system based on the surplus of energy available at each junction; this is after demands and

pressures have been satisfied. Second, the NRI considers the combined e

ffect of both surplus power

/var/www/pavco.com.co/public/site/pdftohtml/5e81077962f627b67949947c14bfa4ef/index-html.html
background image

Water 2020, 12, 1037

21 of 26

and reliable loops in the network as a measurement of resilience. NRI addressed the reliable loops as a
consequence of having uniform diameters connected to a node.

In this analysis, Hanoi, Balerma, Taichung and Pescara networks were considered, leading to the

results shown in Figure

13

. As can be seen, in the case of Hanoi and Balerma, the least-cost design is

obtained by using a sag of 0.25. In Taichung, the least-cost design is reached by using a sag value of 0.1
to 0.2, having an increase in the cost of the network of 0.88% when the sag value is 0.15. In Pescara,
the near-optimal design is reached by using a sag value between 0.1 and 0.25. However, when a
sag of 0.20 is used, an increase of 2.05% to the cost is obtained. Hence, it can be concluded that the
least-cost designs are obtained when between 0.10 and 0.25 are used. Accordingly, the behavior shown
by Pescara and Taichung may be a result of the rounding process, where small changes in the optimal

power surface lead to some increase in pipe diameters.

Water 202012, x FOR PEER REVIEW 

21 of 26 

 

4.5. Sensitivity Analysis OPUS Methodology

Sag  

As described in the Optimal Power Use Surface (OPUS) step (Section 3.2) of this methodology, 

the optimal HGL is predefined using the sag, extending Wu (1975) criterion. In this context, the sag 
is a positive value between 0.0 and 0.50 that represents the difference, in percentage, between a 
straight HGL and the parabolic HGL that represents the minimum cost design.  

To analyze the sensitivity of the resulting OPUS designs regarding the selected sag, several tests 

were developed considering values between 0.0 and 0.5, with increments of 0.05. In order to assess 
the performance of the methodology, two indexes were considered in addition to the cost (in millions) 
of the resulting designs for measuring the reliability of the system: the Resilience Index (RI) proposed 
by Todini [51], and the Network Resilience Index (NRI) proposed by Prasad and Park [30]. In first 
place, the RI measures the capacity to be prepared against the occurrence of sudden failures in the 
supply system based on the surplus of energy available at each junction; this is after demands and 
pressures have been satisfied. Second, the NRI considers the combined effect of both surplus power 
and reliable loops in the network as a measurement of resilience. NRI addressed the reliable loops as 
a consequence of having uniform diameters connected to a node. 

In this analysis, Hanoi, Balerma, Taichung and Pescara networks were considered, leading to 

the results shown in Figure 13. As can be seen, in the case of Hanoi and Balerma, the least-cost design 
is obtained by using a sag of 0.25. In Taichung, the least-cost design is reached by using a sag value 
of 0.1 to 0.2, having an increase in the cost of the network of 0.88% when the sag value is 0.15. In 
Pescara, the near-optimal design is reached by using a sag value between 0.1 and 0.25. However, 
when a sag of 0.20 is used, an increase of 2.05% to the cost is obtained. Hence, it can be concluded 
that the least-cost designs are obtained when between 0.10 and 0.25 are used. Accordingly, the 
behavior shown by Pescara and Taichung may be a result of the rounding process, where small 
changes in the optimal power surface lead to some increase in pipe diameters. 

 

Figure 13. Sensitivity Results for Cost and Network Resilience Index (NRI) for Different Values of Sag 
(a) Hanoi Network. (b) Balerma Network. (c) Taichung Network (d) Pescara Network. 

In addition, for sag values higher than 0.25, the total costs of the network will be greater. The 

reason for this was first presented by Wu (1975), as for bigger sag values there will be a lower energy 
slope upstream (negative sag) or downstream (positive sag). Thus, cost will increase because of the 

Figure 13.

Sensitivity Results for Cost and Network Resilience Index (NRI) for Di

fferent Values of Sag

(a) Hanoi Network. (b) Balerma Network. (c) Taichung Network (d) Pescara Network.

In addition, for sag values higher than 0.25, the total costs of the network will be greater. The

reason for this was first presented by Wu (1975), as for bigger sag values there will be a lower energy
slope upstream (negative sag) or downstream (positive sag). Thus, cost will increase because of the
need for use of larger sizes of pipe in those sections. Therefore, the costs curve will be convex, with the
lower values being between the aforementioned range for sag values (0.10 and 0.25 depending on the
network), and increasing above this threshold.

Regarding the comparison between the Resilience Index and the Network Resilience Index, it can

be seen that they both decrease as the sag values increase. Given this, it can be inferred that for sag

values above 0.25, the obtained results will be more expensive and less resilient. In addition, based on

the topology of the networks, it can be seen that if the system tends to be more branched, the NRI

will approach the RI. On the other hand, if the network tends to be more looped, the NRI will di

ffer

considerably from the RI.

The results presented for the Hanoi, Balerma and Pescara networks will be compared with

some results available in literature to validate the performance of OPUS methodology in terms of the
resilience of the network. As is presented in Figure

13

, the Hanoi network has values ranging from 0.19

/var/www/pavco.com.co/public/site/pdftohtml/5e81077962f627b67949947c14bfa4ef/index-html.html
background image

Water 2020, 12, 1037

22 of 26

to 0.27 for the RI and 0.18–0.26 for the NRI. Similarly, the Balerma network has a RI between 0.30–0.44,
and an NRI between 0.28 and 0.42. Finally, Pescara network has a RI ranging between 0.39 to 0.52,
and an NRI between 0.32–0.43.

Several optimization approaches have been used for the design of water distribution networks.

Reca, Martínez, Baños and Gil [

52

] used genetic algorithms to minimize the network investment cost,

determining an optimal cost of approximately $6.1 million for the Hanoi network with an RI of 0.15.
Reca et al. [

52

] also developed five di

fferent multi-objective optimization models for the design of

water distribution networks, which are based on minimizing the investment cost and maximizing the

RI. These methods proved to reach a lower cost for the same resilience index in comparison with other
approaches. The resilience index interval found for the Hanoi network varies from 0.15 to 0.35 with
cost ranging from around $5.6 million to $8 million [

52

]. On the other hand, Monsef (2019) reported an

RI between 0.15 and 0.25, and an NRI ranging from 0.174 to 0.331 for the Hanoi network [

53

].

Regarding the Balerma network, the RI varies between 0.45–0.85 with cost ranging from around

€2.9 million to €9.5 million [

52

]. The behavior of the pareto front of the Balerma network found by

Reca et al. is similar to the results presented by other authors [

54

,

55

].

Furthermore, Wang, Savi´c and Kapelan [

54

] used di

fferent multi-objective metaheuristics to

minimize the total cost and maximize the NRI of Pescara network resulting in an NRI varying from 0.4
to 0.99 with a cost ranging from 2 million to 18 million. On the other hand, multi-objective optimization
based on decomposition presented values for NRI that ranged from 0.32 to 0.96, corresponding to costs
of between 2 and 3.9 million [

55

].

Hence, through these comparisons, it can be concluded that OPUS is able to reach near-optimal

designs that are comparable to those obtained by other authors in terms of costs, but also in terms of
reliability and resilience. Finally, for reaching least cost designs, it is recommended to use sag values
between 0.10 and 0.25, depending the complexity of the network, considering it in terms of its size,
and also on whether the system tends to be more branched or more looped. As shown in this sensitivity
analysis, using a sag value in this range ensures achieving near the minimum cost for the network.

5. Conclusions

The presented methodology was tested on four benchmark systems, previously mentioned in

current literature: Hanoi, Balerma, Taichung and Pescara. In all cases, the OPUS algorithm reached
near optimal results with the number of iterations up to five orders of magnitude smaller than other
methodologies [

26

]. These results are noteworthy when considering that the methodology introduced

here is entirely deterministic and thus, any user will obtain the same design, as opposed to metaheuristic
methodologies that consider a stochastic component. OPUS algorithm takes into account an optimal
distribution pattern of flow in the network as a way of spending energy properly. This approach
di

fferentiates it from metaheuristic methodologies that do not take into account hydraulic principles to

search the solution space [

28

], proving that it is feasible to develop an entirely hydraulic method that

provides near-optimal results.

In the case of the Hanoi network, OPUS achieved a di

fference of only of 4.81% higher than the

lowest reported cost, requiring a number of iterations three orders of magnitude smaller. Regarding
this network, the benchmark value in the literature was $6.081, given that lower solutions did not meet

pressure constraints using EPANET2 [

27

]. Moreover, for establishing a comparison, the near-optimal

solution reached by OPUS is located within a frequency of 2% of the most frequent solutions shown
by Mora-Meliá et al. in the cited study. Hence, it can be seen that several authors report their best
results in the literature, regardless of whether or not it is a result that can be obtained consistently. As
for the Balerma network, it presented a di

fference of 3.86% in comparison to the lowest cost reported

in previous works, but executed a number of iterations five orders of magnitude below. The Pescara
network presented a di

fference of 18.73% with respect to the results reported previously, considering

that this design presented a violation of velocity constraint in 2 pipes (ID 84 and ID 103). Finally,
the Taichung network presented a di

fference of 1.59% when compared to the result reported in literature.

/var/www/pavco.com.co/public/site/pdftohtml/5e81077962f627b67949947c14bfa4ef/index-html.html
background image

Water 2020, 12, 1037

23 of 26

Even though these results did not report the number of iterations required to reach the solution, with
OPUS, we only required a total for 47 iterations to reach this near-optimal result, which does not
demand a significant computational e

ffort.

Regarding the sag value, it was demonstrated that least cost designs are obtained when this

value is between 0.10 and 0.25. This can be compared to the Wu criterion [

16

], which uses a sag

of 0.15 for series’ of irrigation pipes. A comparison using the Network Resilience Index (NRI) and
the Resilience Index (RI) was done, resulting in a decreasing behavior for higher values of the sag.
Similarly, the NRI and RI results were compared with values available in literature, finding that the
results of the OPUS methodology were comparable to them. Hence, it can be concluded that for higher
sag values, the resulting designs will be more expensive and less reliable. Meanwhile, both resilience
indexes showed a similar behavior among them.

The OPUS methodology showed a good performance on real-size networks, leading to the

minimization of constructive costs in a substantially low number of iterations in comparison to other
methodologies. Conversely, because of the di

fferences in hydraulic behavior between a tree-like

network and a real-life looped network, when OPUS is used in small WDSs, the results are not as good
as expected because that di

fference in hydraulics is higher in smaller networks. In addition, Giudicianni

(2018) found that there are several topological di

fferences between small, large and synthetic WDNs

(Water Distribution Networks) [

51

,

56

], thus it is possible to state that one methodology fits for some

kind of networks but not for others. In this case, the result was highly influenced by the optimization

process and not by the steps based on hydraulic principles, as a result of the fact that the design is

executed based on the branched network (spanning tree).

The algorithm described herein proves to be e

fficient and to achieve near-optimal solutions that

resembles the costs and resilience of other evolutionary algorithm methodologies, but with running a
lower number of iterations, based on understanding and applying the hydraulic principles that govern

WDSs. Some of the results presented correspond to theoretical networks with some constraints that

limited the design possibilities, as in the case of Hanoi. Therefore, it is highly recommended to test this
methodology on real WDSs and include additional objectives in the procedure, such as maximizing
the network reliability or minimizing leakages. Furthermore, energy approaches can be used for the
calibration of models and design optimal system operation. The results obtained through the OPUS
algorithm makes it a promising methodology to be used as a “hot start” for di

fferent metaheuristics.

Namely, the pipe configuration obtained as a result of OPUS could be the initial configuration used
by a stochastic methodology as a hot start to obtain even better results still with a small number
of iterations.

Finally, despite the fact that the OPUS methodology was originally made to be applied on

gravity-driven WDSs only, the algorithm could also be easily applied to the case of pump-driven
networks. OPUS has not been applied to these cases, but one can speculate that the inclusion of pumps
could be carried out with some minor adjustments to the algorithm. In those networks in which the

water is pumped to an initial reservoir, the pumping e

ffect could be easily introduced increasing the

HGL at the entrance of the distribution system. A di

fferent case is when the pump is located in some

node of the WDS downstream of the initial node. In this case, all nodes downstream of the pump must
be changed to take into account the e

ffect of a sudden increase in pressure. A decrease in the HGL

of those downstream nodes should be included so that a simulation of a topological descent of the
network can mimic the action of gravity. This means that the power surface must be changed because
of the pressure addition at specific points of the network [

26

].

Supplementary Materials:

The following are available online at

http:

//www.mdpi.com/2073-4441/12/4/1037/s1

.

Optimal Design for Hanoi using OPUS (Hanoi.inp), Optimal Design for Balerma using OPUS (Balerma.inp),
Optimal Design for Pescara using OPUS (Pescara.inp), Optimal Design for Taichung using OPUS (Taichung.inp).
REDES 2018 software.

/var/www/pavco.com.co/public/site/pdftohtml/5e81077962f627b67949947c14bfa4ef/index-html.html
background image

Water 2020, 12, 1037

24 of 26

Author Contributions:

J.S. conceived the algorithm and both guided and supervised the research process; J.S.,

C.S., D.P., D.C., P.C. and L.L.L. wrote the paper; L.L.L. developed implemented the algorithm in REDES software;
C.S., L.L.L., P.C. and N.L. performed the simulations and analyzed the results; D.P. explained the proposed
Optimal HGL Criteria; C.S. performed the sensitivity analysis for the optimal sag with OPUS. All authors have
read and agree to the published version of the manuscript.

Funding:

This research received no external funding.

Conflicts of Interest:

The authors declare no conflict of interest.

References

1.

Ilaya-Ayza, A.E.; Martins, C.; Campbell, E.; Izquierdo, J. Implementation of DMAs in Intermittentwater
supply networks based on equity criteria. Water 2017, 9, 851. [

CrossRef

]

2.

Giudicianni, C.; Herrera, M.; di Nardo, A.; Carravetta, A.; Ramos, H.M.; Adeyeye, K. Zero-net energy
management for the monitoring and control of dynamically-partitioned smart water systems. J. Clean. Prod.
2020

, 252, 119745. [

CrossRef

]

3.

Central Public Health and Environmental Engineering Organisation; Ministry of Urban Development; World
Health Organisation. Manual on Operation and Maintenance of Water Supply Systems; International Water

Association: New Delhi, India, 2005.

4.

Yates, D.F.; Templeman, A.B.; Bo

ffey, T.B. The computational complexity of the problem of determining least

capital cost designs for water supply networks. Eng. Optim. 1984, 7, 143–155. [

CrossRef

]

5.

Wu, Z.Y.; Walski, T. Self-adaptive penalty approach compared with other constraint-handling techniques for

pipeline optimization. J. Water Resour. Plan. Manag. 2005, 131, 181–192. [

CrossRef

]

6.

Montalvo, I.; Izquierdó, J.; Peerez-Garcia, R.; Herrera, M. Improved performance of PSO with self-adaptive

parameters for computing the optimal design of Water Supply Systems. Eng. Appl. Artif. Intell. 2010, 23,

727–735. [

CrossRef

]

7.

Gessler, J. Pipe network optimization by enumeration. In Proceedings of the Specialty Conference on
Computer Applications in Water, Bu

ffalo, NY, USA, 10–12 June 1985.

8.

Alperovits, E.; Shamir, U. Design of optimal water distribution systems. Water Resour. Res. 1977, 13, 885–900.
[

CrossRef

]

9.

Storn, R.; Price, K. Di

fferential Evolution: A Simple and Efficient Adaptive Scheme for Global Optimization

Over Continuous Spaces. J. Glob. Optim. 1997, 11, 341–359. [

CrossRef

]

10.

Savic, D.A.; Walters, G.A. Genetic Algorithms for Least-Cost Design of Water Distribution Networks. J. Water

Resour. Plan. Manag. 1997, 123, 67–77. [

CrossRef

]

11.

Cunha, M.C.; Sousa, J. Water Distribution Network Design Optimization: Simulated Annealing Approach.
J. Water Resour. Plan. Manag. 1999, 125, 215–221. [

CrossRef

]

12.

Reca, J.; Martínez, J.; Gil, C.; Baños, R. Application of several meta-heuristic techniques to the optimization
of real looped water distribution networks. Water Resour. Manag. 2007, 22, 1367–1379. [

CrossRef

]

13.

Geem, Z.W. Optimal cost design of water distribution networks using harmony search. Eng. Optim. 2006, 38,
259–277. [

CrossRef

]

14.

Lin, M.-D.; Liu, Y.H.; Liu, G.F.; Chu, C.W. Scatter search heuristic for least-cost design of water distribution
networks. Eng. Optim. 2007, 39, 857–876. [

CrossRef

]

15.

Perelman, L.; Ostfeld, A. An adaptive heuristic cross-entropy algorithm for optimal design of water
distribution systems. Eng. Optim. 2007, 39, 413–428. [

CrossRef

]

16.

Geem, Z.W. Particle-swarm harmony search for water network design. Eng. Optim. 2009, 41, 297–311.
[

CrossRef

]

17.

Zheng, F.; Simpson, A.R.; Zecchin, A.C. Coupled binary linear programming-di

fferential evolution algorithm

approach for water distribution system optimization. J. Water Resour. Plan. Manag. 2014, 140, 585–597.
[

CrossRef

]

18.

Reca, J.; Martínez, J.; López, R. A Hybrid Water Distribution Networks Design Optimization Method Based
on a Search Space Reduction Approach and a Genetic Algorithm. Water 2017, 9, 845. [

CrossRef

]

19.

Zheng, F.; Zecchin, A.C.; Simpson, A.R. Investigating the run-time searching behavior of the di

fferential

evolution algorithm applied to water distribution system optimization. Environ. Model. Softw. 2015, 69,
292–307. [

CrossRef

]

/var/www/pavco.com.co/public/site/pdftohtml/5e81077962f627b67949947c14bfa4ef/index-html.html
background image

Water 2020, 12, 1037

25 of 26

20.

Rossman, L.A. Epanet 2 User’s Manual; Environmental Protection Agency: Boston, MA, USA, 2000.

21.

Wu, I. Design of drip irrigation main lines. J. Irrig. Drain. Div. 1975, 101, 265–278.

22.

Featherstone, R.E.; El-Jumaily, K.K. Optimal Diameter Selection for Pipe Networks. J. Hydraul. Eng. 1983,

109, 221–234. [

CrossRef

]

23.

Takahashi, S.; Saldarriaga, J.; Hernández, F.; Díaz, D.; Ochoa, S. An energy methodology for the design of

water distribution systems. In Proceedings of the World Environmental and Water Resources Congress
(EWRI) 2010, Providence, RI, USA, 16–20 May 2010.

24.

Saldarriaga, J.; Takahashi, S.; Hernández, F.; Escovar, M. Predetermining pressure surfaces in water
distribution system design. In Proceedings of the World Environmental and Water Resources Congress 2011:
Bearing Knowledge for Sustainability, Palm Springs, CA, USA, 22–26 May 2011.

25.

Páez, D.; Saldarriaga, J.; López, L.; Salcedo, C. Optimal design of water distribution systems with pressure
driven demands. Procedia Eng. 2014, 89, 839–847. [

CrossRef

]

26.

Saldarriaga, J.; Páez, D.; León, N.; López, L.; Cuero, P. Power use methods for optimal design of WDS:
History and their use as post-optimization warm starts. J. Hydroinf. 2015, 17, 404–421. [

CrossRef

]

27.

Prim, R.C. Shortest Connection Networks and Some Generalizations. Bell Syst. Tech. J. 1957, 36, 1389–1401.
[

CrossRef

]

28.

Saldarriaga, J.; Páez, D.; Cuero, P.; Leön, N. Optimal design of water distribution networks using mock open
tree topology. In Proceedings of the World Environmental and Water Resources Congress 2013: Showcasing
the Future, Cincinnati, OH, USA, 19–23 May 2013.

29.

Bragalli, C.; D’Ambrosio, C.; Lee, J.; Lodi, A.; Toth, P. On the optimal design of water distribution networks:

A practical MINLP approach. Optim. Eng. 2012, 13, 219–246. [

CrossRef

]

30.

Prasad, T.D.; Park, N.S. Multiobjective genetic algorithms for design of water distribution networks. J. Water

Resour. Plan. Manag. 2004, 130, 73–82. [

CrossRef

]

31.

Fujiwara, O.; Khang, D.B. A two-phase decomposition method for optimal design of looped water distribution
networks. Water Resour. Res. 1990, 26, 539–549. [

CrossRef

]

32.

Reca, J.; Martínez, J. Genetic algorithms for the design of looped irrigation water distribution networks.

Water Resour. Res. 2006, 42, 42. [

CrossRef

]

33.

Sung, Y.-H.; Lin, M.-D.; Lin, Y.-H.; Liu, Y.-L. Tabu search solution of water distribution network optimization.
J. Environ. Eng. Manag. 2007, 17, 177.

34.

Mora-Melia, D.; Iglesias-Rey, P.L.; Martinez-Solano, F.J.; Fuertes-Miquel, V.S. Design of Water Distribution
Networks using a Pseudo-Genetic Algorithm and Sensitivity of Genetic Operators. Water Resour. Manag.
2013

, 27, 4149–4162. [

CrossRef

]

35.

Geem, Z.W.; Kim, J.H.; Loganathan, G.V. Harmony search optimization: Application to pipe network design.
Int. J. Model. Simul. 2002, 22, 125–133. [

CrossRef

]

36.

Eusu

ff, M.M.; Lansey, K.E. Optimization of water distribution network design using the shuffled frog leaping

algorithm. J. Water Resour. Plan. Manag. 2003, 129, 210–225. [

CrossRef

]

37.

Liong, S.-Y.; Atiquzzaman, M. Optimal design of water distribution network using shu

ffled complex

evolution. J. Inst. Eng. 2004, 44, 93–107.

38.

Vairavamoorthy, K.; Ali, M. Pipe index vector: A method to improve genetic-algorithm-based pipe
optimization. J. Hydraul. Eng. 2005, 131, 1117–1125. [

CrossRef

]

39.

Zecchin, A.C.; Simpson, A.R.; Maier, H.R.; Leonard, M.; Roberts, A.J.; Berrisford, M.J. Application of two ant
colony optimisation algorithms to water distribution system optimisation. Math. Comput. Model. 2006, 44,
451–468. [

CrossRef

]

40.

Suribabu, C.R.; Neelakantan, T.R. Design of water distribution networks using particle swarm optimization.

Urban Water J. 2006, 3, 111–120. [

CrossRef

]

41.

Kadu, M.S.; Gupta, R.; Bhave, P.R. Optimal design of water networks using a modified genetic algorithm

with reduction in search space. J. Water Resour. Plan. Manag. 2008, 134, 147–160. [

CrossRef

]

42.

Mohan, S.; Babu, K.S.J. Water distribution network design using heuristics-based algorithm. J. Comput.

Civ. Eng. 2009, 23, 249–257. [

CrossRef

]

43.

Suribabu, C.R. Di

fferential evolution algorithm for optimal design of water distribution networks. J. Hydroinf.

2010

, 12, 66–82. [

CrossRef

]

44.

Mohan, S.; Babu, K.S.J. Optimal water distribution network design with honey-bee mating optimization.
J. Comput. Civ. Eng. 2010, 24, 117–126. [

CrossRef

]

/var/www/pavco.com.co/public/site/pdftohtml/5e81077962f627b67949947c14bfa4ef/index-html.html
background image

Water 2020, 12, 1037

26 of 26

45.

Suribabu, C.R. Heuristic-based pipe dimensioning model for water distribution networks. J. Pipeline Syst.

Eng. Pract. 2012, 3, 115–124. [

CrossRef

]

46.

Zheng, F.; Zecchin, A.C.; Simpson, A.R. Self-adaptive di

fferential evolution algorithm applied to water

distribution system optimization. J. Comput. Civ. Eng. 2013, 27, 148–158. [

CrossRef

]

47.

Ochoa, S. Optimal Design of Water Distribution Systems Based on the Optimal Hydraulic Gradient Surface
Concept. Master’s degree Thesis, Universidad de los Andes, Bogotá, Colombia, 2009. Unpublished

(In Spanish).

48.

Tolson, B.A.; Asadzadeh, M.; Maier, H.R.; Zecchin, A. Hybrid discrete dynamically dimensioned search

(HD-DDS) algorithm for water distribution system design optimization. Water Resour. Res. 2009. [

CrossRef

]

49.

Baños, R.; Gil, C.; Reca, J.; Montoya, F.G. A memetic algorithm applied to the design of water distribution
networks. Appl. Soft Comput. J. 2010, 10, 261–266. [

CrossRef

]

50.

Bolognesi, A.; Bragalli, C.; Marchi, A.; Artina, S. Genetic Heritage Evolution by Stochastic Transmission in
the optimal design of water distribution networks. Adv. Eng. Softw. 2010, 41, 792–801. [

CrossRef

]

51.

Todini, E. Looped water distribution networks design using a resilience index based heuristic approach.

Urban Water 2000, 2, 115–122. [

CrossRef

]

52.

Reca, J.; Martinez, J.; Baños, R.; Gil, C. Optimal Design of Gravity-Fed Looped Water Distribution Networks
Considering the Resilience Index. J. Water Resour. Plan. Manag. 2008, 134, 234–238. [

CrossRef

]

53.

Monsef, H.; Naghashzadegan, M.; Farmani, R.; Jamali, A. Deficiency of Reliability Indicators in Water
Distribution Networks. J. Water Resour. Plan. Manag. 2019, 145, 04019022. [

CrossRef

]

54.

Wang, Q.; Savi´c, D.A.; Kapelan, Z. Hybrid metaheuristics for multi-objective design of water distribution
systems. J. Hydroinf. 2013, 16, 165–177. [

CrossRef

]

55.

Yazdi, J. Decomposition based Multi Objective Evolutionary Algorithms for Design of Large-Scale Water
Distribution Networks. Water Resour. Manag. 2016, 30, 2749–2766. [

CrossRef

]

56.

Giudicianni, C.; Di Nardo, A.; Di Natale, M.; Greco, R.; Santonastaso, G.F.; Scala, A. Topological taxonomy of

water distribution networks. Water 2018, 10, 444. [

CrossRef

]

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access

article distributed under the terms and conditions of the Creative Commons Attribution

(CC BY) license (http:

//creativecommons.org/licenses/by/4.0/).

¿Quiere saber más? Contáctenos

Declaro haber leído y aceptado la Política de Privacidad
contactenos