## Abstract

Recently, calibration based methods are considered to determine the quantities and locations of all leakages in a water distribution network (WDN), simultaneously. In this paper, the accuracy of detecting leakages using Ant Colony Optimization (ACO) in two networks including a hypothetical and a laboratory network is investigated. A novel method is introduced to analyze the effects of measurement errors on demand calibration in an under-determined problem. The results confirm the ability of the optimization-based method in detecting the exact locations and values of the leakages in a WDN. Also, it is emphasized that merely suitable fitnesses cannot provide sufficient confidence in result accuracy. To overcome this difficulty, the correctness of leakage detection can be verified in a sampling design problem. Moreover, to check the reliability of leakage detection an alternative method based on various results is introduced in this research.

## HIGHLIGHTS

A new method is presented.

Both experimental and analytical methods are investigated.

## INTRODUCTION

Leakage reduction in a WDN is essentially needed to have an economic, reliable, and safe network. Low quality of pipeline and junction materials, design errors, improper maintenance, high pressure and also random failures are the main reasons causing serious damage within a network (Lay-Ekuakille *et al.* 2009). Currently many methods such as acoustic equipment, mass balance technique, thermography, ground-penetrating radar, tracer gas, and video inspection are used for detecting and locating leakages within networks (Covas & Ramos 2010). Although these approaches provide reliable and precise measurements in some conditions, they are economically expensive, labor intensive and often inaccurate. Amongst these approaches, the equipped acoustic is extensively used by the water industry, which properly measures leakages for metal pipes but gives inaccurate results for other pipe materials such as plastic or PVC (Hunaidi *et al.* 1998, 1999).

Recently model base methods such as transient-based (Liggett & Chen 1994; Brunone 1999; Covas & Ramos 1999; Brunone & Ferrante 2001; Covas *et al.* 2005) and steady state techniques (Carpentier & Cohen 1991; Pudar & Liggett 1992; Andersen & Powell 2000; Poulakis *et al.* 2003; Abhulimen & Susu 2004; Cheng & He 2011) have been considered by many researchers to cope with the shortcomings caused by the conventional detection techniques concisely outlined above. However, transient based methods can find leakage location more accurate than steady based methods. Notice that it has some drawbacks like generation of negative waves which make this method difficult for complex networks. Also, it cannot be used for a large network, simultaneously.

Pudar & Ligget (1992) used the inverse steady state analysis method for detecting leakages in pipe networks (Pudar & Liggett 1992). Almandoz *et al.* (2005) and Walski *et al.* (2006, 2008) developed an approach for detecting leakages by network hydraulic modeling (Almandoz *et al.* 2005; Walski *et al.* 2006). Wu *et al.* (2010) proposed an optimization-based method for locating and sizing the water losses obtained for the WDN through the calibration process (Sage & Wu 2006). In this method, the demand multiplier factors have been adjusted for node groups using a genetic algorithm optimization method. In this leakage detection method, leakage is modeled as additional demands in the nodes, and then total nodal outflows are adjusted for all nodes. The base demand is subtracted from the adjusted demand and the leakage is defined. Based on this research, one logger is needed for each 200 houses in this method. To cope with this problem, the next investigations are focused on pressure-dependent analysis for leak detection (Giustolisi *et al.* 2008; Wu *et al.* 2010). In this case, the nodal demands can be estimated using demand patterns and the leakage value can be computed based on the orifice area and nodal pressure by the use of specific equations. In this method, a SCADA system provide a large number of observations using a few loggers while the pressure dependent analysis prevents increasing the number of unknowns. Todini (1999) developed a Kalman filter approach by converting the WDN model formulation to a linear estimation problem (Todini 1999). Based on his work, in a looped system, the network observability cannot be guaranteed using a unique set of steady-state data even if all the nodal heads and demands are measured.

Often there are a few leaks in WDS and many of the nodes there do not leak. Using this fact can minimise the space search. Wu *et al.* (2010) assumed that the number of leakages are limited to a specified value (N). Nasirian *et al.* (2011) limited the space search using a novel algorithm (Nasirian *et al.* 2011).

Most previous works in network calibration and leakage detection used the Genetic Algorithm (GA) (Savic & Walters 1995; Sage & Wu 2006; Jamasb *et al.* 2008; Behzadian *et al.* 2009; Pérez *et al.* 2011). GA is a revolution base algorithm that used previous answer results to find the correct direction for the next generations. Calibration problem in steady state is inherently an ill-posed problem (Giustolisi & Berardi 2010; Ostfeld *et al.* 2012) and because of low sensitivity of the nodal pressure to demand variation, the results in the first analysis with relatively high fitness cannot guide to the correct direction. The Wu *et al.* (2010) assumption could help to get better answers.

ACO have the especial parameter named ‘pheromone’ which use previous answers. Additionally, it can be guided by another factor named ‘heuristic guidance’. This parameter can be used by the analyzer to lead to correct answers. Using this parameter to apply the number of potential leaks can decrease the time of program execution and accuracy of the results. Measured pressures and flows in WDS are different with models. Some of this discrepancy may be caused by leakage. Other sources of this error include measurement errors, modeling errors and search space discretization errors. Pressure, flow rate, pipe diameter, pipe material and node elevation are examples of measurement errors. Some assumptions in WDS modeling enter errors in modeling. It is assumed that consumption and demand happen in nodes. Minor loss in the WDS component applies in the friction coefficient.

Most of the search methods used for WDS separate the search space. When the solution space is discretized, it is possible that the exact answer does not exist in the answer's set. The effects of these errors appear in the measured pressure and flow. Therefore, it is difficult to detect a small leak in WDS. On the other hand, the inverse steady state method can detect and locate areas with high leakage. Also, it is important to investigate the effect of measurement errors on accuracy of results. When the relation between two parameters is determined, it is simple to estimate the error of one parameter due to the error of the other parameter. Most of the sensitivity analysis methods are developed for even determined state (Kang *et al.* 2009; Kang & Lansey 2010, 2011; Cheng & He 2011; Cheng *et al.* 2014). This paper presents a method to investigate the existence of unique answers in different cases. However, applying the method in a real situation is the best way to evaluate the model. Then, a laboratory scale network can clarify the ability of the method in a real situation. It provides an opportunity to check the location and value of detected leakages with real leaks. It is clear that a real network is more complicated than a laboratory one; however, in a real network it is not possible to check the results exactly. Up to now, laboratory networks have been tested for network calibration; for example, Walski *et al.* (2006, 2008) or in leakage detection using transient analysis; for example, Covas & Ramos (2010) and Lay-Ekuakille *et al.* (2009) and there is no laboratory examination for leakage detection in steady state analysis.

The method is examined for two networks. The first network is considered as a hypothetical network selected from well-known papers in the literature and the second network comprises the laboratory-scale network built for this work. Additionally, a novel sensitivity analysis technique is employed to exactly evaluate the relationship between the amount and the location of the leakage and the number and accuracy of sensors used for the pressure measurements. An in-house computer program in MATLAB has been written for the optimization part. This computer program is coupled with EPANET 2 (Rossman 2000) to compute the discharges within the pipes and the pressure heads for the nodes.

## METHODS

Leakage detection using network calibration is investigated with some assumptions: (1) a roughly calibrated network is provided before demand calibration, (2) a few simultaneous leakages occur in the network, (3) base demands are detected for all nodes, and (4) uncertainty analysis is performed in the presence of one leakage in the network.

### Calibration problem

The WDN calibration is generally used to determine the different model parameters needed to provide a reasonable accord between the measured and calculated pressures and/or flows in a network. Parameters of calibration usually comprise nodal demand and pipe roughness coefficient. Nodal demand includes nodal based demand, hourly and daily multipliers, and nodal leakage value. Nodal base demand can be recorded by customer consumption. In this work, the unknown parameters are nodal leakages (). The formulation of the objective function of the WDN hydraulic model is defined based on a weighted least squares problem.

**x**is the vector of leakages at various nodes, is the fitness function that should be minimized, and are the observed and simulated total heads at node j, respectively, N is the number of nodes in the network and and are the upper and lower bounds of the possible nodal leakages. Usually, the lower bound is zero, which means no leakage happens, and the upper bound can be estimated by the analyzer. Maximum of can be the total leakage in the network or DMA. If the total inflow to the network is defined, by subtracting the total customer consumption is detected as a constraint. It is expected that if the calculated fitness for an answer (a set of demands for all nodes) is smaller than the calculated one for another answer, the detected leakage will be closer to the other. To evaluate and compare various results obtained in the optimization process, the sum of absolute differences between the real and the simulated leakages is calculated. It is divided by the total inflow in the network to be obtained as a unitless parameter:where is the sum of absolute difference between the simulated and observed leakage in iteration t, is the real leakage value at node j, is the simulated leakage in iteration t at node j, N is the node number in the network and is the total inflow in the network. If the leakage quantities and locations are detected exactly, is zero. If the leakage locations and/or quantities are not exact, it takes a nonzero value. Equation (4) can be evaluated when is defined.

### Uncertainty analysis

Measurement devices are not perfect and the field measurements of pressures and flows in the WDN suffer from some uncertainties. Measurement errors are propagated to the determined values of the calibration parameters. Because of the fact that WDN simulation results depend on calibration parameters, error in the measured values influences the model predictions dominantly. In the previous research works usually the First-Order Second Momentum (FOSM) estimation method calculated based on the Taylor series expansion is employed for uncertainties prediction (Kapelan *et al.* 2003, 2005; Behzadian *et al.* 2009; Jung *et al.* 2014, 2016). This method can be used for the even-determined problems. A particular type of uncertainty predictor for under-determined problems is introduced in this paper. For an even-determined problem, uncertainty in any unknown parameter can be obtained as a scalar value based on the measurement errors but if the problem has some freedom degrees, it has no unique answer and a solution space would be obtained. When measurement error is added to the problem, the problem will be more complicated.

*et al.*2001; Kapelan

*et al.*2003, 2005; Behzadian

*et al.*2009):where is the calculated pressure with respect to the change in the demand value , shows the pressure calculation with an assumed demand , and denotes the deviation value added to . The Jacobian matrix elements can be computed based on the following steps: (1) the values of the base demand are assumed and the hydraulic model is built, (2) the hydraulic model is simulated by adding some perturbed values to the base demand for , (3) the derivatives are computed by applying Equation (6), which provides the first row of the Jacobian matrix; (4) steps 2 and 3 are repeated for (Behzadian

*et al.*2009).

The Jacobian matrix can be inversed if the matrix is squared; in other words, the number of measurements should be equal to the number of unknown parameters. For the systems posed for the WDN, the number of equations is less than the number of unknown parameters, which leads to under-determined systems. To cope with this difficulty, some researchers use pipes and nodes grouping to reduce the number of unknowns and/or run the studied network for the different load conditions (Kapelan *et al.* 2005; Behzadian *et al.* 2009).

*i*due to the leakage at node

*j*. and are the upper and lower bounds of the variation of the pressures respectively due to measurement errors, ( ). The accuracy of pressure measurement depends on the accuracy of sensors and the precision topographic mapping and for each WDS should be estimated by analyzer. is one row of Jacobian matrix that shows the relationship between variation in pressure of node

*i*due to changing in the demand of node () where N is the node number in the network. is the vector of demand changing because of the difference in pressure of node

*i*and are the perturbation in upper and lower bounds due to the existence of error with the amount of in pressure measurement.

*k*is limited. Equation (9) computes this bound.where are the upper and lower bounds for the demand variation in the junction

*k*because of error in the pressure measurement at nodes , respectively. H is the Heaviside step function that is zero when the value is negative or zero, and one if it is positive (Abramowitz 1972). The space solution for the leakage is located between upper and lower bound of demand perturbation. For example, if the demand perturbation is caused by a leakage, the interval between can be the correct options as the leakage value and it can cause the observed pressure in the measurement points.

### Optimization methodology

The optimization methodologies can be classified into three different categories: mathematical, statistical and stochastic methods. Among the stochastic methods, heuristic and meta-heuristics methods are used in several fields of science as well as problems related to the WDN (Sage & Wu 2006; Kapelan *et al.* 2007; Behzadian *et al.* 2009; Rezaei *et al.* 2014; Yazdi *et al.* 2014). The ACO, as one of the heuristic methods, is introduced in 1992 by Dorigo (1992). Simpson et al*.* (2001) discussed the parameters used in the ACO algorithm for the network optimization problems. More recently, Maier *et al.* (2003) developed an ACO method for an optimal design in two case studies and compared the results with GA. Their results indicated more preference of ACO than GA in two benchmarks of WDN. Afshar (2007, 2010) proposed a new transition rule for ACO algorithms using elitist strategies and applied the method to network optimization problems. Nasirian *et al.* (2013) compared ACO with Darwin calibrator WaterGems3 in leakage detection. Their results have shown that ACO provides more accurate results with rapid convergence as well as significant improvement in terms of time and accuracy when compared to the GA.

In this paper, calibration of demands is carried out in the following steps: (1) the hydraulic model of the network is simulated in EPANET; (2) the interval between minimum () and maximum () that is the solution space for unknown leakage as a continuous region is discretized in a number of with a number of . On the other hand, is the number of options that the leakage value gets in each node. In other words, for each node is selected from the set of and . This set is fixed for all nodes and the base demand for each node is added to the set to obtain a set of nodal demands . This procedure will be performed for each node. Then ACO is implemented to select one of as the correct demand at each node, (3) ACO parameters are adjusted, (4) the program is run. In this step, each artificial ant selects a demand for each node. Pressures of all nodes using these demands are computed by EPANET, and then the fitness value is calculated using Equation (2). (5) Pheromone for different paths based on these fitnesses is updated. It means that the pheromone of the demand option of each node that is placed in a patch with better fitness is increased, while the pheromone of the other demand options in that node are decreased because of evaporation parameters. (6) Steps 4 and 5 are repeated until the interest fitness is achieved and (7) adjusted demands are manually compared with the base demands to identify leakage locations and values.

*k*as demand at node

*j*in the iteration. is the amount of pheromone for the nodal demand at node in the iteration, which is updated after each iteration, is the heuristic factor favoring option called ‘heuristic guidance’, which is a constant value allocated based on the researcher estimation for demand option at node

*j*. For instance, a few leakage numbers usually occur in a network simultaneously (Wu

*et al.*2010) and most nodes show no leakage; then, the selection possibility of the no leakage option can be multiplied.

*α*and

*β*are two parameters that regulate the comparative weight of pheromone series and heuristic guidance, respectively.

*t*and

*t*

*+*1, is the updating value which is added to the pheromone value of the local best answer, is the evaporation rate of the pheromone, which is ranged between zero and one and decreases the pheromone value of all solution space. These two parameters take constant values defined by trial and error by the analyzer. In the current paper, and are 0.1 and 0.98, respectively. It should be noted that pheromones for all demand options in every node are chosen as equal to one in the first analysis. Then the ACO is run a few times with this regulation to acquire a number of fitnesses. Next, the obtained results are sorted. Finally, is added to the selected value of demand, causing the best achieving fitness; meanwhile, for other demand options in each node, the pheromones are fixed. However, because of the evaporation function, the corresponding pheromones are diminished. In other words, in each pheromone updating, the pheromone of answer with good fitness is increased while the pheromone of other answers is decreased.

## RESULTS AND DISCUSSION

In this study, as the first and second networks are benchmark and experimental, the consumption pattern is not considered and the analysis was carried out in steady condition. It should be noticed that this method is extensible for EPS condition.

### Case 1

To examine the reliability of the calibration method, the Anytown network, which is known as a small network, is created to investigate the leakage detection based on ACO and uncertainty methodology (Figure 1). It should be noted that this network was widely used in the literature to specify the optimal sensor placement (Ferreri *et al.* 1994; Kapelan *et al.* 2003, 2005, 2007; Behzadian *et al.* 2009). The optimal places for installation of the pressure measurement devices are proposed at nodes 90, 120, 110 and 40 (Kapelan *et al.* 2005; Behzadian *et al.* 2009). Also, the optimal number of pressure measurements are suggested as 7 and 6 by Behzadian *et al.* (2009) and Kapelan *et al.* (2005), respectively, when the accuracy of pressure measurement is equal to 0.1 m. Both of them classified the demands and load conditions within five groups. The relative accuracy of calibration with various numbers of pressure measurements for 2, 4 and 8 pressure observations are 0.2, 0.6 and 0.85, respectively. The fitness function in Behzadian *et al.* (2009) and Kapelan *et al.* (2005) is weight mean square error (Kapelan *et al.* 2005; Behzadian *et al.* 2009). The relative accuracy means the accuracy if N nodes pressure are observed relative to the case that all nodes pressures are measured (Kapelan *et al.* 2003). It should be noted that attaining the accuracy of 0.1 m is difficult in some aspects in real networks. A pressure sensor has about 1% error in the best condition. Topography map error is added to the pressure sensor error as well. The HGL^{1} almost has 0.1–1 m accuracy. In this benchmark case, the assumptions of previous studies are implemented. The other error which is added to the measurements is the inflow parameter. The flowmeters, which measure the inflow rate, usually have an accuracy about 0.05–0.1 l/s. In this case, 0.05 l/s was assumed for the flowmeter.

The hypothetical network has 16 nodes, 34 pipes, 3 pumps as well as 2 tanks and 1 reservoir with the elevation of 71.6 and 3.04 m. Length, diameter and Hazen-Williams (H-W) coefficient are denoted above each pipe and elevation and consumption of each node is represented in a box next to the node. The other properties can be found in mentioned references. It should be mentioned that the hydraulic analysis is carried out in EPANET 2. and the Hazen-Williams formula is used in the analysis.

In this research, the network analysis is carried out without demand grouping and also no-load condition (demand grouping decreases the unknown parameters. Using several load conditions increases the number of observations). The demands are ranged from zero to 50 l/s with an increment of 10 l/s and the maximum iteration is 1,000. The heuristic guidance values based on statistics and trial and error are selected. They are for one and two leakages considered as 30 and 20, respectively. Table 1 shows 8 different scenarios for calibration.

1 . | 2 . | 3 . | 4 . | 5 . | 6 . | 7 . | 8 . | ||
---|---|---|---|---|---|---|---|---|---|

Real leakage . | Calculated leak . | ||||||||

Analysis case . | Measurement error (m) . | Position . | Value (l/s) . | Observation . | Position . | Value (l/s) . | Fitness (m^{2})
. | Iteration . | Time (sec) . |

Scenario 1 | 0 | 70 | 30 | P90 inflow | 20 40 70 | 10 10 10 | 2.9E-06 | 43 | 6 |

Scenario 2 | 0 | 70 | 30 | P90 P120 | 70 | 30 | 1.8e-005 | 532 | 86 |

Scenario 3 | 0 | 70 140 | 30 20 | P90 P120 | 110 | 20 | 0.007547 | 157 | 35 |

Scenario 4 | 0 | 70 140 | 30 20 | P90 P120 P110 inflow | 70 140 | 30 20 | 3.46E-05 | 63 | 26 |

Scenario 5 | 0.1 | 70 | 30 | P90 P120 | 70 | 20 | 0.004699 | 92 | 42 |

Scenario 6 | 0.1 | 70 140 | 30 20 | P90 P120 P110 inflow | 100 140 160 | 20 20 10 | 1.42E-02 | 871 | 181 |

Scenario 7 | 0.1 | 70 140 | 30 20 | P90 P120 P110 P40 P170 inflow | 60 70 140 | 10 20 20 | 0.001969 | 914 | 325 |

Scenario 8 | 0.1 | 70 140 | 60 40 | P90 P120 P110 inflow | 160 70 140 | 10 50 40 | 0.00087 | 694 | 289 |

1 . | 2 . | 3 . | 4 . | 5 . | 6 . | 7 . | 8 . | ||
---|---|---|---|---|---|---|---|---|---|

Real leakage . | Calculated leak . | ||||||||

Analysis case . | Measurement error (m) . | Position . | Value (l/s) . | Observation . | Position . | Value (l/s) . | Fitness (m^{2})
. | Iteration . | Time (sec) . |

Scenario 1 | 0 | 70 | 30 | P90 inflow | 20 40 70 | 10 10 10 | 2.9E-06 | 43 | 6 |

Scenario 2 | 0 | 70 | 30 | P90 P120 | 70 | 30 | 1.8e-005 | 532 | 86 |

Scenario 3 | 0 | 70 140 | 30 20 | P90 P120 | 110 | 20 | 0.007547 | 157 | 35 |

Scenario 4 | 0 | 70 140 | 30 20 | P90 P120 P110 inflow | 70 140 | 30 20 | 3.46E-05 | 63 | 26 |

Scenario 5 | 0.1 | 70 | 30 | P90 P120 | 70 | 20 | 0.004699 | 92 | 42 |

Scenario 6 | 0.1 | 70 140 | 30 20 | P90 P120 P110 inflow | 100 140 160 | 20 20 10 | 1.42E-02 | 871 | 181 |

Scenario 7 | 0.1 | 70 140 | 30 20 | P90 P120 P110 P40 P170 inflow | 60 70 140 | 10 20 20 | 0.001969 | 914 | 325 |

Scenario 8 | 0.1 | 70 140 | 60 40 | P90 P120 P110 inflow | 160 70 140 | 10 50 40 | 0.00087 | 694 | 289 |

The network is analyzed for different scenarios applying the pressure measurements and total inflow to the network. The pressure measurements are needed to calculate the fitness function, while the total inflow is applied as a physical constraint. Table 1 presents the results of 8 calibration scenarios. According to Table 1 column 2, nodal pressures are measured 100% accurately in the first 4 analyses. However, the pressure measurement perturbations are assumed as 0.1 m in the next 4 analyses. Also, the next two columns represent the locations and values of leakages located at nodes 70 and 140 in various scenarios. The next five columns show the observational states of pressures and flow, the locations and sizes of the leakages obtained in calibration, the fitness value and finally the total iterations, respectively. The final column presents the run-time for each scenario.

Based on Table 1, it is clearly observed that the calibrations in scenarios 1 and 3 that tried to find one and two leakages using one and two pressure measurements do not attain the correct results. However, scenario 2, which used 2 pressure measurements for finding one leakage has reached perfect answers. Also, scenario 4, which utilized 3 pressure and inflow measurements, attained the correct answer. The observations in scenario 5 are subjected to some errors. In this case, although the leakage location is found exactly, some errors can be observed in the quantity of leakage. Two simultaneous leakages are investigated in scenario 6. Although an excellent fitness is achieved, the consequence demonstrates significant variation. Likewise, scenario 7 can lead to successful outcomes with additional observations. The last analysis is performed by duplicating the quantities of leakages in scenario 8. The results of this case show a dramatic improvement in the leakage detection in comparison to scenario 7. To investigate the relationship between the number of observations, measurement error and the accuracy of the obtained results in the different scenarios, the uncertainty analysis is carried out.

In order to clarify the behavior of the proposed method in detecting the leakage, the results of all scenarios are investigated. Based on the results of scenario 4, which correspond to the exact solution and also its number of iterations, this scenario is chosen and its results are presented in Figure 2. As shown in this figure, the difference between the leakages in the real and calculated answers for 63 checking answers is plotted against the fitness .

This figure shows a positive correlation between and . By decreasing , is also decreased, which shows a kind of improvement to the real answer. Most of the optimization approaches use fitness as a criterion to check the correct trend of progress towards the final best results. When the correlation between and is weak, the optimization methods will not be able to find the correct answer. In ACO, the pheromone updating is applied to the fitness of previous answers to select a new direction towards the correct answer. This property is similar to the evolutionary based optimization methods. However, ACO is equipped with another device named heuristic guidance that helps in finding the correct direction towards the real answer. In this paper, leaks are defined at a few unidentified nodes. Then, naturally, the existence of leakage probability for each node is decreased, which is fulfilled by the heuristic guidance. By the use of this technique, it will be possible to find the correct answer among 19,344 possible solutions by performing only 63 iterations in this scenario.

### Uncertainty analysis

It is assumed that one pressure sensor is placed at node 90 and one leakage with the amount of 30 l/s is introduced at node 70. The magnitude of pressure decrement at node 90 due to this leakage is 0.33 m. This reduction may happen from some other leakages in other nodes. Figure 3 shows some other results that only with one leaking node can reduce the pressure of node 90 to an amount of 0.33 m. For instance, a flow rate of 10 l/s at node 90 can cause a similar reduction. This means there are several choices to set the true value in the pressure of this node.

The leakage in this case is about 7.5 percent of total consumption. If the leakage is decreased to 10 l/s, the pressure reduction in the observed node is 0.11 m. With pressure measurement errors, this method will be inaccurate for small leaks.

With some techniques it is possible to attain better results. For example, in the minimum night flow the consumption is significantly low and both pressure and leakage over consumption are high. To ascertain better results, it is proposed to stop the pressure management in the measurement night. Also, this method can help to detect the zone with high leakage or pipe burst. Using EPS data can increase the accuracy of the methods.

To attain the unique answer, the observational pressure at node 120 is added to the pressure measurements set. When the pressure at nodes 120 and 90 has been measured accurately, the only choice which can adjust the pressures at two junctions is the leakage at node 70 of 30 l/s (see Figure 4 and Table 2).

In the Anytown network, the number of unknowns were considered 16 node consumption and the number of pressure observations considered were among 2–6.

The common method in real WDNs, which is used extensively in studies, is classifying the consumption in nodes. For instance, the network could be divided into 4 or 5 parts. With this technique, the number of unknowns decrease; however, the accuracy of finding the location reduces.

When the measurements are associated with some uncertainties, the points shown in Figure 3 are depicted by vertical lines. Figure 5 shows the leakage range in various nodes that cause 0.33 m reduction in the pressure of node 90. Also, Figure 6 has shown the different leakage values that possibility caused the pressure with perturbation amount at nodes 90 and 120. These two charts indicate that the 16 possible nodes in Figure 4 are limited to 9 nodes in Figure 5 that can be considered as the leakage location. By adding the total inflow measurement as a constraint, the candidate nodes are limited to 3 nodes. Likewise, the pressures of nodes 110 and 40 are joined to the observational set. Figure 7 has clearly shown the analysis has reached an exact value. In Figures 6 and 7, the flowmeter error (0.05 l/s) is shown.

Following issues can be concluded from these investigations. (1) It is clearly shown that the locations and values of leakages can be obtained by a steady state analysis and an optimization method. (2) ACO can find optimum answers of a problem with a few iterations. It should be noted that in this example, if the inflow is not used as observation (e.g. scenario 2, 3 and 5) there are 6^{16} = 2,821,109,907,456 different options and if the total leakage is 50 l/s and used as a constraint, there are 19,344 various cases. Different scenarios can find the correct answer in 14 to 914 iterations. In the worst case, about 5% of the solution space is checked to find the exact answer. (3) The leakage detection occasionally has improved with increasing in the leakage amount, number of observations and their accuracies. (4) Reaching the supreme fitness cannot guarantee either location or quantity of leakages; hence, an optimal sensor placement analysis should be accomplished for each network to understand the capability of the number and quality of the observational set in a WDN as well as the calibration methodology. In other words, the problem is to find the leakage with the available equipment. (5) In the previous researches (Kapelan *et al.* 2003, 2005; Behzadian *et al.* 2009) it has been emphasized that for a relative accuracy of about 0.85 in the calibration of a network, the pressures at 8 nodes among 16 nodes should be measured. Five load conditions should be implemented and the number of unknown demands should be limited to 5 groups. In the other word, to find only 5 unknown parameters with relative accuracy of 0.85, 40 observations are needed. Therefore, it seems that by the use of this method, it is almost impossible to detect the leakage in a real network.

Therefore, using this method seems almost impossible for a real network. However, in this research it is clearly shown that with a few measurements, a suitable result can be gained while neither node/pipe grouping nor load conditions are used. Actually, by implementing these two facilities better results can be achieved. As a matter of fact, it should be said that by the use of these techniques results with higher accuracy can be achieved.

### Case 2: experimental model

To investigate the above-mentioned technique, a laboratory model was set up at the hydraulic laboratory of Ferdowsi University of Mashhad, Iran, for evaluation of the explained method. Figure 8 shows some pictures and Figure 9 indicates the EPANET model of the network. The height and width of the network are 3.70 m and 10.30 m, respectively consisting of eight consumption taps at nodes 1 to 8 for demand simulation and 31 pressure measurement situations within six loops. Moreover, as shown in Figure 8, to calculate the water in the network, water returning pipes, a pump and a storage tank are used. Additionally, a volumetric measurement tank is used to measure the demand flow rate. In this model, electronic pressure transducers with an accuracy of 1% in predetermined spots were used to measure the pressure and then data are transferred to a PC. The pipes were made of Poly Propylene Random Co-polymer (PPRC) with an inner diameter of 16 and 10 mm for loops and returning pipes. A full port ball valve is installed in each loop of pipe to change the configuration of the network. Additionally, an electromagnetic flowmeter was installed in the network entrance to measure the total inflow of the network with an accuracy of 0.03 l/s. As can be seen in Figure 8(d) we have 8 valves. Each valve receives two pipes which end in a reservoir. Near the reservoir we have scaled storage that help us to measure the flowrate.

The pump that is used in the WDN is chosen in a way to respond to the necessary flows and pressures in the network. The power, maximum head and flow of the pump are 4 HP, 52 m and 5 l/s respectively. Also, the pump characteristic curve is shown in Figure 10.

### Calibration of the network

The network is modeled in EPANET with the specification of exact pipes, pump characteristic curve, nodal elevations, friction factors and minor losses of the pipes. The hydraulic head loss and minor loss coefficients are determined through several experimental runs. Since in this paper only the demand is considered as an unknown parameter, other parameters of the model should be calibrated first. The pipe roughness is the most important uncertain parameter in a network. Nasirian *et al.* (2011) have considered the best equation to estimate the longitudinal loss as well as the pipe roughness and minor loss coefficients in the experimental network. H-W coefficient (Walski *et al.* 2002) is calculated for about 150 pipes (Nasirian *et al.* 2011), which are confirmed by Braud & Soom (1981) and Wu & Gitlin (1984). Additionally, they have proved the appropriate accuracy of Darcy-Weisbach (D-W) and H-W's equations in this network, and then the H-W equation is used in the analysis. The other testing shows that by measuring the pressure at the pump and reservoir, the H-W coefficient was measured at about 130 (Nasirian *et al.* 2011).

The head loss between pump and reservoir is calculated using the H-W and D-W equation, and also it is measured experimentally as shown in Figure 11.

Figure 12 shows the head difference between the Darcy-Weisbach method with Colebrok-White roughness in comparison with other methods. It can be attributed to the effect of minor head loss. A good agreement between head loss accounting for minor head loss and real head loss can be observed. When the Hazen-Williams coefficient is 130, the difference between observed hydraulic head loss and the computed head loss relates to the minor head losses. The minor head loss for the whole path is computed and it is 8.6 m. This value is the summation of all minor head losses in each connection and valve that exist in the path. With calculating the minor head loss and the linear head loss, the total head loss is derived. This value corresponds to the laboratory result. Therefore, ignoring the minor head losses leads to significant errors. So, all the minor head losses must be computed separately. Typically, minor losses and friction head loss are related to the velocity with different exponents. The flow velocity doesn't exceed 1.5 l/s. It has a value lower than 1.5 l/s. Using an equivalent pipe friction factor instead of head loss coefficient causes an error that for an urban network with long pipes is negligible, but in a laboratory scale model, it causes a significant error. To prevent propagating this error in the demand calibration, the minor loss coefficient for each pipe is experimentally examined. Measured pipes roughness and minor loss coefficients are used in the model and then the pressures are calculated and compared with the observational pressures.

### One leakage in the network

In this case, a leakage of 0.262 l/s is introduced at node 9 while the main pipes are opened. Figure 9 shows the leakage location. Also, the location of pressure sensors used in various scenarios has been shown as hollow circles. Nodal pressures of nodes 17 and 24 are measured. The demand will be searched within a range of a minimum of zero to a maximum of 1 l/s with an incremental step of 0.05 l/s in all cases.

Table 3 shows the computational and observational pressures in the considered nodes. The program reaches the true answer in the 102th iteration. The results show that the leakage with a quantity of 0.25 l/s happens at node 9, which is the best possible answer.

Node No. . | Observation pressure (m) . | ACO results . | |
---|---|---|---|

Calculated demand (l/s) . | Simulated pressure (m) . | ||

9 | – | 0.25 | 24.35 |

17 | 29.96 | 0 | 30.30 |

24 | 28.72 | 0 | 29.07 |

Node No. . | Observation pressure (m) . | ACO results . | |
---|---|---|---|

Calculated demand (l/s) . | Simulated pressure (m) . | ||

9 | – | 0.25 | 24.35 |

17 | 29.96 | 0 | 30.30 |

24 | 28.72 | 0 | 29.07 |

In this case, an uncertainty analysis is accomplished. Figure 13 represents the solution space for various pressure observations. Figure 13(a) indicates the solution space when the pressure at node 24 is measured. Figure 13(b) show the solution space if the pressure at node 17 in measured. Figure 13(c) shows the answer spaces when the pressures of these two nodes are measured simultaneously. Figure 13(a) and 13(b) illustrate 3 and 4 answers for pressure measurements at nodes 24 and 17, respectively, while the solution space using 2 pressure measurements include 2 answers. In other words, 0.3 l/s of leakages at nodes 5 and 9 can be potentially found as the answer, while as mentioned above the leakage location is found in the calibration process at node 9. This analysis clarified that, although the calibration results in Table 3 have reached the true results, another potential answer exists. Figure 13(d) presents the solution space after adding the pressures at nodes 24, 17 and 25 to previous observations. However, no unique value is observed. It should be noted that more analyses with various observational cases clarified that reaching a unique answer in this case is impossible. Nevertheless, the results are considered very satisfactory and promising because the solution space is closed to the real leaky point.

### Two leakages in a sub-network

In the other case, the network is reconfigured by closing the valves in pipes 8, 12 16, 19, and 22. In other words, the middle pipes are closed and only the peripheral pipes are opened. The installed taps at nodes 5 and 6 are opened with the measured leakages of 0.09 and 0.33 l/s, respectively. Also, the pressures of nodes 17, 24 and 27 are measured. The answers are found at the 195th iteration and the leakages at nodes 5 and 6 are computed as 0.1 and 0.3 l/s. The final results are shown in the Table 4.

Parameter . | Node No. . | Observation . | ACO results . |
---|---|---|---|

Demand (l/s) | 5 | 0.09 | 0.1 |

6 | 0.33 | 0.3 | |

Pressure (m) | 27 | 23.96 | 24.41 |

24 | 22.65 | 22.98 | |

17 | 25.35 | 25.87 |

Parameter . | Node No. . | Observation . | ACO results . |
---|---|---|---|

Demand (l/s) | 5 | 0.09 | 0.1 |

6 | 0.33 | 0.3 | |

Pressure (m) | 27 | 23.96 | 24.41 |

24 | 22.65 | 22.98 | |

17 | 25.35 | 25.87 |

### Two leakages in main network

As the final case, leakages at nodes 8 and 6 with the amount of 0.37, 0.526 l/s are introduced to the nodal while all of the pipes in the network are opened. The pressures at five nodes 20, 33, 25, 14 and 45 are measured. The network is calibrated using two observational cases. The pressures of nodes, 33, 20 and 25 are used in the first case, while all of the five pressure measurements are applied in the second case. The leakages are found in the range of 0 to 0.7 l/s and incremental steps are 0.05 l/s.

To evaluate the errors in the measured data, nodal pressures are calculated using EPANET model applying measured flows in the leaky nodes. A comparison between the observational and computational pressures shows an average discrepancy of about 0.41 m. Therefore, if calibration attains the exact value of leakages, the obtained fitness will not be less than 0.98. Table 5 presents the measured pressures, EPANET results as well as, the ACO optimization results. Table 5 column 5 shows the results using 3 pressure observations. Although the obtained fitness, 0.0185, is great, the leakages are detected incorrectly and only one of the leakages is partly defined. The final examination using five pressure observations is presented in column 6 of Table 5. The results indicate reasonable accuracies in the obtained results.

1 . | 2 . | 3 . | 4 . | 5 . | 6 . |
---|---|---|---|---|---|

Parameter . | Node No. . | Observation . | EPANET results . | ACO results with 3 observations . | ACO results with 5 observations . |

Demand (l/s) | 3 | 0 | 0 | 0 | 0 |

4 | 0 | 0 | 0 | 0 | |

8 | 0.37 | 0.37 | 0 | 0.25 | |

6 | 0.526 | 0.526 | 0.4 | 0.55 | |

9 | 0 | 0 | 0.45 | 0.1 | |

7 | 0 | 0 | 0.05 | 0 | |

5 | 0 | 0 | 0 | 0 | |

10 | 0 | 0 | 0 | 0 | |

Pressure (m) | 20 | 34.9 | 35.54 | 34.82 | 35.19 |

33 | 32.4 | 32.84 | 32.51 | 32.63 | |

25 | 34.4 | 34.93 | 34.4 | 34.52 | |

14 | 37.8 | 37.54 | – | 37.45 | |

45 | 36.1 | 35.92 | – | 35.54 | |

Fitness | – | – | 0.98 | 0.0185 | 0.587 |

1 . | 2 . | 3 . | 4 . | 5 . | 6 . |
---|---|---|---|---|---|

Parameter . | Node No. . | Observation . | EPANET results . | ACO results with 3 observations . | ACO results with 5 observations . |

Demand (l/s) | 3 | 0 | 0 | 0 | 0 |

4 | 0 | 0 | 0 | 0 | |

8 | 0.37 | 0.37 | 0 | 0.25 | |

6 | 0.526 | 0.526 | 0.4 | 0.55 | |

9 | 0 | 0 | 0.45 | 0.1 | |

7 | 0 | 0 | 0.05 | 0 | |

5 | 0 | 0 | 0 | 0 | |

10 | 0 | 0 | 0 | 0 | |

Pressure (m) | 20 | 34.9 | 35.54 | 34.82 | 35.19 |

33 | 32.4 | 32.84 | 32.51 | 32.63 | |

25 | 34.4 | 34.93 | 34.4 | 34.52 | |

14 | 37.8 | 37.54 | – | 37.45 | |

45 | 36.1 | 35.92 | – | 35.54 | |

Fitness | – | – | 0.98 | 0.0185 | 0.587 |

Due to the fact that the best expected fitness is 0.98, there are several other answers that have shown relatively good fitnesses. Table 6 displays the best 10 calibrated results for the first observational case. The first column expresses the best results in this case. Meanwhile, the other columns show that a number of results are approaching the best one. The best answer shows the leakages at nodes 6, 7 and 9 while the second answer shows leakages at two different nodes 4 and 10 accompanied with a significant reduction in the leakage quantities at node 9. Likewise, the third answer is completely different from the best answer. In other words, the results have shown no convergence or correlation between various results attaining the final solution.

Fitness . | 0.019 . | 0.030 . | 0.098 . | 0.099 . | 0.174 . | 0.432 . | 0.541 . | 1.175 . | 1.533 . | 1.857 . |
---|---|---|---|---|---|---|---|---|---|---|

Node No. . | ||||||||||

9 | 0.45 | 0.15 | 0.55 | 0.6 | 0 | 0 | 0 | 0 | 0 | 0 |

8 | 0 | 0 | 0 | 0 | 0.15 | 0 | 0.45 | 0 | 0 | 0.55 |

7 | 0.05 | 0.05 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |

3 | 0 | 0 | 0 | 0 | 0 | 0.4 | 0.55 | 0 | 0.15 | 0 |

4 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |

5 | 0 | 0.2 | 0 | 0 | 0.1 | 0 | 0 | 0 | 0.7 | 0.4 |

6 | 0.4 | 0.4 | 0 | 0 | 0.65 | 0.55 | 0 | 0.5 | 0.1 | 0 |

Fitness . | 0.019 . | 0.030 . | 0.098 . | 0.099 . | 0.174 . | 0.432 . | 0.541 . | 1.175 . | 1.533 . | 1.857 . |
---|---|---|---|---|---|---|---|---|---|---|

Node No. . | ||||||||||

9 | 0.45 | 0.15 | 0.55 | 0.6 | 0 | 0 | 0 | 0 | 0 | 0 |

8 | 0 | 0 | 0 | 0 | 0.15 | 0 | 0.45 | 0 | 0 | 0.55 |

7 | 0.05 | 0.05 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |

3 | 0 | 0 | 0 | 0 | 0 | 0.4 | 0.55 | 0 | 0.15 | 0 |

4 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |

5 | 0 | 0.2 | 0 | 0 | 0.1 | 0 | 0 | 0 | 0.7 | 0.4 |

6 | 0.4 | 0.4 | 0 | 0 | 0.65 | 0.55 | 0 | 0.5 | 0.1 | 0 |

Table 7 shows a number of answers of the second case using five pressure measurements. It is clearly observed that the calculated answer has greatly approached the final results. The results elucidate a good convergence between the three final answers.

Node No. . | 0.59 . | 0.61 . | 0.61 . | 0.73 . | 1.04 . | 1.08 . | 1.88 . | 2.32 . | 2.69 . | 2.69 . |
---|---|---|---|---|---|---|---|---|---|---|

10 | 0 | 0 | 0 | 0 | 0.35 | 0 | 0.5 | 0 | 0 | 0 |

9 | 0.1 | 0.05 | 0 | 0 | 0.15 | 0 | 0.15 | 0 | 0.35 | 0.35 |

8 | 0.25 | 0.3 | 0.3 | 0 | 0.1 | 0.15 | 0 | 0 | 0 | 0 |

7 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |

3 | 0 | 0 | 0 | 0 | 0 | 0 | 0.05 | 0.4 | 0 | 0 |

4 | 0 | 0 | 0 | 0.35 | 0 | 0 | 0 | 0 | 0.55 | 0.55 |

5 | 0 | 0 | 0 | 0 | 0.25 | 0.1 | 0 | 0 | 0 | 0 |

6 | 0.55 | 0.55 | 0.6 | 0.55 | 0.05 | 0.65 | 0.2 | 0.55 | 0 | 0 |

Node No. . | 0.59 . | 0.61 . | 0.61 . | 0.73 . | 1.04 . | 1.08 . | 1.88 . | 2.32 . | 2.69 . | 2.69 . |
---|---|---|---|---|---|---|---|---|---|---|

10 | 0 | 0 | 0 | 0 | 0.35 | 0 | 0.5 | 0 | 0 | 0 |

9 | 0.1 | 0.05 | 0 | 0 | 0.15 | 0 | 0.15 | 0 | 0.35 | 0.35 |

8 | 0.25 | 0.3 | 0.3 | 0 | 0.1 | 0.15 | 0 | 0 | 0 | 0 |

7 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |

3 | 0 | 0 | 0 | 0 | 0 | 0 | 0.05 | 0.4 | 0 | 0 |

4 | 0 | 0 | 0 | 0.35 | 0 | 0 | 0 | 0 | 0.55 | 0.55 |

5 | 0 | 0 | 0 | 0 | 0.25 | 0.1 | 0 | 0 | 0 | 0 |

6 | 0.55 | 0.55 | 0.6 | 0.55 | 0.05 | 0.65 | 0.2 | 0.55 | 0 | 0 |

Three case studies of leakage detections have shown that increasing the number of leakages or complexity of a network increases the difficulty of leakage detection. Also, considering two cases of 3 and 5 observations have determined the leakage at node 6, which is defined in both cases, while leakage at node 8 is only defined after using five observations. This confirms a relationship between the value of leakage and the number of observations. In other words, the requested number of pressure measurements is increased with increasing size and number of leakages as well as the complexity of the network.

Investigating the best 10 results in Tables 6 and 7 has illustrated that the leakage at node 6 is correctly detected in 8 answers as shown in Table 7, while it is obtained only in 5 answers as shown in Table 6. In the other word, although the leakages of node 6 in both observational cases are correctly determined, leakage detection frequency in Table 7 is higher than Table 6. This comparison between leakage detection frequency at nodes 8 and 6 in each of the observational cases confirms a relationship between frequency of leakage detection and reliability of answers.

The experimental investigations in two observational cases have illustrated the fact that a good fitness cannot guarantee attaining a good answer. Based on the results, it is recommended that some good answers are checked out instead of only on the best result. If the other answers with proximate comparatively good fitness present similar leakages, the results are reliable. Otherwise, an uncertainty analysis can help to clarify that with a number of recorded data (obtained from pressure/flow measurement devices) if there is a unique answer set to the problem or the feasible solution space contains a large number of answer sets.

## CONCLUSIONS

In this paper, the leakage detection in a WDN based on calibration and ACO is investigated. Three bases of finding leakage in a network are as the following: gathering network data, hydraulic analysis of the network, and the optimization algorithm. In the current paper, as a main objective it is focused on the improvement of the optimization method. However, in the future works the model will be applied to real networks. It is believed that the most advanced methods in network analysis such as pressure-dependent analysis and EPS should be used. In addition, it is mandatory to use pipe and node groupings in large networks, which make it possible to solve the problem. The method is applied to a hypothetical and a laboratory network without using pipe and/or node grouping for unknown reduction or implementation at different load conditions to increase the observational numbers. So, the analyses are accomplished in a fully under-determined state. Also, a novel approach is used to investigate the uncertainty in an under-determined problem. This method is successfully applied to both networks. The solution space is illustrated for a measurement scenario. Also, the effects of measurement errors on final answers are investigated. Utilization of this method for two networks proves the accuracy of ACO to find a WDN leakage with a few observations. Based on the current investigations, leakage detection in the studied networks can attain appropriate results with a few pressure and flow observations. Also, this method succeeded in recognizing the sizes and locations of the leaks in the laboratory network with actual errors. The experimental and computational investigations in two networks have illustrated that a good fitness cannot guarantee a reliable answer. To cope with this difficulty, it is recommended to check some answers instead of just looking for the best result. If the answers with comparatively good fitness present similar leakages, the results are reliable; otherwise, an uncertainty analysis can be used to understand the feasibility of attaining answers with available pressure/flow measurements. Other investigations have shown that the reliability of the technique has a direct proportionality with the number of correct answers and the total number of answers. The case studies of leakage detections have shown that increasing the number of leakages in a WDN or complicating the network needs a larger number of observations.

## DATA AVAILABILITY STATEMENT

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

Hydraulic Gradient Line.

## REFERENCES

*Optimization, Learning and Natural Algorithms*

*Ph.D. thesis*