Leakage detection and calibration of pipe networks by the inverse transient analysis modified by Gaussian functions for leakage simulation

Fast and reliable leakage detection of a pipe network has significant economic and environmental benefits. This issue has attracted the attention of many engineers and researchers worldwide to develop leakage detection methods and devices. In the past three decades, various leakage detection methods based on hydraulic principles have been introduced.

The numerical modelling of pipe flow is highly sensitive to leakage specifications and pipe friction factors. This fact can be efficiently used for leakage detection and calibration of the system. The aim of model-based methods is to determine the number, location and size of leakages in the system as well as pipe friction factors. The model-based methods may be developed based on steady or unsteady (transient) flow hydraulics. A transient pipe flow results in pressure waves that periodically propagate through the system and are affected by the system features and abnormalities, e.g. leakages or blockages. The transient pressure signals at a few measurement sites can be used to assess the entire network for leak detection and calibration of pipe friction factors.

In general, the transient-based models can be classified into two categories: the frequency-domain techniques (Suo & Wylie 1989; Mpesha et al. 2001; Lee et al. 2005; Sattar & Chaudhry 2008; Ghazali et al. 2012; Duan & Lee 2015) and time-domain techniques (Liggett & Chen 1994; Brunone 1999; Vítkovský et al. 2000; Covas & Ramos 2001; Wang et al. 2002; Kapelan et al. 2003; Nixon & Ghidaoui 2007; Haghighi et al. 2012a, 2012b; Meniconi et al. 2018). While the former category is mostly applicable to single and simple pipelines, the latter category is capable of handling complicated pipe systems with various nonlinear boundary conditions. Inverse Transient Analysis (ITA) is a time-domain leakage detection method that is mathematically applicable to any piping system, from single pipelines to complex and large-scale pipe networks.

Liggett & Chen (1994) first introduced the ITA method. The unknowns of an ITA problem are the leakage parameters as well as pipe friction factors. Since the number of leakages is a priori unknown, the problem dimension (size) is unknown in advance. To make the problem mathematically solvable, it is assumed that all junctions of the network have a candidate leakage with unknown area size. By this trick, the number and location of candidate leakages are determined in the beginning and consequently removed from the unknown parameters. Accordingly, a standard ITA for leakage detection and calibration of a pipe network may be outlined as comprising the following main steps:

1. Transient flow is generated in the network by exciting a nodal consumption.

2. The transient pressure heads at some measurement sites are measured.

3. A hydraulic simulation model is developed for the network at hand. This model predicts the nodal pressure heads at the measurement sites as a function of the unknown parameters.

4. A mathematical programming problem is introduced to match the measured and model-predicted pressure heads at the measurement sites, while the decision variables are the leakage parameters and pipe friction factors.

5. A nonlinear optimization solver is applied to solve the above programming problem.

It is worth noting that the application of transient-based fault detection methods to real piping systems is quite difficult and constrained by several field considerations and deficits. Producing a useful and informative transient flow as well as reliable data acquisition in complex pipe networks are difficult tasks which are well addressed in Meniconi et al. (2015, 2018). While the criticalities coming from the execution of transient tests in practice are acknowledged, the present study intends to focus on the numerical aspects of the ITA method. If the measured data are not well processed and reliably analyzed the entire ITA fails. In this study, the importance of mathematical programming of the ITA and the necessity of management of the problem solution are emphasized.

A literature review shows that the computational cost and the result accuracy are the main challenges of ITA addressed by the previous studies. Liggett & Chen (1994) used the mathematical optimization method of Levenberg–Marquardt (LM) to solve the ITA problem. LM is a local optimization solver that works based on the first-order (Jacobian) and second-order (Hessian) partial derivatives of the objective function. Although they successfully applied the LM method to a benchmark example, Vítkovský et al. (2000) showed that the ITA's search space is multimodal, and hence local solvers may be easily trapped in local solutions. To cope with this problem, they suggested a continuous genetic algorithm (GA) to solve the benchmark example and indicated that the GA can properly overcome the multimodality of the problem. Nevertheless, the GA, similar to the other stochastic evolutionary methods, is computationally expensive and slow to converge. The transient analysis is itself time-consuming. When it is merged with the GA evolution process the entire ITA becomes still more time-consuming and computationally inefficient. As a remedy, Kapelan et al. (2003) hybridized the GA with the LM to take advantage of both methods. By applying the new GA–LM method to the benchmark example, they concluded that the hybridization of a fast local gradient-based method (LM) and a global metaheuristic method (GA) would result in a faster and more reliable ITA. Vítkovský et al. (2007) investigated the ITA experimentally through a laboratory pipeline. They addressed several issues affecting the ITA accuracy. They also introduced a model parsimony to reduce the number of the ITA's decision variables and its complexity. Their experimental observations confirmed that the ITA is successful in detection and quantification of leakages in the system. Shamloo & Haghighi (2010) investigated whether if other steps of the ITA including the transient generation and measurement sites are optimized, the ITA decision space becomes near convex and can be still solved by fast mathematical methods. On this basis, they successfully used the mathematical solver of Sequential Quadratic Programming (SQP) to solve the benchmark example. Later, Haghighi & Ramos (2012) used a new metaheuristic method, namely the Central Force Optimization (CFO), for ITA problems. The CFO, originally introduced by Formato (2007), is a deterministic nature-inspired metaheuristic based on the metaphor of gravitational kinematics. By solving the benchmark example, they found that the CFO has the advantages of both gradient-based and evolutionary optimization methods so that its performance is somehow similar to the hybridized optimization models. Recently, Huang et al. (2015) developed an ITA model using the Simulated Annealing (SA) method. They applied the SA to a real pipeline and a synthetic pipe network from the literature and discussed their method's capabilities.

Since the number of leakages is a priori unknown the ideal ITA approach is intractable to solve. Also, the ITA may be carried out in several stages so that first the leaky junctions are detected and then the pipes connected to them are surveyed one by one. In fact, for detection of a small number of real leakages the traditional ITA must investigate a large number of hypothetical leakages. Despite many studies so far done to develop leakage detection techniques based on hydraulic modeling and inverse problems, several challenges and limitations related to the problem dimension and ill conditions, the process of solving and optimization, the accuracy and reliability of results and handling uncertainties still remain unanswered. To improve the traditional ITA in terms of dimensionality and handling leakage location and distribution in the system, this study proposes a new style of leakage simulation using Gaussian functions. Also, some of the above issues are addressed by introducing an algorithm for reducing the problem search space and increasing the chance of reaching the global solutions by managing the ITA implementation.

The most popular method to solve the above equations is the Method of Characteristics (MOC) that is adopted here to develop the simulation model based on Chaudhry (2014) formulations.

In addition to the friction loss term, another nonlinearity of the MOC is the leakage relationship (Equation (3)). In such nodes the leakage discharge is added to the flow unknowns. In this condition, Equation (3) is required to be simultaneously solved with the characteristic Equations (5) and (6) while where superscripts L and R respectively indicate the pipe discharge at the left and right hand sides of the leakage. Since Equation (3) is nonlinear, a Newton–Raphson scheme is used to determine the flow unknowns at each time step.

Since the number of leakages is unknown in advance, the problem dimension is unknown too. As a remedy, the traditional ITA supposes that all junctions have a candidate leakage with unknown size. As a result, the decision variables of the above objective function become the area sizes of the candidate leakages as well as pipe friction factors. As discussed earlier, this trick increases the problem complexity and multimodality. If one decides to detect the potential leakages in the pipes in addition to the main junctions, the number of decision variables will significantly increase, causing the ITA to become more burdensome and time-consuming to solve. Furthermore, the traditional ITA supposes that leakages only exist on the characteristic nodes. Therefore, to increase ITA accuracy, the network hydraulic simulation domain (time–space grid) must be discretized with a small mesh size. This issue also increases the load of ITA computations. In addition, because of possible uncertainties and errors in the measurements and modelling, the leakage parameters cannot be precisely determined in practice.

Experimental and field investigations on the ITA have shown that under best conditions the model-based methods are able to identify a leakage neighborhood for a true leakage. A leakage neighborhood is an area with a larger leakage at the centre and some smaller leakages around that as shown in Figure 1(a) from Covas & Ramos (2010) and Haghighi et al. (2012a, 2012b). To cover the above problems with the traditional ITA this study proposes a new style of leakage simulation in pressurized pipes using Gaussian functions. The idea arises from observations of the ITA results in real applications that show that a true leakage is practically detected with a quasi-normal distribution of smaller leakages (Haghighi et al. 2012a, 2012b).

To incorporate the Gaussian functions for leakage simulation into the ITA algorithm as well as to manage the problem decision space, a modified version of the ITA is introduced as in the following.

1. In the third step of the ITA, the hydraulic simulation model is extended so that not only the main junctions but also all pipes' characteristic nodes can have a potential leakage.

2. To accelerate the ITA as well as to avoid the above-mentioned problems with large-scale pipe systems the decision variables are managed as follows:

   • Similar to the traditional ITA all pipes of the network have unknown friction factors.

   • Instead of supposing many candidate leakages at the main junctions and the characteristic nodes, a few S errant leakages are introduced to the network.

   • Unlike the traditional ITA, the errant leakages are neither fixed at the junctions nor concentrated on certain points. In the new version, all characteristic nodes have a zero-value leakage except those having an errant leakage (Figure 1(b)). Each errant leakage has a quasi-normal distribution around a peak according to the following Gaussian function.

     

     (5)

     where b, c and > 0 are the Gaussian function constants and e is Napier's number. The plot of a Gaussian function is a symmetric bell-curve shape that quickly falls off towards plus/minus infinity. The parameter b is the height of the curve's peak that justifies the leakage area size; c is the position of the centre of the peak that indicates the location of leakage; and d controls the width of the bell that determines how a leakage is distributed around its centre. The parameter d is useful for handling the uncertainties of measurements and simulations. In addition, each errant leakage needs a pipe identifier to indicate in which pipe of the network the leakage is located. Accordingly, each errant leakage introduces four decision variables to the ITA, including b, c and d. It is also worth noting that x (between zero and pipe length ⁠) is the distance of any location in the pipe from the upstream end node from which the Gaussian function peak location c is estimated (Figure 1(b)). After the Gaussian function of a leakage was obtained, the contribution of each characteristic node in the leakage area size is determined substituting the node distance x into Equation (5).

3. To solve the problem, the objective function (Equation (4)) with the above decision variables is minimized. For this purpose, in the beginning an integer value (greater than one) is supposed for the number of errant leakages in the system. Accordingly, the problem dimension size would be (leakage parameters) (pipe friction factors), thus incorporating the ITA decision variables.

4. A real genetic algorithm adopted from Goldberg (1989) is developed and coupled to the hydraulic simulation model to minimize the objective function. The optimization finally yields one of the following results:

   • at least one of the errant leakages has almost zero area size ⁠;

   • at least two errant leakages have the same pipe identifier and Gaussian function;

   • all errant leakages have different pipe identifiers or Gaussian functions.

5. Cases (i) and (ii) mean that the number of true leakages in the system is less than the initially assumed errant leakages and hence the optimization makes the extra leakages either zero (case (i)) or exactly equal to the other detected leakages (case (ii)). If so, the ITA is done and the optimum decision variables represent the leakage parameters and friction factors in the network. Otherwise, if a third case applies, it means that the initially supposed number of errant leakages has been insufficient and the system may have more leakages. If so, to continue the ITA the number of leakages is updated as where is an integer. Then, the ITA is repeated from step (4) and the algorithm is continued until either of the result cases (i) or (ii) happens.

It is worth mentioning that there is a trade-off between the number of errant leakages and the ITA computations in each step. The decision on this issue is case-dependent and needs engineering judgment. While it is generally suggested to start with a small number for S for better management of the ITA solution, one may initially consider a large number of errant leakages, e.g. ⁠, in the beginning. If so, the entire ITA is very likely terminated in one single run however, with the ITA computations becoming more expensive as well as the multimodality of the search space increasing. Even if 10 errant leakages are assumed in advance, the total number of leakage parameters is at most 40, which is significantly less than the total number of network characteristic nodes. Also, with the aid of the Gaussian function the number of decision variables becomes totally independent of network size as well as the number and length of pipes. Furthermore, because of the reduction in the search space the convexity of the objective function is better preserved. As a result, the applied optimization solver would find the global solution faster and more reliably.

To investigate the proposed scheme as compared with the traditional ITA, a pipe network example shown in Figure 2 from Pudar & Liggett (1992) is taken into account. As a benchmark case study this pipe network has been already considered for leakage detection in several studies cited in the introduction. The pipe network consists of eleven pipes, seven main junctions and a single reservoir with 30 m constant head at junction 1 that gravitationally supplies the system. The system also has a 20 l/s inflow to junction 7 and a 58 l/s outflow from junction 4. All pipes have the same length (762 m), diameter (254 mm) and wave speed (1,316 m/s). The system has a leakage with 1 cm² area at the midpoint (node 48) of pipe 5. The original values of pipe friction factors are also presented in Figure 2. In this case study all friction factors as well as the leakage parameters are considered unknown. To solve the problem inversely, both traditional and modified versions of the ITA method are applied to the problem. The obtained results are then compared.

First, the network is mathematically simulated using the MOC. For this purpose, all pipes are discretized with 10 reaches of 76.2 m length. Consequently, the system computational time step is 0.058 s. In total there exist 106 characteristic nodes (seven main junctions and 99 in-pipe characteristic nodes) in the simulation model. All nodes but the reservoir one (node 1) have candidate leakages resulting in 105 unknown leakage area sizes. Also, there exist 11 unknown pipe friction factors. Therefore, the problem has a total of 116 unknown parameters to be determined using the ITA.

To start the ITA a transient flow is produced in the system and the transient pressures in some sensitive locations are measured. Through an optimization model, Shamloo & Haghighi (2010) proposed a transient excitation scenario for this problem in a manner that the outflow from junction 4 is linearly decreased from 58 l/s to zero in 13.25 s and increased in the same time duration to 58 l/s. Through a sensitivity analysis they also introduced junction 4 as the best location, with priority one, for the measurements.

After the transient excitation and measurement in the system, the objective function (Equation (4)) is defined according to the measurements and the network simulation model. To solve the problem inversely the objective function is minimized by calibrating the total 116 unknown decision variables. To this end, two approaches of the traditional and modified ITA are applied and investigated as in the following. For both models, the same genetic algorithm is used; however, the GA parameters and population size in each approach are tuned according to the problem's dimension size.

1. First, the network is solved using the traditional ITA through which all 116 decision variables are calibrated together to minimize the objective function. The GA's population size is decided as 120 chromosomes through some initial trial-and-error runs. Then, the GA was run several times. Almost all runs converged to the same best solution in about 1,500 to 2,000 generations. For each run of the traditional ITA, on average about 180,000 objective function evaluations corresponding to 180,000 pipe network transient simulations are required. At the end of the optimization, the nodal leakage area sizes as well as pipe friction factors are obtained as illustrated in Figure 3(a) and 3(b), respectively. Also, the model-predicted and measured pressure heads at junction 4 exhibit an excellent agreement as shown in Figure 4(a). However, in spite of that and expensive computations for the traditional ITA, the decision variables are not in good agreement with their true values as depicted in Figure 3(a) and 3(b). In fact, the ITA performance can be checked on two criteria: first the minimum value of the objective function and second the errors associated with the predicted variables. Although the objective function has been minimized very well, the parameter estimation errors are not satisfactory. Figure 3(a) and 3(b) show that the traditional ITA has distributed the system's leakage over the other nodes. Also, the pipe friction factors are not well calibrated. For better understanding of the model performance it is worth noting that the calibrated friction factors include about 25.72% relative error, and the detected leakage area sizes on average have about 4.35 cm² absolute error. It is also concluded that the current ITA problem has a relatively flat decision space about the global minima so that the problem has no strict extremum. Consequently, a wide range of solutions may minimize the objective function and make it close to zero. This issue introduces an ill condition to the ITA optimization and threatens its reliability in practice. The example reveals that although the ITA with any dimension size is mathematically solvable, the obtained results could be unreliable and misleading.

2. Using the modified scheme of ITA, the problem's dimensionality as well as the optimization search are systematically managed. For this purpose, the system is considered to have only two candidate errant leakages in advance. Each candidate leakage is simulated by the Gaussian function introducing the four unknown parameters of ⁠, b, c and d to the problem, as explained. Accordingly, the ITA objective function is minimized by calibrating 19 decision variables in total, consisting of 8 leakage parameters and 11 friction factors in the first step. In this condition, the number of decision variables is about 6 times lower than the number of traditional ITA variables. This issue significantly results in the reduction of the problem search space and consequently increases the chance of reaching the global optima faster and more reliably. The GA with 40 population size is applied to solve the problem. Similarly to the previous approach, the GA was run several times. Almost all runs successfully converged to the best solution in about 100 to 150 generations. Table 1 and Figure 5(b) respectively present the optimum leakage parameters and pipe friction factors at the end of the first step optimization. The results are required to be investigated to decide on terminating or continuing the ITA. As seen in Table 1, between the two initially assumed errant leakages one has a zero-valued peak ⁠. Therefore, we can terminate the ITA and state that the system has only one leakage. As shown in Table 1 the detected leakage has an area of 0.978 cm² (2.2% relative error in size) which is located in pipe 5 at 373.4 m distance from the upstream end. Also, the effect of parameter d = 26.67 m on the detected leakage in pipe 5 is illustrated in Figure 5(a). As seen in the figure, the leakage has mostly concentrated on the characteristic node 48 that is the true location of the leakage and has no significant effect on its vicinity nodes 47 and 49. As a result, from the leakage detection point of view the modified ITA outperforms the traditional scheme remarkably. According to Figure 5(b), the method is also successful in calibrating the pipe friction factors, although the objective function is much less sensitive to them than to the leakage parameters. Figure 4(b) also demonstrates the calibrated and measured transient pressure heads at the measurement site of node 4. Comparing Figure 4(a) and 4(b) reveals that both traditional and modified versions of the ITA have precisely matched the simulation with the measured data. Nevertheless, only the modified scheme has been successful in identifying the true leakage as well as pipe friction factors. In fact, the modified ITA considerably lessens the burden of unnecessary computations for the optimization, transient analysis and Newton–Raphson iterations for handling the nonlinearity of the leakage and friction loss equations. Also, it improves the problem's accuracy and chance of escaping local optima or fake solutions.

Finally, Table 2 presents a general comparison between the above approaches on the computational speed and accuracy.

The ITA method is a well known model-based approach for leakage detection and calibration of pipe networks. Since the number of unknowns in an ITA problem is much less than the sampled pressure heads, the problem is mathematically underdetermined and needs optimization to solve inversely. A traditional ITA problem is generally high dimensional and has a multimodal search space which may be ill conditioned particularly when the ratio of sampling sites to the unknowns is very small. For large-scale pipe networks with a big number of decision variables, although the ITA is mathematically solvable, the results may include significant errors. It is found that large dimensionality of the problem would significantly deform the ITA's search space so that the problem is no longer strictly convex around the global solution. Accordingly, a wide range of different solutions may minimize the objective function. Another problem with the traditional ITA is that the candidate leakages are introduced only at the characteristic nodes, while in reality a leakage may happen anywhere in a pipe. To address the above issues, some modifications were introduced here to improve the ITA reliability and efficiency. Based on the benchmark example, it is concluded that the assumption of quasi-normal distribution for leakages and simulating that by Gaussian functions result in the following. First, it significantly reduces the number of unknowns and causes the number of decision variables to be independent of the size of the network; second, it makes it possible to allocate candidate leakages everywhere in the network; and third, it helps to handle uncertainties of modeling and measurements as well as limitations of sampling rate. Also, the proposed step-by-step approach for implementing the ITA is found very useful and promising for managing the inverse solution. This approach not only overcomes the curse of dimensionality and multimodality of the search space but also significantly avoids unnecessary computations for optimization and transient analysis. Nevertheless, the proposed modifications were only applied to a theoretical example. For more investigations on the methods and its capability and limitations, more applications and experiments are highly required.