CN112906317B - Robust dynamic state estimation method for natural gas pipe network - Google Patents
Robust dynamic state estimation method for natural gas pipe network Download PDFInfo
- Publication number
- CN112906317B CN112906317B CN202110177597.7A CN202110177597A CN112906317B CN 112906317 B CN112906317 B CN 112906317B CN 202110177597 A CN202110177597 A CN 202110177597A CN 112906317 B CN112906317 B CN 112906317B
- Authority
- CN
- China
- Prior art keywords
- node
- natural gas
- time
- obtaining
- matrix
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- VNWKTOKETHGBQD-UHFFFAOYSA-N methane Chemical compound C VNWKTOKETHGBQD-UHFFFAOYSA-N 0.000 title claims abstract description 128
- 239000003345 natural gas Substances 0.000 title claims abstract description 64
- 238000000034 method Methods 0.000 title claims abstract description 25
- 238000005259 measurement Methods 0.000 claims abstract description 30
- 238000013178 mathematical model Methods 0.000 claims abstract description 17
- 238000001914 filtration Methods 0.000 claims abstract description 12
- 239000011159 matrix material Substances 0.000 claims description 52
- 239000007789 gas Substances 0.000 claims description 48
- 239000013598 vector Substances 0.000 claims description 27
- 238000002347 injection Methods 0.000 claims description 4
- 239000007924 injection Substances 0.000 claims description 4
- 238000012937 correction Methods 0.000 claims description 3
- 230000005484 gravity Effects 0.000 claims description 2
- 230000009897 systematic effect Effects 0.000 claims description 2
- 230000005540 biological transmission Effects 0.000 description 4
- 238000009530 blood pressure measurement Methods 0.000 description 2
- 238000012544 monitoring process Methods 0.000 description 2
- 238000013341 scale-up Methods 0.000 description 2
- 238000004088 simulation Methods 0.000 description 2
- 238000012360 testing method Methods 0.000 description 2
- 230000001052 transient effect Effects 0.000 description 2
- 238000013459 approach Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000006835 compression Effects 0.000 description 1
- 238000007906 compression Methods 0.000 description 1
- 238000013523 data management Methods 0.000 description 1
- 230000007423 decrease Effects 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 230000007613 environmental effect Effects 0.000 description 1
- 238000012806 monitoring device Methods 0.000 description 1
- 238000010248 power generation Methods 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/28—Design optimisation, verification or simulation using fluid dynamics, e.g. using Navier-Stokes equations or computational fluid dynamics [CFD]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2113/00—Details relating to the application field
- G06F2113/08—Fluids
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
- G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
- Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation
Abstract
The invention discloses a robust dynamic state estimation method of a natural gas pipe network, which comprises the following steps: the method comprises the steps of (1) obtaining a dynamic mathematical model of a natural gas pipeline; (2) obtaining natural gas pipe network boundary conditions; (3) obtaining a model of natural gas pipe network dynamic state estimation; (4) And obtaining a state quantity estimated value of the natural gas pipe network by adopting a robust Kalman filtering algorithm. The invention can accurately obtain the accurate state of the large natural gas pipeline network under the condition of bad data and colored noise; measurement redundancy is improved; the state quantity does not increase sharply with the increase of the network scale, and the efficiency is improved.
Description
Technical Field
The present invention relates to a dynamic state estimation method, and in particular, to a robust dynamic state estimation method for a natural gas pipe network.
Background
With the increasing severity of environmental problems, renewable energy power generation systems continue to scale up. Due to uncertainty of renewable energy sources, the consumption of natural gas varies randomly to a certain extent, which affects the safe and stable operation of the system. Therefore, the accurate dynamic characteristics such as the pressure and the flow of the natural gas pipe network are obtained, and the method is a key for ensuring the safety and the optimized operation of the comprehensive energy system containing the electric power and the natural gas. To obtain these dynamic states, the network must be equipped with a large number of monitoring devices, which requires a large investment, and even so it is not possible to install sensors on each node of the network and obtain all the states. On the other hand, the measurement device inevitably encounters random errors and bad data, so that state estimation of the collected data is required to be able to apply to leak detection and control problems.
The existing natural gas pipe network dynamic state estimation is mainly based on a Kalman filtering gas transmission pipe network state estimation method, namely a transient flow equation set is solved through a numerical approximation technology called finite element. These approaches have serious limitations with scale up, such as scalability, unacceptable computational load, and lack of robustness in the event of sensor failure. In practical systems, random errors and bad data are inevitably present in the monitoring and data management systems due to sensor errors and electromagnetic interference. The existing state estimation method based on Kalman filtering reduces random errors to a certain extent, but is easily affected by bad data. Therefore, the enhancement of the robustness of the system in the process of accurately obtaining the transient state of the large natural gas pipeline network under the conditions of severe data and colored noise is more important for the safe and stable operation and the optimal control of the natural gas transmission pipeline network.
Disclosure of Invention
The invention aims to: the invention aims to solve the defects in the prior art, provides a robust dynamic state estimation method for a natural gas pipe network, and solves the problem of the robustness of the state estimation of the natural gas pipe network under the conditions of severe data and colored noise.
The technical scheme is as follows: the invention relates to a robust dynamic state estimation method of a natural gas pipe network, which comprises the following steps:
(1) Obtaining a natural gas pipeline dynamic mathematical model;
(2) Obtaining a natural gas pipe network boundary condition;
(3) Obtaining a model of natural gas pipe network dynamic state estimation;
(4) And obtaining a state quantity estimated value of the natural gas pipe network by adopting a robust Kalman filtering algorithm.
The step (1) comprises the following steps:
(11) Obtaining partial differential mathematical model of natural gas pipeline
Wherein ρ is the density of the gas in the pipeline, p is the pressure of the gas in the pipeline, v G The gas velocity in the pipeline is t, time, zeta is a space coordinate, f is a friction factor, G is a gravity constant, theta is a pipeline inclination angle, and d is a pipeline diameter;
(12) Obtaining a simplified partial differential mathematical model of a natural gas pipeline
wherein ,p＝c2 ρ，c 2 ＝ZRT，/>A is the flow rate of a pipeline, a is the cross-sectional area of the pipeline, Z is the gas compressibility factor, R is the gas constant, T is the gas temperature, c is the sound velocity, and if the temperature is unchanged for natural gas, c 2 Constant (I)>Is the average value of the gas flow rate;
(13) Obtaining a natural gas pipeline differential mathematical model
Wherein i, j, k are three nodes of the pipe network, and k<i<j, gas flows into the large node from the smaller node, and the subscript ij is the number of the head and the tail nodes of the pipeline, L ij For the length of pipe i-j ρ i,t and ρi,t+1 The densities of the gas at the node i at the time t and the time t+1 are respectively, ρ j,t and ρj,t+1 The density of the gas at the node j at the time t and t+1, a ij For the cross-sectional area of conduit i-j, d ij For the diameter of the pipe i-j, Δt is the time step, and />The flow from node i to node j at times t and t+1, respectively, +.> and />And respectively defining positive directions as outflow and negative flow values of source nodes for flows flowing from the node j to the node i at the times t and t+1.
The step (2) comprises the following steps:
(21) Boundary conditions for sink node flow balancing
wherein ,representing the flow from node i to node k at time t, k e i representing that node k is connected to node i through pipe i-k; i epsilon N Sink Indicating i is a sink node, N Sink Indicates the number of bus nodes, ">For the traffic of node i at time t, and />Representing the sum of traffic flowing from and into node i.
(22) Boundary condition of source node gas density
ρ g,t ＝ ρ Sg,0 , g∈N Source (8)
wherein ,ρSg,0 The gas density of the source node g at the initial moment; n (N) Source Represents the number of source nodes, g.epsilon.N Source G is denoted as source node.
The step (3) comprises the following steps:
(31) Obtaining system equation x t+1 ＝Fx t +u t+1 (21)
wherein ,xt and xt+1 State vectors, x, at times t and t+1, respectively t ＝[x r,t ,x m,t ] T ， i<j，x r,t Is the state vector of the gas density, x m,t Is the state vector of the flow rate, ρ 1,t ,ρ 2,t ,/>Nodes 1,2 and n at time t respectively N Density, x m,t Is 2n L ，n N and nL Representing the number of nodes and the number of pipes, a system equation matrix +.>Controlled variable n in and nout The number of injection nodes and the number of outflow nodes are respectively represented; a is that 11 ，A 12 ，A 21 ，A 22 ，B 11 ，B 22 A submatrix, a zero matrix and a matrix dimension as subscripts; u (u) r,t+1 ，u m,t+1 Density and flow control amount, respectively; />(l, i) represents the first row and the first column; i.alpha.l indicates that node i is one end of pipe l, +.>As a set of real numbers, α 1 ，α 2 ，α l ，/>as matrix A 12 Middle element, beta l As matrix A 21 Middle element, gamma l As matrix A 22 The element j-l represents that the node j is one end of the pipeline l, (s, g) represents the s-th row and the g-th column, (j, 2 l) represents the j-th row and the 2 l-th column, (i, 2 l-1) represents the i-th row and the 2 l-1-th column;
(32) Obtaining measurement equations
z t+1 ＝ H t+1 x t+1 (22)
wherein ,zt+1 Measuring vector, H for time t+1 t+1 Measuring a matrix for the time t+1; i is an identity matrix; z r1,t+1 ,z r2,t+1 ,/>Is a measurement of pressure at time t+1, z m1,t+1 ,z m2,t+1 ,…,/>Is a measurement of flow at time t+1;
(33) Model for obtaining natural gas pipe network dynamic state estimation
wherein ,vt Is the systematic error vector, w k+1 Is the vector of the measurement error and,u k+1 for the control vector, H is the measurement equation matrix, v t+1 Variance matrix is Q t+1 ，w t+1 Variance matrix R t+1 。
The step (4) comprises the following steps:
(41) According to the state quantity estimated value at the time tCovariance matrix P t|t Obtaining state quantity predictive value ∈ ->And a prediction covariance matrix P t+1|t
P t+1|t ＝ FP t|t F T + Q t (26)
wherein ,predicted value of state quantity at time t+1, < >>The state quantity estimated value at the moment t; p (P) t+1|t and Pt|t The prediction error variance and the estimation error variance, Q t Measuring an error equation array for the time t;
(42) Correction state quantity predicted valueUpdating covariance matrix P t+1|t+1
K k+1 ＝ P t+1|t H T (P e,t+1 ) -1 (27)
P t+1|t+1 ＝ P t+1|t -K k+1 HP t+1|t (29)
wherein ,Kk+1 Is a gain matrix; p (P) e,t+1 Is an innovation matrix;the state quantity estimated value at time t+1.
P in step (42) e,t+1 ＝HP t+1|t H+μ t+1 R t+1 Time-varying scale factor For the innovation vector, m W Is the window length.
P in step (42) e,t+1＝ HP t+1|t H T +μ t ′ +1 R t+1 ，μ i ′＝max(1,μ t+1,ii ),i＝1,2,…,2n N ,μ t+1,ii Is mu t+1 Is the i-th diagonal element of (c).
The beneficial effects are that: compared with the prior art, the invention has the remarkable advantages that: (1) And adjusting the measurement variance matrix according to the information, so that the state quantity predicted value can be accurately corrected by measurement, the robustness of the Kalman filter to bad data and colored noise is improved, and the accurate state of the large natural gas pipeline network under the condition of containing the bad data and the colored noise is obtained. (2) And zero injection flow at the confluence node is used as virtual measurement, so that the measurement redundancy is improved. (3) The average gas flow rate is used to approximate a second order term, and the pressure and the flow of nodes at two ends of the pipeline are used as state quantities, so that the state quantities cannot be increased sharply along with the increase of the network scale, and the efficiency is improved.
Drawings
FIG. 1 is a mathematical model of a natural gas pipeline of the present invention;
FIG. 2 is a schematic diagram of a test system according to the present invention;
FIGS. 3 (a) -3 (b) are pressure and flow estimates for bad data conditions of the present invention;
FIGS. 4 (a) -4 (b) are time-varying scale factors for pressure and flow under bad data conditions of the present invention;
FIGS. 5 (a) -5 (b) are the results of pressure and flow estimation under colored noise conditions of the present invention;
fig. 6 (a) -6 (b) are time-varying scale factors for pressure and flow under colored noise conditions of the present invention.
Detailed Description
The technical scheme of the invention is further described below with reference to the accompanying drawings.
The invention relates to a robust dynamic estimation method of a natural gas pipe network, which comprises the following steps:
(1) Establishing a natural gas pipeline dynamic mathematical model;
(2) Establishing a natural gas pipe network boundary condition;
(3) Establishing a model of natural gas pipe network dynamic state estimation;
(4) And obtaining a state quantity estimated value of the natural gas pipe network by adopting a robust Kalman filtering algorithm.
The step (1) comprises the following steps:
(11) Obtaining partial differential mathematical model of natural gas pipeline
Wherein: ρ is the gas density in the pipeline, p is the gas pressure in the pipeline, v G The gas velocity in the pipeline, t is time, ζ is space coordinate, and f is friction factor; g is a gravitational constant; θ is the inclination of the pipeline; d is the diameter of the pipe.
(12) The simplified partial differential mathematical model of the natural gas transmission pipeline is obtained by:
a function between pressure, flow and gas density at a gas pressure of less than 10bar and a temperature of less than 20 ℃):
p ＝ c 2 ρ (2)
in the formula ：c2 ＝ZRT，And a is the flow and cross-sectional area of the pipe, respectively; z, R and T are the gas compression coefficient, gas constant and gas temperature, respectively. c is the sound velocity, for natural gas, if the temperature is unchanged, c 2 Is constant. Here the absolute value of the mean value of the gas flow rate is used +.>Instead of |u|.
Since the natural gas flow rate is much less than the speed of sound and assuming that the pipeline is horizontal, the simplified partial differential mathematical model of the natural gas pipeline becomes:
wherein ,is the average value of the gas flow rate.
(13) The obtained natural gas pipeline differential mathematical model is as follows:
the finite differential process of partial differential equation discretization is:
wherein: x represents a generalized variable; Δt and Δx are the time and space steps, respectively; subscripts s and t are node number and time, respectively; x is X S+1,t+1 and XS+1,t Variable values at time t+1 and time t S+1, respectively; x is X S,t+1 and XS,t The variable values at time t+1 and time t S, respectively. Equation (4) can be converted into a natural gas transmission pipeline differential mathematical model:
wherein i, j, k are three nodes of the pipe network, and k<i<j, gas flows from the smaller node into the larger node. Subscript ij is the number of the head and the tail nodes of the pipeline; l (L) ij For the length of pipe i-j ρ i,t and ρi,t+1 The densities of the gas at the node i at the time t and the time t+1 are respectively, ρ j,t and ρj,t+1 The density of the gas at the node j at the time t and t+1, a ij For the cross-sectional area of conduit i-j, d ij Representing the diameter of the pipe i-j, deltat is the time step, and />Respectively representing the traffic flowing from node i to node j at times t and t +1, and />Respectively representing t and t+1 time slave nodesj flows to node i define the positive direction as the outflow, so the flow value of the source node is negative.
Step (2) comprises the following steps:
(21) Boundary condition representation of sink node flow balance
The boundary condition is that the pressure of the air source node is constant and the sum of the flow of the confluence nodes is 0. The boundary conditions for flow balancing at the sink node are expressed as:
in the formula ：representing the flow from node i to node k at time t, k e i representing that node k is connected to node i through pipe i-k; i epsilon N Sink Representing i as a source node, N Sink Indicates the number of bus nodes, ">Is the traffic of node i at time t. and />The sum of the traffic of the outgoing node i and the incoming node i is represented, respectively.
(22) Boundary condition of source node gas density
For the source node, the gas volume is assumed to be large enough so that the gas density can remain constant during the dynamic process, creating the following boundary conditions:
ρ g,t ＝ ρ Sg,0 , g∈N Source (8)
in the formula ：ρSg,0 The gas density of the source node g at the initial moment; n (N) Source Represents the number of source nodes, g.epsilon.N Source G is denoted as source node.
The step (3) comprises the following steps:
(31) Obtaining system equation x t+1 ＝Fx t +u t+1
The dynamic state estimation equation takes the gas density at the node and the flow at the two ends of the pipeline as state quantities. For having n N Node and n L Natural gas pipe network of pipeline, state quantity number is n N +2n L 。x t and xt+1 State vectors, x, at times t and t+1, respectively t ＝[x r,t ,x m,t ] T ，x r,t Is the state vector of the gas density at the time t; x is x m,t Is the state vector of the flow at the moment t, ρ 1,t ,ρ 2,t ,/>Respectively t time nodes 1,2 and n N Density, x m,t Is 2n L ，n N and nL Representing the number of nodes and the number of pipes. According to the above state quantity definition, the formula (6) can be rewritten as a matrix form as follows:
wherein ：is a real number set, x r,t+1 Is the state vector of the gas density at time t+1; x is x m,t+1 Is the state vector of the flow at the time t+1, A 11 ，A 12 ，A 21 ，A 22 Is a submatrix; alpha 1 ，α 2 ，α l ，/>As matrix A 12 Middle element, beta l As matrix A 21 Middle element, gamma l As matrix A 22 The element, (l, i) represents the first row and the i-th column; i.oc l, j.oc l indicates that node i and node j are one end of pipe l. For a natural gas pipeline network, the total number of state quantities is n N +2n L 。
The number of equations is smaller than the number of state quantities, and boundary conditions are needed as supplementary constraints. Formulas (7) and (8) can be written in the form of a matrix as follows:
wherein ,B11 ，B 22 A submatrix, a zero matrix and a matrix dimension as subscripts; u (u) r,t+1 ，u m,t+1 Density and flow control amount, respectively; n is n in and nout The number of injection nodes and the number of outflow nodes are respectively represented; (s, g) represents the s-th row and the g-th column, (j, 2 l) represents the j-th row and the 2 l-th column, (i, 2 l-1) represents the i-th row and the 2l-1 column; .
Formulas (9), (17) may be written in the following form:
in the formula ：
two sides of (20) are simultaneously multiplied byThen equation (20) becomes:
x t+1 ＝ Fx t + u t+1 (21)
wherein: system equation matrixControl vector->
Equation (21) is a system equation of the natural gas pipeline network, and is used for representing a dynamic process.
(32) Obtaining measurement equation z t+1 ＝H t+1 x t+1
Pressure and flow measurement information from the monitoring control and data acquisition system provides redundancy constraints for dynamic state estimation. Assuming that all nodes are equipped with flow meters and barometers, and taking these information as measurement vectors, the measurement equation is:
z t+1 ＝ H t+1 x t+1 (22)
in the formula ：zt+1 and Ht+1 Respectively measuring vectors and measuring matrixes at the time t+1; i is an identity matrix; z r1,t+1 ,z r2,t+1 ,/>Is a measurement of pressure at time t+1, z m1,t+1 ,z m2,t+1 ,…,/>Is a measurement of flow at time t+1; 0 is zero matrix, and subscript is zero matrix dimension; n is n N and nL The number of nodes and the number of pipes, respectively.
(33) Obtaining a natural gas pipe network dynamic state estimation model
in the formula ：vt and wk+1 Other than the system and the measurement error vector,u k+1 is a control vector. v t+1 and wk+1 The variance matrix is Q t+1 and Rt+1 。
The step (4) comprises the following steps:
the state quantity predicted valueAnd a prediction covariance matrix P t+1|t The method comprises the following steps:
P t+1|t ＝ FP t|t F T + Q t (26)
in the formula , and />Respectively a state quantity predicted value and an estimated value; p (P) t+1|t and Pt|t The prediction error variance and the estimation error variance, respectively.
The pair of state quantity predictorsCorrecting and updating covariance matrix P t+1|t+1 The method comprises the following steps:
K k+1 ＝ P t+1|t H T (HP t+1|t H T + R t+1 ) -1 (27)
P t+1|t+1 ＝ P t+1|t -K k+1 HP t+1|t (29)
in the formula ：Kk+1 Is a gain matrix; p (P) e,t+1 ＝HP t+1|t H T +R t+1 Is an innovation covariance matrix, which represents the error between a measured value and its predicted value,an estimate of time t+1.
time-varying scale factor mu at time t+1 t+1 The following conditions are satisfied:
P e,t+1 ＝HP t+1|t H+μ t+1 R t+1 (30)
estimating an innovation covariance matrix by adopting a sliding window method:
wherein ：is an innovation vector; m is m W Is the window length. Mu from (31) t+1
Mu obtained t+1 May be off-diagonal, for which a diagonal matrix μ is defined t ′ +1 ：
wherein ：μi ′＝max(1,μ t+1,ii )，i＝1,2,…,2n N ,μ t+1,ii Mu is t+1 The i-th diagonal element. Formula (27) becomes:
K k+1 ＝ P t+1|t H T (HP t+1|t H T + μ t ′ +1 R t+1 ) -1 (34)
the dynamic state estimation method is applied to the 30-node natural gas pipe network shown in fig. 2. And (5) calculating the dynamic process of the node pressure and flow by using the formulas (6) to (8). And taking the simulation result as a true value, and adding a random error to the true value to obtain a measured value. The state estimation is carried out on the natural gas pipe network under the three measurement conditions of normal, bad data and biased noise by adopting the traditional Kalman filtering and the robust Kalman filtering algorithm.
To further illustrate the present invention, the above schemes are simulated by using specific examples, and the parameters of the test system are shown in table 1. The pressure was constant at source nodes 1 and 2, 27.8bar and 28.5bar, respectively. The simulation time was 24 hours with 15 minutes intervals. Standard deviations of the pressure and flow measurements were 0.01bar and 2%, respectively.
Table 1 natural gas pipeline data
The estimation result is evaluated by the following filter coefficients:
wherein epsilon is a filter coefficient;is an estimated value; />The true value of the measurement quantity.
The filter coefficients of the two DSE methods are shown in table 2. Node 1 and node 2 are source nodes whose pressures are constant and filter coefficients are not calculated. As can be seen from table 2, the filter coefficients are all less than 1, indicating that the robust kalman filter algorithm based on kalman filter and described in the present invention is effective. In addition, most of the filter coefficients of the robust Kalman filter are smaller than those of the Kalman filter, which means that the robust Kalman filter has better performance.
Table 2 filter coefficients
Bad data is artificially added to the measurement vector based on normal measurement conditions. The pressure measurements at node 30 were set at 12, 10.7 and 13.8bar at times 5, 5.25 and 5.5 hours, respectively; set to 13, 15.5 and 23bar at times 13.25, 13.5 and 13.75 hours, respectively; the flow measurement at node 11 was set to 3, 2.1, 3 and 2.2kg/s at times 7.5, 7.75, 8 and 8.25 hours, respectively; reset to 3, 2.1 and 1.7kg/s at times 15.75, 16 and 16.25 hours, respectively.
The estimation results of the Kalman filtering and the robust Kalman filtering algorithm are shown in figure 3, and the time-varying scale factors are shown in figure 4. It can be seen that the kalman filter curve deviates from the true value when bad data occurs. This is because the measurement error variance matrix does not conform to the actual situation under poor data conditions. However, the curve and the true-value curve of the robust kalman filter always agree well. When bad data appear, the scale factor increases sharply, the correction effect of the measured value decreases, and the estimated value is not affected by the bad data.
For biased noise, constant deviations of 0.2bar and 0.1kg/s were added to the pressure measurement and flow measurement for 10 to 19.75 hours and 5 to 12.5 hours, respectively. The estimation results are shown in fig. 5 and 6. It can be seen that under biased noise conditions, the estimated curve of the kalman filter has a significant deviation from the true curve. When the measured mean value deviates from zero, the scale factor is obviously increased, and an estimated curve based on the robust Kalman filtering is always consistent with a true value.
Claims (3)
1. The robust dynamic state estimation method of the natural gas pipe network is characterized by comprising the following steps of:
(1) Obtaining a natural gas pipeline dynamic mathematical model;
(2) Obtaining a natural gas pipe network boundary condition;
(3) Obtaining a model of natural gas pipe network dynamic state estimation;
(4) Obtaining a state quantity estimated value of a natural gas pipe network by adopting a robust Kalman filtering algorithm;
the step (1) of obtaining the natural gas pipeline dynamic mathematical model comprises the following steps:
(11) Obtaining partial differential mathematical model of natural gas pipeline
Wherein ρ is the density of the gas in the pipeline, p is the pressure of the gas in the pipeline, v G The gas velocity in the pipeline is t, time, zeta is a space coordinate, f is a friction factor, G is a gravity constant, theta is a pipeline inclination angle, and d is a pipeline diameter;
(12) Obtaining a simplified partial differential mathematical model of a natural gas pipeline
wherein ,p＝c2 ρ，c 2 ＝ZRT，/>Is the flow of the pipeline, a is the cross-sectional area of the pipeline, Z is the gas compressibility factor, R is the gas constant, and T isThe gas temperature, c is the sound velocity, c for natural gas if the temperature is unchanged 2 Constant (I)>Is the average value of the gas flow rate;
(13) Obtaining a natural gas pipeline differential mathematical model
Wherein i, j, k are three nodes of the pipe network, and k<i<j, gas flows into the large node from the smaller node, and the subscript ij is the number of the head and the tail nodes of the pipeline, L ij For the length of pipe i-j ρ i,t and ρi,t+1 The densities of the gas at the node i at the time t and the time t+1 are respectively, ρ j,t and ρj,t+1 The density of the gas at the node j at the time t and t+1, a ij For the cross-sectional area of conduit i-j, d ij For the diameter of the pipe i-j, Δt is the time step, and />The flow from node i to node j at times t and t+1, respectively, +.>Andrespectively defining positive directions as outflow and negative flow values of source nodes for flows flowing from the node j to the node i at times t and t+1;
the step (2) of obtaining the natural gas pipe network boundary condition comprises the following steps:
(21) Boundary conditions for sink node flow balancing
wherein ,representing the flow from node i to node k at time t, k e i representing that node k is connected to node i through pipe i-k; i epsilon N Sink Indicating i is a sink node, N Sink Indicates the number of bus nodes, ">For the traffic of node i at time t, +.> and />Representing the sum of traffic flowing from and into node i;
representing the sum of traffic flowing from and into node i;
(22) Boundary condition of source node gas density
ρ g,t ＝ρ Sg,0 ,g∈N Source (8)
wherein ,ρSg,0 The gas density of the source node g at the initial moment; n (N) Source Represents the number of source nodes, g.epsilon.N Source G is represented as a source node;
the step (3) of obtaining the model of the natural gas pipe network dynamic state estimation comprises the following steps:
(31) Obtaining system equation x t+1 ＝Fx t +u t+1 (21)
wherein ,xt and xt+1 State vectors, x, at times t and t+1, respectively t ＝[x r,t ,x m,t ] T ， x r,t Is the state vector of the gas density, x m,t Is the state vector of the flow rate, ρ 1,t ,ρ 2,t ,/>Nodes 1,2 and n at time t respectively N Density, x m,t Is 2n L ，n N and nL Representing the number of nodes and the number of pipes, a system equation matrix +.>Control amount->
n in and nout The number of injection nodes and the number of outflow nodes are respectively represented; a is that 11 ，A 12 ，A 21 ，A 22 ，B 11 ，B 22 A submatrix, a zero matrix and a matrix dimension as subscripts; u (u) r,t+1 ，u m,t+1 Density and flow control amount, respectively;(l, i) represents the first row and the first column; i.alpha.l indicates that node i is one end of pipe l, +.>As a real number set，
α 1 ，α 2 ，α l ，/>As matrix A 12 Middle element, beta l As matrix A 21 Middle element, gamma l As matrix A 22 The element j-l represents that the node j is one end of the pipeline l, (s, g) represents the s-th row and the g-th column, (j, 2 l) represents the j-th row and the 2 l-th column, (i, 2 l-1) represents the i-th row and the 2 l-1-th column;
(32) Obtaining measurement equations
z t+1 ＝H t+1 x t+1 (22)
wherein ,zt+1 Measuring vector, H for time t+1 t+1 Measuring a matrix for the time t+1;
i is an identity matrix; z r1,t+1 ,z r2,t+1 ,/>Is a measurement of pressure at time t+1, z m1，t+1 ，Is a measurement of flow at time t+1;
(33) Model for obtaining natural gas pipe network dynamic state estimation
wherein ,vt Is the systematic error vector, w k+1 Is the vector of the measurement error and,u k+1 for the control vector, H is the measurement equation matrix, v t+1 Variance matrix is Q t+1 ，w t+1 Variance matrix R t+1 ；
The step (4) of obtaining the state quantity estimated value of the natural gas pipe network by adopting a robust Kalman filtering algorithm comprises the following steps:
(41) According to the state quantity estimated value at the time tCovariance matrix P t|t Obtaining state quantity predictive value ∈ ->And a prediction covariance matrix P t+1|t
P t+1|t ＝FP t|t F T +Q t (26)
wherein ,predicted value of state quantity at time t+1, < >>The state quantity estimated value at the moment t; p (P) t+1|t and Pt|t The prediction error variance and the estimation error variance, Q t Measuring an error equation array for the time t;
(42) Correction state quantity predicted valueUpdating covariance matrix P t+1|t+1
K k+1 ＝P t+1|t H T (P e,t+1 ) -1 (27)
P t+1|t+1 ＝P t+1|t -K k+1 HP t+1|t (29)
wherein ,Kk+1 Is a gain matrix; p (P) e,t+1 Is an innovation matrix;the state quantity estimated value at time t+1.
2. The method of robust dynamic state estimation for a natural gas network according to claim 1, wherein in step (42)P e，t+1 ＝HP t+1|t H+μ t+1 R t+1 Time-varying scale factor For the innovation vector, m W Is the window length.
3. The method of robust dynamic state estimation for a natural gas network according to claim 1, wherein P in step (42) e,t+1 ＝HP t+1|t H T +μ′ t+1 R t+1 ，μ′ i ＝max(1，μ t+1，ii )，i＝1，2，…，2n N ,μ t+1,ii Is mu t+1 Is the i-th diagonal element of (c).
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110177597.7A CN112906317B (en) | 2021-02-09 | 2021-02-09 | Robust dynamic state estimation method for natural gas pipe network |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110177597.7A CN112906317B (en) | 2021-02-09 | 2021-02-09 | Robust dynamic state estimation method for natural gas pipe network |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN112906317A CN112906317A (en) | 2021-06-04 |
| CN112906317B true CN112906317B (en) | 2023-08-22 |
Family
ID=76122985
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202110177597.7A Active CN112906317B (en) | 2021-02-09 | 2021-02-09 | Robust dynamic state estimation method for natural gas pipe network |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN112906317B (en) |
Citations (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO1999010783A1 (en) * | 1997-08-22 | 1999-03-04 | Voyan Technology | A method for real-time nonlinear system state estimation and control |
| JP2014215822A (en) * | 2013-04-25 | 2014-11-17 | 日本電信電話株式会社 | State estimating apparatus, method, and program |
| CN107590317A (en) * | 2017-08-17 | 2018-01-16 | 河海大学 | A kind of generator method for dynamic estimation of meter and model parameter uncertainty |
| CN109376910A (en) * | 2018-09-28 | 2019-02-22 | 河海大学 | A kind of power distribution network dynamic state estimator method based on historical data driving |
| CN110222309A (en) * | 2019-05-06 | 2019-09-10 | 河海大学 | A kind of generator method for dynamic estimation based on robust volume Kalman filtering |
| CN110532517A (en) * | 2019-08-28 | 2019-12-03 | 北京理工大学 | Gas pipeline method for parameter estimation based on improved ARUKF |
| CN110619487A (en) * | 2019-10-12 | 2019-12-27 | 东北大学 | Electric-gas-thermal coupling network dynamic state estimation method based on Kalman filtering |
| CN111898816A (en) * | 2020-07-21 | 2020-11-06 | 南京信息工程大学 | Dynamic state estimation method for comprehensive energy system |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US8706424B2 (en) * | 2009-06-30 | 2014-04-22 | H2Scan Corporation | System for estimating a gas concentration in a mixed atmosphere |
-
2021
- 2021-02-09 CN CN202110177597.7A patent/CN112906317B/en active Active
Patent Citations (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO1999010783A1 (en) * | 1997-08-22 | 1999-03-04 | Voyan Technology | A method for real-time nonlinear system state estimation and control |
| JP2014215822A (en) * | 2013-04-25 | 2014-11-17 | 日本電信電話株式会社 | State estimating apparatus, method, and program |
| CN107590317A (en) * | 2017-08-17 | 2018-01-16 | 河海大学 | A kind of generator method for dynamic estimation of meter and model parameter uncertainty |
| CN109376910A (en) * | 2018-09-28 | 2019-02-22 | 河海大学 | A kind of power distribution network dynamic state estimator method based on historical data driving |
| CN110222309A (en) * | 2019-05-06 | 2019-09-10 | 河海大学 | A kind of generator method for dynamic estimation based on robust volume Kalman filtering |
| CN110532517A (en) * | 2019-08-28 | 2019-12-03 | 北京理工大学 | Gas pipeline method for parameter estimation based on improved ARUKF |
| CN110619487A (en) * | 2019-10-12 | 2019-12-27 | 东北大学 | Electric-gas-thermal coupling network dynamic state estimation method based on Kalman filtering |
| CN111898816A (en) * | 2020-07-21 | 2020-11-06 | 南京信息工程大学 | Dynamic state estimation method for comprehensive energy system |
Non-Patent Citations (1)
| Title |
|---|
| "Trajectory tracking method of natural gas, district heating and power systems";Liang Chen等;Energy Conversion and Management;第259卷;正文第1-12页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN112906317A (en) | 2021-06-04 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| US11922335B2 (en) | Method and system for evaluating macro resilience of offshore oil well control equipment | |
| CN102735386B (en) | Bending stiffness-considered numerical computation method for stay cable forces | |
| CA3060419C (en) | Detecting and correcting for discrepancy events in fluid pipelines | |
| US7213454B2 (en) | Method and apparatus for obtaining improved accuracy and range for air data parameters inferred from independent measurements of interdependent pressures | |
| CN109827629B (en) | Distributed reliability estimation method for urban river water level | |
| CN111852463B (en) | Gas well productivity evaluation method and equipment | |
| CN108873091B (en) | The full tensor of Satellite gravity field restores the determination method and system of earth gravitational field | |
| Moutsoglou | An inverse convection problem | |
| CN113551744A (en) | Ultrasonic flowmeter performance online monitoring method and system | |
| Ishido et al. | A new indicator for real-time leak detection in water distribution networks: design and simulation validation | |
| CA3019887A1 (en) | Method and apparatus for estimating down-hole process variables of gas lift system | |
| Chen et al. | Stress influence line identification of long suspension bridges installed with structural health monitoring systems | |
| CN112906317B (en) | Robust dynamic state estimation method for natural gas pipe network | |
| Zhao et al. | Evaluation and early warning of vortex-induced vibration of existed long-span suspension bridge using multisource monitoring data | |
| Zhu et al. | Numerical study on the displacement reconstruction of subsea pipelines using the improved inverse finite element method | |
| KR101348817B1 (en) | Method and system for estimation of error or omitted measured values using artificial neural network | |
| So et al. | Modeling Reynolds-number effects in wall-bounded turbulent flows | |
| Chu et al. | Life-cycle assessment of long-span bridge’s wind resistant performance considering multisource time-variant effects and uncertainties | |
| Sollund et al. | Effects of seabed topography on modal analyses of free spanning pipelines | |
| CN112925206B (en) | Distributed robust fault diagnosis method for nonlinear multi-inverted pendulum interconnection system | |
| KR20180001617A (en) | Displacement and inclination data fusion method for estimating structural deformation | |
| CN114063456B (en) | Fault prediction and early warning method using autoregressive model and Kalman filtering algorithm | |
| CN115047213B (en) | Method for improving long-term stability of MEMS accelerometer | |
| McNeill et al. | A Method for Estimating Quasi-Static Riser Deformation and Applied Forces From Sparse Riser Inclination Measurements | |
| CN113916181B (en) | Data processing method of surface-internal integrated deformation monitoring device |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |