Detection of Multiple Leaks in a Natural Gas Pipeline Using

Subscriber access provided by Nanyang Technological Univ

Article

Detection of Multiple Leaks in a Natural Gas Pipeline using Observer and Mixed-Integer Partial Differential Equation-Constrained Optimization Xinghua Pan, Weiting Tang, and Muhammad Karim Ind. Eng. Chem. Res., Just Accepted Manuscript • DOI: 10.1021/acs.iecr.7b00527 • Publication Date (Web): 18 Sep 2017 Downloaded from http://pubs.acs.org on September 24, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Industrial & Engineering Chemistry Research is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Detection of Multiple Leaks in a Natural Gas Pipeline using Observer and MixedInteger Partial Differential Equation-Constrained Optimization Xinghua Pana, Weiting Tangb, M. Nazmul Karima* a: Artie McFerrin Department of Chemical Engineering, Texas A&M University, 3122 TAMU, College Station, Texas 77843, United States b: Department of Chemical Engineering, Texas Tech University, 6th Street and Canton Ave, Lubbock, Texas 79409, United States * Corresponding Author: [email protected]; Phone: +1 979 845 9806; Fax: +1 979 845 3266

ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 34

ABSTRACT: Software-based methods are common approaches for detecting faults in chemical processes. In this paper, new software-based methodologies were developed for locating multiple leaks in a natural gas pipeline. Two types of multiple leaks, subsequent and simultaneous multiple leaks, from a natural gas pipeline were studied, separately. For both subsequent and simultaneous leaks, case studies with two leak occurrences were demonstrated using MATLAB® simulation. For detecting and locating subsequent multiple leaks, an unknown input observer was designed and applied, which was modified from our previous study. New optimization methodology for locating simultaneous multiple leaks was demonstrated. Leak locations were estimated by solving a nonlinear global optimization problem. The global optimization problem contained constraints of partial differential equations, integer, and continuous variables.

A new

discretization approach was proposed and demonstrated, which required significant less computation effect comparing to conventional global optimization algorithms. Key Words: Multiple fault detection; Unknown input observer; Global optimization; Partial differential equation


Page 3 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1. INTRODUCTION Algorithms for detection and location estimation of a single leak in a pipeline have been studied extensively, using both isothermal and non-isothermal modeling for liquid/gas transportation.1–4 Detection of multiple leaks from a natural gas pipeline have not been demonstrated. Multiple leaks from a pipeline can be categorized into subsequent multiple leaks and simultaneous multiple leaks. The term ‘simultaneous’ indicated that multiple leaks can occur at the same time. The term ‘subsequent’ indicated that leaks can occur one after another. Detection of subsequent two leaks from a water pipeline has been demonstrated using multiple observers by Verde et al.5,6 Major difference between a liquid and a gas pipeline is the change of gas phase with thermal properties. Process models for isothermal (liquid) and non-isothermal (gas) transportation have different structures due to the presence of thermal-related parameters.7–10 We showed the effect of temperature and pressure change on the flow rate of natural gas in a pipeline in our previous study.7 Thus, detection of multiple leaks from a natural gas pipeline needs to consider the effect of thermal properties. In this paper, we extended the work in our previous study to detect subsequent multiple leaks, which can reduce the effect of thermal properties on leak detection.7 The above-mentioned or other methods for detection of a single leak or subsequent multiple leaks cannot be applied directly to detect simultaneous multiple leaks because they use steady-state value of the boundary flow rate.4,5,11–13 For a pipeline with simultaneous multiple leaks, multiple leaks can be mistakenly interpreted as a single leak at a different location when the steady-state value of the boundary flow rate is used.14 Thus, it is necessary to investigate the dynamic response of the flow rate to locate simultaneous multiple leaks.15 Detection of simultaneous multiple leaks from a water pipeline has been showed by Verde et al.14,16 However,



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

this method cannot be applied directly to a non-isothermal natural gas pipeline due to the high complexity of the non-isothermal natural gas flow model. The isothermal model cannot reproduce the non-isothermal dynamic response of boundary flow rates under leaks, which is due to the effect of thermal properties in a natural gas pipeline. This indicates that detection of multiple simultaneous leaks from a natural gas pipeline requires application of a non-isothermal model. The estimation of the location of simultaneous multiple leaks in a natural gas pipeline can be solved as a nonlinear mixed-integer partial differential equation-constrained global optimization problem with both continuous and integer variables. The continuous variables are sizes of the leaks. The integer variables are the leak locations, which are sections from discretization of a pipeline. The objective function of the global optimization is to minimize the error between measurement and estimation of the boundary flow rates. The existence of the integer variables, which are the leak location, hinders the application of gradient-based optimization algorithm. Partial differential equation-constrained optimization has been studied by many researchers.17–19 However, due to the existence of integer variables, current global partial differential equation-constrained optimization algorithms cannot be directly applied to solve the problem. Besides, approach such as linearization of the nonlinearity in this problem will introduce optimization error, which will reduce the accuracy of the result. Available derivative-free global optimization algorithms such as genetic algorithm20, Monte Carlo simulation21,22, or particle swarm23 need extensive computation iterations and time, which are not suitable for an urgent task for reducing financial lost and potential fire hazard.24,25 The mixed-integer PDE constrained nonlinear optimization is an ongoing research topic. The theoretical and computational challenges such as functional analysis and large number of


Page 4 of 34

Page 5 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


variables remain as issues for the optimization methods. Different strategies have been reported to solve the mixed-integer PDE constrained global optimization problems.26 Total discretization, such as discretization of the PDE with small discretization steps, can lead to large number of equations for the system. 27,28 For the integer variables, nonlinear branch and bound method for global optimization can lead to very large branch and bound trees. Another strategy is to reformulate the PDE constrained or ODE constrained problems into other problem with available solvers.29,30 However, some of the results indicated the need for a large and numerically difficult reformulation. Other algorithms such as relaxation methods have also been studied to change the optimization problem with integer variables into a continuous problem.31,32

However, the

solution assumes certain degree of the smoothness of the system. Certain system with inherent discontinuity may not be applied.31 In this paper, to solve this optimization problem in a time-efficient manner, we improved the total discretization method by a new global optimization algorithm with sequential zooming discretization. The new algorithm can significantly reduce the size for the optimization problems. Briefly, the algorithm first discretizes the pipeline with large mesh size (the integer variable), and all possible pairs of leak locations are screened using Newton’s method. Leak locations with smallest objective function value are selected. The pipeline is further discretized around the selected leak locations and Newton’s method is repeated. After repeated discretization around the selected locations, more and more accurate location approximations can be achieved. Due to noticeable effect of environment temperature on boundary flow rate, prior knowledge of temperature, either from measurement or estimation, was required for the



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 34

optimization algorithm to locate multiple simultaneous leaks. In the study of locating simultaneous multiple leaks, no process disturbance was considered. The proposed optimization algorithm can also be applied to solve other engineering problems with mixed-integer partial differential equation constraints.

2. METHODOLOGIES 2.1 Detection and location estimation of subsequent multiple leaks: design of observer Basing on the non-isothermal and isothermal models introduced in our previous study, a same procedure can be applied for designing an observer for detecting subsequent multiple leaks.7 The procedure is briefly summarized here. The first step is to reduce the isothermal flow model to a linear model as shown in Equation (1). qin (t)=qin.st. + ∑ni=0 ai qusage (t-i)

(1a)

qout (t)=qout. st. + ∑ni=0 bi qusage (t-i)

(1b)

qusage represents the consumer usage in the pipeline. For a straight pipeline without any consumer usage, the parameter associated with the consumer usage becomes zero. qin (t) and qout t are the time-variant inlet and outlet flow rate due to the existence of consumer usage. Parameter i represents different time point. qin.st. and qout.st. are the inlet and outlet flow rate at steady state.


Page 7 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The second step is to design an unknown input observer and solve the parameters of the observer. The detailed development of the unknown input observer and model reduction are explained in our previous study.7 The reduced linear equation can be generalized as the following equation. xt= ∑ni=0 Ai xt-τi + ∑ni=0 Bi ut-τi +Wwt +

(2a)

yt=Cxt+Md(t)

(2b)

in which, xt is the state, yt is the measurement, ut-τi is the input (consumer usage), d(t) is measurement noise which is assumed white noise, and wt is the unknown input. τi is constant time intervals. W is the parameter matrix for the unknown input, which was determined by the simulation study on effect of the unknown inputs in the system. An unknown input observer can be designed as shown below. z(t)= ∑ni=0 Fi zt-τi + ∑ni=0 TBi ut-τi + ∑ni=0 Gi yt-τi x t=zt+Ny(t)

(3)

in which, F, T, B, and G are process parameters for the unknown input observer that need to be determined. x t is the estimation of the state and y(t) is the measurement from the system. Process disturbances such as temperature change and pressure oscillation at pump station were treated as unknown inputs. The effects of process disturbance were studied in our previous study and applied for constructing the observer. 7 The third step is to apply the unknown input observer to estimate the flow rate in a pipeline, and compare the estimated flow rate with the measurement. Detection of leaks from a pipeline is based on the analysis of a residual value. A residual value is the difference in flow



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 34

rate between measurement and estimation from the observer, as shown in our previous study.7 A leak is identified when the residual value exceeds a preset threshold, according to the noise level of the system (0.5% in our case). Multiple leaks can be detected subsequently according to the residual value which will be shown in the results section. The fourth step involves identifying the leak from the residual signal and locating the leak in a pipeline. Equation (4) is applied to calculate the leak location.12

XL =L/(1-

qin.mea -qin

(qout.mea -qout)

)

(4a)

X is the estimated leak location and L is the total length of the pipe. qin.mea and qout.mea are the measured inlet and outlet flow rate. qin and qout are inlet and outlet flow rate from the reduced model in Equation (1). The convergence criterion for the estimation of the location of a leak can be calculated using Equation (4b). In the equation, is the estimate of the location of the leak at the time

point i, is the estimate of the location of the leak at time point n+1. The estimation of the

location of the leak converged when s is smaller than a threshold. In this case, the threshold is set as 2% of the length of the pipe. n was set at 10.

s = ∑| − |

(4b)

After a leak was identified basing on the residual value, the steady state values qin.st. and

qout. st. in Equation (1) were updated, so qin (t) and qout t in Equation (1) will match the flow rate measurement at steady state. The updated reduced linear model with qin.st. and qout. st. can also account for any temperature change and pressure change. The updated equations were used for detecting possible more leaks and estimating the leak location.


Page 9 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


2.2 Location estimation of simultaneous multiple leaks: development of optimization algorithm New global optimization algorithm is proposed to locate multiple simultaneous leaks in a natural gas pipeline. The objective function of the optimization problem is to minimize the estimation error from dynamic response of the boundary flow rate. The optimization problem assumes a prior knowledge of temperature and does not consider oscillation of pump pressure. There are both integer variables (leak location) and continuous variables (leak size) in the optimization problem. The optimization problem is formatted in Equation (5).

minu (q ) ∑N -y ) Q y-y i=0 (y i i i i T

L

∂P

s.t.

∂t x=0

=0,

∂P ∂t x=L

=0; qL1x , qL2x ∈ Θ; qL1s , qL2s ∈Γ

q yi = q1 ; u=qL1x qL2x

qL1s

0

1 ∂q

∂P ∂t

∂q ∂t

=

4UT"Tg q dT T∂Z q dP ACP dx + (Z∂T + 1)APZRT dx# D

A∂x A∆xqL + ZC ∂T + TC ! P P 2DA3 P2

=−A

1 ∂Z

1

1

fq3 z2 R2 T2

!ZRT 2 ∂P $∂T% 2 ZC ∂T TC # Z RT Z CP P P P ∂Z

1

∂P ∂x

−

AP ZRT

gsinθ −

fq2 2DAP

∂Z 2

T

ZRT −

2 ∂Z

1

1 qL q $ % ZRT A ∆x P



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

ZRT 2 ∂P ∂t Z RT 1

∂T ∂t

=

P ∂Z ∂P

2 ∂T + Z RT ZRT2 P ∂Z

P

+

1 ∂q A∂x P ∂Z P 2 ∂T + Z RT ZRT2

+

qL

1

A∆x P ∂Z + P Z2 RT∂T ZRT2

Page 10 of 34

(5)

yi is the measured boundary flow rate at time step i, and y i is the estimated boundary flow rate at time step i. N refers to the total time of dynamic response. Boundary conditions (fixed boundary pressure) are applied to the optimization problem and assumed validated during the dynamic response. qL1x and qL2x refer to the two leak locations, which are constrained integer

variables from discretization of a pipeline. qL1s represents the size of leak one. The other leak

size can be calculated by subtracting q ' from the total leak size (considering a two-leak case). The non-isothermal equations in the constraints were explained in our previous study.7 To avoid massive computation time and improve the calculation efficiency, a new sequential zooming discretization method was developed and applied in our solution. The optimization problem was solved using Hessian-based Newton's method. Penalty method was applied for the boundary of the continuous variables in the constraints. The goal of the optimization problem is to locate multiple leaks and their corresponding leak sizes. The accuracy of location estimation depends on the mesh size of the discretization. The sequential zooming discretization method is illustrated in Figure 1. qL1x and qL2x are integer variables (leak location) that need to be searched as accurate as possible. A process for detecting two simultaneous multiple leaks is demonstrated below.


Page 11 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 1. Illustration of the optimization algorithm

As shown in Figure 1, a middle point was first calculated to locate two leaks in two sections along the pipe using Equation (6), in which qin.st. and qout.st. are inlet and outlet flow rate before leak occurrence at steady state. qin and qout are inlet and outlet flow rate after leak occurrence. This is a procedure to reduce the computation time for two simultaneous leaks, which is not necessarily applied to more simultaneous leaks.

Lmiddle =(1-

qin.st -qin

(qout.st -qout)

)

(6)

The second step is to discretize the pipeline into several large sections and perform a global search of leak locations from all possible combinations of discretized sections as shown in Figure 1. Two possible discretization sections (leak locations) from both sides (separated by the calculated middle point) were chosen as a leak location pair (for the two-leak case). In the second step, for each leak location pair, Newton’s method was used to search for leak sizes (continuous variables) to minimize the objective function. All the possible leak location pairs were screened and compared using the minimized objective function value. A first approximate



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 34

leak location pair was selected with the smallest numerical value of the objective function, as labeled using stars in Figure 1. Because the leak locations can lie near the boundary of each discretization, four pairs of possible leak locations with smallest objective function values were selected after screening all the possible combinations of the leak locations instead of one pair. In step three, further discretization around the chosen pairs of approximate leak locations is performed and the same procedure is repeated to calculate more accurate leak locations. The third step is repeated until satisfied estimation accuracy of leak locations is achieved basing on the mesh size of the discretization. The optimization was terminated when the step decrease of the objective function value is less than 0.1 for the first iteration of optimization. The termination criteria for other iterations is decreased accordingly. In the application of Newton’s method, only a couple optimization steps were needed, which were sufficient enough to select the most possible candidates, due to the relative small size of the leaks. We demonstrated two examples on location estimation of simultaneous multiple leaks with different leak locations and leak sizes. Both leaks are smaller than 5%. Newton’s method was applied to find the minimal value of the objective function in all these possible leak location pairs. Hessian matrix was calculated numerically in Newton’s method.

2.3 Simulation study For subsequent multiple leaks, artificial measurement data was generated using the nonisothermal model, with 0.5% process noise, as described in our previous paper.7 Two subsequent leaks were introduced into the system with both sizes of 1.5% of the stable flow rate, at the locations of 0.7 and 0.45 of the pipeline. The two subsequent leaks were introduced at 500 min and 650 min, separately.


Page 13 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


For simultaneous multiple leaks, two leaks were introduced simultaneously at different locations. The details of the leaks are described in Table 1. No process noise was introduced into the system for the simultaneous case study. ‘Artificial measurement’ was generated from the nonisothermal simulation with 74 discretization nodes.

Table 1. Case study of two simultaneous leaks Leak I location

Leak I size

( % of Length)

Leak II location

Leak II size

(% of Length)

Case I

0.736

4.0%

0.277

1.6%

Case II

0.418

3.6%

0.797

1.2%

3. RESULTS AND DISCUSSION 3.1 Detection of subsequent multiple leaks Several research papers have proposed methods to detect subsequent leaks from a water pipeline Due to the complicity of the natural gas flow model, the method proposed for water pipeline cannot been applied for detecting a natural gas pipeline when considering process noise in terms of temperature change or pressure drop. Furthermore, our observer method can be directly applied to detect three or more subsequent leaks without modifying the algorithm.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 34

Figure 2. Comparison of boundary flow rate between measurement and observer estimation

Figures 2 to 4 demonstrate the identification of subsequent leaks and estimation of the location of the leaks. The estimation for subsequent leaks starts with a model reduction method. Figure 2 shows the estimation of flow rate from observer compared to the original nonisothermal model. Figure 3 shows the value of the residual from state estimation, which is an indicator of leak occurrence. The jump of the residual value indicates a new leak. Figure 4 shows the estimation of locations of two subsequent leaks.

Figure 2 shows the comparison of observer estimation and measurement of flow rate from the simulation when the consumer usage is zero, without leak. The observer estimation starts from 150 minutes. As can be seen in the figure, the estimation from the unknown input observer overlaps with the measurement which indicates a good estimation of the boundary flow rate. Current available methods for process monitoring under process disturbances perform real-time modeling/simulation when the pressure and temperature information is available.


Page 15 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


However, our observer method provides easier way to reduce the effect of pressure and temperature drop without the measurements of pressure and temperature, which reduces the requirement of extra instrumentation.7 The residual value of two subsequent leaks is shown in Figure 3, which shows the detection of subsequent leaks. As can be seen in the figure, subsequent two increase of the residual value indicates two leaks at different time. Without changing any structure and parameter of the observer, the observer is able to identify more leaks using the residual value. Estimation of leak locations for the two leaks is shown in Figure 4. As shown in the figure, the estimation of leak location can converge within 20 minutes. Oscillation of the location estimation is due to the process noise.

Figure 3. Residual value of the observer for two subsequent leaks



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 34

Figure 4. Estimation of leak location for two subsequent leaks It is easy to extend the same principle to detect a third or fourth leak. Further examples of more leaks are not demonstrated here.

3.2 Detection and location estimation for simultaneous multiple leaks Two leaks with both sizes less than 5% were introduced into a straight pipeline at different locations as shown in Table 1. For both cases, the pipeline was discretized into 8, 16, 32, and 64 sections in the process of optimization. Four pairs of leak locations with the smallest numerical values of the objective function from each iteration of optimization are demonstrated in Table 2 and Table 3. The optimization problem contains an integer variable, which is a location-related variable for the discretization of the pipeline. Calculation of gradient for the integer variable is


Page 17 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


not feasible for this problem. To search for possible leak locations across the whole pipeline, a global optimization algorism is required.33–35 Global derivative-free optimization methods requires many sampling points, which will significantly delay the leak detection process.23,36 It is desired to have a faster approach to search for the integer variable in a time-efficient manner rather than the global derivative-free optimization. Studies have been reported for PDEconstrained mixed integer for linear programming.37 Adaptive discretization approach for PDEconstrained optimization have also been reported by Ziem et al.38

Application of PDE-

constrained optimization for transient optimization of gas network has also been reported by Steinbach.39 Following the procedure in Equation (6), a middle point was calculated using the steady state flow rate after leak occurrence. The calculation of middle point was only for the case of two leaks. The calculated middle point using Equation (6) was 0.42, so section 4 was chosen as a middle point. For the first iteration of optimization in case I, the pipeline was discretized into 8 sections, and 20 optimizations were performed to search for the leak location and sizes. The first iteration of optimization involved optimizations with two fixed integer variables containing the discretized pipeline section pairs of (1,4), (1,5), (1,6), (1,7), and (1,8). The first iteration of optimization also contained leak location pairs starting with section 2: (2,4), (2,5), (2,6), (2,7), (2.8); section 3: (3,4), (3,5),(3,6), (3,7), and (3,8); and section 4: (4, 4), (4,5), (4,6), (4,7), (4,8). Each optimization used Newton’s method to search two possible leak sizes (continuous variable) to minimize the objective function. To shorten the calculation time, initial guess of the continuous variable (one of the leak size) was set half of the total leak size, which can be measured by subtracting inlet flowrate using outlet flowrate.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 34

Table 2 and Table 3 show the optimization results for case I and case II. As shown in Table 2, section pairs of (3, 6), (3, 7), (4, 6) and (4, 7) were selected which has the smallest objective function value for the first iteration of optimization. For the second iteration of optimization, the pipeline was further divided into 16 sections, and 16 optimizations were performed near the selected sections in iteration 1. In iteration 3 and iteration 4, the pipeline was discretized into 32 and 64 sections, and 16 optimizations were performed for each iteration. In our study, we performed 4 optimization iterations, and further iterations of optimization with more discretization sections can be continued. Table 3 shows the optimization result for location estimation for case II, which followed the exact same procedure. For each iteration of optimization, after the integer variable (discretization section) was fixed, Newton’s method was applied to calculate the corresponding leak size. It is worth noticing that the partial differential equation-constraint in Equation (5) is not second-order differentiable, because the higher order term of the compressibility factor is not available. Both first and second order derivative were calculated numerically. Four lowest pairs of possible of leak locations were further discretized instead of one possible pair of leak location, which counts for the calculation error during the optimization procedure. This procedure will require more computation time, however, it will deliver more accurate results. The number of discretization nodes will affect the total number of optimization iterations. More discretization nodes lead to more accurate optimization results, but increase the computational time. When the discretization node is made smaller, more accurate estimation of the location leak can be obtained through the optimization process. The final leak location estimation is claimed to be inside the finest discretization node, after the sequential discretization steps. As long as a discretization node contains the leak, the objective function value from that discretization node is


Page 19 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


smaller than other cells without the leak. However, the specific location of a leak inside the discretization node is unknown at first. It can only be searched through optimization at a different level. Overall, it is a trade-off, and balance between accuracy of the estimation of the location of a leak and the computational time. The values of the minimal objective functions at different levels are different as shown in Table 2 and Table 3. Smaller discretization nodes will lead to smaller objective function value. This indicates that large discretization nodes constitute a starting search point for the estimation of the location of the leak. With more steps of discretization, the location of the leak can be estimated more accurately. In our demonstration, the pipeline is first divided into 8 nodes. After the first iteration of optimization, 4 of them were chosen for further discretization into 8. For every iteration, the chosen discretization node was divided into two nodes in next level of iteration.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 34

Table 2. Optimization results of leak detection for simultaneous leaks case I Iteration/ Discretization number I/8

II/16

III/32

IV/64

Number of optimizations

20

16

16

16

Pipeline section number value of objective function

Pipeline section number Section number

6

7

3

1.6519

0.6916

4

1.9201

0.9056

Section number

12

13

5

0.3156

0.6648

6

0.3415

0.4841

Section number

23

8

0.1806

0.1171

9

0.1481

0.0974

Section number

48

17

0.0763

18

0.0766


24

Page 21 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 3. Optimization results of leak detection for simultaneous leaks case II Iteration/ Discretization number I/8

II/16

III/32

IV/64

Number of optimizations

20

16

16

16

Pipeline section number value of objective function

Pipeline section number Section number

6

7

3

0.6543

0.6953

4

1.1743

1.0227

Section number

13

14

7

0.2207

0.1579

8

0.1289

0.2193

Section number

25

26

13

0.1010

0.1099

14

0.0782

0.1006

Section number

52

27

0.0600

28

0.0564

Table 4 summarizes the comparison between simulated leak location and location estimation from optimization. The estimated location was an average of the calculated location as shown in Table 2 and Table 3. The location of the leak is inside the final discretization node, but its exact location inside that discretization node is unknown. The final location of the leak location is given by the middle point of the discretization node, which is an estimated average. As can been seen in the table, after 4 iterations of optimizations, the estimated leak location is very close to the real location. With the assumption of a pipeline between two stations is 10 km, the average estimation error is 125 m. As demonstrated in case I, a total of 68 optimizations were



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 34

performed to search the leak locations as accurate as 1/64 of the pipeline length. In the nonisothermal model, the second order derivate is calculated numerically from the first order derivate. A total number of 1396 first order gradient calculation was performed in the system. The proposed work in the manuscript is to improve the ‘total discretization’ approach. The total discretization method plus branch and bound method can lead to a global optimization result, but with increased computational time. For example, to achieve the accuracy of 64 discretization nodes of a pipeline, two of the 64 discretization nodes have to be tested, which will lead to a total of 2016 optimizations calculations. After applying the middle line to reduce the computation time, a total of 975 optimizations are required to achieve the same estimation accuracy without the sequential zooming discretization method. Table 4. Comparison between simulated leak location and optimization result Location of Leak 1 ‘simulated leak / optimization result

Leak size

Location of Leak 2 ‘simulated leak / optimization result

Leak size

Case I

0.736 / 0.750

3.9%

0.277 / 0.273

1.7%

Case II

0.797 / 0.813

3.5%

0.418 / 0.429

1.3%

Detection of simultaneous multiple leaks is generally more difficult than detecting subsequent leaks. Due to the fact that boundary flow rates are the only online measurements, no prior knowledge can be obtained about the number of simultaneous leaks in the system. In this paper, we demonstrate a methodology to detect two simultaneous leaks.


Page 23 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


For detection of multiple subsequent leaks, the same observer method can be directly applied without modification. For example, for three subsequent leaks, we will see similar transient responses as seen in Figure 3 for the first and the second leak; it will be delayed by the time of occurrence. The detection of more subsequent leaks can be identified both spatially and temporally. For detection of more simultaneous leaks, the same optimization method can still be applied, although it is a rare situation. However, the computational time will be increased dramatically. The first step for detecting multiple simultaneous leaks is to assume a certain number of the leaks.

For example, for detection of three leaks, a ()* number of initial

optimization calculations need to be performed for eight discretization node, similar to Table 2. Further iterations with more numbers of optimization will need to be performed. The optimization is based on the dynamic response of the flow rate due to the occurrence of leaks. However, for the cases with more simultaneous leaks, resolution and accuracy of the flow rate measurements need to be improved dramatically, which might be a challenge in a real-world application.

4. CONCLUSIONS Detection and location estimation of subsequent and simultaneous multiple leaks from a natural gas pipeline was studied. For multiple subsequent leaks, extension of an unknown input observer with modification from our previous research was demonstrated. Global optimization method was proposed for detection of simultaneous multiple leaks. The global optimization problem was solved using a sequential zooming discretization approach, under the constraints of partial differential equations, integer variable, and continuous variable. The sequential zooming discretization approach is effective to search for leak locations after several iterations of optimizations. However, the optimization method is not able to handle process disturbance from



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 24 of 34

temperature or pressure change, which needs further investigation. Detection of multiple subsequent leaks can be directly applied based on the two-leak case.

Detection of more

simultaneous leaks will require much more computation effort and sensor accuracy and resolution.

ACKNOWLEDGEMENTS This work is supported by T. Michael O’Connor Chair II Endowment Fund and Energy Institute Fellowship from Texas A&M University.

REFERENCES (1)

Murvay, P. S.; Silea, I. A Survey on Gas Leak Detection and Localization Techniques. J. Loss Prev. Process Ind. 2012, 25, 966.

(2)

Hauge, E.; Aamo, O.; Godhavn, J. M. Model-Based Monitoring and Leak Detection in Oil and Gas Pipelines. SPE Proj. Facil. Constr. 2009, 4, 53.

(3)

Benkherouf, A.; Allidina, A. Y. Leak Detection and Location in Gas Pipelines. IEE Proc. D Control Theory Appl. 1988, 135, 142.

(4)

Reddy, H. P.; Narasimhan, S.; Bhallamudi, S. M.; Bairagi, S. Leak Detection in Gas Pipeline Networks Using an Efficient State Estimator. Part-I: Theory and Simulations. Comput. Chem. Eng. 2011, 35, 651.

(5)

Verde, C.; Visairo, N. Bank of Nonlinear Observers for the Detection of Multiple Leaks in a Pipeline. In Proceedings of the 2001 IEEE International Conference on Control Applications; 2001; pp 714–719.


Page 25 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(6)

Verde, C.; Morales-Menendez, R.; Garza, L. E.; Vargas, A.; Velasquez-Roug, P.; Rea, C.; Aparicio, C. T.; De La Fuente, J. O. Multi-Leak Diagnosis in Pipelines - A Comparison of Approaches. In 7th Mexican International Conference on Artificial Intelligence Proceedings of the Special Session, MICAI 2008; 2008; pp 352–357.

(7)

Pan, X.; Tang, W.; Raftery, J. P.; Karim, M. N. Design of an Unknown Input Observer for Leak Detection under Process Disturbances. Ind. Eng. Chem. Res. 2017, 56, 989.

(8)

Osiadacz, A. J.; Chaczykowski, M. Comparison of Isothermal and Non-Isothermal Pipeline Gas Flow Models. Chem. Eng. J. 2001, 81, 41.

(9)

Brouwer, J.; Gasser, I.; Herty, M. Gas Pipeline Models Revisited: Model Hierarchies, Nonisothermal Models, and Simulations of Networks. Multiscale Model. Simul. 2011, 9, 601.

(10)

Abbaspour, M.; Chapman, K. S. Nonisothermal Transient Flow in Natural Gas Pipeline. J. Appl. Mech. 2008, 75, 31018.

(11)

Wang, S.; Carroll, J. Leak Detection for Gas and Liquid Pipelines by Online Modeling. Energy Process. 2007, 39, 21.

(12)

Verde, C. Multi-Leak Detection and Isolation in Fluid Pipelines. Control Eng. Pract. 2001, 9, 673.

(13)

Liu, M.; Zang, S.; Zhou, D. Fast Leak Detection and Location of Gas Pipelines Based on an Adaptive Particle Filter. Int. J. Appl. Math. Comput. Sci. 2005, 15, 541.

(14)

Verde, C.; Visairo, N.; Gentil, S. Two Leaks Isolation in a Pipeline by Transient Response. Adv. Water Resour. 2007, 30, 1711.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(15)

Colombo, A. F.; Lee, P.; Karney, B. W. A Selective Literature Review of Transient-Based Leak Detection Methods. J. Hydro-Environment Res. 2009, 2, 212.

(16)

Verde, C.; Molina, L.; Torres, L. Parameterized Transient Model of a Pipeline for Multiple Leaks Location. J. Loss Prev. Process Ind. 2014, 29, 177.

(17)

Biegler, L.; Ghattas, O. Large Scale PDE Constrained Optimization: An Introduction; Springer-Verlag Berlin Heidelberg: Berlin, 2003.

(18)

Protas, B. Adjoint-Based Optimization of PDE Systems with Alternative Gradients. J. Comput. Phys. 2008, 227, 6490.

(19)

Cioaca, A.; Alexe, M.; Sandu, A. Second-Order Adjoints for Solving PDE-Constrained Optimization Problems. Optim. Methods Softw. 2012, 27, 625.

(20)

Storn, R.; Price, K. Differential Evolution -- A Simple and Efficient Heuristic for Global Optimization over Continuous Spaces. J. Glob. Optim. 1997, 11, 341.

(21)

Hesselbo, B.; Stinchcombe, R. B. Monte Carlo Simulation and Global Optimization without Parameters. Phys. Rev. Lett. 1995, 74, 2151.

(22)

Vanderbilt, D.; Louie, S. G. A Monte Carlo Simulated Annealing Approach to Optimization over Continuous Variables. J. Comput. Phys. 1984, 56, 259.

(23)

Rios, L. M.; Sahinidis, N. V. Derivative-Free Optimization: A Review of Algorithms and Comparison of Software Implementations. J. Glob. Optim. 2013, 56, 1247.

(24)

Rees, T.; Dollar, H. S.; Wathen, a. J. Optimal Solvers for PDE-Constrained Optimization. SIAM J. Sci. Stat. Comput. 2010, 32, 271.


Page 26 of 34

Page 27 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(25)

Hinze, M.; Pinnau, R.; Ulbrich, M.; Ulbrich, S. Optimization with PDE Constraints; Springer Science & Business Media: New York, 2009.

(26)

Steinbach M. C., On PDE solution in transient optimization of gas networks, Comput. Appl. Math. 2007, 203, 345.

(27) Geißler B.; Morsi A.; Schewe L. A new algorithm for MINLP applied to gas transport energy cost minimization. In Facets of combinatorial optimization, pages 321–353. Springer, Heidelberg, 2013. (28) Göttlich S.; Herty M.; Ziegler U. Modeling and optimizing traffic light settings in road networks. Comput. Oper. Res. 2015, 55, 36. (29) Gerdts M. A variable time transformation method for mixed-integer optimal control problems. Optim. Control Appl. Meth. 2006, 27, 169. (30) Baumrucker B.T.; Biegler L.T. MPEC strategies for cost optimization of pipeline operations. Comput. Chem. Eng. 2010, 34, 900. (31) Hante F.M.; Sager S. Relaxation methods for mixed-integer optimal control of partial differential equations. Comput. Optim. Appl. 2013, 55, 197. (32) Hante F.M. Relaxation methods for hyperbolic PDE mixed-integer optimal control problems. Optim. Control Appl. Meth. 2017, 1 (33)

Boukouvala, F.; Misener, R.; Floudas, C. A. Global Optimization Advances in MixedInteger Nonlinear Programming, MINLP, and Constrained Derivative-Free Optimization, CDFO. Eur. J. Oper. Res. 2016, 252, 701.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(34)

Romeijn, H. E.; Schoen, F. Handbook of Global Optimization; Kluwer Academic Publishers: Boston, 2002; Vol. 62.

(35)

Floudas, C. A.; Gounaris, C. E. A Review of Recent Advances in Global Optimization. J. Glob. Optim. 2009, 45, 3.

(36)

Pham, N.; Malinowski, A.; Bartczak, T. Comparative Study of Derivative Free Optimization Algorithms. IEEE Trans. Ind. Informatics 2011, 7, 592.

(37)

Armin Fügenschuh,Stefan Vigerske,Henning Homfeld, H. S. Mixed-Integer Nonlinear Problems in Transportation Applications. In 2nd International Conference on Engineering Optimization; 2010; pp 1–10.

(38)

Ziems, J. C.; Ulbrich, S. Adaptive Multilevel Inexact SQP Methods for PDE-Constrained Optimization. SIAM J. Optim. 2011, 21, 1.

(39)

Steinbach, M. C. On PDE Solution in Transient Optimization of Gas Networks. J. Comput. Appl. Math. 2007, 203, 345.


Page 28 of 34

Page 29 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


TOC



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 1. Illustration of the optimization algorithm


Page 30 of 34

Page 31 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 2. Comparison of boundary flow rate between measurement and observer estimation 677x244mm (96 x 96 DPI)



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 3. Residual value of the observer for two subsequent leaks 677x244mm (96 x 96 DPI)


Page 32 of 34

Page 33 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 4. Estimation of leak location for two subsequent leaks 99x92mm (300 x 300 DPI)



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

TOC


Page 34 of 34

Detection of Multiple Leaks in a Natural Gas Pipeline Using

Recommend Documents