Robust Kalman filter-based dynamic state estimation of natural gas pipeline networks

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 networks 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.


INTRODUCTION
In comparison with the traditional coal-fired power units, gas-fired electric generators can response 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][21][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 establish the state estimation model applying the PDEs directly, and solve 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 is 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 is can be seen that the robustness of the Kalman filter is a big challenge for its practical applications, attacking 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 the 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 the 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, making that 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. Section 2 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, Section5 concludes this paper.

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.

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 (PDE). These PDEs are derived from the momentum conservation principles and materiel balances, which is given as the following state space form [11,4,27]: 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; 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 [11,4]: where, c 2 = ZRT; m and a are the mass flow rate and cross section area of pipelines respectively; Z, R and T are the gas compressibility factor, the specific gas constant and the gas temperature, respectively. For the natural gas, c 2 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, equation (1) can be rewritten as the following equation [4]: The simplified partial differential equation (4) is able to describe the practical natural gas dynamic process in continuous form, which should be transformed to difference equation before numerical solving.

Discretization of PDE
A generalized model of natural gas pipelines is shown in Fig. 1. Assuming 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 ij m and ji m , respectively. The mass flow rate of gas load at node iis i m .
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 Fig. 1, and satisfy the PDEs (4). One of the common methods for solving PDEs is the Euler finite difference technique, which is applied to solve (4) in our works. The differencing scheme is [28]: where, X represents the generalized stats; t and x are the time and spatial step widths respectively; the subscripts s and t represent the sth node and time step t, respectively. Equation (4) can be transformed to (6) according to the differencing scheme (5) for pipeline i-j.
where, the subscript ij represents the pipeline i-j;Lij is the length of pipeline i-j. In addition to the difference equations, some boundary conditions should be satisfied.

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: where, ki  means the node k connected to node i with pipeline i-k; Sink iN represents i is the source node. The first and second terms of (7) 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: where, Si,0 is the gas density at source node i at initial time; Source iN  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.

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 with nN nodes and nL pipelines, the number of states is nN + 2nL.
where, (l, i) represents the row l and column i; il 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 (9) is 2nL. Equation (9) cannot be solved because the states are more than the equations. The boundary conditions can be applied as the supplementary constraints. Similarly, (7) and (8) 22 1, if is a outlet node, ( Equation (9), (17) can be written as the following concentrated form: where, A=

Premultiplied by the inverse matrix of the left matrix, (20) becomes
Equation (21) 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 zt+1 = Ht+1xt+1 (22)  I is an identity matrix; 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 Hx w (24) where, vt and wk+1 are the system and measurement error vectors respectively, . 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 experience 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 (24) should be estimated by the DSE algorithm based on Kalman filter.

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

Kalman filter
Kalman filter algorithm [29] includes two basic steps: prediction and filtering. Prediction step. Given the initial estimated state 0 x and its covariance matrix Pt|t, the state prediction x Fx u (25) (28) Pt+1|t+1 = Pt+1|tKk+1HPt+1|t (29) The right part of (27), Pe,t+1 = HPt+1|tH T + 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.

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 (30).
To obtain the value of t+1, the sliding window method is used to estimate the innovation covariance matrix: where, Hx is the innovation vector; mW is the window length. By solving (31), t+1 can be obtained as The obtained t+1 by (32) may be non-diagonal, as a result, the inversion of the innovation covariance matrix is singular probably. Aiming at this problem, a diagonal matrix where, 1, max (1, ) i t ii , i = 1, 2, …, 2nN, t+1,ii is the ith diagonal element. Equation (27) becomes 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.

Case study
In this section, the proposed DSE method is applied to a 30-node natural gas pipeline network shown in Fig. 2. 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 (6)~(8). 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.

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 Tab.1. The friction factor f and the gas speed c 2 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.
The 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.

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].  (35) where, t z is the estimated value of measurements, which can be calculated by the estimated states; t z is the true value of measurements without errors.
The filter coefficients of the two kinds of DSE methods are shown in Tab.2. Nodes 1 and 2 are source nodes, the pressures of which are constant, so the filter coefficients are not computed. All results in Tab. 2 are smaller than1, meaning that both the DSE 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 (35) is approximately equal to the numerator, as a result, the filter coefficients is 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.

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  Fig. 3, and the time-varying scalars are shown in Fig. 4. It is 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 appeared, 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.

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 hour and 5~12.5 hour, respectively.
The estimating results are shown in Fig. 5 and Fig. 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.

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]. Estimating results of the pressure and mass flow rate undernon-zero mean noises condition.