Robust Data-Driven Leak Localization in Water Distribution Networks Using Pressure Measurements and Topological Information

Abstract

This article presents a new data-driven method for locating leaks in water distribution networks (WDNs). It is triggered after a leak has been detected in the WDN. The proposed approach is based on the use of inlet pressure and flow measurements, other pressure measurements available at some selected inner nodes of the WDN, and the topological information of the network. A reduced-order model structure is used to calculate non-leak pressure estimations at sensed inner nodes. Residuals are generated using the comparison between these estimations and leak pressure measurements. In a leak scenario, it is possible to determine the relative incidence of a leak in a node by using the network topology and what it means to correlate the probable leaking nodes with the available residual information. Topological information and residual information can be integrated into a likelihood index used to determine the most probable leak node in the WDN at a given instant k or, through applying the Bayes’ rule, in a time horizon. The likelihood index is based on a new incidence factor that considers the most probable path of water from reservoirs to pressure sensors and potential leak nodes. In addition, a pressure sensor validation method based on pressure residuals that allows the detection of sensor faults is proposed.

1. Introduction

Water distribution networks are complex systems that are difficult to manage and monitor with extreme importance nowadays. The detection and location of leaks have become crucial for water distribution because when there are bursts or leaks, this can generate not only economic losses but also an environmental issue and represents a potential risk to public health with contaminated water [1]. Another concern is the scarcity of water that can occur in 2025, which may affect half the world’s population that will not have access to safe and accessible water for their basic needs [2]. However, with all these risks, currently, this infrastructure does not perform satisfactorily in practice. According to [3], a global volume of water loss called Non-Revenue Water (NRW) has been calculated at 346 million cubic meters per day or 126 billion cubic meters per year.

The infrastructure in a medium-sized city can have pipes that span hundreds of kilometers connected to hundreds of nodes (pipe junctions or customers that connect to the network). Therefore, several factors can generate water loss during transport between the treatment plants and the reservoir for consumers, usually attributed to several causes, including leaks, measurement errors, and theft. Water loss can be divided into two terms, “real losses” and “apparent losses”. Apparent losses are constituted by badly read measurements, data handling errors, and illegal water tapping. In contrast, the real losses comprise leakage from all system parts and overflow at storage tanks. Real losses are divided into “background leakage” made up of small undetectable and into detectable leaks relevant for detection as they represent significant losses for the water distribution company.

Effective leak management is vital for all of the factors mentioned above to save financial resources and water. The methods of leak localization can be classified into two categories: Hardware-based system and Software-based.

The Hardware-based utilizes hardware sensors to detect a leak directly and help the localization of the leak. As there are various types of sensors and instruments available, they can be further subclassified as: acoustic [4,5] and non-acoustic detection methods [6].

Software-based methods generally rely on an algorithm or model for detecting leaks. Unlike hardware-based methods, these methods do not seek to locate the leak point accurately but minimize possible leakage areas. Since these methods are based on information, such as the pressure of the pipe network, flow data, and so forth, they work well on any type of pipe. These methods can be divided into physical modeling methods and data-driven methods. The physical modeling methods or model-based methods identify the leak using a numerical model and compare the results with the field data, for example, Ref. [7] which uses pressure sensitivity analysis, Ref. [8] uses leak signature space, Ref. [9] analyzes the sensitivity matrix and residuals, and [10] uses pressure and flow measurements to perform leakage detection through model-invalidation. On the other hand, data-driven methods analyze the monitoring data, combining tools such as artificial intelligence (e.g., classifiers [11,12,13,14] or artificial neural networks [15,16]). Thus, it is possible to identify potential areas of the leak based on certain rules or principles without resorting to the simulation of the physical model results [17]. However, these methods need, in general, an important number of non-leak and leak data scenarios in the training process to obtain reasonable results. As an exhaustive amount of leak scenarios are not available in general, a hydraulic simulator can be used to generate leak data. This work deals with the problem of leak localization and it is assumed that it is available a leak detection method that determines if a leak is present or not in the WDN. In particular, a non-numerical localization method, focused on a data-driven approach, is proposed.

Like other recent works [18,19], it requires only topological information of the network and historical data without leakage of the available measurements. In this work, the topological information provides the most probable paths for extra flows produced by leaks. A new incidence factor from every combination of nodes and sensors is computed with this information. Every incidence factor determines how a leak in a particular node affects a specific pressure sensor. On the other hand, historical data are used to calculate non-leak pressure estimations at sensed inner nodes. Residuals are generated using the comparison between these estimations and leak pressure measurements. Incidence factors are integrated with residuals in likelihood indexes to give the most probable leak node in a leak scenario. In addition, pressure residuals are used to detect sensor faults by means of a novel sensor validation algorithm.

The remainder of this paper is organized as follows: Section 2 presents the theory of graphs applied to WDN and explains the structure of the reduced-order model used in this work. The developed leak localization has been elaborated in Section 3. In Section 4 a sensor validation method that allows the detection of pressure sensor faults is presented. Section 5 introduces the case studies of Hanoi and Modena’s WDNs. Section 6 presents the conclusions and future scope of the research work.

2. Water Distribution Networks

2.1. Preliminaries

A water distribution network is composed of m pipes, n internal consumer nodes and can be described by a directed graph $\mathcal{G}=\{\mathcal{V},\mathcal{E}\}$,[20], with $\mathcal{V}=\{v_1,\dots ,v_n\}$ is the set of vertices that represent connections between the components of the network, additionally the last $\{v_{n-n_I+1},\dots ,v_n\}$, represent the vertices of the system’s input, $n_I$ being the number of the inlets, with $n_I\ge 1$. The elements of the set $\mathcal{E}=\{e_1,\dots ,e_m\}$ are the edges, which represent the m pipes in the network.

The graph $\mathcal{G}$ can be represented by the incidence matrix $\mathbf{H}=\left[h_{ij}\right]$, in which the elements $h_{ij}$ are defined as:

$$h_{ij}=\left\{\begin{array}{ll}-1& \mathrm{if}\mathrm{the}j\mathrm{th}\mathrm{edge}\mathrm{is}\mathrm{entering}i\mathrm{th}\mathrm{vertex}.\\ 0& \mathrm{if}\mathrm{the}j\mathrm{th}\mathrm{edge}\mathrm{is}\mathrm{not}\mathrm{connected}\mathrm{to}\\ & \mathrm{the}i\mathrm{th}\mathrm{vertex}.\\ 1& \mathrm{if}\mathrm{the}j\mathrm{th}\mathrm{edge}\mathrm{is}\mathrm{leaving}i\mathrm{th}\mathrm{vertex}.\end{array}\right.$$

The direction of the edge represents a reference direction for the flow in the corresponding pipe. The incidence matrix is composed of $H\in {\{-1,0,1\}}^{n\times m}$ with each row corresponding to a node and column corresponding to a pipe.

The WDN must fulfill mass conservation law, which expresses the conservation of mass in each vertex, described by:

$$\mathbf{H}·\mathbf{q}=\mathbf{d},$$

(1)

where $\mathbf{d}\in {\mathbb{R}}^n$ is the vector of nodal demands, with ${\mathbf{d}}_i$ > 0 when the flow is into the node i, and $\mathbf{q}\in {\mathbb{R}}^m$ is the vector of flows in the edges. By virtue of the mass conservation, it is possible to have only $n-1$ independent nodal demand, ${\sum }_{i=1}^n{d}_i=0$, therefore the supply flow must equal the end-user demands as there is no storage in the network.

Let $\mathbf{p}$ be the vector of absolute pressures at the nodes and $\Delta \mathbf{p}$ be the vector of differential pressures across the pipes, both in meters of water column $\left[\mathrm{mwc}\right]$, then the energy law for water networks gives:

$$\begin{array}{c}\hfill \Delta \mathbf{p}={\mathbf{H}}^T\mathbf{p}=f\left(\mathbf{q}\right)-{\mathbf{H}}^T\mathbf{h},\end{array}$$

(2)

where $\mathbf{p}\in {\mathbb{R}}^n$, and $f:{\mathbb{R}}^m\to {\mathbb{R}}^m$, $f\left(\mathbf{q}\right)=(f_1\left(q_1\right),\dots ,f_m\left(q_m\right))$. The function $f_j(·)$ describes the flow dependent pressure drop due to the hydraulic resistance in the jth edge. The relationship between pipe flow and energy loss caused by friction in individual pipes can be computed using the Hazen–Williams formula [21] for expression $f_j(·)$:

$$f_j\left(q_j\right)=\frac{10.7·L_j}{{\rho }_j^{1.852}·D_j^{4.87}}·q_j^{1.852},$$

(3)

where $L_j$ is the length of the pipe and $D_j$ is the diameter of the pipe, both in meters $\left[\mathrm{m}\right]$, $q_j$ is the pipe flow in m³/s and ${\rho }_j$ is the pipe roughness coefficient.

The term ${\mathbf{H}}^T\mathbf{h}$ is the pressure drop across the pipes due to the difference in geodesic level (i.e., elevation) in meters $\left[\mathrm{m}\right]$ between the ends of the pipes with $\mathbf{h}\in {\mathbb{R}}^n$ the vector of geodesic levels at each vertex.

2.2. Structure of the Reduced Order Model

The reduced-order network model is used in this paper to calculate the nominal pressure at the measured internal nodes. The model uses the pressure dependence of the network’s internal nodes with the pressure and flow measurements of the inlets. The details of the model derivation can be found in [22,23].

A network can be divided into nodes connected with reservoirs (the inlets nodes) and internal nodes that compose the system. To facilitate the explanation in this work, the information regarding inlet nodes will be represented by $\left(r\right)$ superscript and those of the internal nodes, which will be expressed by the $\left(i_n\right)$ superscript. In particular, vector ${\mathbf{p}}^{\left(in\right)}$ will contain pressure node values $p_1,\dots ,p_{n-n_I}$ and ${\mathbf{p}}^{\left(r\right)}$ inlet pressure values $p_{n+1-n_I},\dots ,p_n$.

The network needs to fulfill some conditions for using the reduced model proposed:

Condition 1: corresponds to the demands of the internal nodes of the system, where Equation (1) can be defined as:

$$\mathbf{d}\left(k\right)=-\mathbf{v}\left(k\right)\sigma \left(k\right),$$

(4)

where $\sigma \left(k\right)$ denotes the total inlet flow into the network at time instant k, the vector $\mathbf{v}\left(k\right)$ defines the distribution of the total demand in the internal nodes at every time k, with the property ${\sum }_i^n{v}_i\left(k\right)=1$. Notice that if all consumers are residential, all nodes demand have the same consumption profile, in consequence, the  $\mathbf{v}\left(k\right)$ will be constant $\mathbf{v}\left(k\right)=\mathbf{v}$.

Condition 2: is a particularly case when the vector ${\mathbf{p}}^{\left(r\right)}$ of control inputs fulfill the following case,

$${\mathbf{p}}^{\left(r\right)}\left(k\right)+{\mathbf{h}}^{\left(r\right)}=\kappa \left(k\right)\mathrm{𝟙},$$

(5)

for some $\kappa \in \mathbb{R}$, which is the total head at the inlets in $\left[\mathrm{mwc}\right]$ and where $\mathrm{𝟙}$ denote the vector consisting of ones. In [23], there is a discussion on this definition’s feasibility where the controllers should satisfy this premise at least in networks with the low total consumption.

If these two conditions are fulfilled, the pressure at the ith internal node can be expressed by:

$$p_i^{\left(i_n\right)}\left(k\right)={\alpha }_i{\sigma }^2\left(k\right)+\sum _{j=1}^{n_I}{\beta }_{ij}\left(k\right)p_j^{\left(r\right)}\left(k\right),$$

(6)

where ${\alpha }_i$ is parameter dependent on the network topology and the distribution of demands in the network, and ${\beta }_{ij}$ is dependent on the network topology with $j=1,\dots ,n_I$. The total inlet flow $\sigma$ is typically well-known since inlet flows are measured. ${\alpha }_i$ is a parameter dependent on the network topology and the distribution of demands in the network, and ${\beta }_{ij}$ is dependent on the network topology with $j=1,\dots ,n_I$. The total inlet flow $\sigma$ is typically well-known since inlet flows are measured.

Some methods of identifying parameters can be used to identify parameters ${\alpha }_i$ and ${\beta }_{ij}$ since model (6) of $p_i^{\left(i_n\right)}$ is linear [24], using the measures of $\sigma$, ${\mathbf{p}}^{\left(r\right)}$ and $p_i^{\left(i_n\right)}$ with nodes that contain pressure sensors that will be denoted as $p_{s_i}$$\forall i=1,\dots ,n_s$ in the following, where $n_s$ is the number of sensors installed in the inner nodes.

Once inner pressure model (6) has been calibrated, the accuracy of the model can be assessed by applying the computation of the model error or pressure residual defined by:

$$r_{s_i}={\widehat{p}}_{s_i}\left(c\right)-p_{s_i}\left(c\right),\forall i=1,\dots ,n_s,$$

(7)

where c denotes the boundary conditions (heads and inflows in inlets) necessary to compute pressure estimation by means of (6). For example, minimum and maximum residual bounds $\underline{{\sigma }_i}$ and ${\overline{\sigma }}_i$, considering the available data, can be computed for every sensor $i=1,\dots ,n_s$ to obtain an idea of the accuracy of model (6). Sensor noises and error models can produce residual errors. If big values of residual bounds are obtained, improvements in model (6) should be considered. For example, the assumption that all the nodes have the same consumption profile can lead to a big error in some networks. In this case, the error could be decreased if model (6) is calibrated only using data from the same hour but on different days. It would be assumed that different users can have different profiles at a given hour, but a particular user will have the same profile at a particular hour for all the different days. This possible improvement will imply the calibration of 24 different models (6) (one for each hour) and will require more historical data to obtain good accuracy. Another method to obtain an estimation of the pressure in inner nodes ${\widehat{p}}_{s_i}\left(c\right)$ is to use historical data directly as a lookup table, as was proposed in [18]. That is, given particular operating conditions, c provides the inner pressures from historical data that had the closest operating conditions $\widehat{c}$ to c. Residuals (7) considering leak pressure measurements will be used in the leak localization, as will be explained in detail in the next section.

3. Leak Localization

The location of the leak on the WDN is typically divided into two steps: leak detection and leak localization [25]. The focus in this work is the leak localization assuming that the detection has already been effectuated. In addition, it is assumed that the leaks can only happen in the nodes of the network (as considered in [7,26], or [8]), making the number of nodes equal to the number of potential leaks. The nodes correspond to water users, pipe junctions, and other structures such as hydrants. However, if the number of nodes will not provide a representative discretization of the network, some artificial nodes could be considered.

In this Section, two leak localization methods will be proposed. The first one will only use available measurements, and its diagnosis will point to one of the inner pressure sensors installed in the WDN. Therefore, the detected leak should be in an area around this sensor (cluster). The second method will combine the information of the first method with the topological information: characteristics of the pipes and connections between the nodes of the WDN, in a likelihood index that will allow the leak localization at the node level.

3.1. Leak Localization at Cluster Level

As stated before, the proposed leak localization is applied after the leak detection. In addition to inlet pressure and flow sensors, it is assumed that $n_s$ pressure sensors are installed at different inner nodes. Consider a leak $l_j$ acting on the node j of the network, and the used measurements are assumed to be captured under a leaky situation. Additionally, admitting leak-free historical data of all the sensors are available. The residual pressure in internal nodes that contain a sensor, defined in (7), can be computed as:

$$r_{s_i}={\widehat{p}}_{s_i}\left(c\right)-p_{s_i}\left(c^{l_j}\right),\forall i=1,\dots ,n_s,$$

(8)

where ${\widehat{p}}_{s_i}\left(c\right)$ is the pressure estimation considering boundary conditions c in a leak-free scenario. On the other hand, $p_{s_i}\left(c^{l_j}\right)$ is the pressure value measured by the inner pressure sensor i under boundary conditions $c^{l_j}$ (the same heads and inflows in inlets as in c but with a leak in node j).

Following the ideas in [18], positive residuals can be obtained from the following transformation:

$${\overline{r}}_{s_i}=r_{s_i}-\min (r_{s_1},\dots ,r_{s_{ns}})\forall i=1,\dots ,n_s.$$

(9)

Then, as the leak localization can be achieved by determining the residual pressure component with maximum size (see [22,27]), leak localization can be formulated as:

$$\widehat{J}=\underset{i\in \{1,\dots ,n_s\}}{arg\; max}\left\{{\overline{r}}_{s_i}\right\}.$$

(10)

Notice that the result of the leak localization method (10) is one of the $n_s$ pressure sensor locations.

Then, the leak localization results in $\widehat{J}$ point not only to sensor location $s_j$ but also to the nodes that produce a higher incidence for this sensor than the other sensors (cluster j).

3.2. Leak Localization at Node Level

Considering the Hazen–Willians Equation (3) for every pipe (edge $e_z$) a resistance $R_z$ can be defined:

$$R_z=\frac{10.7·L_z}{{\rho }_z^{1.852}·D_z^{4.87}}.$$

(11)

Among the multiple pipe paths that can connect every pair of nodes $ij$, a path ${\mathcal{P}}_{ij}^{min}$ with a minimum total resistance $R_{ij}$ can be computed by means of:

$${\mathcal{P}}_{ij}^{min}=\underset{{\mathcal{P}}_{ij}^{\left(k\right)}\in {\mathcal{P}}_{ij}}{arg\; min}\sum _{e_z\in {\mathcal{P}}_{ij}^{\left(k\right)}}{R}_z,$$

(12)

where ${\mathcal{P}}_{ij}=\{{\mathcal{P}}_{ij}^{\left(r\right)},\dots ,{\mathcal{P}}_{ij}^{\left(e\right)}\}$ denotes the set of paths connecting nodes i and j.

On the other hand, the minimum path from the $n_I$ inlets to a node j, ${\mathcal{I}}_j^{min}$, can be obtained by applying the computation of the minimum paths from the $n_I$ inlets to node j by means of (12) and determine which is the one with the minimum resistance among the $n_I$ paths.

When a leak is produced in node j, ${\mathcal{I}}_j^{min}$ is the most probable path for the extra flow produced by the leak. So the effect of a leak in node j to sensor $s_i$ depends on the intersection of the paths from inlets to node j and the node where the sensor is located $s_i$: ${\mathcal{I}}_j^{min}$ and ${\mathcal{I}}_{s_i}^{min}$. To quantify the degree of incidence of the leak to the sensor, an incidence factor $g_{j,s_i}$ is defined as:

$$g_{j,s_i}=R_{j,s_i}^c{\overline{g}}_{j,s_i},$$

(13)

where $R_{j,s_i}^c$ is the resistance of the path defined by ${\mathcal{I}}_j^{min}\cap {\mathcal{I}}_{s_i}^{min}$, the superscript c refers to the common path between node and sensors, and ${\overline{g}}_{j,s_i}$ is a normalization factor that takes into account the inverse of the resistance from the node j to the different sensors:

$${\overline{g}}_{j,s_i}=\left\{\begin{array}{ll}\frac{\frac{1}{R_{j,s_i}}}{{\sum }_{l=1}^{n_s}\frac{1}{R_{j,s_i}}}& \mathrm{if}j\ne s_i\\ 1& \mathrm{if}j=s_i.\end{array}\right.$$

The $n_s$ incidence factors associated to a leak in node j, $g_{j,s_i}$$i=1,\dots ,n_s$ can be normalized:

$${\lambda }_{j,s_i}=\frac{g_{j,s_i}}{{\sum }_{l=1}^{n_s}{g}_{j,s_l}},$$

(14)

where coefficient ${\lambda }_{j,s_i}$ determines the relative incidence of a leak in node j to sensor $s_i$ regarding all the $n_s$ sensors and the need to fulfill:

$$\sum _{i=1}^{n_s}{\lambda }_{j,s_i}=1\forall j=1,\dots ,n-n_I.$$

(15)

For every node $j=1,\dots ,n-n_I$, the most sensitive sensor to a leak in this node can be computed as:

$$\widehat{J}=\underset{i\in \{1,\dots ,n_s\}}{arg\; max}\left\{{\lambda }_{j,s_i}\right\}.$$

(16)

The $n_s$ clusters used in the leak localization defined in (10) can be computed using the set of nodes that provide the same value of $\widehat{J}$. The following equation is the definition of the cluster associated with the sensed node l:

$$C_l=\{v_j\in V|\underset{i\in 1,\dots ,n_s}{arg\; max}\{{\lambda }_{j,s_i}\}=l\},$$

(17)

where $l=1,\dots ,n_s$. The topological information of ${\lambda }_{j,s_i}$ and the measurement information of residuals ${\overline{r}}_{s_i}$ can be integrated in a parameter ${\theta }_j$ defined as:

$${\theta }_j=\frac{1}{\overline{\theta }}\sum _{i=1}^{n_s}{\lambda }_{j,s_i}{\overline{r}}_{s_i}$$

(18)

where $\overline{\theta }$ is a normalization factor. Then, ${\theta }_j$ can be interpreted as a likelihood index, and the leak localization at cluster level defined in (10) can be formulated at node level as:

$$\widehat{J}=\underset{j\in \{1,\dots ,n-n_I\}}{arg\; max}\left\{{\theta }_j\right\}.$$

(19)

In order to improve the performance of the leak localization method, the information of the residuals at different time instants k can be taken into account applying the Bayes’ rule as:

$$P_j\left(k\right)=\frac{P_j(k-1){\theta }_j\left(k\right)}{{\sum }_{l=1}^{n-n_I}{P}_l(k-1){\theta }_l\left(k\right)},$$

(20)

where $P_j(k-1)$ is the prior probability whose initial value $P_j(k-1)$ has to be determined (for example $P_j\left(0\right)=1/(n-n_I)$). Then, the leak node localization can be estimated by using posterior leak probabilities by:

$$\widehat{J}\left(k\right)=\underset{j\in \{1,\dots ,n-n_I\}}{arg\; max}\left\{P_j\left(k\right)\right\}.$$

(21)

4. Sensor Validation

When a leak is not detected by the leak detection method, anomalous values of pressure residuals $r_{s_i}\left(k\right)$$i=1,\dots ,n_s$ defined in (7) can be used to detect sensor faults. In the same operating conditions, the historical data of inner pressure sensors (leak-free data or data for a particular leak scenario) can be used first to calibrate a pressure estimation model as described in Section 2.2. Secondly, to determine residual bounds $\underline{{\sigma }_i}$ and ${\overline{\sigma }}_i$ that allows the implementation of pressure sensor fault detection through checking:

$$\left\{\begin{array}{c}{r}_{s_i}\left(k\right)\in \left[\underline{{\sigma }_i},{\overline{\sigma }}_i\right]\Rightarrow \mathrm{No}\mathrm{Fault}\left({\varphi }_i\left(k\right)=0\right)\hfill \\ r_{s_i}\left(k\right)\notin \left[\underline{{\sigma }_i},{\overline{\sigma }}_i\right]\Rightarrow \mathrm{Fault}\mathrm{in}\mathrm{sensor}{s}_i\left({\varphi }_i\left(k\right)=1\right).\hfill \end{array}\right.$$

(22)

The accuracy of this fault detection method depends on the length of residual bounds $\underline{{\sigma }_i}$ and ${\overline{\sigma }}_i$ and, therefore, on the accuracy of pressure estimation (6). In order to increase the accuracy of the fault detection method, spatial residuals [28] between pressure residuals (7) can be computed

$$S{r}_{s_i,s_j}\left(k\right)=r_{s_i}\left(k\right)-r_{s_j}\left(k\right)\forall i=1,\dots ,n_s-1\mathrm{and}j=i+1,\dots ,n_s.$$

(23)

In the same way as the pressure residuals, spatial residual bounds ${\underline{\epsilon }}_{i,j}$ and ${\overline{\epsilon }}_{i,j}$ can be computed using leak-free data, and the fault detection can be implemented as follows:

$$\left\{\begin{array}{c}S{r}_{s_i,s_j}\left(k\right)\in \left[{\underline{\epsilon }}_{i,j},{\overline{\epsilon }}_{i,j}\right]\Rightarrow \mathrm{No}\mathrm{Fault}({\mathsf{\Phi }}_{i,j}\left(k\right)=0)\hfill \\ S{r}_{s_i,s_j}\left(k\right)\notin \left[{\underline{\epsilon }}_{i,j},{\overline{\epsilon }}_{i,j}\right]\Rightarrow \mathrm{Fault}({\mathsf{\Phi }}_{i,j}\left(k\right)=1).\hfill \end{array}\right.$$

(24)

As model errors will affect in a similar way as close pressure sensors, it is expected that some spatial residual bounds will be smaller than pressure residual bounds. Therefore, fault detection defined by (24) will be more sensitive to pressure sensor faults than the one defined by (22). The accuracy of the sensor fault detection can be increased by means of average computing residuals in a time window leading to smaller residual bounds.

Once a residual has been violated, that is, at least one of the sensor faulty signals ${\varphi }_i\left(k\right)$$i=1,\dots ,n_s$ or spatial faulty signals ${\mathsf{\Phi }}_{i,j}\left(k\right)$$i=1,\dots ,n_s-1\mathrm{and}j=i+1,\dots ,n_s$ is equal to one, the sensor fault isolation can be implemented in two stages as described in Algorithm 1:

Algorithm 1 Sensor validation search for sensor fault

Stage 1: In the case of the activation of one or more sensor faulty signals ${\varphi }_i\left(k\right)$$i=1,\dots ,n_s$, as these signals are uniquely related to sensors $s_i$$i=1,\dots ,n_s$, the isolation is trivial and faulty sensors must be discarded for future leak localization, and the number of available healthy sensors $n_s$ should be updated. Stage 2: Only considers Spatial faulty signals ${\mathsf{\Phi }}_{i,j}\left(k\right)$ of the $n_s$ non-faulty sensors from $Stage1$. As these fault signals are potentially affected by two possible sensor faults $s_i$ and $s_j$, the fault isolation can be implemented iteratively by the following steps:

1: for $i\leftarrow 1,n_s-1$ do

2:  for $j\leftarrow i+1,n_s$ do

3:   for ${\mathsf{\Phi }}_{i,j}\left(k\right)==1$ then

4: $$\widehat{ı}=\underset{i\in \{1,\dots ,n_s\}}{arg\; max}\{\sum _{j=i+1}^{n_s-1}{\mathsf{\Phi }}_{i,j}\left(k\right)+\sum _{j=1}^{i-1}{\mathsf{\Phi }}_{j,i}\left(k\right)\}.$$ (25)

5:    Discard sensor $s_{\widehat{ı}}$, eliminate faulty signals related to this sensors, update $n_s$.

6:   end if

7:  end for

8: end for

In the case that two or more sensors obtain the same cost function in (25) and less than the maximum possible value $n_s-1$, the computation of (25) should be performed in a time window until new Spatial faulty signals are activated.

5. Case Study

5.1. Hanoi WDN

The Network used for the case study is a reduced city’s network model from Hanoi’s WDN (Vietnam). It is composed of one inlet (reservoir), 34 pipes, and 31 nodes, represented by Figure 1.

To analyze the performance of the proposed approach, data with different conditions have been generated artificially using the EPANET hydraulic simulator [29]. In order to consider realistic scenarios, some uncertainty has been added to the data [30]: the magnitude of the leak is random with a range of 25 to 75 [l/s], that is, between 1% and 2.5% of the average inlet flow of the WDN. Furthermore, white noise has been combined to emulate the noise present in real measurements, and uncertainty of 10% (uniform distribution) was added in the nominal demand value.

The daily water consumption pattern used for the calibration of Equation (6) is shown in Figure 2, having four days of operation.

The sample rate is 10 min, but average hourly measurements are calculated to reduce uncertainties on the diagnostic stage.

Results

The evaluation of the performance of the proposed leak localization method at node level defined in Equation (21) will be analyzed utilizing Average Topological Distance (ATD) [11]. The ATD represents the node’s distance between the node predicted as leaking and the actual node with the leak. To calculate the ATD, it is first necessary to create a matrix containing the minimum topological distance (in nodes or meters), $A\in {\mathbb{R}}^{n-n_I\times n-n_I}$.

Finally, the confusion matrix ${\Gamma }_{i,j}\left(n-n_I\times n-n_I\right)$ defined in [18] and depicted in

Table 1. Confusion matrix $\mathsf{\Gamma }$.

	${\widehat{\mathit{l}}}_1$	⋯	${\widehat{\mathit{l}}}_{\mathit{i}}$	⋯	${\widehat{\mathit{l}}}_{\mathit{n}-{\mathit{n}}_{\mathit{I}}}$
$l_1$	${\mathsf{\Gamma }}_{1,1}$	⋯	${\mathsf{\Gamma }}_{1,i}$	⋯	${\mathsf{\Gamma }}_{1,n-n_I}$
⋮	⋮	⋮	⋮	⋮	⋮
$l_i$	${\mathsf{\Gamma }}_{i,1}$	⋯	${\mathsf{\Gamma }}_{i,i}$	⋯	${\mathsf{\Gamma }}_{i,n-n_I}$
⋮	⋮	⋮	⋮	⋮	⋮
$l_{n-n_I}$	${\mathsf{\Gamma }}_{n-n_I,1}$	⋯	${\mathsf{\Gamma }}_{n-n_I,i}$	⋯	${\mathsf{\Gamma }}_{n-n_I,n-n_I}$

is used to assess the performance of Equation (21). The rows of this matrix correspond to the leak scenario and the columns to where the leak is located ($\widehat{l}$) by the leak localization method.

Considering the confusion matrix $\mathsf{\Gamma }$, the ATD can be computed as follows:

$$ATD=\frac{{\sum }_{i=1}^{n-n_I}{\sum }_{j=1}^{n-n_I}{\mathsf{\Gamma }}_{i,j}{A}_{i,j}}{{\sum }_{i=1}^{n-n_I}{\sum }_{j=1}^{n-n_I}{\mathsf{\Gamma }}_{i,j}}.$$

(26)

Four cases have been considered with different quantities of sensors in the network to analyze how this affects the final result.

Table 2. Nodes with sensors.

Case	Nodes with Sensors
1	12, 17, 23, 29
2	6, 12, 17, 23, 29, 21
3	6, 12, 15, 17, 23, 21, 27, 30
4	6, 9, 12, 15, 17, 24, 21, 22, 28, 29, 31

presents the distribution of the selected nodes to contain a sensor. As seen in [31], the positioning of the sensors produces different results. As this work did not discuss the adequate sensors’ arrangement, they were chosen to consider an improvement in the results regarding the number of sensors.

Using the inlet flow data and non-leak historical pressure measurements of the selected sensors, the ${\beta }_i$ and the ${\alpha }_i$ with $i=1,\dots ,31$ in (6) have been identified (notice that $n_I=1$). With these parameters, the pressure estimations under a non-leak condition in the network can be calculated considering inlet measurements using Equation (6) and posteriorly applied to calculate the residuals (9) with measured pressures in leak scenarios. In addition, non-leak pressure measurements and estimations are used to generate fault-free pressure residuals $r_{s_i}\left(k\right)$ and bounds $\underline{{\sigma }_i}$, ${\overline{\sigma }}_i$$i=1,\dots ,n_s$, as well as spatial residuals $S{r}_{s_i,s_j}\left(k\right)$ and bounds ${\underline{\epsilon }}_{i,j}$, ${\overline{\epsilon }}_{i,j}$$\forall i=1,\dots ,n_s-1\mathrm{and}j=i+1,\dots ,n_s$.

For every sensor configuration, normalized incidence factors (14) have been computed with topological information: node connections and pipe characteristics (length, diameter, and roughness). Figure 3 has the objective of comparing the information of the incidence of single leaks to pressure sensors obtained by a hydraulic model with the one obtained by means of topological information. The nodes selected to have sensors are the ones defined in the first case in Table 2. In particular, Figure 3c shows the clustering that groups the nodes that produce the highest effect in a specific pressure sensor. Nodes in yellow define the cluster of nodes where a leak produces a maximum pressure deviation from the non-leak scenario in the sensor installed in node 12, and the same is true for nodes in violet, red and green regarding pressure sensors in nodes 17, 23, and 29, respectively. Finally, nodes in black are nodes that produce a similar variation of pressure (difference of variation less than 0.1 [mwc]) in at least two different pressure sensors. In order to obtain this information, a hydraulic model to compute the difference of non-leak and leak pressures in all the nodes for the different leaks is required. On the other hand, Figure 3a shows the clustering that takes into account the shortest weighted pipe length (hydraulic distance), that is, the sum of $(L_z/D_z^{4.87})$ for all edge $e_z$ in the path to the sensors, being the smallest one used to define the most resemblance to the sensor, used in Ref. [18]. Finally, the clustering depicted in Figure 3b is defined by Equation (17), which is based on the common resistance path explained in Section 3. These two last clusters that only require topological information could be used in the leak localization method at the cluster level defined in Equation (10). It is important to emphasize that the clustering based on the resistance common path, proposed in this paper and depicted in Figure 3b, resembles the clustering based on the actual leak effect in the network (given by the model) depicted in Figure 3c much more than the clustering based in the hydraulic distance depicted in Figure 3a. Therefore, the clustering proposed in this paper provides more accurate information for leak localization purposes than that based on the hydraulic distance. For example, as shown in Figure 3c, when a leak is present in nodes 3, 4, 5, 6, 7, 8 or 9, the sensor most affected by the leak is the sensor in node 12. This information is the same as the one provided by the clustering depicted in Figure 3b and is computed only with topological information. However, using the clustering of Figure 3a based on the hydraulic distance between nodes and sensors, the closest sensor to these nodes is the sensor in node 17.

Figure 4 shows the correlation analyses of the relative incidence index ${\lambda }_{j,s_i}$ defined in Equation (14) for all the nodes $j=1,\dots ,31$ depicted in every subplot for every sensor $s_i$$i=1,\dots ,4$. As this index is normalized, its values are in the range [0,1). The nodes with the higher index (more brown color) are those that produce a higher effect in the pressure sensor $s_i$.

Figure 5a displays the evolution of the ATD (in nodes) obtained by the leak localization method based on the Kriging spatial interpolation methodology presented in [18] with the time horizon (in hours) used recursively by the Bayes’ rule in (20). Four different sensor configurations are considered with 4, 6, 8, and 10 sensors placed optimally in order to maximize the performance of the leak localization proposed [18]. The performance can be compared with the one obtained by the new leak localization method proposed in this paper at node level defined in Equation (19) with the same dataset and the same sensor configurations as in [18], depicted in Figure 5b and with the sensor configurations shown in Table 2, depicted in Figure 5c.

Figure 5a shows that the leak detection performance of the Kriging method improves significantly from four to eight sensors and more moderately compared to ten sensors, still having a good result, even with noise data managing to reach an ATD equal to 2.5 node. When compared to the newly proposed leak localization method, as can be seen in Figure 5b,c, the new leak localization method always outperforms the Krigring method, even in the case of using the sensor configurations proposed in [18] that were computed to optimize the performance of the Kriging method. Figure 5c shows that the sensor configurations proposed in [18] are not optimal for the proposed method but the performance can be improved by changing the sensor configurations, in this case manually.

In order to illustrate the performance of the proposed sensor validation method, Case 1 (four sensors) will be considered. The four-sensor residuals computed by Equation (7) have been considered in a time window of 24 h using leak-free data leading to upper residual bounds equal to:

$$\left[{\overline{\sigma }}_1,{\overline{\sigma }}_2,{\overline{\sigma }}_3,{\overline{\sigma }}_4\right]=\left[0.11,0.06,0.09,0.11\right]$$

and the lower residual bounds equal to:

$$\left[{\underline{\sigma }}_1,{\underline{\sigma }}_2,{\underline{\sigma }}_3,{\underline{\sigma }}_4\right]=\left[-0.14,-0.10,-0.10,-0.08\right].$$

In the same way, the six spatial residuals defined by (23) have been computed in the same conditions as sensor residuals leading to spatial residual bounds:

$$\left[{\overline{\epsilon }}_{1,2},{\overline{\epsilon }}_{1,3},{\overline{\epsilon }}_{1,4},{\overline{\epsilon }}_{2,3},{\overline{\epsilon }}_{2,4},{\overline{\epsilon }}_{3,4}\right]=\left[0.06,0.07,0.08,0.04,0.05,0.03\right]$$

and

$$\left[{\underline{\epsilon }}_{1,2},{\underline{\epsilon }}_{1,3},{\underline{\epsilon }}_{1,4},{\underline{\epsilon }}_{2,3},{\underline{\epsilon }}_{2,4},{\underline{\epsilon }}_{3,4}\right]=\left[-0.04,-0.06,-0.08,-0.06,-0.09,-0.03\right].$$

Figure 6 and Figure 7 depict the evolution of sensor and spatial residuals with their respective residual bounds in a fault scenario of sensor 1 that corresponds to the pressure sensor in node 12. The fault is a drift of $0.1\left[mcw\right]$ that starts on the 5th day. As shown in Figure 6, by applying (22) to sensor residuals it is impossible to detect the fault until the end of the day 9 (i.e., 4 days later) when residual sensor 1 violates the bounds. However, by applying (24) to spatial residuals it is possible to detect the fault in 10 h: $S{r}_{s_1,s_2}$ violates its bounds in 10 h, and $S{r}_{s_1,s_3},S{r}_{s_1,s_4}$ violate their bounds in 16 h and 22 h, respectively.

5.2. Modena WDN

The second case study selected to test the performance is the reduced model of the real water distribution network of the Italian city Modena. This large-scale network is comprised of 268 junctions (nodes) connected through 317 pipes and served by four reservoirs. There are no pumps in the network since it is entirely gravity-fed [32,33].

The EPANET hydraulic simulator was used to generate artificial data to analyze the performance of the proposed method. The following simulation conditions were considered:

• The leak scenario consists of data samples collected every 10 min and filtered to hourly values to reduce the uncertainty in the data;
• The uncertainty of demand is considered by introducing the uncertainty of 10 [%] (normal distribution) of the nominal demand value. In addition, white noise is deemed to emulate the noise in the measurements;
• The leak size is randomly selected with a range of 3 to 6 [l/s], representing 1% to 2.5% of the network consumption.

The sensor bias, sensor drift, and abrupt sensor failure of sensor faults were proposed to analyze the sensor validation method. The sensor bias fault was simulated as a step change, and the drift fault was given as a time-varying ramp signal. In both cases, the fault magnitude was randomly chosen with a range of $0.1$ to $0.2$ [mwc]. The last fault was simulated by turning the sensor output to zero.

Results

As applied in the previous case study, the Average Topological Distance (ATD) was used to assess the performance of the proposed leak localization method at the node level defined in (19). Two scenarios have been considered with five and ten pressure sensors that are presented in Figure 8a,b respectively. As emphasized in the last section, performance in the leak localization task is highly dependent on the number of sensors installed in the network [34,35,36].

Figure 9 shows the result of ATD evolution as defined in (26) applied with Bayes’ posterior time reasoning (20) to represent the leak location performance of the proposed method. This figure shows that the leak localization performance reached an ATD of 8 and 5.5 nodes with 5 and 10 inner pressure sensors installed in the network respectively. Considering that the proposed leak localization method only requires topological information and non-leak historical data in available measurements, the obtained performance is reasonably good.

On the other hand, a total of 6000 scenarios were simulated with 10 days each to evaluate the sensor validation method for the five sensor configurations depicted in Figure 8a. Thus, 1000 scenarios were generated for each sensor with sensor bias, sensor drift, and abrupt sensor failure applied randomly, and the remainder 1000 without faults.

To calculate the residual and spatial residuals bounds, a 6-month leak-free scenario was generated. The five sensor residuals computed by Equation (8) considering the time window of 24 h and increasing $24\%$ observed bounds. Leading to upper residual bounds equal to:

$$\left[{\overline{\sigma }}_1,{\overline{\sigma }}_2,{\overline{\sigma }}_3,{\overline{\sigma }}_4,{\overline{\sigma }}_5\right]=\left[0.10,0.06,0.04,0.01,0.04\right]$$

and to lower residual bounds equal to:

$$\left[{\underline{\sigma }}_1,{\underline{\sigma }}_2,{\underline{\sigma }}_3,{\underline{\sigma }}_4,{\underline{\sigma }}_5\right]=\left[-0.08,-0.05,-0.03,-0.01,-0.06\right].$$

Following, the ten spatial residuals defined by (23) were computed in the same conditions as sensor residuals leading to spatial residual bounds:

$$\begin{array}{cc}[{\overline{\epsilon }}_{1,2},\hfill & {\overline{\epsilon }}_{1,3},{\overline{\epsilon }}_{1,4},{\overline{\epsilon }}_{1,5},{\overline{\epsilon }}_{2,3},{\overline{\epsilon }}_{2,4},{\overline{\epsilon }}_{2,5},{\overline{\epsilon }}_{3,4},{\overline{\epsilon }}_{3,5},{\overline{\epsilon }}_{4,5}]=\hfill \\ \hfill & \left[0.07,0.08,0.08,0.10,0.06,0.05,0.06,0.03,0.06,0.05\right]\hfill \end{array}$$

and

$$\begin{array}{cc}[{\underline{\epsilon }}_{1,2},\hfill & {\underline{\epsilon }}_{1,3},{\underline{\epsilon }}_{1,4},{\underline{\epsilon }}_{1,5},{\underline{\epsilon }}_{2,3},{\underline{\epsilon }}_{2,4},{\underline{\epsilon }}_{2,5},{\underline{\epsilon }}_{3,4},{\underline{\epsilon }}_{3,5},{\underline{\epsilon }}_{4,5}]=\hfill \\ \hfill & \left[-0.07,-0.09,-0.08,-0.09,-0.05,-0.04,-0.06,-0.03,-0.04,-0.04\right].\hfill \end{array}$$

For this study, the evaluation metric applied was classification accuracy. To this purpose, the confusion matrix was used, which presents the classification accuracy and the misclassification error, and the horizontal axis of the confusion matrix describes the predicted labels of samples, while the longitudinal axis depicts the true labels of samples. The right side shows the percentages of correctly and incorrectly classified observations for each true class.

Figure 10 illustrates the result for the confusion matrix for all scenarios generated, and depicts that the accuracy of detecting faults in the sensor is very high, where the lowest accuracy is presented in fault sensor number five with an accuracy of $95.4\%$ and the highest in fault sensor number three with $100\%$ accuracy. Regarding the accuracy of the scenario with no-fault, eight of the 1000 fault free scenarios presented one false alarm among the 240 samples of the scenario; therefore providing an average interval between false detections of 240,000/8 = 30,000 h.

6. Conclusions

A new data-driven method for leak localization in WDN based on historical non-leak data and the topological information of the network is proposed. The proposed method is triggered when a leak is detected, and it is based on the evaluation of residuals generated by leak pressure measurements in some inner nodes and the estimation of leak-free pressures in these nodes utilizing a reduced-order model and historical data. Topological information is used to compute a new incidence factor that considers the most probable path of water from reservoirs to pressure sensors and potential leak nodes. The proposed incidence factor combined with residual information generates a likelihood index that allows leak localization at the node level. In addition, a sensor validation method based on the sensor pressure residuals, which is able to detect and isolate pressure sensor faults, is proposed.

The proposed method’s general performance for leak location and sensor validation is evaluated in reduced models of the Hanoi and Modena water distribution networks. The results of the leak localization are compared to another technique published with satisfactory results. Future works can be developed to improve the leak localization and sensor validation performances, with a study of an algorithm to automatically determine the optimal sensors required to maximize the leak localization performance.