Comparison of Algorithms for the Optimal Location of Control Valves for Leakage Reduction in WDNs

: The paper presents the comparison of two different algorithms for the optimal location of control valves for leakage reduction in water distribution networks (WDNs). The former is based on the sequential addition (SA) of control valves. At the generic step N val of SA, the search for the optimal combination of N val valves is carried out, while containing the optimal combination of N val − 1 valves found at the previous step. Therefore, only one new valve location is searched for at each step of SA, among all the remaining available locations. The latter algorithm consists of a multi-objective genetic algorithm (GA), in which valve locations are encoded inside individual genes. For the sake of consistency, the same embedded algorithm, based on iterated linear programming (LP), was used inside SA and GA, to search for the optimal valve settings at various time slots in the day. The results of applications to two WDNs show that SA and GA yield identical results for small values of N val . When this number grows, the limitations of SA, related to its reduced exploration of the research space, emerge. In fact, for higher values of N val , SA tends to produce less beneﬁcial valve locations in terms of leakage abatement. However, the smaller computation time of SA may make this algorithm preferable in the case of large WDNs, for which the application of GA would be overly burdensome.


Introduction
The pursuit of benefits in terms of leakage and pipe burst reduction as well as of infrastructure life extension has spurred water utilities to perform the active control of service pressure in water distribution networks (WDNs) [1][2][3]. In fact, if service pressure is reduced without violating pressure constraints for water supply, positive effects arise without the service effectiveness being affected. The active control of service pressure can also be performed through devices that enables electric power generation [4,5].
When WDN sources have pressure surplus compared to demanding nodes, the active control can be carried out by inserting control valves in key locations of the network. Two main lines of research are currently underway in this context: the former, which started with the paper by Jowitt and Xu [6], concerns the optimal location and control of valves, while the latter, which started with the paper by Campisano et al. [7], concerns the real time regulation of these devices.
In the former line of research, which is the topic of the present paper, some initial works were dedicated to setting up algorithms for the optimization of control valve settings to minimize daily leakage, while meeting minimum pressure requirements at WDN demanding nodes. In this context, Jowitt and Xu [6] proposed an algorithm based on iterated linear programming (LP) to minimize leakage at each time slot of the day. Vairavamoorthy and Lumbers [8], instead, proposed the use of sequential quadratic programming, with an objective function that allows minor violations in the targeted pressure requirements. Other papers [9][10][11][12][13][14] addressed both the issues of valve location and control. Reis et al. [9] tackled the optimization of valve locations and settings by using genetic algorithms (GA) and LP, respectively. In detail, LP was embedded in the genetic algorithm to search for the optimal valve settings for each location of control valves proposed as a solution by the genetic algorithm. Araujo et al. [10] and Cozzolino et al. [13] made use of genetic algorithms for both issues. Liberatore and Sechi [11] proposed a two-step procedure for the optimal valve location and control. In the former step, candidate sets for the location of valves are restricted to pipes defined based on hydraulic analysis. The meta-heuristic Scatter Search routines were used in the latter step to identify the best solution in the location and control problems by optimizing a weighted multi-objective function that considers the cost of inserting valves and the penalty for node pressures that do not meet the requirements. Ali [12] made exclusive use of genetic algorithms, in which he incorporated the physical knowledge of the WDN to improve the algorithm efficiency. De Paola et al. [14] made use of the harmony search approach to optimize both control valve locations and settings.
Whereas all the algorithms cited above are based on the single-objective approach, other algorithms [15][16][17][18] were conceived using the multi-objective approach, to construct Pareto fronts of optimal solutions in the tradeoff between number of control valves, as a surrogate for the installation cost, and daily leakage. Pezzinga and Gueli [15] proposed, for optimal valve location, a fully deterministic procedure, based on the sequential addition (SA) of beneficial valves up to a maximum number of valves installable in the WDN has been fixed. LP is used, instead, for optimal valve control. Nicolini and Zovatto [16] used a multi-objective GA to optimize both valve locations and control settings. Creaco and Pezzinga [17,18] used a GA to search for the optimal valve locations, and LP to optimize the control valve settings for each solution proposed by the GA.
Despite the plethora of works available on the subject, there are only a few comparative analyses. In fact, authors have often proposed new algorithms without testing them against those already available in the scientific literature. Furthermore, the merits and demerits of deterministic and probabilistic multi-objective approaches have never been compared. The lack of this kind of comparison in the scientific literature has motivated the present work. The algorithms of Pezzinga and Gueli [15] and of Creaco and Pezzinga [17,18] were chosen as representative for the comparison. The former was chosen because of its fully deterministic approach, which confers significant computational lightness on the methodology. The choice of the latter is due, instead, to its computational effectiveness, which was acknowledged [19] during the battle of background leakage assessment for water networks (BBLAWN) [20].
In the remainder of the paper, algorithms are first described and then applied to two WDNs.

Materials and Methods
In the present section, first the LP used in both SA [15] and GA [17,18] for the optimization of control valve settings is described (Section 2.1). Sections 2.2 and 2.3 are dedicated to describing the different approaches used in SA and GA, respectively, to tackle optimal valve locations.

Fitness Evaluation of the Generic Location of Control Valves
The generic location of control valves, proposed as solution by whatever algorithm (e.g., the algorithms described in Sections 2.2 and 2.3), can be evaluated in terms of installation cost and daily leakage volume W L . The former assessment comes straight away after the control valves have been installed in selected pipes of the WDN model. In the applications of this work, the total number of control valves is considered as a surrogate for the cost. The latter assessment requires determination of the daily pattern of optimal valve settings in the context of WDN extended period simulation (EPS).
Let us assume a WDN with n p pipes and n n nodes = n demanding nodes + n 0 sources. WDN operation can be represented through a succession of N ∆t time slots, all featuring the same duration ∆t. At each time slot, in which source heads and nodal demands are assigned, pipe water discharges and nodal heads can be derived by solving the following equations, derived from previous works [17,18,21]: A 11 MQ + A 12 H = −A 10 H 0 momentum equation where H (n × 1) and Q (n p × 1) are the unknown vectors of nodal heads at the demanding nodes and of pipe water discharges, respectively. H 0 (n 0 × 1) and d (n × 1) are the assigned vectors of source heads and nodal demands, respectively. Matrixes A 12 and A 10 come from the incidence topological matrix A, with size n p × n n . In the generic row of A, associated with the generic network pipe, the generic element can take on the values 0, −1 or 1, whether the node corresponding to the matrix element is not at the end of the pipe, it is the initial node of the pipe, or it is the final node of the pipe, respectively. Starting from A, matrix A 12 (n p × n) is obtained by considering the columns corresponding to the n network nodes with unknown head. A 21 (n × n p ) is the transpose matrix of A 12 . Matrix A 10 (n p × n 0 ) is obtained by considering the columns corresponding to the n 0 nodes with fixed head. A 11 (n p × n p ) is a diagonal matrix, the elements of which identify the resistances of the n p network pipes through the following relationship: where D i is the diameter of the i-th pipe and where roughness coefficient k i , coefficient b i and exponents α, β and γ depend on the formula used to express pipe head losses. As an example, α = 2, β = 5.33 and γ = 2 in the Strickler formula. Matrix M in Equation (1) is used to increase the resistance of the N val pipes fitted with the control valve. It is a diagonal matrix (n p × n p ), in which the diagonal elements corresponding to the N val pipes fitted with control valve are equal to V k (k = 1, . . . , N val ), whereas those corresponding to the n p -N val pipes without valve are equal to 1. Indeed, V k is the valve setting ranging from 0 to 1. For the generic pipe subject to a certain head loss, it represents the ratio of the pipe water discharge in the presence of the regulated valve to the water discharge in its absence. The extreme values of the range represent the fully closed and fully open valve, respectively. In Equation (1), the vector q L (n × 1) of leakage allocated to network demanding nodes can be calculated starting from the following Equation (3): where the elements of the vector Q L (n p × 1) of leakage outflows from WDN pipes can be assessed through the following relationship [6]: where, for the i-th pipe, h i and L i are the mean pressure head (ratio of pressure to specific weight of water) and the length, respectively. C L,i and n leak are empirical coefficients. The mean pressure head h i can be calculated as: where H i,1 , H i,2 and z i,1 , z i,2 are the heads and elevations, respectively, for the end nodes of the pipe. At each time slot of WDN operation, the vector V of valve settings V k (k = 1, . . . , N v ) can be optimized to minimize the total leakage volume W L,j from the network, while meeting the minimum pressure head requirement h des at the demanding nodes: This can be done using iterated linear programming (LP). This algorithm was first proposed by Jowitt and Xu [6] and then refined by Pezzinga and Gueli [15] and by Creaco and Pezzinga [17,18] through the implementation of an underrelaxation method for setting update throughout the iterations.
The total daily leakage volume can be obtained as the sum of the leakage volumes W L,j at each time slot:

Optimal Location through Sequential Addition of Valves (SA)
The algorithm based on the sequential addition of valves (SA) was proposed by Pezzinga and Gueli [15]. In other words, SA is based on the optimal location of each single valve one after the other, independently from the potential benefit of the remaining ones. First, a maximum number N max of valves installable in the network has to be fixed. The total number of steps in the algorithm is then N max + 1. At Step 0 of this algorithm, the WDN has no control valves. At Step 1, the optimal location of 1 valve is searched for among the available locations in the WDN, i.e., all the n p pipes. To this end, the control valve is placed inside the WDN model at one potential location at a time. The algorithm described in Section 2.1 is applied, enabling assessment of the installation cost of the system (or rather the total number of control valves as a surrogate, banally coinciding with the step of SA) and assessment of W L . Then, the most beneficial valve, which yields the largest W L reduction compared to the no valve scenario, is detected. At Step 2, the optimal location of two control valves is searched for, while keeping the first optimal control valve obtained from Step 1 installed in the WDN. The second valve for the optimal combination of two valves is searched for among the available locations, i.e., the n p − 1 remaining pipes. The algorithm of Section 2.1 is applied considering n p − 1 combinations of two control valves and the most beneficial one is detected in terms of W L reduction compared to the 1 valve scenario. Other steps of SA can be carried out up to N max , always considering the following basic assumption: at the generic Step N val , the optimal combination of N val valves is searched for, while containing the optimal combination of N val − 1 valves detected at the previous step. Considering the number N max of control valves installable in the WDN, SA would require the following number C * v of locations of control valves to be evaluated with the algorithm described in Section 2.1: C * v is much smaller than the total enumeration C v with no repetitions, which is given by: Therefore, SA consists in a deterministically driven exploration of the research space of possible locations of control valves. The exploration of SA is expected to be very effective in the cases where the sequentially added valves do not affect each other. Otherwise, in those cases where the non-linearities of the problem are high, the effectiveness of SA is expected to diminish.
By plotting the W L values obtained through SA as a function of N val , with 0 ≤ N val ≤ N max , a Pareto front of optimal solutions with 0 ≤ N val ≤ N max is obtained.

Optimal Location through Multi-Objective Genetic Algorithm (GA)
The multi-objective Non-Dominated Sorted Genetic Algorithm II (NSGA II) [22] was used to optimize the optimal location of control valves. In this GA, the adaptive mutation operator described by Carvalho and Araujo [23] was implemented. The decisional variables are encoded in individuals (i.e., solutions), with several genes equal to the number n p of potential locations of control valves. Each gene can take on two possible values, 0 and 1, which indicate the absence and presence of the control valve at the generic pipe, respectively. The number N val of valves for the generic solution is banally obtained by summing up the gene values.
In GA, the number of population individuals and the total number of generations must be fixed in light of the trade-off between computation time (which is related to the total number of objective function evaluations and is then a function of the available computation capabilities) and accuracy of the results. Incidentally, the total number of function evaluations is given by the total number of population individuals times the total number of generations.
The various individuals of the GA, which represent various locations of control valves installed in the network, are compared based on their fitness, composed of two objective functions evaluated through the algorithm described in Section 2.1, i.e., cost, or rather number of control valves, and W L . Both the objective functions need to be minimized inside the GA. Similar to the algorithm for SA, the results of the multi-objective optimization are Pareto fronts. Unlike SA, in which the Pareto fronts feature a maximum number of valves equal to N max depending on the number of SA steps, the GA Pareto fronts are made up of solutions with N val ranging from 0 to n p .
GA carries out a probabilistically driven exploration of the research space. Unlike SA, GA is theoretically able to explore all the possible locations.

Case Studies
The application of this work concerned two case studies (Figure 1a,b, respectively). The first is the skeletonized WDN of Santa Maria di Licodia, a small town in Sicily, Italy [24]. The use of a skeletonized WDN is not expected to undermine the validity of the results since the pipes with larger diameter, which compose the WDN skeleton, are typically the best candidates for the installation of control valves.
The network layout is made up of n n = 34 nodes (of which n = 32 with unknown head and n 0 = 2 source nodes with fixed head, i.e., Nodes 33 and 34) and of n p = 41 pipes. The features of the WDN demanding nodes in terms of ground elevation z, mean daily demand d, and desired pressure head h des for full demand satisfaction, are reported in Table 1. At the generic demanding node, the latter variable was set at the minimum value between 15 m and the daily lowest pressure head under uncontrolled conditions. This was done to avoid any pressure deficit increase at nodes that were pressure deficient ab initio. Table 2 reports the features of the pipes in terms of length L, diameter D and Strickler roughness coefficient k. The two source nodes, i.e., Nodes 33 and 34, have a ground elevation of 477.5 m a.s.l.
N ∆t = 12 2-h-long time slots were used to represent the WDN daily operation. The patterns of the heads at the source nodes and of the hourly multiplying coefficient C h for nodal demands are reported in Table 3.
As to leakage, two conditions were considered. In the former, values of the coefficient C L = 2.79 × 10 −8 m 0.82 /s and of the exponent n leak = 1.18 were assumed for all the pipes of the network. The values of C L and n leak lead to a daily leakage volume W L = 1242.6 m 3 , which is about 44% of the water volume leaving the source nodes, as evaluated in the real network. In the latter, the coefficient C L was reduced to 0.85 × 10 −8 m 0.82 /s, obtaining a daily leakage volume W L = 397.7 m 3 (about 20% of the water volume leaving the source nodes) without control valves.
The second case study is a district of the pipe network model used as benchmark in the Battles of Water Networks of the last WDSA conferences [19,20]. This district is made up of n n = 46 nodes (of which n = 45 with unknown head and n 0 = 1 source node with assigned head, i.e., tank Node 46) and n p = 52 pipes. The features of the district demanding nodes and pipes are provided by Creaco and Pezzinga [16], who also reported the values of the leakage exponent n leak = 0.9 and of the leakage coefficient C L ranging from 0.2 to 1 m 1.1 /s in the various pipes. In this district, the leakage volume adds up to 118.7 m 3 , which is 6.1% of the water volume entering the district, in the no valve scenario.   250  75  22  110  125  65  2  314  175  65  23  214  150  65  3  1100  125  75  24  85  100  65  4  350  100  65  25  398  100  65  5  96  250  75  26  242  100  65  6  282  100  75  27  118  175  65  7  148  250  75  28  324  175  65  8  256  250  75  29  140  125  65  9  192  100  75  30  206  125  65  10  58  100  75  31  70 125 In the application, values equal to 105.8 m a.s.l. and 2.4 m, respectively, were assigned to the elevation and to the average water level of the source node. An inflow (that is a negative demand) takes place in correspondence to Node 1, at certain hours of the day, due to the activation of a pumping system that takes water from another district of the network. The detailed operation of this system, in terms of relationship between water discharge and head gain across the pump, is neglected in this work.
The district operation can be summarized in three time slots, for each of which the head of the source node, the inflow from the pump and the hourly demand coefficient for the demanding nodes were specified by Creaco and Pezzinga [18]. No valve installation is allowed in the main line (Figure 1b) that connects the pump with the tank node, to prevent any interferences with the filling/emptying process of the tank.

Results
In the application with SA [15] to the first case study, the optimal locations of up to N max = 10 control valves were searched for. This required 365 objective function evaluations.
In the application with GA [17,18], a population of 50 individuals and a total number of 50 generations, corresponding to total number of 2500 objective function evaluations, were used. A healthy initial population was generated in GA to guarantee high computation efficiency at the end of the optimization.
The applications of this work were performed in the Matlab(R) 2011b environment by making use of a single processor of an Intel(R) Core(TM) i7-7700 3.60 GHz unit.
The tradeoff curves of W L as a function of N val obtained with SA and GA in the first leakage condition are reported in Figure 2. As for GA, two Pareto fronts are shown, the final one (after 50 generations) and that after seven generations. The latter is associated with 350 objective function evaluations (very close to SA). Though GA explored solutions with N val values up to n p (see Section 2.3), only the solutions with N val ≤ 10 were reported in the graph for a consistent comparison with SA. This graph shows that the curve of SA and the final curve of GA have identical trend up to N val = 3. To the right of N val = 3, the curve obtained with GA dominates that of SA, in that it provides lower values of W L for each value of N val . Furthermore, the distance between the two curves increases as N val grew. This seems to suggest that, for low number of N val , SA can provide identical or close results to those of GA. Conversely, for high number of N val , the better performance of GA stands out, due to the wider exploration of the research space. However, this comes at an about seven times larger computational cost. As was expected, the curve of GA after seven generations features worse solutions than the final curve of GA to the right of N val = 3. The curve of GA after seven generations enables the results of GA to be compared with those of SA, given the same computational budget. Interestingly, neither curve prevails in absolute terms. In fact, identical solutions are observed up to N val = 3. GA after seven generations is slightly better for N val = 7, 9 and 10, while not offering any solutions for N val = 8. SA, instead, is better for N val = 5 and 6. Summing up the results in Figure 2, the better performance of GA stands out only with a higher computational budget. An insight into the different results obtainable with SA and GA is provided in Table 4, which reports the optimal valve locations and the W L values for the two algorithms. Whereas the valve positions suggested by SA and GA are the same up to N val = 4, the two algorithms behave differently starting from the optimal location of five valves. In fact, at Step 5, SA suggests adding valve in Pipe 33 to the optimal location of N val = 4 valves, in which the valve-fitted pipes are 3, 7, 14 and 27. Instead, for the sake of optimality, GA proposes valve elimination in Pipe 14 and valve insertion in Pipes 25 and 26, while moving from N val = 4 to N val = 5. The locations of SA and GA for N val = 5 valves are shown in Figure 3. As Table 4 shows, this yields a 4.36% benefit of GA compared to SA, in terms of W L . The percentage benefit of GA tends to grow up to almost 10% for N val = 10.  Another example of the results obtained is provided in the graphs in Figure 4, associated with the optimal location of three valves, installed in Pipes 3, 7 and 27, respectively, for both SA and GA. Figure 4a shows the temporal pattern of V, optimized through the iterated LP at each time slot. Figure 4b shows the pressure head h down at the downstream node, which can be used as the time varying setting to be prescribed if the water utility chooses to perform the pressure regulation by means of pressure reducing valves (PRVs). In fact, PRVs has the downstream pressure head as setting. This variable is always equal to h des (that is 15 m) at the valve in Pipe 3 because the downstream node of Pipe 3, namely Node 32, is a terminal node of the network. It is always equal to h des also for the valve in Pipe 7, though the downstream node of this pipe, namely Node 3, is not terminal. This happens because the descending topography downstream of Pipe 7 facilitates the meeting of the pressure requirements. The case of the valve in Pipe 27, which has Node 9 as downstream node, is similar. However, this valve cannot reduce h down to h des at nighttime, even if it is almost fully closed (V ≈ 0). An analysis should also be carried out concerning the possibility to convert the control valves proposed by the optimizer to isolation valves. This could be done when the optimal settings V are very close to 0 throughout the day. Furthermore, isolation valves could be closed in some pipes in parallel to the installed control valves if water flow is remarked to bypass the latter. However, neither situation was remarked to occur in the present case study. Furthermore, it must be noted that closing an isolation valve may decrease WDN redundancy and reliability.
The tradeoff curves of W L as a function of N val obtained with SA and GA in the second leakage condition are reported in Figure 5. This figure leads to similar remarks as Figure 2, as far as the comparison of the two algorithms is concerned.  Table 5 gives some insight into optimal valve locations, W L values and benefits of GA compared to SA. The comparison of Table 5 (low starting leakage) and Table 4 (high starting leakage rate) points out that the reduction in the starting leakage has the following minor effects:

•
The benefits of GA in Table 4 tend to grow with N val increasing, whereas the values in Table 5 tend to stabilize around 6.5%. The applications to the second case study were carried out similarly to the first case study. The results are reported in Figure 6 and Table 6. The former reports the Pareto fronts obtained through SA and GA, whereas the latter reports optimal valve locations and W L provided by the two algorithms, as well as the benefits of GA. Similar to the first case study, these benefits stand out only for N val ≥ 5. However, they are smaller (always below 4.32%), due to the different structure of the WDN. In fact, while the network of Santa Maria di Licodia is very interconnected, the network of the second case study is made up of quite independent branch structures fed by the main line, along which valves cannot be installed. This attenuates the non-linearities pertinent to optimal valve location.

Conclusions
In this work, two algorithms for the optimal location of control valves, aimed at reducing service pressure and leakage in water distribution networks (WDNs), were compared. The former algorithm is deterministic and explores the research space of possible locations by adding one control valve at each step (sequential addition of valves SA). This exploration results in a significant reduction in the research space. The latter algorithm is a multi-objective genetic algorithm (GA), which explores, using a probability driven approach, the whole research space without restraints. The applications to skeletonized WDNs show that the two algorithms give identical results for a small number of installed control valves. However, GA outperforms SA when the number of installed valves increases. This may be because the non-linear effects of the optimization problem, totally neglected by SA, become more and more predominant as the number of installed valves grows. The benefits of GA are larger in the very interconnected WDNs, in which the non-linearities are emphasized. However, it must be noted that SA has a much smaller computational overhead and may then be preferable in the case of overly large networks. Furthermore, valve cost increase typically discourages water utilities from installing too many valves in the network. However, an advantage of GA lies in the possibility of accounting for other issues in the evaluation of the fitness of the generic solution. As an example, different and more complete objective functions could be easily and more rigorously considered inside GA, such as a relationship to express valve cost as a function of the diameter of the valve-fitted pipe. Furthermore, the closure of isolation valves, which are already available in WDNs to isolate segments [25], can be accounted for in GA to improve leakage reduction, as was shown in the works [17,18].