(I)

CrossMark

Available online at www.sciencedirect.com

ScienceDirect

Procedía Engineering 154 (2016) 19 - 26

Procedía Engineering

www.elsevier.com/locate/procedia

12th International Conference on Hydroinformatics, HIC 2016

A graph theoretical sectorization approach for energy reduction in

water distribution networks

M.E. Castro-Gamaa*, Q. Pana, A. Jonoskia, D. Solomatinea.

aUNESCO-IHE, Institute for Water Education, Westvest 7, 2611AX, Delft, The Netherlands.

Abstract

Most water supply utilities in the world need to deal with daily operations of their systems while at the same time they are required to increase their efficiency. One way to increase the efficiency of the system is to divide the water distribution network into smaller subsystems called sectors. This article presents a new sectorization approach, which has been developed as part of the European FP7 research project ICeWater, focused on increasing energy efficiency. The approach is based on the combination of graph theory and multi-objective optimization applied to the system of Milan (Italy). Most existing approaches for sectorization deal with gravity sources. The methodology presented here has been applied to a system with 26 pump stations, which is an innovation and a step forward for other utilities facing similar challenges. The application of sectorization to the network of Milan has the advantage of allowing to fix the number of sectors and enabling the exploration of multiple scenarios of system operation.

© 2016Publishedby ElsevierLtd. This isanopenaccess article under the CC BY-NC-ND license

(http://creativecommons.Org/licenses/by-nc-nd/4.0/).

Peer-review under responsibility of the organizing committee of HIC 2016

Keywords: Sectorization, graphs, water distribution networks, Milan

1. Introduction

In the past, many researchers have proposed solutions to Single Objective Optimization (SOO) and Multi-Objective Optimization problems (MOO) formulated for reducing energy consumption in Water Distribution Networks (WDN), such as, for example, pump scheduling problems. However, in most cases the number of decision

* Corresponding author. Tel.: +31-6-8281-4926. E-mail address: m.castrogama@unesco-ihe.org

1877-7058 © 2016 Published by Elsevier Ltd. This is an open access article under the CC BY-NC-ND license

(http://creativecommons.Org/licenses/by-nc-nd/4.0/).

Peer-review under responsibility of the organizing committee of HIC 2016

doi: 10.1016/j.proeng.2016.07.414

variables used for optimization is low as compared to the number of decision variables required in real systems [1, 2].

Nomenclature

A s^pn incidency matrix of the Global Gradient Algorithm

D degree matrix of a graph

DMA District Metering Areas

Ex a set of graph edges

Gx a graph containing a set of vertex (Vx) and a set of edges (Ex)

k number of sectors

L Lapacian matrix of a graph

MOO Multi-Objective Optimization problems

n number of subsectors

nc number of configurations of sectorization to be generated for a particular number of sectors k

neig number of non-zero eigenvalues of matrix L

nn is the number of nodes, tanks and reservoirs in a network

np number of pipes, valves and pumps in a network

r number of remaining subsectors for allocation

Sk A sectorization configuration with k sectors.

Si a subsector of a list of n available

SOO Single Objective Optimization

WNP Water Network Partitioning

WNS Water Network Sectorization

WNSeg Water Network Segmentation

WDN Water Distribution Networks

Vx a set of graph vertices

It is necessary to find new ways to approach the problem of energy minimization for real networks of large size. Such approaches include: Water Network Partitioning (WNP) which corresponds to dividing a network into a number of district metering areas (DMA) [3, 4]; Water Network Segmentation (WNSeg) [5, 6], in which the goal is to identify locations for pressure regulating valves; Water Network Sectorization (WNS) when a part of a WDN, named 'sector' is supplied only by a single source and the DMAs are completely isolated [7].

All of the approaches mentioned above are intended for application with gravity sources, and in small to medium size WDNs. For that reason these approaches are not suitable for most large WDNs.

In the case of Milan, the system has a total of 26 pumping stations. The system contains 103 booster pumps, 41 tanks, 143,684 junctions, 113,784 pipes and 35,338 valves. Currently the system is operated as a single Pressure Management Zone (one large sector), which leads to excessive energy use, bringing costs to the utility of nearly €16m per year [8]. Recently, several different methodologies have been applied for pump scheduling in a subsector of Milan known as Abbiategrasso [9] and a WNS approach was proposed for the northern area of Milan [10]. In another approach, a WNS of the full system of Milan was proposed [11], making the boundaries between isolated subsectors part of the optimization process, but lacking the possibility to fix the number of sectors required.

One of the main problems with sectorization approaches is that the number of sectors can not be fixed as desired. The only attainable way to overcome this is to pose a constraint to the number of sectors or partitions and include this in the MOO. This is a known issue even for most complicated approaches using modularity analysis for WNSeg of WDN [6, 7].

The approach presented here allows us to deal with two main issues of WNS: 1) the sources of each sector created can be pump stations and not only gravity sources, 2) the number of possible sectors can be fixed for a particular problem setup.

The methodology is applied to the case of Milan WDN. The new sectorization approach enables splitting of the network in such a way that a large WDN can be converted in a set of smaller isolated networks of a manageable size.

This article presents only the structure of the sectorization algorithm, and does not deal with the pump scheduling problem which is inherent part of the whole approach.

The outline of the paper is as follows: in section 2 specific definitions of graph theory and its connection to WDN is presented, in section 3 a discussion on advantages and limitations of the proposed sectorization algorithm are presented. Finally, we present our conclusions and future work.

2. Graph theory and WDN

The main algorithm for hydraulic simulations of WDN is known as the Global Gradient Algorithm (GGA), as it is implemented in EPANET2.0 [12]. Inside this algorithm a definition for a topological incidence matrix is elaborated. This matrix, known as incidence matrix Apn has dimensions of [np x nn], where np is the number of pipes (including pipes, valves and pumps) and nn is the number of nodes in the network (including demand nodes and water sources such as tanks and reservoirs). Each row contains only two non-zero values which correspond to the initial and ending node of the pipe. The values are {1,-1} depending on the positive direction of the assumed flow in the pipe.

By analogy, this matrix can be interpreted as the graph incidence matrix where nodes are graph vertices (Vj) and pipes are graph edges (Ej). By definition a graph Gj contains a set of GJ(VJ,EJ) elements, so it is possible to analyze WDN using same principles of graph theory [13, 14]. A WDN can always be expressed as a planar graph, because all of its elements are located in 2D and its connectivity and layout allows for its study differently than complex networks [13].

2.1. Algorithm for sectorization

The algorithm for sectorization of a WDN is composed of two different steps. This first step is initialization and the second one is sectorization. Initialization deals with the identification of elements of the WDN which may isolate smaller subsystems supplied by a source(s). As a result of initialization, smaller subsectors are identified. The number of subsectors is established as the maximum number of sectors for a given WDN. If all subsectors are connected, then a WDN operates as a single sector, while if all subsectors are isolated the maximum isolation possible has been reached. As mentioned earlier, currently the WDN of Milan operates as a single connected system.

It is also possible to create larger sectors from the subsectors depending on their neighbors. This allows for the possibility to create a diverse number of sectors of different size (that we may call 'sectorization configurations'), where each sector can be supplied with single or multiple sources. This process of splitting a graph is known as a k-way graph partitioning [15, 16].

2.1.1. Initialization

Previously, in [11], the authors presented an algorithm for initialization based on graph theory using Shortest Dissipated Power Path (SDPP) [3]. It corresponds to a shortest path analysis of a WDN, in which the distance measured along each vertex is the head loss from a hydraulic simulation. However, this algorithm was presented in [3] only for use when the sources of each subsector are gravity sources.

In the case of gravity sources the energy supplied, or the maximum distance reached from any source, is the result of the head losses of pipes connected to a certain source. In the case of pump stations as sources, the limitation is not only the head loss of the pipes added to each source, but also the maximum available discharge which may be supplied by the pumps. This is overcome by taking into account the maximum extent or SDPP to which a pump station can provide flow to satisfy the demands.

By applying this methodology, the system of Milan can be split into a total of 18 subsectors (sj to sj8) identified in Table 1, and a total of 38 cut-sets (explained below and presented in Fig. 1), composed of a total of 346 isolation valves. It is also possible to identify that each subsector contains a different number of pump stations (sources).

Previously the initialization algorithm presented in [11] identified the boundaries between subsectors, named cutsets (sets of isolation valves). However, in that approach, once the isolation valves were identified, these cut-sets

were included as decision variables in a MOO problem. Given that the opening or closing of cut-sets was driven only by the MOO, it was not possible to fix the number of sectors during optimization.

2.1.2. Sectorization

To overcome the issue of variable number of sectors it was required to think of the WDN with an additional abstraction. After performing initialization, a graph of the network is used G2(V2,E2), but in this case the set V2 corresponds to the subsectors and the set E2 corresponds to the cut-sets obtained during initialization.

Table 1, Subsectors identified in WDN Milan by initialization.

Subsector si S2 S3 S4 S5 S6

Pump Stations Salemi Vialba Suzzani Gorla Chusabella Cimabue

n-pumps 4 4 5 3 3 4

Subsector S7 S8 S9 Sl3 Sl4 Sl5

Pump Stations Comasina Novara San Siro Cantore Este Anfossi

n-pumps 4 4 4 3 2 2

Subsector sio Sil S12 Si6 Sl7 Sl8

Pump Stations Armi Feltre Tonnezza Ovidio Abbiategrasso Martini

Italia Crescenzago Baggio Linate

Parco Lambro Assiano

Padova

n-pumps 10 17 12 7 4 3

Fig. 1. (a) Real network showing pressure variation in the system, (b) Subsectors abstraction.

Each cut-set becomes a boundary between subsectors which is known in advance, making it possible to associate subsectors to create larger sectors. For the system of Milan the adjacency matrix of G2 is known and with it is possible to identify the neighbors of the each subsector. The values of the adjacency matrix in each row are equal to 1 if the subsector x is adjacent to subsector y and 0 otherwise.

A comparison between the real system (Fig. 1a) and the graph abstraction (Fig. 1b) is presented in Fig. 1. Lines in Fig. 1b represent cut-sets (E2) between subsectors. As all subsectors in Fig. 1b are connected it represents the current system of Milan with all the valves open among subsectors.

Once the subsectors are identified, it is possible to create different sectorization configurations with different preset, desired number of sectors - k. Here k ranges between 2 and 18 and it is guaranteed that k sectors are created or a k-way graph partitioning of the system is performed. The flowchart of the algorithm is presented in Fig. 2.

Fig. 2. Flowchart of the sectorization algorithm.

One particular sectorization configuration will be named Sk, which is a list of the subsectors sn which belongs to each sector k. A random combination of the n subsectors (18 for Milan) is performed to choose k subsectors as the nuclei of sectorization. These k nuclei are marked and removed from the list of available sectors. Each nuclei serves as initial subsector of a larger sector.

The total number of possible initial nucleus in Milan is C(n,k), where C, is the number of combinations. For example, for a total of n = 18, and k = 2 sectors there are 153 initial nucleus in Milan. Following the same calculation for k = {2,3,4,...17}, the total number of initializations for the case of Milan is 262,124, plus one configuration where all sectors are isolated and one configuration in which all sectors are connected (Fig. 1b).

In addition, we will call r the set of remaining subsectors not belonging to the nucleus. The total number of allocations that the set r can have into k sectors is kr. Following the same example for n = 18, and k = 2, k = 65,536 possible allocations for each initialization. However, not all possible allocations are feasible, depending on the neighborhood of the subsectors r ^k. For this a different approach for subsectors allocation must be used, based on the initial k nuclei.

The algorithm proceeds to identify the adjacent subsectors to each sector k. In the first iteration, only one subsector (the nuclei) belongs to each k sector. Only subsectors belonging to r that are adjacent to each of the k sectors can be attached to it. For this, the algorithm allocates one by one the remaining subsectors r to any of the k sectors if and only if they share a boundary. This is done until all subsectors are allocated (r = 0). Each time a subsector in r is allocated to a sector k, its cut-set is attached to k, and deleted form the set r. This last step guarantees that subsectors can be allocated also to previously allocated subsectors now belonging to k.

For example, a valid sectorization configuration S2 for Milan is {{s¡, s2, s3, s4, s5, s6, s7, sn, s16}{s«, s9, s10, s12, s13, s14, s ¡s, s17, s«}}, while an invalid sectorization S2 is {{s¡, s2, s4, s5, s6, s7, s¡¡, s16}{s3, s8, s9, s10, s12, s13, s14, s1s, s17, s18}}. The reason for the latter to be invalid is that s3 has no connecting cut-sets to any of the subsectors of the second sector, making this sectorization not feasible, so the algorithm would never generate such sectorization, while a simple random permutation of subsectors may generate such configuration.

Third, once all the subsectors have been allocated (r = 0), it is possible to open the valves which are contained in the same k sector. This is done just by storing a list of the valves of all the cut-sets in a particular sector. The remaining valves remain closed as they were set at the beginning, allowing the isolation into different k sectors.

Finally, many different initializations can converge to the same final sectorization. For that reason, once a new sectorization configuration is created an iterative procedure must be performed to compare it with other configurations stored previously. If the new sectorization is different to all other stored configurations then it is attached at the end of the list of sectorization configurations, until the total required number of configurations nc is reached. This requires a one by one comparison, which has a computational complexity comparable to a sorting algorithm of an order nclog (nc).

In Fig. 3, after applying the algorithm with k = 2, 3, 4, 18, one may obtain different sectorization configurations with nc = 5. This represents only a small subset of sectorization configurations, although for k = 18 only one possible configuration exists. All sectorization configurations are different which is seen by the comparison with previously

obtained sectorization configurations. Different configurations can be collected for future analysis of the WDN of Milan.

Fig. 3. Sample of sectorization configurations in Milan WDN, for k = 2, 3, 4 and 18.

In order to estimate which sectorization is better for Milan a subjacent pump scheduling for each of the sectorization configurations obtained must be performed because each sector contains a number of pumps (Table 1) which may be operated in different ways. However, the application of the subsequent pump scheduling that is part of the overall approach is not covered in this article.

3. Discussion on advantages and limitations of the sectorization algorithm

As it was evidenced by the results presented in Fig. 3, it is possible to fix the number of sectors which compose a WDN, which overcomes one of the issues in current literature of WNS and WNP. The algorithm also has the possibility to obtain multiple configurations of sectorization starting from diverse initializations of the nucleus.

If a large number of configurations must be studied (nc >> 1,000), the final step of comparison of a new one with previously obtained configurations can become computationally expensive.

The algorithm can be extended also to other WDN where the sources are not pump stations but rather gravity sources. This is possible as long as the cut-sets between subsectors have been identified previously and G2 graph can be obtained.

One of the largest limitations of the algorithm originates from the random nature of the selection of subsectors during allocation. This limitation, balances the sizes of the number of subsectors per sector, but impedes the algorithm to easily generate a sectorization where there are large variations of the sizes of sectors. For example, if an operator in Milan wishes to obtain 2 sectors (k = 2), a valid sectorization is S2 = {{si, £2, £3, £4, S6, £7, Ss, £9, s1(), su, S12, £13, S14, S15, Si6, S17, sis},{s5}}. Notice how S5 is isolated from the rest of the WDN. This corresponds to generating the nuclei (Sx, £5), and then allocating the r = 16 subsectors in a row to the same nucleus Sx, or to a probability of 0.516 —0.0000152, which means that on average one must generate 65,536 different sectorization configurations to obtain such a specific one.

On the assumption that, for other WDN's, G2 is available, it is possible to generate all possible sectorization configurations by simple permutations (not combinations). This becomes a computationally expensive problem (N-P complexity) but feasible with current machines for a small number of subsectors (n < 15).

Fig. 4. Cases of sectorization for a WDN composed of n = 4, and Von Neumann neighborhood.

Fig. 5. Eigenvalues of the Laplacian matrix for the first configuration of Fig 4 with k = 1, 2, 3, 4.

One way to analyze whether or not each of the sectorization configurations obtained with a random permutation is valid sectorization or not, is through the analysis of the Laplacian matrix (L) of G2. Matrix L is calculated as D -A, where D is the graph degree (number of edges connected to node a in the graph) as a diagonal matrix, and A is the graph adjacency matrix (Fig. 5). The number of eigenvalues of L equal to zero (neig) corresponds to the number of connected blocks or components of A [17, 18]. In general, a particular value of neig, will define the number of sectors k obtained with a random permutation for G2 if total enumeration is executed. For example, if a WDN is composed of n = 4 subsectors with a Von-Neumann neighborhood (Fig. 4), the total number of random permutations with repetitions is 256, while the real number of sectorization configurations is easily obtained by using the zero multiplicity property of eigenvalues of L. as presented in (Fig. 4), where if k=1, nc=1; k=2, nc=6, k=3, nc=4; and k=4, nc=1. As a proof of concept of Laplacian matrix identification of sectorization, the first sectorization configuration for each k in Fig. 4 is analyzed separately in Fig. 5. The estimation of L is performed. Values marked

in red are equal to zero corresponding to neig which is equal to k for each particular configuration. In general, for each n and k, nc < S(n,k), where S(n,k) are the Stirling numbers of the first kind.

4. Conclusions

We have presented a new algorithm for sectorization of a large WDN, which aims to solve two different open issues in the subject of WNS and WNP: (1) it guarantees that the obtained solution is with a fixed number of sectors, and (2) allows for at least one source to provide water for a specific sector, or multiple sources if aggregation of subsectors is performed. Also a procedure for sectorization verification is introduced requiring only the adjacency matrix of the sectorization configuration.

Subsectors created with this algorithm guarantee that areas inside cut-sets satisfy operational condition in the system (e.g. pressure, demand). If an aggregation of subsectors is performed, larger sectors could be used to address optimization for energy consumption reduction in large WDN by using pump scheduling.

Acknowledgements

The authors would like to acknowledge Metropolitana Milanese S.p.A. for sharing the data of the topology of their network, to SurfSARA, for providing access to their High Performance Computing (HPC). Finally the authors would like to acknowledge EU-FP7 ICeWater project which funded part of this research.

References

[1] von Lucken C.; Baran B. & Sotelo, A. Pump scheduling optimization using asynchronous parallel evolutionary algorithms CLEI, Electronic

Journal, 2004, 7, 1-20

[2] Kurek W. & Ostfeld, A. Multiobjective Water Distribution Systems Control of Pumping Cost, Water Quality, and Storage-Reliability

Constraints Journal of Water Resources Planning And Management, 2014, 140, 184-193

[3] Di Nardo, A. & Di Natale, M. A heuristic design support methodology based on graph theory for district metering of water supply networks.

Engineering Optimization, 2011, 43, 193-211.

[4] Di Nardo, A.; Di Natale, M.; Santonastaso, G.; Tzatchkov, V.G.. & Alcocer-Yamanaka, V.H.. Performance indices for water network

partitioning and sectorization Water Science and Technology: Water Supply, 2014

[5] Giustolisi O.; Ridolfi L. & Berardi, L. General metrics for segmenting infrastructure networks Journal of Hydroinformatics, 2015, 17, 505-

[6] Giustolisi O. & Ridolfi, L. A novel infrastructure modularity index for the segmentation of WDN. Water Resources Research, 2014, 50, 7648-

[7] Di Nardo, A.; Di Natale, M.; Santonastaso, G.; Tzatchkov, V. & Alcocer-Yamanaka, V. Water network sectorization based on a genetic

algorithm and minimum dissipated power paths. Water Science and Technology: Water Supply, 2013, 13, 951-957

[8] Gama, M.C.; Lanfranchi, E. A.; Pan, Q. & Jonoski, A. Water Distribution Network Model Building, Case Study: Milano, Italy. Procedia

Engineering , 2015, 119, 573 - 582

[9] Castro Gama M.E.; Pan Q.; Salman, M. A. & Jonoski, A. Multivariate optimization to decrease total energy consumption in the water supply

system of Abbiategrasso (Milan, Italy). Environmental Engineering & Management Journal, 2015, 14, 2019-2029

[10] Scarpa F.; Lobba A. & Becciu, G. Elementary DMA Design of Looped Water Distribution Networks with Multiple Sources. Journal of Water Resources Planning And Management, 2016. 10.1061/(ASCE)WR.1943-5452.0000639 , 04016011

[11] Castro Gama, M.E. et al., 2015. Model-Based Sectorization Of Water Distribution Networks For Increased Energy Efficiency. In Proceedings of the HIC2014, NY USA. Paper 233.

[12] Rossman, L. EPANET 2 User Manual 2000. U.S. Environmental Protection Agency.

[13] Yazdani A. & Jeffrey, P. WDS vulnerability analysis using weighted and directed network models. Water Resources Research, 2012, 48, W06517, DOI:10.1029/2012WR011897.

[14] Yazdani, A. & Jeffrey, P. Complex network analysis of water distribution systems. Chaos, 2011, 21, 1-11.

[15] Kernighan, B.W. & Lin, S. An efficient heuristic procedure for partitioning of graphs. The Bell System Technical Journal. 1970, February, 291-307.

[16] Karypis, G. & Kumar, V. Parallel Multilevel series k-Way Partitioning Scheme for Irregular Graphs. SIAM Rev, 1999, 41(2)278-300, DOI:10.1137/S0036144598334138.

[17] Buekenhout, F. & Parker, M. The Number of Nets of the Regular Convex Polytopes in Dimension <=4. Disc. Math. 186, 69-94, 1998.

[18] Skiena, S. Implementing Discrete Mathematics: Combinatorics and Graph Theory with Mathematica. Reading, MA: Addison-Wesley, p. 235, 1990.