Optimal selection of number and location of pressure sensors in water distribution systems using geostatistical tools coupled with genetic algorithm

Network design has been widely developed and used in various branches of engineering science in the last three decades or so. Application of network design can be classified into two categories: (1) service networks (e.g., bank and/or airport networks) and (2) monitoring or data collection networks. In water resources planning and management, the monitoring networks are broadly designed and applied in different disciplines including precipitation monitoring (Bastin et al. 1984; Pardo-Igúquiza 1998; Shahidi & Abedini 2018a), groundwater monitoring (Loaiciga et al. 1992; Li & Chan Hilton 2005), river quality monitoring (Karamouz et al. 2009), discharge monitoring (Alfonso et al. 2013), sewer system monitoring (Banik et al. 2017), wastewater treatment plant monitoring (Villez & Corominas 2016), and water distribution systems (WDS) monitoring (Kapelan & Savic 2005; Romano et al. 2013). As a result, a logical and time-efficient approach in delineating monitoring instruments in space, time, and/or space–time domain is considered as an essential component for decision-makers.

Monitoring of WDS is usually required for a variety of purposes. Model calibration, leak detection, and pollutant delineation are typical examples requiring monitoring of state variables, hence sensor placement. In WDS, the sensors are generally divided into: (1) hydraulic sensors for monitoring the nodal pressure or pipe flow rate and (2) water quality sensors to detect contamination events posing a growing threat to public health. Pressure sensors are more frequently used compared to flow meter sensors as they are cheaper and easier to install and collect nodal pressure data. The pressure sensors give instantaneous readings, whereas the flow ones do not react instantaneously to changes in flow rate (de Schaetzen et al. 2000).

Monitoring of pressure field is considered to be a fundamental requirement for proper management of the system. The knowledge of pressure values and its subsequent management generally drives and controls operational actions such as available pressure for various usage, leakage, and demand control. An urban WDS usually contains hundreds or thousands of nodes and pipes, but the budget constraints do not allow the deployment of pressure sensors at every potential measurement location. Hence, one has to decide where to optimally place the finite number of pressure sensors, out of a vast number of potential locations. No doubt, the accuracy of monitoring network depends on the quality and quantity of the collected data. As a result, the measurement layout should be optimal in the sense that it should maximize the information content of the system in an efficient and cost-effective way. Therefore, the selection of an appropriate collection location, called the sampling design or sensor placement design, is important and has long been a challenge for researchers and practitioners in the water industry (Kapelan & Savic 2005).

In the design of sensor placement in WDS, a number of interrelated factors might affect the final network arrangement. These factors are considered to be the overall objective of designing a sensor network, the demand condition (i.e., pressure-driven versus demand-driven approaches), system dynamic (i.e., steady, extended period simulation (EPS) and/or unsteady), attribute under consideration (pressure or flow), topography of the system (flat or mountainous), the nature of objective function (e.g., variance-based, entropy-based, fractal-based, and distance-based techniques) and the optimization algorithm used for minimization or maximization purposes (e.g., exhaustive search or random search). In the present study, the overall objective is to make pressure field monitoring – in particular, the average pressure over the entire network – more accurate within the WDS. As a consequence, the attribute for monitoring will be nodal pressure. The study area itself will dictate the nature of topographic variation. In this study, the newly proposed approach for pressure sensor placement design will be implemented on two case studies.

In this paper, both the accuracy of the mean pressure estimation in a direct way (as expressed by the variance of residuals) and the economic cost of the data collection network in an indirect way will be the sole objective. To solve this problem, a well-known geostatistical variance-based reduction method, the so-called variance of residuals computed via block ordinary kriging (BOK), is successfully used. This variance-based procedure has been widely used for the design of surface water facilities such as rain-gauge networks (Pardo-Igúquiza 1998; Shaghaghian & Abedini 2013; Attar et al. 2018), groundwater monitoring networks (Ben-Jemma et al. 1994; Yeh et al. 2006), etc.

To address a typical sensor placement problem, there are three different approaches: (1) simplified approach; (2) exhaustive search approach; and (3) random search approach. The exhaustive search surveys all possible combinations independently which somehow creates a problem when the number of nodes within the system is quite large, running into the curse of dimensionality issue. In other words, for a fixed number of nodes, N, the objective function for small and large values of n sensors in is quite small which can be examined and quantified in an efficient way. However, for intermediate values of n in ⁠, this exhaustive and comprehensive examination of each combination is not generally possible. Researchers surmount this combinatorial, NP-hard problem using Bellman's principle of optimality via some simplified approaches (Bastin et al. 1984; Kassim & Kottegoda 1991; Chen et al. 2008) or using random search optimization such as artificial bee colony (ABC) (Attar et al. 2018) or genetic algorithm (GA) (Shahidi & Abedini 2018b).

The current study intends to present a new technique whereby a variance-based objective function will be implemented to place pressure sensors via simplified, exhaustive, and random search procedures (to be described in more detail later). The paper is organized as follows. A critical review of the literature on sensor placement is briefly presented. The sensor placement problem is then formulated followed by introducing the proposed methodology along with a geostatistical framework. In subsequent sections, the proposed methodology is tested and verified on a literature case study of Anytown and a real-life WDS of C-Town consisting of five district metered areas (DMAs). Finally, conclusions which can be drawn from this study are summarized in the last section.

During the past three decades or so, numerous investigators have addressed the problem of optimal sensor placement layout. Walski (1983) was one of the first to determine where to measure pressures and flows in a WDS to calibrate a model. He proposed monitoring stations to be placed near the high-demand zones and on the perimeter of the system, away from water resources. Later on, most of the sampling design approaches have been based on the sensitivity criteria of measurement points with respect to calibration parameters and these are generally grouped into two categories: (1) ranking potential locations in reference to a sensitivity-based criteria (Ferreri et al. 1994; Bush & Uber 1998; Piller et al. 1999) and (2) developing an optimization problem (Lee & Deininger 1992; Yu & Powell 1994; de Schaetzen et al. 2000; Meier & Barkdoll 2000; Kapelan et al. 2003; Vitkovsky et al. 2003). Generally speaking, some unique invariant features of either Jacobian matrix or variance-covariance matrix are used to identify the most sensitive locations for monitoring purposes. Additionally, most approaches have solved the associated optimization problem by means of GAs (de Schaetzen et al. 2000; Meier & Barkdoll 2000; Kapelan et al. 2003; Behzadian et al. 2009).

Kapelan et al. (2003) introduced a deterministic, multi-objective genetic algorithm for sampling design. The two objective functions proposed were to maximize network accuracy and to minimize the total cost. The first-order second-moment (FOSM) uncertainty model, previously introduced by Lansey et al. (2001), was used to assess both parameters and model prediction uncertainties. Later, Behzadian et al. (2009) used Kapelan et al. (2003) theoretical development by considering a stochastic approach rather than a deterministic way to design a pressure sensor layout for model calibration. Finally, Simone et al. (2016) proposed a completely different methodology for designing a pressure sampling design based on the WDS topological analysis. For water quality monitoring, optimal sensor placement has also been considered with the aim of identifying accidental or intentional contamination sources in a timely manner (Kessler et al. 1998; Ostfeld & Salomons 2004; Berry et al. 2006; Aral et al. 2010; Rathi & Gupta 2014). The cited scholars introduced various methods to minimize timely detection of contamination and the risk incurred on the population via optimal placement of sensors in WDS.

These aforementioned methods have several limitations, as follows. (1) Some methods solve the sensor placement problem via ranking nodes which rely solely on the fact that the optimal set for n measurement locations is always a superset of the optimal set for n-1 locations. Needless to say, even though the ranking type approach is computationally superior and easy to set up, it does not solve an optimization problem and they might fail to identify the optimal solution. (2) The approach adopted by most studies has been deterministic not giving appropriate weight to the most important state variable within the system, pressure. (3) Most (if not all) of the studies have addressed pressure sensor placement as an optimization problem based on the sensitivity criteria of measurement locations with respect to predefined variations in unknown parameters of interest such as hydraulic parameters (e.g., nodal demands, pipes' roughness coefficients, etc.). These studies used FOSM uncertainty model based on Jacobian and variance-covariance matrices to find the most sensitive locations for monitoring. FOSM approach is based on some assumptions that may not always be true in the WDS modeling (Kapelan et al. 2007). In addition, the conventional paradigm adapted in sensor placement design also requires calculation of derivatives of model-dependent variables with respect to calibration parameters that may be computationally demanding and prone to possible numerical errors and discontinuity issues. The aim of the current study is to develop and apply a totally different approach using geostatistical jargon to overcome some of the above limitations for optimally placing pressure sensors in a given WDS.

In this paper, the major objective is to propose a methodology whereby a combination of geostatistical tools as an estimator of pressure field and an algorithm of optimization is effectively utilized to design an optimal pressure sensor network. In the current study, nodal pressure is considered as a regionalized variable having spatial statistical structure of its own. The sensor network design is implemented using a variance-based approach whereby the objective function is the variance of residuals over a block as large as the WDS. The proposed algorithm was first tested and verified on a benchmark, widely used synthetic system of Anytown WDS (Kapelan et al. 2003) and then applied to C-Town WDS (Ostfeld et al. 2012). EPANET 2.0 (Rossman 2000) is used as the hydraulic simulator of the corresponding WDS. Nodal pressures, obtained through running the hydraulic solver, are used to delineate and quantify the model of spatial variability, i.e., the variogram function.

Formulating and solving the sensor placement optimization problem presented here involves identifying a model of spatial variability, discretizing the study area, developing an optimization algorithm to delineate sensor spatial location(s) for various values of n in through the various approaches, including simplified, exhaustive, and random search procedures.

As the proposed algorithm is coined with geostatistical jargon, for the sake of paper integrity, a brief account of geostatistical concepts will be provided first. Interested readers may want to consult numerous textbooks written on the subject elsewhere (Isaaks & Srivastava 1989; Cressie 1993; Chilés & Delfiner 1999).

The suggested methodology is first applied and verified on a hypothetical case study known as ‘Anytown’. This system (Figure 1) was first set up by Walski et al. (1987) as a benchmark to test different types of calibration models, and contains many aspects of real-world systems. It resembles reality in terms of topology layout and nature of pipes including city center area with older pipes and boundary, and residential areas with newer pipes (Kapelan et al. 2007). Anytown WDS consists of 34 pipes, 16 nodes, two elevated storage tanks, and three parallel pumps that transfer water from a reservoir to the system. The Anytown configuration and data used here are taken from Kapelan et al. (2003). The steady-state model is only used for average demand condition. A total of N= 16 nodes are considered to be possible potential locations for pressure monitoring. Anytown WDS, with the above specifications, is modeled and then nodal pressure is determined. In order to assign the x-y coordinate to each node, a local Cartesian coordinate system is defined whereby its origin is on Node 70. The coordinate of each node along with pressure values are given in Table 1.

In this work, pressure field is considered as a two-dimensional random process which could vary in time and space. However, it is assumed that the hydraulic model of Anytown WDS is steady-state, so that nodal pressure is invariant with respect to time, and the process only varies with space. Assuming intrinsic hypothesis and considering no variation with direction (i.e., isotropic process), the variogram will be solely a function of separation distance.

It should be noted that the experimental omnidirectional variogram cannot be used directly in the kriging system as positive definiteness of variance computation cannot be guaranteed. Christakos (1984) touched on admissibility of variogram and covariance further, and can be consulted for more detailed explanations. The variogram model should have a mathematical expression that can describe the variance of random process with changing distance. This step is the most challenging task, as selection of model parameters has direct impact on the weight coefficients of BOK and the variance of residuals. For steady-state condition, nodal pressures do not change with time and the associated variogram model could be considered as time-invariant. The experimental variogram is constructed with nodal pressure data given in Table 1. Furthermore, some of the best-fitted variogram models used are tabulated in Table 2. The parameters of variogram models are determined based on nonlinear least-square fit. The selected theoretical model is based on smallest residual sum of squares (RSS) consisting of a spherical structure as shown in Figure 2. In the current study, the proposed technique is implemented and rationalized on Anytown WDS via two different approaches as follows.

In this approach, Bellman's principle of optimality (Bellman 1957) is assumed. According to this principle, an optimal combination at a particular stage considers the earlier chosen sensors to be in place and to build a new combination in reference to the earlier one. In a sense, this principle converts an objective function with n decision variables into n objective functions with only one decision variable.

In the simplified approach, for various values of n out of N nodal points, one has to compute the variance of residuals with increasing number of sensors sequentially. In Table 3, the values of variance of residuals are determined for various numbers and optimal placement of pressure sensors. As can be seen, for a one-sensor case, Node 90 is found to produce minimum variance of residuals. It is the nearest node to the center of the system and has high base demand as well.

Figure 3 shows the plan of spatial locations for sensors delineated based on the simplified approach (solid circles). For five-sensor case, Nodes 90, 130, 30, 100, and 80 are determined, which are generally near the center or the perimeter of WDS, away from water sources.

The minimum variance of residuals versus the number of pressure sensors is depicted in Figure 4. As can be seen, a comparatively small number of nodes needs to be monitored to achieve minimum variance of residuals. Further increase in the number of sensors would not affect the variance of residual dramatically. For example, if five out of 16 system nodes are monitored, the variance of residuals is approximately only 9% more than the case whereby all nodes are monitored.

Figure 5 plots the numerical values of against n for Anytown WDS. As can be seen, for Anytown WDS, with only 16 nodes, the curse of dimensionality is not considered to be a big issue. Therefore, the exhaustive search could be effectively used as an applicable tool to reach the global optimum solutions.

The optimal monitoring locations via exhaustive search approach can also be seen in Figure 3 (solid triangles). It can be noted that the best pressure monitoring locations include Nodes 40, 70, 80, 160, and 100, which are totally different compared to the simplified approach result.

Figure 6 demonstrates the variation of the minimum variance of residuals against number of sensors via the exhaustive search approach. As can be observed, like the simplified approach, only a limited number of nodes is required to significantly reduce the variance of residuals.

It can be noted from Table 4 that the optimal set of locations for n sensor is not necessarily a superset of the optimal set for n-1 sensor placement. For instance, the set of four pressure sensor locations contains Nodes 40, 70, 80, and 160, while the set of three optimal pressure sensor contains Nodes 50, 70, and 160. Therefore, sensor placement design based on the simplified approach or any other methodology in which the optimal set of n sensor locations is derived from the set of n-1 optimal sensors may fail to recognize a truly optimal solution. This issue has previously been proved by Kapelan et al. (2003) as well.

Finally, in reference to the proposed approach, in Table 5, a comparison of the presented method for designing the monitoring network in the Anytown WDS with other studies is given. Whereas the optimum number of pressure sensors is five in the current study, in similar studies with totally different objective function, the total number of optimum sensors was found to be five (Kapelan & Savic 2005) and six (de Schaetzen et al. 2000). Furthermore, for the optimum number of five, in both simplified and exhaustive search approaches, four selected nodes are similar to other studies which shows the capability of the proposed methodology.

With respect to the findings obtained, a critical discussion is comparatively summarized as follows. (1) For the case study analyzed, the graphs clearly demonstrate that reduction of the estimation variance is not uniform for both simplified and exhaustive search approaches. In other words, at an early stage of selection process, the rate of variance reduction is quite significant. However, later on, the rate of variance reduction diminishes significantly (approximately 2%). (2) For both approaches, the first and last selection of pressure sensors is identical. This means when n equals one, the center of the study area will be filled first (i.e., Node 90). This deployment layout results in prevent­ing the accumulation of the measurement station within a small area of the WDS. (3) In addition, the variance of residuals for ultimate values of n (15, 16) in all approaches is identical over Anytown. This means for a certain number of sensors, the precision will not change even though the total costs increase. (4) While the optimal number of pressure sensors from both approaches would be four, variance of residuals in the exhaustive search approach will behave as lower limit compared to the simplified approach. In a sense, rate of variance reduction in the exhaustive search approach is more rapid compared to the simplified viewpoint. As a consequence, the exhaustive search approach tends to reach the minimum variance of residuals much earlier.

The second case study presented here is the C-Town WDS, made up of one reservoir, seven tanks, 388 nodes, 432 pipes, 11 pumps grouped into five pumping stations, and four valves. Water is supplied to the system by a large reservoir with a constant head and seven water tanks. The system is divided into five DMAs, each of which features a pumping system and a system of tanks as shown in Figure 7. Water demands at each node, hourly tank levels and pumping station flows are available for a period of 168 h (1 week) in C-Town.inp and C-Town.xls files (Ostfeld et al. 2012).

For the C-Town case study, a local coordinate system is considered for assigning the x-y coordinate to each node which is specified with an ID number. The first 24 h data are divided into four intervals of 6-h long and then nodal pressure is calculated for each interval resulting in four pressure realization patterns. For each DMA, experimental variogram is calculated by the aforementioned method. As shown in Figure 8(a)–8(e), the cluster of points corresponding to each time interval imitates a similar trend. Although this trend in Figure 8(d) and 8(e) are a little different, the trend still exists. Also, the points are so close that they can hardly be distinguished. For modeling the variogram, one possible approach would be to model the theoretical variogram corresponding to each realization of the random field in time. The best-fit of variogram model for each DMA is specified as the exponential model.

After assigning the model to the experimental variogram for each interval, model parameters (‘c’ as sill and ‘a’ as range) are determined by the least-squares fit. After comparing model parameters illustrated in Table 6, it can be seen that the variation in the variogram function resulting from changing sill value (c) is more significant at each interval for DMA1. Survey of model parameters for other DMAs confirmed the above-mentioned outcome. As a result, if variogram function is considered as γ[h; c(t),a], it can be written in the form of c(t); γ*(h; a) where c(t) is called a scaling factor multiplied by γ*(h; a) as a shaping factor or the spatial autocorrelation factor. Index t shows that c factor is dependent on time steps (Bastin et al. 1984). It is a kind of separation of variables used extensively in analytical solution of partial differential equations. Indeed, the variogram model is separated into two factors: c(t), which is time-varying but space-invariant, and γ*(h; a) which is time-invariant but space-dependent. Thus, it can be pointed out that it is possible to apply only an exponential model as a mean time-invariant variogram over four intervals as a unique representative model γ_rep presented for each DMA in Figure 8(a)–8(e). It will be proved that the results would be true for time steps less than 6-h long and sensor placement outcomes are independent of time variation of nodal pressure. As demonstrated for the previous case, the pressure sensors' placement problem will be typically solved via two different approaches, discussed below.

Principally, the DMA design is a typical way to manage monitoring stations in WDS. For this purpose, the simplified approach is independently applied for each DMA on C-Town WDS. The optimal nodes for pressure sensor placement for five DMAs are depicted with solid circles in Figure 9(a)–9(e). The results show that the sensors have a probable tendency to be homogeneously distributed on all parts of the DMAs which can be considered to be effective owing to the increase in the reliability level of the sensor network. In addition, nodes which are the final receivers of water (away from the tanks) have the best chance to be the first candidates in each DMA. For example, for a single sensor scenario, Node 149 in DMA1, Node 260 in DMA2, Node 231 in DMA3, Node 185 in DMA4, and Node 90 in DMA5 attain the highest priority for sensor placement. This confirms that the optimal sensor placement satisfies suggestions presented by Walski (1983) and verified by Kapelan et al. (2003). Moreover, it can be noted that the optimal number of sensors is determined when the rate of reduction of the objective function tends to be nearly invariant; where it is six for DMA1 (with 134 nodes), five for DMA2 (with 110 nodes), and four for each of DMA3 (with 39 nodes), DMA4 (with 56 nodes), and DMA5 (with 49 nodes). The set of critical nodes are identified as follows: 149, 343, 312, 377, 137, and 107 for DMA1; 260, 28, 248, 283, and 11 for DMA2; 231, 218, 228, and 238 for DMA3; 185, 159, 200, and 299 for DMA4; 90, 82, 100 and 89 for DMA5.

Table 7 demonstrates time variation of objective function for the optimal number of sensors for DMA1. As can be seen, the deployment of the best sensors is quite identical for all time steps, and only the numerical value of the objective function (variance of residuals) varies from one time step to another. The values of objective functions are only scaled by the c factor, i.e., the objective function (Equation (4)) could also be expressed as σ²_BOK(t) =c(t).V*_BOK where c is the time-variant part and V* as a shaping factor is time invariant which needed to be computed once for all time steps. Under this condition, the outcomes of best placement of pressure sensors depend purely upon shaping factor which is time-invariant. It will be illustrated that even if time interval changes, this stability in the outcome will persist. This also indicates robustness of the proposed algorithm in the real-world WDS.

In this work, for C-Town WDS, with 388 nodes, the curse of dimensionality imposes limitation on exhaustive search algorithm. As stated earlier, for such systems, as all possible combinations cannot be searched out, a random search optimization algorithm should be undoubtedly implemented. For example, for this system with 388 nodes, for selection of five-sensor case, this gives nearly 72 billion possible combinations. In recent years, many heuristic algorithms have been developed, which try to find optimal solutions of these problems through a random or unstructured search procedure. Since the decision variables are integers (ID of the nodes), the use of a discrete genetic algorithm (GA) can be efficiently applicable. This algorithm is best suited to solve combinatorial optimization problems that cannot be solved using more conventional methods. It has been used to reduce the global search space as well as the overall computational cost for the optimal sensor placement in WDS. During the optimization process, the maximum number of candidate nodes for each DMA was selected as 20. The following GA parameters are used based on a number of trial runs: population size of 50, single-point crossover with probability 0.80, and random-by-gene mutation with probability 0.20. All runs were performed for 50 generations. For the first several generations, the objective function value shows a significant improvement and then the convergence rate decreases dramatically.

After running the code, the optimal sensor placements for each of the five DMAs are marked with solid triangles in Figure 9(a)–9(e). As can be seen, from the random search approach, the best nodes for sensors are placed either in the same nodes or in the proximity of selected nodes corresponding to the simplified approach. The result of this procedure, along with the one obtained from the simplified approach, is plotted in Figure 10 for DMA5 and graphically plots the minimum variance of residuals as a function of the number of sensors. It can be observed from Figure 10 that both graphs match reasonably well. Very similar results were obtained for other DMAs.

In Table 8, results of sensor placement obtained through the random search approach are further compared with the simplified approach for all DMAs. It can be noted from Table 8 that optimal sensor placement via the simplified approach is very close compared to near-optimal solution via the random search approach. Although the value of variance of residuals for any optimal number of sensors for random search approach has to be less than the simplified approach, the difference is not quite significant.

This paper addressed the optimization problem and the solution process for design of pressure sensor placement in order to make the estimation of the pressure field more accurate and efficient in the WDS. The objective is to find the best measurement locations in the system to collect the required data. We tend to design the pressure sensor network which maximizes accuracy with the minimum total costs associated with it. Obviously, these two factors depend on the number and spatial location of the pressure sensors in the WDS. The current study intends to use geostatistical tools to define the single-objective function. Therefore, for finding the best combinations of pressure sensors, a measure of fitness which minimizes the variance of residuals is determined. The developed objective function is subsequently used to extract the optimum pattern of pressure sensors.

The proposed methodology is first implemented on a literature WDS of Anytown to design a pressure sensor network. After coupling a variance-based objective function with an exhaustive search optimizer as an algorithm of minimization, ranking of pressure sensors at potential nodes is achieved. The proposed methodology is compared and contrasted with conventional paradigms to indicate the advantages of the proposed scheme. In particular, the proposed scheme outperforms other approaches in terms of mathematical labor (e.g., no need for Jacobian matrix, variance-covariance matrix and its subsequent creation via numerical differentiation, etc.), accuracy, robustness, and computational time. After testing and verifying the proposed approach, the novel technique has been applied to a real-world case study of C-Town. However, in this case, to overcome the curse of dimensionality, a random search algorithm (GA) is used. The proposed methodology benefits from the following advantages:

• 1.

  A geostatistical tool (BOK) coupled with an optimization algorithm (exhaustive and random search) could clearly be an appropriate, robust, and efficient approach to extract the decaying function of block variance of residuals versus the number of pressure sensors.

• 2.

  The proposed method gives planners, water utilities, and many others involved in decision-making an opportunity to better invest and manage limited budgets by considering the cumulative costs in designing an optimum network implicitly if not explicitly.

• 3.

  Note that, even though sensor placement methodologies would be computationally expensive, especially, when applied to real-case scenarios, the method presented here clearly demonstrates that a system with DMAs would help to decrease the required computational operations. In reference to the results obtained and the peculiarity associated with variogram modeling for the entire system, one could argue that, for a system containing DMAs, delineation of optimal sensor placement for each district is more logical and defendable compared to such delineation taking the WDS model as a single unit. For large-scale systems without separate districts, any type of clustering can be applied for simplification.

• 4.

  After delineating the optimal arrangement of sensors for a given WDS, collected data can be used as a tool for real-time estimation of the pressure field and then the results can be considered as an input for purposes such as model calibration, leak detection and localization.

The recommendations for further work on the proposed methodology involve testing and validation on more complex real-life systems with thousands of nodes for real-time estimation of the pressure field where the demand at each node will be governed by corresponding pressure. Moreover, the capability of this technique on contamination detection should also be investigated. These suggestions would pave the way to investigating the use of geostatistical tools for designing monitoring network in the large-scale real-world WDS.