Improved community detection for WDN partitioning in the dual topology based on segments and valves

One of the most explored topics in the recent scientific literature of water distribution networks (WDNs) is partitioning, which enables subdividing large networks into sufficiently small parts, typically called district-metered areas (DMAs), which can be easily monitored and managed by water utilities in terms of users' consumption and leakage (Morrison et al. 2007). This also results in fast detection of anomalies, e.g., unexpectedly high consumption and pipe bursts, and in more effective maintenance campaigns. WDN partitioning has also been implemented in commercial software over the last decade (e.g., WaterGEMS).

In most cases, the partitioning is carried out in two phases, i.e., clustering and dividing (e.g., Di Nardo et al. 2018; Giudicianni et al. 2020). The clustering is typically performed by graph theory techniques and aims to identify the size and extension of DMAs, as well as the location of boundary pipes. As the final aim of partitioning is the real-time knowledge of water exchanges between DMAs to facilitate water mass balance at the WDN level, the dividing follows to identify the boundaries to be closed or fitted with a flow meter, to have known water exchanges equal to 0 or measured water exchanges, respectively.

A widely adopted concept in the context of WDN partitioning is modularity, which was first defined by Clauset et al. (2004) in the context of complex network theory to assess the goodness of network subdivision into communities, associated with balanced cluster sizes and few boundary links. While the concepts of modularity and community detection have been adopted in numerous algorithms for WDN partitioning (Campbell et al. 2016; Zhang et al. 2017; Brentan et al. 2018; Rahmani et al. 2018; Creaco & Haidar 2019; Creaco et al. 2019; Yao et al. 2021; Sharma et al. 2022), Giustolisi & Ridolfi (2014a, 2014b) pointed out that the original formulation of modularity falls short of the representation of hydraulic properties for WDNs, e.g., water demand, pressure, pipe length, and background leakage, thus proposing a revisited one. Another formulation of modularity tailored to WDNs has recently been proposed by Shao et al. (2022), in which the modularity component associated with boundary pipes is independent from that associated with the water demand to be made uniform across DMAs. Though being significant contributions to this research field, the revisited formulations of modularity proposed so far fail to consider the actual distribution of isolation valves, which subdivides the WDN into segments. Incidentally, segments are the smallest WDN parts that can be disconnected for maintenance purposes. As Santonastaso et al. (2019) pointed out, DMAs must originate from the merging of WDN segments and therefore WDN partitioning algorithms must be formulated based on the dual topology (Walski 1993) based on segments and valves present in the WDN.

To bridge this research gap, the present work aims to propose a new formulation of modularity specifically tailored to WDNs and based on the dual segments/valves topology of the WDN. As an additional improvement from the previous works, the new formulation can also deal with properties to be made uniform inside DMAs, rather than only across DMAs. This is the case with nodal ground elevations or pipe age, the uniformity of which inside each DMA is a desired condition to facilitate pressure regulation/leakage mitigation or pipe maintenance, respectively. Furthermore, an optimization algorithm inspired by the simulated annealing (Laarhoven & van Peter 1987) approach is proposed for the clustering of WDN segments into the desired number of DMAs.

The following sections describe the methodology and its applications to two case studies. The paper ends with the discussion and the conclusions.

The new clustering algorithm is based on the dual topology based on segments and isolation valves (Walski 1993), instead of the traditional topology based on nodes and pipes. The disconnection of the generic segment can be obtained by closing the isolation valves that separate it from its adjacent segments. In this context, the reason for the choice of the dual topology is that its adoption guarantees that isolation valves are present in the boundaries between DMAs after the clustering, as was originally proven by Santonastaso et al. (2019) and later considered by other authors, such as Creaco et al. (2022). The use of the dual topology brings about the two following advantages, which make the clustering solutions obtained more feasible:

• – In the boundaries to be closed during the dividing phase, there is no need for installing additional isolation valves.

• – In the boundaries to be left open during the dividing phase, the valve vault can sometimes be used to host the flow meter, which must be installed for monitoring the inter-DMA exchange of water, entailing other potential money savings.

In Equation (3), U_i is the property value in the generic i-th DMA, obtained as the sum of the segment values inside, and is the total value of the property in the WDN. As an example, by considering the DMA demand D_i and the total WDN demand D_tot in place of U_i and U_tot, respectively, the property to be made uniform is the demand. The lower the difference in the demand between the DMAs, the lower the value of H₂. When all the M DMAs feature the same value of demand, the value of H₂ is the smallest. It grows, instead, when the values of DMA demands are dissimilar.

The three contributions subtracted from 1 in Equation (1) are all positive. Therefore, index Q cannot exceed 1. It will be close to 1 in the case of good clustering, in which the three contributions H are small. Otherwise, it will decrease from 1. Finally, it must be remarked that coefficients α play the role of weights, to give leverage to the three elements of modularity in the context of the clustering. Indeed, though all values can be assigned to the weights, the maximization of the modularity in an optimization framework is only affected by their mutual relationship. Hereinafter, weight coefficient values with a sum equal to 2 are considered. Since the generic weight can then range from 0 to 2, a low value close to 0 will be assigned to decrease the leverage of the corresponding contribution of the modularity. A high value close to 2 will, instead, be used to increase its leverage. As general rules for assigning these values do not exist, a sensitivity analysis must be carried out to analyze the clustering solutions obtained in each case study for different triplets of weights. Even if the number M of DMAs can be preventively assumed, based on the desired partitioning of the whole demand, the sensitivity analysis can also be extended to assess the impact of this parameter. In comparison with most methodologies present in the scientific literature, which are based on fixed parameters, the variable parameters considered in the present methodology can give flexibility to the clustering, while accounting for various issues such as those related to topology, demands and ground elevations. The various clustering configurations obtainable as the values of weights α and parameter M change must be assessed based on topological, hydraulic and engineering judgment criteria, leading to the choice of the ultimate solution. While the results of the clustering can be of assistance to water utility managers, the ultimate solution may be edited to better match the presence of geographic/natural constraints related to, as an example, water bodies, major limited access highways, railroad tracks and parks/golf courses present in the territory.

The clustering algorithm proposed in the present work enables refining a clustering solution of the first attempt featuring the desired number M of DMAs. A good solution of the first attempt can be obtained by applying the FGA of Newman & Girvan (2004) to the dual segment/valve topology (e.g., see Figure 2(b)). Initially, the number of DMAs is therefore equal to the number N_s of segments, in other words, each segment represents a DMA of its own. Then, all possible merging actions between adjacent DMAs are considered in step 1. These actions are easily identified based on the topological connections in the dual segment/valve topology. As an example, Figure 2(b) shows that segments S1 and S2 can be merged because they are adjacent. S1 cannot be merged, instead, with segment S8 because no valve connects these two segments. For each possible merging in the network, the associated variation ΔQ is calculated. After testing all possible merging actions, the one yielding the largest positive value of ΔQ is applied at the end of step 1, producing a new WDN clustering into a number of DMAs equal to N_s − 1. Starting from this new clustering configuration, step 2 explores all potential merging actions, to identify that yielding the largest positive variation ΔQ, which will be the basis for step 3, and so forth. The steps are repeated N_s − M times until the number of DMAs equals M. Although being very effective, FGA does not reach the optimal WDN clustering into M DMAs in terms of Q. In fact, considering the largest ΔQ value at each step does not guarantee that the highest possible value of Q is reached at the (N_s − M)-th step. As there is no guarantee that a sequence of optimal steps results in the optimal final value of Q, the solution obtained with FGA lends itself to being refined, as is explained below in Section 2.4.

The optimization algorithm considered in the present work is inspired by the simulated annealing (Laarhoven & van Peter 1987) approach and operates by refining the WDN clustering solution obtained with FGA, or alternatively a generic configuration of the first attempt indicated as C₀, in which the constraints (6) and (7) hold. In the context of the present methodology, a boundary segment is a segment connected to a valve lying in the boundary between two adjacent DMAs. At the generic iteration, the algorithm i) explores the potential exchanges of boundary segments between adjacent DMAs, while ii) preventing the transfer of segments from DMAs made up of only one segment. The actions (i) and (ii) enable constraints (6) and (7) to hold and the number M of DMAs to be unaltered during the whole optimization process.

The following subsections explain the operation of the optimization algorithm at the generic iteration and the framework of iterations, respectively.

The boundary segments to be exchanged can be easily identified in the dual segment/valve topology by looking at the isolation valves actively separating the DMAs. As in Figure 4, Table 3 links active isolation valves to the potential segment exchanges between the DMAs.

In most cases, the transfer of one segment from a giving DMA does not create any disconnections in this DMA, which then stays feasible. However, in some cases, the transfer of one segment from a giving DMA may cause this DMA to remain subdivided into n_r > 1 remaining parts. As considering these n_r parts as DMAs of their own would cause the increase in M beyond the prefixed value, only the largest of the n_r parts in terms of number of segments is kept in the giving DMA; the other n_r − 1 parts are transferred, instead, to the receiving DMA.

The generic k-th of the N_e possible segment exchanges will be associated with a variation ΔQ_r,k in the modularity Q evaluated with Equation (1), with the subscript r standing for the word refinement.

Considering the set of all N_e possible variations ΔQ_r,k, the selection of the segment exchange to consider for the n-th iteration is made in a probabilistic way. Instead of deterministically selecting the segment exchange yielding the largest positive ΔQ_r,k, the potential variations ΔQ_r,k are sorted in ascending order and their cumulative frequency F_k for the selection is calculated. Then, a random number rand is generated between 0 and 1 and the segment exchange corresponding to the ΔQ_r,k value with the smallest cumulative frequency F larger than rand is selected, leading to the new WDN clustering configuration C_n. In comparison with C_n−1, C_n may be better or worse, as variations ΔQ_r,k may be either positive or negative.

Summing up, the calculation of F_k based on Equation (8) can modulate the probabilistic selection of segment exchanges. Indeed, for low values of k_val, the approach is fully probabilistic with no privilege for the good segment exchange solutions. In this case, a worse WDN clustering solution C_n than C_n−1 in terms of Q is likely to be obtained. The privilege for good segment exchanges grows with k_val increasing. For high values of k_val, the probability of obtaining a worse C_n than C_n−1 decreases. For k_val = 99, the exchange becomes fully deterministic leading to the choice of the best exchange solution for the growth of Q (steepest ascent).

A maximum number N_max of iterations must be set (e.g., N_max = 2,000).

The overall operation of the algorithm can then be explained as follows. When the algorithm is run, it starts with n = 0 and n_stag = 0. As n grows, k_val grows according to Equation (9). The algorithm then changes smoothly from fully probabilistic to fully deterministic (e.g., see Figure 5). This change in behavior enables first exploring the research space, consisting of modified WDN clustering configurations in comparison with the one of the first attempt, and then searching for a local maximum in a certain research space area. When Q increases no more, n_stag is set to n. Then, starting from the next iteration, i.e., when n = n_stag + 1, the algorithm becomes fully probabilistic again. Again, the growth of n causes the nature of the research space exploration to progressively change from probabilistic to deterministic, until a new local maximum of Q is found, in correspondence to which n_stag is set to n. The iterations proceed up till n = N_max, with the ultimate WDN clustering solution being the one featuring the highest value of Q during all the N_max iterations, i.e., the maximum of all the local maxima of Q encountered during the iterations.

Two case studies are considered in the present work. The first case study is the WDN of Ferrara (Alvisi et al. 2011), with two source nodes with a fixed head of 30 m, 536 demand nodes and 825 pipes. Overall, 969 isolation valves are present in the WDN, plus two isolation valves at the exit of the sources, leading to the formation of 682 segments. The first case study was used to search for WDN clustering configurations featuring uniform demands, as the uniformity of nodal elevations inside each DMA is per se guaranteed by the flat topography of the WDN area.

The second case study is the WDN of Saadnayel in Lebanon (Creaco & Haidar 2019), made up of one source node with a fixed head of 962 m, 152 demand nodes and 186 pipes. Each pipe is assumed here to be equipped with an internal isolation valve, leading to a total of 153 segments, each made up of a single node. Nodal ground elevations lie between 874.3 and 911.8 m above sea level (asl). Due to the larger diversity of ground elevations than the WDN of Ferrara, the second case study was then used to search for WDN clustering configurations featuring uniform demands across the DMAs and uniform ground elevations inside each DMA.

More information about the two case studies can be found in the referenced works. Besides the application of the clustering methodology presented in this paper, the following subsections also report some results of the dividing phase, in which the inter-DMA boundaries are either closed or fitted with a flow meter. In the end, the flow meters present in the WDN enable estimating the overall demand (user consumption + leakage) of each DMA. When all DMAs are supplied in parallel by single flow meters (one for each DMA), each flow meter reading represents the overall demand of one DMA. Otherwise, trivial addition/subtraction calculations must be carried out on the readings to estimate the overall demand for intermediate DMAs.

The applications of this work were carried out in the Matlab^® 2022b environment by using a single unit of a 12th Gen Intel(R) Core(TM) i9-12900H 2.50 GHz processor. To speed up calculations, the algorithm can also be applied with parallel computation to WDNs of any size, with no restraints on the number of nodes and pipes. H₁ in Equation (2) was calculated as n_b/n_v. H₂ and H₃ were calculated considering demands and ground elevations, respectively.

In the first case study, the clustering into M = 8 DMAs was initially considered, like in Creaco et al. (2022). The performance of the novel algorithm in comparison with FGA is evaluated in terms of values of Q, H₁, H₂, number n_b of boundary pipes and coefficient of variation C_v,d for DMA demands. Incidentally, C_v,d is a dimensionless coefficient calculated as the ratio of mean to standard deviation, which expresses the heterogeneity of DMA demands: the larger C_v,d, the more heterogeneous the DMA demands. The choice of reporting this coefficient in the analysis is due to its presence in previous works in the scientific literature (Creaco et al. 2019, 2022).

In the case of M = 8, FGA and the novel optimization algorithm were applied considering six different triplets of values for the weights in Equation (1). In all triplets, α₃ was set to 0 since the uniformity of ground elevations was not investigated in the first case study. While the sum of α₁ and α₂ was kept equal to 2 in all triplets, various combinations of α₁ and α₂ were considered, as reported in Table 4. This was done to analyze the extent to which the WDN clustering configurations react to the increase in the weight α₁ from 0.1 to 1.9 given to the first term of the modularity, associated with the boundary pipes, and then to the simultaneous decrease in the weight α₂ from 1.9 to 0.1 given to the second term of the modularity, associated with the demand. Table 4 shows that, in every triplet of α₁, α₂ and α₃, the optimization algorithm applied to the solution of the first attempt obtained by means of FGA was beneficial in that it yielded increases in Q, as well as decreases in H₁, H₂, n_b and C_v,d. These improvements are the evidence of the suboptimal performance of FGA, which lends itself to be improved by means of the optimization algorithm. Looking at the values of n_b and C_v,d obtained after the optimization algorithm, which are more meaningful from the engineering viewpoint, the increase in α₁ at the expense of α₂ tends to produce WDN clustering configurations with fewer boundary pipes and more heterogeneous DMA demands. This is because the increase in α₁ and the decrease in α₂ give more weight to H₁ and less weight to H₂ in Equation (1).

Though neglected in the present work, water quality implications may arise from closing valves. Additional water quality simulations could be carried out to identify potentially negative issues, suggesting changes in the boundaries to be closed or fitted with a flow meter.

An additional analysis was carried out in the first case study to investigate the extent to which the change in the number of DMAs affects the performance of the clustering optimization. In this context, Table 5 shows the results obtained in terms of Q, H₁, H₂, n_b and C_v,d when the number of DMAs is increased from eight to 13, using the same triplets α₁, α₂, α₃ as in Table 4. The analysis of Table 5 points out for M = 13 similar remarks to M = 8. Once more, the optimization algorithm applied to the starting solution obtained by means of FGA yielded increases in Q, as well as decreases in H₁, H₂, n_b and C_v,d, confirming the suboptimality of FGA and the possibility of improving the results with the optimization algorithm. With α₁ and α₂ increasing and decreasing, respectively, WDN clustering configurations with fewer boundary pipes and more heterogeneous DMA demands are obtained, because of the higher weight given to H₁ than to H₂ in Equation (1). Overall, the comparison of the optimization results in Tables 4 and 5 shows that the increase in M has detrimental effects in n_b and C_v,d. The triplets being the same, configurations with less uniform DMAs and more numerous boundary pipes are obtained when M increases.

In the second case study, the clustering into M = 5 DMAs was considered, like in Creaco & Haidar (2019). Two triplets were considered for the parameters in Equation (1), as is shown in Table 6. In the first triplet, α₁ = 0.2, α₂ = 1.8 and α₃ = 0.0 were considered, i.e., no weight was given to the uniformity of nodal elevations inside each DMA. The second triplet, instead, was set to α₁ = 0.15, α₂ = 0.15 and α₃ = 1.7 to take into account the impact of the uniformity of nodal elevations inside each DMA.

Then, the FGA and the optimization algorithm were applied to the second case study considering both triplets. In the application with the second triplet, the two algorithms were fictitiously applied considering M = 6, as the source, which has a very different elevation from the remainder of the WDN, was always identified as a whole DMA of its own.

In this work, a new formulation of modularity tailored to WDNs was proposed and implemented inside an effective optimization algorithm inspired by the simulated annealing approach, able to refine an initial clustering solution of the first attempt. As an improvement in comparison with the others present in the scientific literature, this formulation has the following two beneficial features for the sake of WDN clustering into DMAs:

• – It is based on the dual topology based on segments and valves.

• – It considers properties to be made uniform both across and inside DMAs.

As proven by the applications to the first case study of this paper, the use of the dual topology in the modularity enables obtaining configurations of DMAs that take advantage of the isolation valves actually present in the WDN. For the boundaries to be closed, there is then no need of installing additional isolation valves. Furthermore, the already available valve-vaults can sometimes be used to host the flow meters in the boundary pipes that are left open after the dividing phase. As proven by the applications to the second case study, the possibility of considering properties to be made uniform inside DMAs can be exploited to obtain clustering configurations with uniform nodal elevations, suitable for carrying out pressure management in the WDN.

The performance of the optimization algorithm for the maximization of the modularity in the clustering phase proved outstanding in the trade-off between numerical effectiveness and computational burden, taking as benchmark results obtained in previous works.

All the clustering solutions obtained perform well from the hydraulic viewpoint during the dividing phase, in which the decision of closing the valve or installing a flow meter must be taken at each boundary.

Furthermore, the methodology developed enables overcoming one of the main difficulties encountered by previous authors with the clustering, that is the possibility of considering nodal (instead of pipe) properties as an objective to pursue. In fact, while weighted coefficients are typically assigned to links, the new framework of this work can assign weights to both nodes and links.

Though the clustering applications of this work considered only the minimization of number of boundary pipes, heterogeneity of demand across DMAs and heterogeneity of nodal ground elevations inside DMAs, the framework proposed is very flexible. Therefore, other properties can be easily included, such as boundary pipe diameter, pipe length across DMAs and pipe age inside DMAs.

As the applications conducted in this paper concerned small to medium-sized networks, future work will be dedicated to extending the methodology to complex and large networks with numerous control devices (e.g., pumps and control valves).

Support from Italian MIUR and University of Pavia is acknowledged within the program Dipartimenti di Eccellenza 2023–2027.

All relevant data are included in the paper or its Supplementary Information.

The authors declare there is no conflict.