Burst detection and localization in water pipelines based on an extended differential evolution algorithm

This paper presents a new burst detection and location technique for pressurized pipelines based on an extension of the Differential Evolution (DE) algorithm. The proposed approach addresses the burst location problem as an optimization task, by considering the dynamic model that describes the behavior of a ﬂ uid through a pipeline and the presence of ﬂ uid losses produced by a burst. The optimization problem relies on ﬁ nding suitable estimations related to the burst parameters, i.e. magnitude, pressure and position of a burst, while a de ﬁ ned cost function is minimized. In order to deal with this problem, three strategies are proposed to extend and adapt the DE algorithm: (i) an informed de ﬁ nition of the physical restrictions of the problem according to the pipeline characteristics; (ii) a training stage of the algorithm that allows to ﬁ nd the appropriate synthesis parameters; (iii) a multi-start structure, in order to track dynamical variations of the problem. Experiments on a pipeline prototype illustrate the results obtained by the proposed algorithm on the estimation of the burst parameters, comparing its performance with an algorithm based on the Extended Kalman Filter, which is widely used in the literature.


INTRODUCTION
Water loss management is one of the most important issues facing water suppliers around the world.In many water pipelines, a significant percentage of the total water volume supplied is lost due to bursts and leaks, resulting in discontinuous, expensive and inefficient services.Furthermore, the presence of a burst can lead to the rapid deterioration of the pipeline, which might allow contaminant ingress and affect the urban infrastructure (Xie et al.As far as we know, the use of heuristic techniques to address the BDL problem has not been explored before. The DE algorithm has been widely applied in the water management and optimal design of water distribution systems (Suribabu ; Zheng et al. ).In this paper, the DE algorithm is proposed to address the burst detection and localization problem as an optimization task that relies on finding suitable estimations related to the magnitude, pressure and position of a burst, while a proposed cost function is minimized.This algorithm addresses the optimization problem through a sequence of evolutive operations (mutation, recombination and selection) in order to find the candidate solution with the best cost function evaluation.
In order to achieve feasible solutions with a low number of iterations, the present paper also proposes a three-

Pipeline model
The equations describing the transient pressure and flow rate in closed conduits are a set of partial differential equations known as Water Hammer Equations (WHE), which are derived applying the conservation laws on a control volume (Roberson et al. ).In general, the transient flow is considered unidirectional since the terms of the convective acceleration are smaller than the axial terms.
Therefore, by neglecting these terms and the spatial density variation, besides considering a fluid slightly compressible, pipes slightly deformable, and a leveled pipeline, the WHE can be expressed as (Chaudhry ): where a is the pressure wave speed in the fluid In order to compute friction effects for transient flow conditions, different methods have been developed.In general, the friction losses are divided into quasi-steady (J s ) and unsteady (J u ) parts, as J w ¼ J s þ J u .The quasi-steady term is usually estimated by the Darcy-Weisbach equation: where the friction factor f can be estimated by the Swamee and Jain equation: Considering a pipeline discretization in two sections, as shown in Figure 1, the WHE can be approximated by the next nonlinear model: where z L is the distance from the beginning of the pipeline to the burst position, u ¼ [H 1 H 3 ] T is the input vector that contains the pressures at the ends of the pipeline, y ¼ [Q 1 Q 3 ] T is the output vector that contains the flow rates and f(Q 1 ) and f(Q 3 ) represent the friction coefficients calculated with the input and output flows, respectively.
In this work, the proposed method only considers a single burst.Equation ( 6) can be written in the classical form of a nonlinear state space representation: To obtain a computationally feasible representation of model ( 6), a temporal discretization based on the Heun's method is performed as in Delgado-Aguiñaga et al. (), as follows: where and k is the index of the discrete time and Δt is the time step.
Model (8) can be written in compact form as follows: where In a practical scenario, the typical instrumentation of a transmission pipeline is composed of pressure and flow rate sensors located at the ends of the pipeline.The main objective of this sensor configuration is to evaluate the supplied volumes to the distribution zones and to monitor the efficiency of the system.This pipeline instrumentation scheme is suitable for the use of model ( 9) in the implementation of a BDL algorithm.Also, it is necessary that the sampling rate satisfies the Courant-Friedrichs-Lewy (CFL) criterion in order to guarantee the numerical stability of model ( 9).The CFL condition can be expressed as: Δz ! a Δt, where Δt is the time step.

Longitudinal adjustment
In general, pipeline systems are composed of different fittings in order to adapt the installation to the topography features.When the fluid flows through a fitting, such as a contraction or an elbow, the viscous effects and the turbulence produced by the shape of these accessories cause an energy loss, known as a local loss, which is added to the friction produced in the straight pipe sections.The WHE are developed considering a straight pipeline, without accessories.Thus, in order to obtain a dynamical model that includes fittings, in terms of the WHE, it is necessary to express the length of the pipeline in an equivalent straight length (L eq ).Hence, each accessory in the studied pipeline will be seen as a virtual straight pipe section that generates the same local loss associated with its corresponding accessory (Mataix ).Following this, L eq is computed by using the Darcy-Weisbach equation: where ΔH is the differential head pressure at the ends of the pipeline in a burst-free condition [m].Finally, by letting L ¼ L eq in Equation ( 8), it is possible to introduce the effects of the fittings into the pipeline model.
Using the longitudinal compensation, the burst is localized in a virtual equivalent straight length, however, in a practical situation, this does not represent a complete solution.In addition to the burst localization, it is necessary to establish a conversion from the obtained virtual position to a real-length value.In this work, we use the procedure described by Badillo-Olvera et al. () in order to perform this conversion.

Differential Evolution algorithm
The DE is a powerful optimization algorithm, which has proved to be very robust in a wide variety of problems This process will be described as follows.
• Mutation: In this step, mutant vectors v i,g ¼ [v  given by fixing F near to 2.
• Recombination: In order to improve the diversity of the population, a recombination process is defined after the mutation step, performing a crossover operation in each individual χ i,g , by using as its partner the mutant vector v i,g .This operation generates a trial vector μ i,g ¼ [μ {i,1,g} , μ {i,2,g} , . . ., μ i,d,g ] T by means of: where C r is the crossover rate, which is a constant fixed in the interval [0, 1], and it controls the admittance of the entries from the mutant vector in the formation of the trial vector.
• Selection: In order to determine the best candidate sol- utions that will pass to the next generation g þ 1, each trial vector μ i,j,g is compared against its corresponding parent, as follows: The processes of mutation, recombination and selection are iteratively executed until the optimum parameters are obtained or a pre-specified stop criterion is satisfied.

BURST PARAMETER IDENTIFICATION STRATEGY
The general formulation of identification techniques is based on the comparison between the response of a real system y (sensor measurements), against the outputs given by a related mathematical model ŷ , while a performance criterion that quantifies the value of the approximation error is used.Therefore, the optimization task relies on finding a set of values of the burst parameters minimizing a criterion function of the estimation error.
In order to evaluate (9), it is necessary to count with the values of z L (k), H 2 (k) and λ(k) since these are unknown and correspond directly to the solution of the BDL problem.In this way, the proposed extended DE algorithm is used to estimate such values, using a population of candidate sol- Now, considering the measurements Q 1 (k) and Q 3 (k), and its estimations c Q 1 (k) and c Q 3 (k) given by ( 9), it is possible to establish a performance measure in order to quantify how close ŷ is from y, for each candidate solution.So, model ( 9) is completed by using the values obtained by the proposed algorithm and then, the performance of each χ(i, g) is evaluated by the cost function ( 14) in order to select the best candidate solution: where order to select the best candidate solution.
The BDL problem can be expressed as an optimization task as: where S is the feasible space and χ up and χ low are the lower and upper bounds for each parameter to be estimated.
The following subsections present three strategies that allow to extend the classical DE algorithm in order to tackle the BDL problem efficiently.

Initial bounds of the search space
Before the population of the DE algorithm can be initialized, both upper χ up and lower χ low bounds for all the candidate solutions must be specified.In this case, for each characteristic of the burst, there is a certain range within which the value of the bounds should be restricted, due to its physical correspondence.A good definition of the feasible search space and their boundaries provides extra knowledge about the optimization problem, making the deployment of the algorithm faster and efficient (avoiding the consideration of unfeasible candidate solutions).
The BDL problem is defined by the unknown values of the burst position, pressure and magnitude of model ( 9), whose physical boundaries are included in Table 1.
In Table 1, ϵ m is selected as 0:001 in order to guarantee the numerical stability of the algorithm.On the other hand, it is clear that the upper bound for the burst position is a value near to the total length of the pipeline and the pressure at the burst point must be between the inlet and outlet pressure.However, in order to define the upper bound of the magnitude factor (λ), it is necessary to determine the maximum tolerable crack size before a total fracture of the pipeline occurs.Under this consideration, the upper bound of the magnitude factor is given by: where e c is the stable grow crack size supported by the pipe and C d is a discharge coefficient that can be selected near to 1.
When a burst localization method is used, it is desirable that the orifice or rupture remains stable in order to provide sufficient time to isolate and repair it.The analysis that describes the necessary conditions for this situation is known as Leak-Before-Break (LBB) criterion.A basic analysis, based on Linear Elastic Fracture Mechanics (LEFM), to establish the LBB criterion is presented by Anderson ().
Considering this, the crack size will remain stable if: where K IC is the fracture toughness of the pipe material

Calibration zone
One of the main difficulties that a user faces when a heuristic algorithm is applied is to select an appropriate set of synthesis parameters (Lobo et where the lower bounds in both cases are 0.0001 in order to guarantee numerical stability and the upper bounds of F and C r are 2 and 1, respectively, in agreement with Price et al. ().The information obtained from the CZ allows a relationship to be established between the calibration and performance of an algorithm, providing a statistical measurement of the probability to obtain a high-performance execution by using a calibration setting that lies in a specific region.Of course, the intention is to use the calibration settings that present the higher probabilities of providing good performance for the algorithm.In this case, the number of iterations needed to reach a specific error bound is defined as a performance measure of the algorithm.

Multi-start scheme
In Equation ( 9), the pressure at the burst point H  Burst position Burst pressure

Prototype description
The BDL methodology based on the extended DE algorithm is evaluated using datasets acquired from the prototype built   where the pressure waves arrive at the pressure sensors located at the ends of the pipeline are registered.Once the time instants are obtained, the wave speed is estimated as: where L r is the distance between pressure sensors, z r is the known valve position, V is the fluid velocity and Δt is the difference of times t 1 and t 2 in which the pressure drop, generated by the burst, arrives at the upstream and downstream sensors, respectively.In order to estimate the arrival time of the pressure wave to the sensors, a Discrete Wavelet Transform (DWT) method is used as in Romero-Delgado & Begovich ().With the only purpose to calibrate the pressure wave speed, a burst is produced by means of the valve 1 (located at 17.045 m) under the experimental conditions of the Table 3.As result of this experiment, a pressure wave speed value 355:914 m=s is obtained.
On the other hand, in order to find the suitable synthesis parameters F and C r of the algorithm to address the BDL problem, the CZ procedure is performed.In this way, the algorithm is executed in an iterative way on the acquired data from the pipeline prototype by producing a burst by the valve 1.To define the CZ, a cost function based on the total mean square error is considered for each execution: where K is the total number of samples in a database.It is expected that the proposed algorithm reaches a good estimation of the burst parameters at the end of each execution, allowing a maximum number of 500 iterations in each execution.In addition to this, an execution is considered successful when the algorithm reaches an error bound defined as f MSE 1 × 10 À12 in less than 20 iterations, as poor execution when the algorithm reaches the error bound in more than 20 iterations, but less than 500 and as failed execution when the algorithm does not reach the error bound, even with 500 iterations.The results obtained after 1000 executions are shown in Figure 3.The circles correspond to the successful executions and the squares the failed executions.The successful executions represent 18% of the total of the executions.As Figure 3 shows, when the parameters F and C r are specified arbitrarily, there is a great probability to obtain a bad performance.On the other hand, it can also be seen that there is an interval of F and

RESULTS AND DISCUSSION
In order to show the performance of the proposed algorithm in the BDL problem, this subsection presents two different experiments using real data from the prototype described above.In addition to this, the proposed algorithm is compared against the EKF, as one of the most widely used algorithms that presents good performance in the BDL problem (Delgado-Aguiñaga et al. ).
In Experiment 1, a single burst is emulated under con-    In order to validate the obtained estimations, Figure 6 shows the comparison between the measured flows and their estimations given by model ( 9), obtained by the incorporation of the best parameters found by the extended DE algorithm at each sample.shows the dataset used during this new experiment and   In order to validate the obtained estimations, Figure 9 shows the comparison between the measured flows and their estimations given by Equation ( 9) considering the best parameters found by the extended DE algorithm at each sample.
Table 7 summarizes the obtained results and their corresponding statistical values, considering a time interval between 130 s and 300 s.
Remember that in Experiment 2, the head pressure H 2 after the burst occurrence is not constant and considering this condition the average value is not expressed.However, in Figure 8(c), it is possible to see that the pressure variations are well estimated.
The results obtained from the performed experiments show that the proposed extended DE algorithm offers similar results than the EKF.Nevertheless, in the first experiment, the burst position presents a lower STD with the EKF, in contrast, in the second experiment, the STD is less with the DE algorithm.This fact shows that the data ).Therefore, in the current literature, different works are focused on developing leak and burst diagnosis algorithms.These algorithms have a common objective to detect and locate leaks or bursts automatically in transmission pipelines, with good precision and minimum invasion of the pipeline (e.g.Verde ; Wang et al. ; Kowalczuk & Gunawickrama ; Torres et al. ; Navarro et al. ; Brunone et al. ; Lizarraga-Raygoza et al. ).In general, there are two widely used approaches to address the burst diagnosis problem, the Transient Test Based (TTB) and the Fault Model Approach (FMA).The idea of the TTB is to exploit the information of the transient produced by burst (Wang et al. ; Brunone et al. , ), where the governing equations are expressed in terms of a Fourier series, and it is distinguished the deviation of the damping characteristics to the burst occurrence from the damping characteristics of the normal operation conditions.On the other hand, the FMA is based on the use of an analytical representation of the dynamic behavior of the flow rate and the pressure head in a pipeline, where the effect of a possible burst that affects the pipeline dynamics is considered within the model.This makes possible a direct estimation of the burst parameters (Santos-Ruiz et al. ).The extended Differential Evolution (DE) algorithm proposed to address the Burst Diagnosis and Location (BDL) problem in the present paper is focused in the FMA.In general, burst detection and localization algorithms based on the FMA are developed under different design considerations and mathematical assumptions.However, in real applications, these considerations are not always fulfilled, causing a reduction in the algorithms reliability and precision.Then, it is convenient to use algorithms developed under different estimation techniques that can work in parallel (analytical redundancy) in order to offer a more informed and reliable diagnosis (Carvajal-Rubio & Begovich ; Santos-Ruiz et al. ).This fact motivates the development of new techniques to address the burst detection and location problem in transmission pipelines.In view of the above, the main contribution of this paper is the development of a new BDL algorithm based on an extended DE algorithm and under FMA, which presents characteristics such as reliability, easy implementation and simple configuration.
step procedure to improve the adaptation of the classical DE algorithm in the BDL problem.These steps are: first, the boundaries of the search space are carefully proposed as a set of restrictions of the optimization problem according to the physical characteristics and mechanical properties of the pipeline.The boundaries of the search space play a fundamental role in the performance of the DE algorithm.Thus, by providing informed bounds, as a priori knowledge, the DE algorithm can reach a good feasible solution in a short time.Second, to select the appropriate synthesis parameters for the DE algorithm, a training step is performed using the method presented by Pérez-González et al. ().This method consists in evaluating the performance of the algorithm in a particular case, under an iterative execution that uses different synthesis parameters.The information acquired from the iterative execution allows to find the values of the parameter settings where few iterations are needed to approximate a solution to the problem.On the other hand, the BDL problem as an optimization task requires a special consideration due to its dynamic characteristics.Therefore, the proposed algorithm includes, as its third strategy, a multistart structure that allows to track the pressure dynamical variations in the burst point.The proposed burst diagnosis method based on the extended DE algorithm is implemented and evaluated in a pipeline prototype installed at Cinvestav, Guadalajara, México.As it will be shown, the obtained results illustrate the effectiveness of the proposed algorithm compared with the Extended Kalman Filter (EKF), which is widely used for the burst diagnosis problem.The paper is organized as follows: Preliminaries section presents the pipeline dynamic model, introduces the BDL problem, and describes the classical DE algorithm.In the Burst parameter identification strategy section, the improved DE algorithm to tackle the BDL problem solution is presented.The Experiments and results section presents and discusses the results obtained.Finally, the respective conclusions are stated.PRELIMINARIES This section introduces the burst detection and localization problem, starting with a brief description of the pipeline fluid dynamic equations, and afterwards the model considered in this work to address the BDL problem is introduced.Also, a brief description of the classical DE algorithm is presented.
[m=s], A is the cross-sectional area of the pipeline [m 2 ], Q is the flow rate [m 3 =s], z is the spatial coordinate along the pipeline [m], g is the gravitational acceleration [m=s 2 ], H represents the pressure head [m], t is the temporal coordinate [s], J w represents the friction factor head losses [dimensionless], and D is the diameter of the pipe [m].The pressure wave speed can be approximated by different formulations; however, the precision of the estimations highly depends on the mechanical properties of the material, which change with the pipeline deterioration.Furthermore, in plastic materials, the pressure wave speed estimations present variations in comparison with the values suggested in classical textbooks (Evangelista et al. ).Therefore, different authors as Meniconi et al. () and Soares et al. () propose experimental procedures in order to obtain more suitable pressure wave speed estimations.

Figure 1 |
Figure 1 | Spatial discretization of a pipeline in two sections.
) where r 1 and r 2 are positive integers from a uniform random distribution in the interval [1, NP], F is the mutation factor, which usually ranges in the interval [0,2].In general, the ability to perform local searches (exploitation) is achieved through values of F near to 0 and, on the contrary, the global search (exploration) is
2 (k) corresponds to a state of the model, which turns the optimization task into a dynamic optimization problem.That is, the proposed algorithm must track the dynamic behavior of H 2 (k) at each instant of the sampling time.In order to do that, the use of a multi-start scheme is proposed (Branke ), that allows consideration of the evaluation of model (9), in each sampling time k, as an independent optimization problem.So, for the k-th sample, there exists a set of values {z L (k), H 2 (k), λ(k)} that minimizes (14), therefore, the DE algorithm estimates these values.Before the dataset that corresponds to the sampling time k þ 1 arrives, the DE scheme as a set of independent optimization tasks.Finally, the proposed DE algorithm to address the BDL problem can be illustrated by the pseudocode presented in the Appendix section (available with the online version of this paper).
coefficient (C d ) of 0.9.Second, the pressure wave speed is determined.Despite the fact that valid results of burst localization have been obtained in publications such as Delgado-Aguiñaga et al. (), Delgado-Aguiñaga & Begovich (), and Lizarraga-Raygoza et al. (), in this work attempting to find improved results, the pressure wave speed is experimentally determined as in Pérez-González et al. ().In the experiment, a burst is produced by a valve with known position and then the difference between instants of time iterations needed to reach the error bound is 14, obtained by using F ¼ 0:53 and C r ¼ 0:85.Once the CZ is determined, the parameters setting for the implementation of the DE algorithm is taken as NP ¼ 30, since for many applications, the population size is taken as NP ¼ 10d, as it is reported in Price ().The control parameters C r ¼ 0:53, F ¼ 0:85 obtained from the CZ are used, and a total of 20 iterations per sample are considered.Finally, in order to detect the presence of a burst, and in consequence, start the algorithm, the burst detection alarm is established when the expression|Q 1 (k) À Q 2 (k)j !0:5 × 10 À4becomes true, in order to avoid the uncertainty in the flow measurement introduced by the noise level of the sensors.It is important to point out that the model-based burst diagnosis methods present high dependence to some pipeline stant flow conditions, by means of the valve 2, located at 49:895 m at time t L ¼ 104 s.The mean temperature during the experiment is 24.31 C and the mean viscosity is 9:0641 [unitless].The pressure head and flow rate at the ends of the pipeline are shown in Figure 4 and the experimental conditions are presented in Table 4.After the burst occurrences, two friction factors and Reynolds numbers are estimated with the input and output flows.In this experiment, f(Q1) ¼ 0:015513 and f(Q2) ¼ 0:015786.On the other hand, we have R e (Q 1 ) ¼ 2:0263 × 10 5 and R e (Q 2 ) ¼ 1:8529 × 10 5 .The results obtained from Experiment 1 with the extended DE algorithm and with the EKF are shown in

Figure 5 :
Figure 5: (a) presents the real burst position z L and its estimation with the proposed DE algorithm and the EKF; the real pressure head H 2 and its estimation are shown in (b), and finally, the magnitude factor λ and its estimation are presented in (c).As it is shown in Figure 5, the burst value estimations obtained from the extended DE algorithm and the EKF are close to the real values.

Figure 3 |
Figure 3 | CZ for the BDL problem.

Figure 4 |
Figure 4 | Flow rate and pressure head at the ends of the pipeline (Experiment 1).
considering a time interval between 130 s and 300 s.A second experiment is performed, in order to highlight the effectiveness of the proposed method.A new burst is emulated using again the valve 2 (49.895 m) at time t L ¼ 40 s, but also, producing variations in the operation point of the pump, by means of the frequency driver.Figure7

Figure 5 |
Figure 5 | Burst values and their estimations in Experiment 1.

Figure 6 |
Figure 6 | Flow measures and their estimations in Experiment 1.

Figure 7 |
Figure 7 | Flow rate and pressure head at the ends of the pipeline (Experiment 2).

Figure 8
Figure 8 shows the estimation results the second experiment.The plot (a) presents the estimation of z L , where the estimated position of the burst point does not change in spite of the variations in the operating point of the pump.In this way, the estimation of z L maintains its behavior in the neighborhood of the real value with a displacement of 1.45%.The same behavior is presented by the magnitude factor, presenting a displacement of 8.23% with respect to the real value of 1:326 × 10 À4 as Figure 8(c) shows.The case of the pressure head at the leak point is different from the previous two, since H 2 changes with respect to the variations of the operation point, despite this, the proposed extended DE algorithm can follow the dynamic variations.

Figure 8 |
Figure 8 | Burst values and their estimations in Experiment 2.

Figure 9 |
Figure 9 | Flow measures and their estimations in Experiment 2.
[χ {i,1,g}, χ {i,2,g} , . . ., χ {i,d,g} ] T , where i ¼ 1, 2, . . ., NP is an integer index that distinguishes each individual from the other members of the NP-size population, g is the generation index and d is the dimension of each individual, given by the number of parameters to be optimized, that is, each (Price et al. ; Coello et al. ); furthermore, it only needs two control synthesis parameters, which makes it easy to configure.In general, the DE algorithm is used to address problems in which the decision variables are real numbers and this feature is suitable to address the BDL problem.In the DE algorithm, the candidate solutions play the element of each individual corresponds to an unknown parameter of the optimization problem.The initial population g ¼ 0 is generated by a uniform random process, settled in the interval {χ low j χ i,j,0 χ j up } for j ¼ 1, 2, . . ., d.After the initialization, the DE process becomes an iterative sequence of evolutionary operations based on mutation, recombination, and selection (Zhang & Sanderson ).

Table 1 |
Unknown burst parameters

Table 2 |
Pipeline physical parameters

Table 3 |
Experimental test conditions to pressure wave speed estimation parameters, such as the friction and diameter, then it is necessary to periodically calibrate the pipeline parameter.

Table 5
presents the average, root mean square error (RMSE) and standard deviation (STD) of each estimated value obtained with the extended DE and EKF algorithms,

Table 4 |
Mean values of the first experiment

Table 5 |
Estimated parameters in Experiment 1

Table 6
presents the experimental scenario.

Table 6 |
Mean values of the second experiment

Table 7 |
Estimated parameters in Experiment 2