# Robust Kalman Filter-Based Dynamic State Estimation of Natural Gas Pipeline Networks

Authors: Liang Chen; Peng Jin; Jing Yang; Yang Li; Yi Song
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

To obtain the accurate transient states of the big scale natural gas pipeline networks under the bad data and non-zero mean noises conditions, a robust Kalman filter-based dynamic state estimation method is proposed using the linearized gas pipeline transient flow equations in this paper. Firstly, the dynamic state estimation model is built. Since the gas pipeline transient flow equations are less than the states, the boundary conditions are used as supplementary constraints to predict the transient states. To increase the measurement redundancy, the zero mass flow rate constraints at the sink nodes are taken as virtual measurements. Secondly, to ensure the stability under bad data condition, the robust Kalman filter algorithm is proposed by introducing a time-varying scalar matrix to regulate the measurement error variances correctly according to the innovation vector at every time step. At last, the proposed method is applied to a 30-node gas pipeline network in several kinds of measurement conditions. The simulation shows that the proposed robust dynamic state estimation can decrease the effects of bad data and achieve better estimating results.


---

## Body

## 1. Introduction

In comparison with the traditional coal-fired power units, gas-fired electric generators can respond to the power load fluctuation rapidly, enhancing the operating flexibility of electrical energy systems [1–3]. This may help to improve security of power system with large-scale renewable energy. The random change of the natural gas consumptions due to the uncertainties of renewable energies makes it essential for obtaining the accurate dynamic states just like the pressures and mass flow rates of the natural gas pipeline networks to ensure the security and optimal operation of the integrated energy system containing electric powers and natural gases [4–7]. To capture the states, the pipeline networks have to be equipped with a mount of measuring devices, which requires heavy investments. Even so, it is impossible to install sensors at every node of the networks and obtain all states. On the other hand, the measuring devices experience random errors and bad data unavoidably, so the measured data cannot be applied directly to the leak detections [8, 9] and control problems [10] before state estimations.In recent years, some research works about state estimations for natural gas pipeline networks have appeared [11–15]. These works are based on the nonlinear partial differential equation (PDE) describing the characteristics of transient gas flow [16, 17]. To linearize the PDEs of gas systems, in [18, 19], the dynamic states are redefined as deviations from steady states for linearization purpose. In [20], an iterative linear approximation is proposed, which is used to solve optimal power flow problem, rather than the linearization of PDEs. In [11], a state estimation method for natural gas pipe lines by using the linearization of PDEs [20–22] is proposed. To improve the estimating performance, the state estimation with a pair of Kalman filter-based estimators running in parallel is carried out. Some researchers established the state estimation model applying the PDEs directly, and solved the model by nonlinear algorithms. In [12], the extended Kalman filter (EKF) is chosen to design an efficient observer for natural gas transmission system, and then an algorithm is proposed to handle the discontinuities that appear in the dynamic model of a gas transmission networks. In [13], a two-step Lax-Wendroff method is used for the discretization of the PDEs to obtain finite-dimensional discrete-time state-space representations, and the particle filter that fits for the nonlinear filter problems is applied to estimate the transient states. These methods can capture the accurate states in the transient processes. However, as the scale expands, the centralized implementation of the Kalman filter has severe limitations such as tuning, scalability, unacceptable calculating loads, and lack of robustness in the case of sensor failures. Aiming at this problem, a strategy for the distributed and decentralized state estimation of state variables in big scale systems is proposed in [14]. In addition, the algorithm for a joint state and parameter estimation problem for large-scale networks of pipelines is presented in [15], and the gradient descent algorithm is applied to solve the optimization problem.In practical systems, the supervisory control and data administration system experience random errors and bad data inevitably due to the sensor error and electromagnetic interference. The existing state estimation methods based on Kalman filter can reduce the random errors to some extent but are vulnerable to bad data. To solve this problem, a variety of improved Kalman filter algorithms are proposed. Based on variational Bayesian technique, an adaptive Kalman filter for linear Gaussian state-space models is proposed in [23, 24], which has better robustness to resist the uncertainties of process and measurement noise covariance matrices, as well as the colored measurement noise. In [25], a robust filter in a batch-mode regression form is developed to process the observations and predictions together, making it very effective in suppressing multiple outliers. To ensure the stability of the unscented Kalman filter, the unknown time-varying matrix is introduced to describe the prediction error of the unscented transformation in [26]. It can be seen that the robustness of the Kalman filter is a big challenge for its practical applications, attracting more and more researchers’ attentions.The existing gas pipeline network state estimation methods based on Kalman filter solve the transient flow equations by using a numerical approximation technique called finite element methods. These estimation methods have the following drawbacks:(1)
The gas pipe lines are divided into several linear elements, and the number of variables increases accordingly. For large-scale networks with long pipes, to ensure the estimating accuracy, many linear elements have to be created, as well as the states, causing low measurement redundancy and heavy computation load.(2)
The algorithm robustness against bad data is not considered. Once the measurements in the practical system experience bad data or non-zero mean noises, the estimating performances cannot be guaranteed.To cope with the above problems, and obtain accurate transient states of large-scale gas pipeline networks under practical operating conditions, this paper focuses on the robust dynamic state estimation method along with the following contributions:(1)
To deal with the problem that the transient equations are less than the states, the boundary conditions are used as supplementary constraints, and the linear process functions are established. The zero mass flow rates at the sink nodes are taken as virtual measurements, improving the measurement redundancy.(2)
A time-varying scalar matrix is proposed to regulate the measurement variance matrix according to the innovations, where the measurements can correct the predicted states accurately, and then increasing the robustness of the Kalman filter against the bad data and non-zero mean noises.(3)
The proposed dynamic state estimation method is tested on a 30-node natural gas pipeline network under the normal measurement condition, bad data condition, and non-zero mean noises condition. The simulation results show that the proposed approach manages to perform robust DSE of natural gas pipeline networks and the performance of the proposed method is better than the traditional Kalman filter.The rest of this paper is organized as follows. Section2 builds the mathematical model of the DSEs for gas pipelines. Section 3 proposes the robust Kalman filter algorithm based on the time-varying scalar matrix. Section 4 shows and analyzes the simulation results under various conditions. Finally, Section 5 concludes this paper.


## 2. DSE Modeling of Gas Pipeline Network

As the basic components of natural gas systems, the pipelines can store a certain amount of natural gases due to the compressibility, which is called the linepack storage. For the natural gas systems with a large number of pipelines, the linepack storage drives the dynamic process of the nodal pressure and mass flow of the natural gas in the pipelines. To describe the above physical dynamic mathematically, the state space representation of transient gas flows is introduced.


### 2.1. State Space Representation of Transient Gas Flow

Under stable operating conditions, the gas states such as the nodal pressure and mass flow rate in pipelines are constant. However, in practical natural gas systems, the stable operating condition can be disrupted easily due to the continuous changes of gas loads, supplies of gas companies, and other operational activities, causing the fluid dynamic process along pipelines, which can be represented by a set of nonlinear partial differential equations (PDEs). These PDEs are derived from the momentum conservation principles and material balances, which is given as the following state space form [4, 11, 27]:(1a)∂ρ∂t+∂ρu∂x=0,(1b)∂ρu∂t+∂ρu2∂x+∂p∂x=−fρuu2d−ρgsinθ,where, ρ, p, and u are the density, pressure, and gas velocity, respectively; t and x are the time and spatial coordinate, respectively; f is the friction factor; g is the gravitational constant; θ is the pipe inclination angle; and d is the pipe diameter. The second term on the left-hand side of (1b) describes the convective effect of the natural gas, which can be omitted when the gas velocity is significantly smaller than the sound speed in practical operation [4, 11]. Further, we assume that the pipelines are horizontal, hence the inclination angle is zero, and the second term on the right-hand side of (1b) can also be omitted. Under the condition that the gas pressure is less than 10 bar and the temperature is lower than 20°C, the functions between the pressure, mass flow rate, and the gas density are as follows [4, 11]:(2)p=c2ρ,m˙=ρua,where c2 = ZRT; m˙ and a are the mass flow rate and cross section area of pipelines, respectively, and Z, R, and T are the gas compressibility factor, the specific gas constant, and the gas temperature, respectively. For the natural gas, c2 is constant if the temperature is not changing. Here, the absolute value of average gas velocity u¯ is used instead of u. Under the aforementioned assumptions, equations (1a) and (1b) can be rewritten as the following equation [4]:(3a)∂ρ∂t+∂m˙a∂x=0,(3b)∂m˙a∂t+c2∂ρ∂x+fm˙u¯2da=0.The simplified partial differential equations (3a) and (3b) are able to describe the practical natural gas dynamic process in continuous form, which should be transformed to difference equation before numerical solving.


### 2.2. Discretization of PDE

A generalized model of natural gas pipelines is shown in Figure1. Assume that the node numbers of the three nodes of pipelines are i, j, and k, respectively, and k < i < j. The gas flow direction is defined as from the smaller node to the bigger node. Thus, the gas flow directions of pipeline k-i and i-j are k ⟶ i and i ⟶ j, respectively. The gas densities are ρk, ρi, and ρj, while the mass flow rates on the two sides of the pipeline i-j are m˙ij and m˙ji, respectively. The mass flow rate of gas load at node i is m˙i. For simplification, the compressor at k end of pipeline ki is modeled as a constant ratio of densities crki. The positive direction is defined as flowing out; thus the mass flow rate values are negative for the source nodes. Every natural gas pipeline can be modeled according to Figure 1 and satisfy the PDEs (3a) and (3b). One of the common methods for solving PDEs is the Euler finite difference technique, which is applied to solve (3a) and (3b) in our works. The differencing scheme is [28](4a)∂X∂t≈12Xs+1,t+1−Xs+1,tΔt+Xs,t+1−Xs,tΔt,(4b)∂X∂x≈12Xs+1,t+1−Xs,t+1Δx+Xs+1,t−Xs,tΔx,(4c)X≈14Xs+1,t+1+Xs,t+1+Xs+1,t+Xs,t,where X represents the generalized stats; Δt and Δx are the time and spatial step widths respectively; and the subscripts s and t represent the sth node and time step t, respectively. Equations (3a) and (3b) can be transformed to (5a) and (5b) according to the differencing scheme (4a)–(4c) for pipeline i-j.(5a)crjiρj,t+1−crjiρj,t+crijρi,t+1−crijρi,t+ΔtLPijaijm˙ji,t+1−m˙ij,t+1+m˙ji,t−m˙ij,t=0,(5b)m˙ji,t+1−m˙ij,t+1+m˙ji,t−m˙ij,t+aijΔtc2LPijcrjiρj,t+1−crijρi,t+1+crjiρj,t−crijρi,t+γv¯GΔt4dijaijm˙ij,t+1+m˙ji,t+1+m˙ij,t+m˙ji,t=0,where the subscript ij represents the pipeline i-j and Lij is the length of pipeline i-j. In addition to the difference equations, some boundary conditions should be satisfied.Figure 1
The general model of natural gas pipelines.


### 2.3. Boundary Conditions

The gas flow dynamic processes and gas density distribution are influenced by the changes of operating conditions at both the source and sink nodes. Under the practical operating condition, the mass flow rates at sink nodes are changing, while the gas density is constant at the source nodes. These constraints are defined as boundary conditions.The boundary condition of the mass flow rate balance at sink nodes is represented as(6)∑k∈i,k<im˙ik,t−∑j∈i,j>im˙ij,t=m˙i,t,i∈NSink,where k∈i means the node k connected to node i with pipeline i-k; i∈NSink represents i is the source node. The first and second terms of (6) represent the sum of mass flow rates from node i and to node i, respectively.For the source nodes, the assumption is made that the gas capacities are big enough so that the gas density is able to maintain constant in the dynamic process. This forms the following boundary condition:(7)ρi,t=ρSi,0,i∈NSource,where ρSi,0 is the gas density at source node i at initial time; i∈NSource represents i is the source node.The boundary conditions supply more constraints, making it possible to build the mathematical model of DSE for the natural gas pipeline networks.


### 2.4. Mathematical Model of DSE

The DSE based on Kalman filter includes two basic steps, the prediction and filtering. In the prediction step, the difference equations and boundary conditions are used to predict the states. In this work, the gas densities at nodes and mass flow rates at the two ends of pipelines are taken as states. Hence, for a natural gas pipeline network withnN nodes and nL pipelines, the number of states is nN + 2nL. The state vector at time step t is xt = [xr,t, xm,t]T, xr,t=ρ1,t,ρ2,t,…,ρnN,tT, xm,t=…,m˙ij,t,m˙ji,t,…T, i < j. xr,t is the state vector of gas densities; xm,t is the state vector of mass flow rates, and the dimension of xm,t is 2nL. Based on the above state definition, (5a) and (5b) can be rewritten as(8)A11A12A21A22xr,t+1xm,t+1=A11−A12−A21−A22xr,txm,t,(9)A11l,i=crij,ifi∝l,j∝l,0,else,A11∈ℝnP×nN,(10)A12=−α1α10−α2α2⋱0−αnLαnL,A12∈ℝnL×2nL,(11)αl=ΔtLijaij,i∝l,j∝l,(12)A21l,i=−crijβl,ifi∝l,j∝l,i<j,crjiβl,ifi∝l,j∝l,i>j,0,else,A21∈ℝnP×nN,(13)βl=aijΔtc2Lij,i∝l,j∝l,(14)A22l,i=γl−1,i∝l,j∝l,i<j,γl+1,i∝l,j∝l,i>j,0,else,A22∈ℝnL×2nL,(15)γl=fu¯Δt4dijaij,i∝l,j∝l,where (l, i) represents the row l and column i; i∝l means node i is one end of pipeline l. For the natural gas pipeline systems, the total number of state variables is nN + 2nL, but the number of equations in (8) is 2nL. Equation (8) cannot be solved because the states are more than the equations. The boundary conditions can be applied as the supplementary constraints. Similarly, (6) and (7) can be rewritten in the matrix form.(16)B110nin×2nL0nout×nNB22xr,t+1xm,t+1=ur,t+1um,t+1,(17)B11s,i=1,i∈NSource,0,else,B11∈ℝnin×nout,(18a)B22j,2l=1,j∈NSink,j∝l,0,else,B22∈ℝnout×2nL,(18b)B22i,2l−1=−1,ifiis a outlet node,i∝l,0,else,where ur,t+1 = […, ρSi,0, …]T, i∈NSource, ur,t+1∈ℝnSource; um,t+1 = […, m˙j,t+1, …]T, i∈NSink, um,t+1∈ℝnSink; nSource and nSink are the numbers of source and sink nodes, respectively; s is the source nodes number, s = 1, 2, …, nSource; 0 is the zero matrix and the subscript of 0 is the dimension.Equations (8) and (16) can be written as the following concentrated form:(19)Axt+1=Bxt+Ut+1,where A = A11A12A21A22B110nin×2nL0nout×nNB22, B = A11−A12−A21−A220nin×nN0nin×2nL0nout×nN0nout×2nL, Ut+1 = 0nL×10nL×1ur,t+1um,t+1.Premultiplied by the inverse matrix of the left matrix, (19) becomes(20)xt+1=Fxt+ut+1,where F = A−1B, ut+1 = A−1Ut+1.Equation (20) is the system equation of the gas pipeline networks, which is used to represent the dynamic processes. Additionally, the measured pressure and mass flow rate information from supervisory control and data administration systems provide redundant constraints for the DSE. In this work, assume that all of the nodes are equipped with flow meters and barometers; hence such information is taken as the measurement vectors. The measurement function is(21)zt+1=Ht+1xt+1,Ht+1=c2InN×nN0nN×2nL0nN×nNH′,where H′j,2l=1,j∝l,0,else,,H′i,2l−1=−1,i∝l,0,else,,H′∈ℝnN×2nL, zt+1 = [zr1,t+1, zr2,t+1, …, zrnN,t+1, zm1,t+1, zm2,t+1, …, zmnN,t+1]T; I is an identity matrix; and zri,t+1 and zmi,t+1 are the note pressure and mass flow rate measurements at time instant t + 1, respectively.Up to now, the DSE model of gas pipeline networks is formed:(22)xt+1=Fxt+ut+1+vt,zt+1=Hxt+1+wt+1,where vt and wk+1 are the system and measurement error vectors, respectively, vt∈ℝnN+2nL, wt+1∈ℝ2nN. The variance matrixes of vt+1 and wk+1 are Qt+1 and Rt+1, respectively. In this paper, the errors are assumed to satisfy the normal distribution. However, the practical measurement experiences bad data inevitably due to the electromagnetic interference or transmission errors. Additionally, in some cases, the mean values of the measurement errors are not zero, so the errors are non-zero mean noises.The states in (22) should be estimated by the DSE algorithm based on Kalman filter.


## 2.1. State Space Representation of Transient Gas Flow

Under stable operating conditions, the gas states such as the nodal pressure and mass flow rate in pipelines are constant. However, in practical natural gas systems, the stable operating condition can be disrupted easily due to the continuous changes of gas loads, supplies of gas companies, and other operational activities, causing the fluid dynamic process along pipelines, which can be represented by a set of nonlinear partial differential equations (PDEs). These PDEs are derived from the momentum conservation principles and material balances, which is given as the following state space form [4, 11, 27]:(1a)∂ρ∂t+∂ρu∂x=0,(1b)∂ρu∂t+∂ρu2∂x+∂p∂x=−fρuu2d−ρgsinθ,where, ρ, p, and u are the density, pressure, and gas velocity, respectively; t and x are the time and spatial coordinate, respectively; f is the friction factor; g is the gravitational constant; θ is the pipe inclination angle; and d is the pipe diameter. The second term on the left-hand side of (1b) describes the convective effect of the natural gas, which can be omitted when the gas velocity is significantly smaller than the sound speed in practical operation [4, 11]. Further, we assume that the pipelines are horizontal, hence the inclination angle is zero, and the second term on the right-hand side of (1b) can also be omitted. Under the condition that the gas pressure is less than 10 bar and the temperature is lower than 20°C, the functions between the pressure, mass flow rate, and the gas density are as follows [4, 11]:(2)p=c2ρ,m˙=ρua,where c2 = ZRT; m˙ and a are the mass flow rate and cross section area of pipelines, respectively, and Z, R, and T are the gas compressibility factor, the specific gas constant, and the gas temperature, respectively. For the natural gas, c2 is constant if the temperature is not changing. Here, the absolute value of average gas velocity u¯ is used instead of u. Under the aforementioned assumptions, equations (1a) and (1b) can be rewritten as the following equation [4]:(3a)∂ρ∂t+∂m˙a∂x=0,(3b)∂m˙a∂t+c2∂ρ∂x+fm˙u¯2da=0.The simplified partial differential equations (3a) and (3b) are able to describe the practical natural gas dynamic process in continuous form, which should be transformed to difference equation before numerical solving.


## 2.2. Discretization of PDE

A generalized model of natural gas pipelines is shown in Figure1. Assume that the node numbers of the three nodes of pipelines are i, j, and k, respectively, and k < i < j. The gas flow direction is defined as from the smaller node to the bigger node. Thus, the gas flow directions of pipeline k-i and i-j are k ⟶ i and i ⟶ j, respectively. The gas densities are ρk, ρi, and ρj, while the mass flow rates on the two sides of the pipeline i-j are m˙ij and m˙ji, respectively. The mass flow rate of gas load at node i is m˙i. For simplification, the compressor at k end of pipeline ki is modeled as a constant ratio of densities crki. The positive direction is defined as flowing out; thus the mass flow rate values are negative for the source nodes. Every natural gas pipeline can be modeled according to Figure 1 and satisfy the PDEs (3a) and (3b). One of the common methods for solving PDEs is the Euler finite difference technique, which is applied to solve (3a) and (3b) in our works. The differencing scheme is [28](4a)∂X∂t≈12Xs+1,t+1−Xs+1,tΔt+Xs,t+1−Xs,tΔt,(4b)∂X∂x≈12Xs+1,t+1−Xs,t+1Δx+Xs+1,t−Xs,tΔx,(4c)X≈14Xs+1,t+1+Xs,t+1+Xs+1,t+Xs,t,where X represents the generalized stats; Δt and Δx are the time and spatial step widths respectively; and the subscripts s and t represent the sth node and time step t, respectively. Equations (3a) and (3b) can be transformed to (5a) and (5b) according to the differencing scheme (4a)–(4c) for pipeline i-j.(5a)crjiρj,t+1−crjiρj,t+crijρi,t+1−crijρi,t+ΔtLPijaijm˙ji,t+1−m˙ij,t+1+m˙ji,t−m˙ij,t=0,(5b)m˙ji,t+1−m˙ij,t+1+m˙ji,t−m˙ij,t+aijΔtc2LPijcrjiρj,t+1−crijρi,t+1+crjiρj,t−crijρi,t+γv¯GΔt4dijaijm˙ij,t+1+m˙ji,t+1+m˙ij,t+m˙ji,t=0,where the subscript ij represents the pipeline i-j and Lij is the length of pipeline i-j. In addition to the difference equations, some boundary conditions should be satisfied.Figure 1
The general model of natural gas pipelines.


## 2.3. Boundary Conditions

The gas flow dynamic processes and gas density distribution are influenced by the changes of operating conditions at both the source and sink nodes. Under the practical operating condition, the mass flow rates at sink nodes are changing, while the gas density is constant at the source nodes. These constraints are defined as boundary conditions.The boundary condition of the mass flow rate balance at sink nodes is represented as(6)∑k∈i,k<im˙ik,t−∑j∈i,j>im˙ij,t=m˙i,t,i∈NSink,where k∈i means the node k connected to node i with pipeline i-k; i∈NSink represents i is the source node. The first and second terms of (6) represent the sum of mass flow rates from node i and to node i, respectively.For the source nodes, the assumption is made that the gas capacities are big enough so that the gas density is able to maintain constant in the dynamic process. This forms the following boundary condition:(7)ρi,t=ρSi,0,i∈NSource,where ρSi,0 is the gas density at source node i at initial time; i∈NSource represents i is the source node.The boundary conditions supply more constraints, making it possible to build the mathematical model of DSE for the natural gas pipeline networks.


## 2.4. Mathematical Model of DSE

The DSE based on Kalman filter includes two basic steps, the prediction and filtering. In the prediction step, the difference equations and boundary conditions are used to predict the states. In this work, the gas densities at nodes and mass flow rates at the two ends of pipelines are taken as states. Hence, for a natural gas pipeline network withnN nodes and nL pipelines, the number of states is nN + 2nL. The state vector at time step t is xt = [xr,t, xm,t]T, xr,t=ρ1,t,ρ2,t,…,ρnN,tT, xm,t=…,m˙ij,t,m˙ji,t,…T, i < j. xr,t is the state vector of gas densities; xm,t is the state vector of mass flow rates, and the dimension of xm,t is 2nL. Based on the above state definition, (5a) and (5b) can be rewritten as(8)A11A12A21A22xr,t+1xm,t+1=A11−A12−A21−A22xr,txm,t,(9)A11l,i=crij,ifi∝l,j∝l,0,else,A11∈ℝnP×nN,(10)A12=−α1α10−α2α2⋱0−αnLαnL,A12∈ℝnL×2nL,(11)αl=ΔtLijaij,i∝l,j∝l,(12)A21l,i=−crijβl,ifi∝l,j∝l,i<j,crjiβl,ifi∝l,j∝l,i>j,0,else,A21∈ℝnP×nN,(13)βl=aijΔtc2Lij,i∝l,j∝l,(14)A22l,i=γl−1,i∝l,j∝l,i<j,γl+1,i∝l,j∝l,i>j,0,else,A22∈ℝnL×2nL,(15)γl=fu¯Δt4dijaij,i∝l,j∝l,where (l, i) represents the row l and column i; i∝l means node i is one end of pipeline l. For the natural gas pipeline systems, the total number of state variables is nN + 2nL, but the number of equations in (8) is 2nL. Equation (8) cannot be solved because the states are more than the equations. The boundary conditions can be applied as the supplementary constraints. Similarly, (6) and (7) can be rewritten in the matrix form.(16)B110nin×2nL0nout×nNB22xr,t+1xm,t+1=ur,t+1um,t+1,(17)B11s,i=1,i∈NSource,0,else,B11∈ℝnin×nout,(18a)B22j,2l=1,j∈NSink,j∝l,0,else,B22∈ℝnout×2nL,(18b)B22i,2l−1=−1,ifiis a outlet node,i∝l,0,else,where ur,t+1 = […, ρSi,0, …]T, i∈NSource, ur,t+1∈ℝnSource; um,t+1 = […, m˙j,t+1, …]T, i∈NSink, um,t+1∈ℝnSink; nSource and nSink are the numbers of source and sink nodes, respectively; s is the source nodes number, s = 1, 2, …, nSource; 0 is the zero matrix and the subscript of 0 is the dimension.Equations (8) and (16) can be written as the following concentrated form:(19)Axt+1=Bxt+Ut+1,where A = A11A12A21A22B110nin×2nL0nout×nNB22, B = A11−A12−A21−A220nin×nN0nin×2nL0nout×nN0nout×2nL, Ut+1 = 0nL×10nL×1ur,t+1um,t+1.Premultiplied by the inverse matrix of the left matrix, (19) becomes(20)xt+1=Fxt+ut+1,where F = A−1B, ut+1 = A−1Ut+1.Equation (20) is the system equation of the gas pipeline networks, which is used to represent the dynamic processes. Additionally, the measured pressure and mass flow rate information from supervisory control and data administration systems provide redundant constraints for the DSE. In this work, assume that all of the nodes are equipped with flow meters and barometers; hence such information is taken as the measurement vectors. The measurement function is(21)zt+1=Ht+1xt+1,Ht+1=c2InN×nN0nN×2nL0nN×nNH′,where H′j,2l=1,j∝l,0,else,,H′i,2l−1=−1,i∝l,0,else,,H′∈ℝnN×2nL, zt+1 = [zr1,t+1, zr2,t+1, …, zrnN,t+1, zm1,t+1, zm2,t+1, …, zmnN,t+1]T; I is an identity matrix; and zri,t+1 and zmi,t+1 are the note pressure and mass flow rate measurements at time instant t + 1, respectively.Up to now, the DSE model of gas pipeline networks is formed:(22)xt+1=Fxt+ut+1+vt,zt+1=Hxt+1+wt+1,where vt and wk+1 are the system and measurement error vectors, respectively, vt∈ℝnN+2nL, wt+1∈ℝ2nN. The variance matrixes of vt+1 and wk+1 are Qt+1 and Rt+1, respectively. In this paper, the errors are assumed to satisfy the normal distribution. However, the practical measurement experiences bad data inevitably due to the electromagnetic interference or transmission errors. Additionally, in some cases, the mean values of the measurement errors are not zero, so the errors are non-zero mean noises.The states in (22) should be estimated by the DSE algorithm based on Kalman filter.


## 3. Robust Dynamic State Estimation Algorithm

One of the commonly used dynamic state estimation algorithms is Kalman filter, which is used in this work.


### 3.1. Kalman Filter

Kalman filter algorithm [29] includes two basic steps: prediction and filtering.


#### 3.1.1. Prediction Step

Given the initial estimated statex⌢0 and its covariance matrix Pt|t, the state prediction x˜t+1 and predicting covariance matrix Pt+1|t can be calculated by(23)x˜t+1=Fx⌢t+ut+1,(24)Pt+1|t=FPt|tFT+Qt.


#### 3.1.2. Filtering Step

The predicted statex˜t+1 should be corrected, and the estimated state x⌢t+1 and covariance matrix Pt+1|t+1 are obtained.(25)Kk+1=Pt+1|tHTHPt+1|tHT+Rt+1−1,(26)x⌢t+1=x˜t+1+Kt+1zt+1−Hx˜t+1,(27)Pt+1|t+1=Pt+1|t−Kk+1HPt+1|t.The right part of (25), Pe,t+1 = HPt+1|tHT + Rt+1, is the innovation covariance matrix, which represents the errors between measurements and their predicting values. When the measurements experience bad data, Pe,t+1 cannot represent the real errors correctly. Thus, the estimating performance of Kalman filter decreases. In this work, we introduce a time-varying scalar to regulate the measurement variance matrix, making Kalman filter robust.


### 3.2. Time-Varying Scalar

Kalman filter estimates the states by trading off the predicted states against the measurements according to the measurement and prediction variances. If the measurements experience bad data, Kalman filter would not correct the predicted states accurately and result in a poor estimation results. Aiming at this problem, to make the algorithm robust to the bad data, a time-varying scalar matrixμt is proposed to regulate the measurement variance matrix. The objective of μt is to fulfill(28)Pe,t+1=HPt+1|tH+μt+1Rt+1.To obtain the value ofμt+1, the sliding window method is used to estimate the innovation covariance matrix:(29)Pe,t+1=1mW∑i=0mW−1et+1−iet+1−iT=HPt+1|tH+μt+1Rt+1,where et+1=zt+1−Hx˜t+1 is the innovation vector and mW is the window length. By solving (29), μt+1 can be obtained as(30)μt+1=1mw∑i=0mw−1et+1−iet+1−iT−HPt+1|tHRt+1−1.The obtainedμt+1 by (30) may be non-diagonal; as a result, the inversion of the innovation covariance matrix is singular probably. Aiming at this problem, a diagonal matrix μt+1′ is defined(31)μt+1′=diagμ1′,μ2′,…,μ2nN′,where μi′=max1,μt+1,ii, i = 1, 2, …, 2nN, μt+1,ii is the ith diagonal element. Equation (25) becomes(32)Kk+1=Pt+1|tHTHPt+1|tHT+μt+1′Rt+1−1.With the help of the scalar, the modified robust Kalman filter can decrease the influence of bad data, and a better performance can be obtained than the traditional Kalman filter.


## 3.1. Kalman Filter

Kalman filter algorithm [29] includes two basic steps: prediction and filtering.


### 3.1.1. Prediction Step

Given the initial estimated statex⌢0 and its covariance matrix Pt|t, the state prediction x˜t+1 and predicting covariance matrix Pt+1|t can be calculated by(23)x˜t+1=Fx⌢t+ut+1,(24)Pt+1|t=FPt|tFT+Qt.


### 3.1.2. Filtering Step

The predicted statex˜t+1 should be corrected, and the estimated state x⌢t+1 and covariance matrix Pt+1|t+1 are obtained.(25)Kk+1=Pt+1|tHTHPt+1|tHT+Rt+1−1,(26)x⌢t+1=x˜t+1+Kt+1zt+1−Hx˜t+1,(27)Pt+1|t+1=Pt+1|t−Kk+1HPt+1|t.The right part of (25), Pe,t+1 = HPt+1|tHT + Rt+1, is the innovation covariance matrix, which represents the errors between measurements and their predicting values. When the measurements experience bad data, Pe,t+1 cannot represent the real errors correctly. Thus, the estimating performance of Kalman filter decreases. In this work, we introduce a time-varying scalar to regulate the measurement variance matrix, making Kalman filter robust.


## 3.1.1. Prediction Step

Given the initial estimated statex⌢0 and its covariance matrix Pt|t, the state prediction x˜t+1 and predicting covariance matrix Pt+1|t can be calculated by(23)x˜t+1=Fx⌢t+ut+1,(24)Pt+1|t=FPt|tFT+Qt.


## 3.1.2. Filtering Step

The predicted statex˜t+1 should be corrected, and the estimated state x⌢t+1 and covariance matrix Pt+1|t+1 are obtained.(25)Kk+1=Pt+1|tHTHPt+1|tHT+Rt+1−1,(26)x⌢t+1=x˜t+1+Kt+1zt+1−Hx˜t+1,(27)Pt+1|t+1=Pt+1|t−Kk+1HPt+1|t.The right part of (25), Pe,t+1 = HPt+1|tHT + Rt+1, is the innovation covariance matrix, which represents the errors between measurements and their predicting values. When the measurements experience bad data, Pe,t+1 cannot represent the real errors correctly. Thus, the estimating performance of Kalman filter decreases. In this work, we introduce a time-varying scalar to regulate the measurement variance matrix, making Kalman filter robust.


## 3.2. Time-Varying Scalar

Kalman filter estimates the states by trading off the predicted states against the measurements according to the measurement and prediction variances. If the measurements experience bad data, Kalman filter would not correct the predicted states accurately and result in a poor estimation results. Aiming at this problem, to make the algorithm robust to the bad data, a time-varying scalar matrixμt is proposed to regulate the measurement variance matrix. The objective of μt is to fulfill(28)Pe,t+1=HPt+1|tH+μt+1Rt+1.To obtain the value ofμt+1, the sliding window method is used to estimate the innovation covariance matrix:(29)Pe,t+1=1mW∑i=0mW−1et+1−iet+1−iT=HPt+1|tH+μt+1Rt+1,where et+1=zt+1−Hx˜t+1 is the innovation vector and mW is the window length. By solving (29), μt+1 can be obtained as(30)μt+1=1mw∑i=0mw−1et+1−iet+1−iT−HPt+1|tHRt+1−1.The obtainedμt+1 by (30) may be non-diagonal; as a result, the inversion of the innovation covariance matrix is singular probably. Aiming at this problem, a diagonal matrix μt+1′ is defined(31)μt+1′=diagμ1′,μ2′,…,μ2nN′,where μi′=max1,μt+1,ii, i = 1, 2, …, 2nN, μt+1,ii is the ith diagonal element. Equation (25) becomes(32)Kk+1=Pt+1|tHTHPt+1|tHT+μt+1′Rt+1−1.With the help of the scalar, the modified robust Kalman filter can decrease the influence of bad data, and a better performance can be obtained than the traditional Kalman filter.


## 4. Case Study

In this section, the proposed DSE method is applied to a 30-node natural gas pipeline network shown in Figure2. The time-varying mass flow rates of gas loads are simulated artificially and the dynamic processes of the node pressures and mass flow rates are calculated by (5a)∼(7). The simulation results are taken as true values, and the measurements are derived by adding random numbers to the real values. The traditional Kalman filter and the proposed algorithm are used to estimate the states of the gas pipeline system. The DSE is carried out in the following three measurement conditions: normal condition, bad data condition, and non-zero mean noises condition.Figure 2
The 30-node test system. The tiny black points are pipeline conjunction nodes. Nodes 1 and 2 are the inlet nodes, and the others are the outlet nodes. The outlet nodes without arrows mean that the gas loads are zero.


### 4.1. Description of the Test System

The test system includes 30 nodes and 29 pipelines, which is used to evaluate the performances of DSE methods. The lengths and cross-sectional diameters are given in Table1. The friction factor f and the gas speed c2 are 0.015 and 340 m/s, respectively. Gas is supplied at the source node 1 and 2, and the pressures of these two nodes are 27.8 bar and 28.5 bar, respectively. We assume that the capacities of gas sources are infinite, which means the pressures are constant in the dynamic processes.Table 1
Parameters of gas pipelines.
PipelineLength (km)Diameter (m)1, 350.63, 430.64, 540.55, 660.56, 770.52, 720.53, 830.48, 950.28, 1070.29, 1150.44, 1240.412, 1380.413, 14100.414, 1590.215, 16100.214, 1740.25, 18100.418, 2030.220, 2170.218, 1920.26, 22100.422, 2360.223, 2470.223, 2540.225, 2690.226, 2740.27, 2820.228, 2970.228, 3050.2The simulation interval is 15 minutes over the time horizon of 24 hours. All nodes are equipped with flow meters and barometers, so the pressures and mass flow rates can be measured. The standard deviations of pressure and mass flow rate measurements are 0.01 bar and 2%, respectively. Note that the mass flow rates at the conjunction nodes are zero, which can be taken as virtual measurements. The errors of virtual measurements are smaller than others, and the errors are set to 0.001 in our work.


### 4.2. Normal Measurement Condition

In this section, the DSE is carried out under the normal measurement condition that the measurement errors obey the Gauss distribution with zero mean. The traditional Kalman filter and robust Kalman filter are used to estimate the states, and the results are compared by the filter coefficientε [30, 31].(33)ε=∑tz⌢t−zt+2∑tzt−zt+2,where z⌢t is the estimated value of measurements, which can be calculated by the estimated states and zt+ is the true value of measurements without errors.The filter coefficients of the two kinds of DSE methods are shown in Table2. Nodes 1 and 2 are source nodes, the pressures of which are constant, so the filter coefficients are not computed. All results in Table 2 are smaller than 1, meaning that both the DSEs based on Kalman filter and robust Kalman filter are effective. It should be noted that the virtual measurement errors of mass flow rates at the sink nodes are very small and the denominator of (33) is approximately equal to the numerator; as a result, the filter coefficients are almost 1. Besides that, most of the coefficients of the robust Kalman filter are smaller than the traditional Kalman filter, meaning that the performance of the robust Kalman filter is better.Table 2
Filter coefficients of DSEs.
NodePressureMass flow rateKFRKFKFRKF1——0.11110.10842——0.48310.449230.00010.00010.98880.988840.00010.00010.98880.988850.00020.00010.98880.988860.00020.00010.98880.988870.00030.00020.98880.988880.00030.00020.98880.988890.00140.00110.56190.5161100.00610.00490.51010.3902110.00110.00080.13940.0706120.00110.00070.43010.3681130.00160.00110.29840.1957140.00270.00210.50320.4625150.01870.02450.64880.6339160.03550.03630.04050.0151170.03080.02060.00610.0021180.00060.00040.98880.9888190.03310.01440.16340.0956200.03880.02930.29030.2361210.06010.04730.20870.1523220.00040.00030.42320.3063230.00640.00420.98880.9888240.03720.02490.06610.0281250.03610.02980.43380.2521260.05010.03630.33760.2185270.06870.03650.13160.0748280.01930.01140.98880.9888290.05510.03180.09440.0606300.05430.03940.10210.0661


### 4.3. Bad Data Condition

In practical systems, the measuring equipments experience bad data inevitably [32], which should be considered in DSEs. In this section, based on the normal measurement condition, the bad data are added to the measurement vectors artificially to test the performances of the proposed DSE for gas pipeline networks. Here, the value of the pressure measurements at node 30 is set to 12, 10.7, and 13.8 bar at time 5, 5.25, and 5.5 hours, and 13, 15.5, and 23 bar at time 13.25, 13.5, and 13.75 hours, respectively. Meanwhile, the mass flow rate measurements at node 11 are reset to 3, 2.1, 3, and 2.2 kg/s at time 7.5, 7.75, 8, and 8.25 hours, and 3, 2.1, and 1.7 kg/s at time 15.75, 16, and 16.25 hours, respectively.The estimating results of Kalman filter and robust Kalman filter are shown in Figure3, and the time-varying scalars are shown in Figure 4. It can be seen that the curve of Kalman filter (the green dotted line) is deviated from the true value (the blue solid line) when the bad data appeared. This is because the measurement error variance matrix does not correspond with the truth in the bad data condition. However, with the help of the time-varying scalar, the curve of robust Kalman filter (the yellow dotted line) and the true value curve are well coincident all the time. The time-varying scalar can regulate the measurement error variance matrix according to the actual errors precisely and the accurate estimating results can be obtained. When the bad data appear, the values of the scalar increase dramatically, causing the reduction of the corresponding elements in Kalman gain matrix Kt+1. As a result, the correcting effects of the measurements decrease and the estimated values are little affected by the bad data.Figure 3
Estimating results of the pressure and mass flow rate under bad data condition.
(a)(b)Figure 4
Time-varying scalars of the pressure and mass flow rate under bad data condition.
(a)(b)


### 4.4. Non-Zero Mean Noises Condition

The errors of measurements in the above two sections are Gaussian white noises, but in practice the non-zero mean noises are more common, which are considered in this section. Based on the condition of section B, we added constant deviations 0.2 bar and 0.1 kg/s to the pressure measurements and mass flow rate measurements from 10∼19.75 hours and 5∼12.5 hours, respectively.The estimating results are shown in Figures5 and 6. It can be seen that the estimating curves of Kalman filter deviate from the true value curves obviously under the non-zero mean noises condition. At the same time, we can notice that the scalars increase obviously while the measurement mean values are deviating from zero, and the estimating curves based on robust Kalman filter fit with the true value all the time. Furthermore, the Kalman filter curves deviate from true value gradually at the beginning of the non-zero mean noises. This is because the negative effects of the non-zero mean noises are accumulating continuously throughout the whole process. If the non-zero mean noises exist all the time, eventually, the estimated value of the Kalman filter will coincide with the deviation.Figure 5
Estimating results of the pressure and mass flow rate under non-zero mean noises condition.
(a)(b)Figure 6
Time-varying scalars of the pressure and mass flow rate under non-zero mean noises condition.
(a)(b)


## 4.1. Description of the Test System

The test system includes 30 nodes and 29 pipelines, which is used to evaluate the performances of DSE methods. The lengths and cross-sectional diameters are given in Table1. The friction factor f and the gas speed c2 are 0.015 and 340 m/s, respectively. Gas is supplied at the source node 1 and 2, and the pressures of these two nodes are 27.8 bar and 28.5 bar, respectively. We assume that the capacities of gas sources are infinite, which means the pressures are constant in the dynamic processes.Table 1
Parameters of gas pipelines.
PipelineLength (km)Diameter (m)1, 350.63, 430.64, 540.55, 660.56, 770.52, 720.53, 830.48, 950.28, 1070.29, 1150.44, 1240.412, 1380.413, 14100.414, 1590.215, 16100.214, 1740.25, 18100.418, 2030.220, 2170.218, 1920.26, 22100.422, 2360.223, 2470.223, 2540.225, 2690.226, 2740.27, 2820.228, 2970.228, 3050.2The simulation interval is 15 minutes over the time horizon of 24 hours. All nodes are equipped with flow meters and barometers, so the pressures and mass flow rates can be measured. The standard deviations of pressure and mass flow rate measurements are 0.01 bar and 2%, respectively. Note that the mass flow rates at the conjunction nodes are zero, which can be taken as virtual measurements. The errors of virtual measurements are smaller than others, and the errors are set to 0.001 in our work.


## 4.2. Normal Measurement Condition

In this section, the DSE is carried out under the normal measurement condition that the measurement errors obey the Gauss distribution with zero mean. The traditional Kalman filter and robust Kalman filter are used to estimate the states, and the results are compared by the filter coefficientε [30, 31].(33)ε=∑tz⌢t−zt+2∑tzt−zt+2,where z⌢t is the estimated value of measurements, which can be calculated by the estimated states and zt+ is the true value of measurements without errors.The filter coefficients of the two kinds of DSE methods are shown in Table2. Nodes 1 and 2 are source nodes, the pressures of which are constant, so the filter coefficients are not computed. All results in Table 2 are smaller than 1, meaning that both the DSEs based on Kalman filter and robust Kalman filter are effective. It should be noted that the virtual measurement errors of mass flow rates at the sink nodes are very small and the denominator of (33) is approximately equal to the numerator; as a result, the filter coefficients are almost 1. Besides that, most of the coefficients of the robust Kalman filter are smaller than the traditional Kalman filter, meaning that the performance of the robust Kalman filter is better.Table 2
Filter coefficients of DSEs.
NodePressureMass flow rateKFRKFKFRKF1——0.11110.10842——0.48310.449230.00010.00010.98880.988840.00010.00010.98880.988850.00020.00010.98880.988860.00020.00010.98880.988870.00030.00020.98880.988880.00030.00020.98880.988890.00140.00110.56190.5161100.00610.00490.51010.3902110.00110.00080.13940.0706120.00110.00070.43010.3681130.00160.00110.29840.1957140.00270.00210.50320.4625150.01870.02450.64880.6339160.03550.03630.04050.0151170.03080.02060.00610.0021180.00060.00040.98880.9888190.03310.01440.16340.0956200.03880.02930.29030.2361210.06010.04730.20870.1523220.00040.00030.42320.3063230.00640.00420.98880.9888240.03720.02490.06610.0281250.03610.02980.43380.2521260.05010.03630.33760.2185270.06870.03650.13160.0748280.01930.01140.98880.9888290.05510.03180.09440.0606300.05430.03940.10210.0661


## 4.3. Bad Data Condition

In practical systems, the measuring equipments experience bad data inevitably [32], which should be considered in DSEs. In this section, based on the normal measurement condition, the bad data are added to the measurement vectors artificially to test the performances of the proposed DSE for gas pipeline networks. Here, the value of the pressure measurements at node 30 is set to 12, 10.7, and 13.8 bar at time 5, 5.25, and 5.5 hours, and 13, 15.5, and 23 bar at time 13.25, 13.5, and 13.75 hours, respectively. Meanwhile, the mass flow rate measurements at node 11 are reset to 3, 2.1, 3, and 2.2 kg/s at time 7.5, 7.75, 8, and 8.25 hours, and 3, 2.1, and 1.7 kg/s at time 15.75, 16, and 16.25 hours, respectively.The estimating results of Kalman filter and robust Kalman filter are shown in Figure3, and the time-varying scalars are shown in Figure 4. It can be seen that the curve of Kalman filter (the green dotted line) is deviated from the true value (the blue solid line) when the bad data appeared. This is because the measurement error variance matrix does not correspond with the truth in the bad data condition. However, with the help of the time-varying scalar, the curve of robust Kalman filter (the yellow dotted line) and the true value curve are well coincident all the time. The time-varying scalar can regulate the measurement error variance matrix according to the actual errors precisely and the accurate estimating results can be obtained. When the bad data appear, the values of the scalar increase dramatically, causing the reduction of the corresponding elements in Kalman gain matrix Kt+1. As a result, the correcting effects of the measurements decrease and the estimated values are little affected by the bad data.Figure 3
Estimating results of the pressure and mass flow rate under bad data condition.
(a)(b)Figure 4
Time-varying scalars of the pressure and mass flow rate under bad data condition.
(a)(b)


## 4.4. Non-Zero Mean Noises Condition

The errors of measurements in the above two sections are Gaussian white noises, but in practice the non-zero mean noises are more common, which are considered in this section. Based on the condition of section B, we added constant deviations 0.2 bar and 0.1 kg/s to the pressure measurements and mass flow rate measurements from 10∼19.75 hours and 5∼12.5 hours, respectively.The estimating results are shown in Figures5 and 6. It can be seen that the estimating curves of Kalman filter deviate from the true value curves obviously under the non-zero mean noises condition. At the same time, we can notice that the scalars increase obviously while the measurement mean values are deviating from zero, and the estimating curves based on robust Kalman filter fit with the true value all the time. Furthermore, the Kalman filter curves deviate from true value gradually at the beginning of the non-zero mean noises. This is because the negative effects of the non-zero mean noises are accumulating continuously throughout the whole process. If the non-zero mean noises exist all the time, eventually, the estimated value of the Kalman filter will coincide with the deviation.Figure 5
Estimating results of the pressure and mass flow rate under non-zero mean noises condition.
(a)(b)Figure 6
Time-varying scalars of the pressure and mass flow rate under non-zero mean noises condition.
(a)(b)


## 5. Conclusion

This paper proposes a robust dynamic state estimation method against bad data for natural gas pipeline networks. The method is applied on a 30-node pipeline network in several conditions, and the filter coefficient is used to evaluate the performances of the dynamic state estimation based on the traditional Kalman filter and the proposed robust method. The results show that most of the filter coefficient values of the proposed method are smaller than the traditional Kalman filter in the normal measurement condition. This means that the filtering efficiency of the robust method is better. Furthermore, the two methods are studied under the bad data and non-zero mean measurement noises conditions, and the results show that the proposed method can decrease the effects of bad data, obtaining accurate estimating states under these two conditions.Future work will focus on the accurate DSE modeling based on the nonlinear transient equations of the pipeline networks with compressors under the operating limitation constraints, and calculating method of predicting error covariance; another interesting topic is the adaptive DSE algorithm that can cope with the model uncertainty due to parameter errors and malicious cyber attacks [33, 34]. Besides, it is interesting to extend the presented approach to the DSE of integrated energy systems with integration of new elements such as cogeneration units [35], electrical vehicles [36], modular multilevel converter based high-voltage direct current systems [37], and uncertain renewable generations [38].



---

*Source: 5590572-2021-03-08.xml*