一种天然气管网的鲁棒动态状态估计方法

1.本发明涉及一种动态状态估计方法，尤其涉及一种天然气管网的鲁棒动态状态估计方法。


背景技术：

2.随着环境问题的日益严峻，可再生能源发电系统的规模不断扩大。由于可再生能源的不确定性，天然气的消耗量随机变化一定程度上影响系统的安全稳定运行。因此获得天然气管网压强和流量等精确的动态特性，是保证含电力和天然气综合能源系统安全、优化运行的关键。为了获取这些动态状态，管网必须配备大量的监测设备，这需要大量的投资，即便如此，也不可能在网络的每个节点上安装传感器并获得所有状态。另一方面，测量装置不可避免地会遇到随机误差和不良数据，因此需要对收集的数据进行状态估计，才能够应用于泄漏检测和控制问题。
3.现有的天然气管网动态状态估计主要是基于卡尔曼滤波的输气管网状态估计方法，即通过一种称为有限元的数值逼近技术来求解瞬态流动方程组。这些方法在随着规模的扩大，存在着严重的局限性，如可调整性、可扩展性、不可接受的计算负载以及在传感器故障时缺乏鲁棒性。在实际系统中，由于传感器误差和电磁干扰，监控和数据管理系统不可避免地会出现随机误差和不良数据。现有的基于卡尔曼滤波的状态估计方法在一定程度上减少了随机误差，但容易受到坏数据的影响。因此在准确地获得大型天然气管网在恶劣数据和有色噪声条件下的瞬态状态的过程中增强系统的鲁棒性对于天然气传输管网的安全稳定运行以及优化控制更加重要。


技术实现要素：

4.发明目的：本发明旨在解决现有技术的上述不足，提供一种天然气管网的鲁棒动态状态估计方法，解决在恶劣数据和有色噪声条件下天然气管网状态估计的鲁棒性问题。
5.技术方案：本发明所述的一种天然气管网的鲁棒动态状态估计方法，包括以下步骤：
6.(1)获得天然气管道动态数学模型；
7.(2)获得天然气管网边界条件；
8.(3)获得天然气管网动态状态估计的模型；
9.(4)采用鲁棒卡尔曼滤波算法获得天然气管网的状态量估计值。
10.步骤(1)中包括以下步骤：
11.(11)获得天然气管道偏微分数学模型
[0012][0013]
[0014]
其中，ρ为管道中气体密度，p为管道中气体压强，v
g
为管道中气体速度，t为时间，ζ为空间坐标，f为摩擦因子，g为重力常数，θ为管道倾角，d为管道直径；
[0015]
(12)获得天然气管道简化的偏微分数学模型
[0016][0017][0018]
其中，p＝c2ρ，c2＝zrt，为管道的流量，a为管道的截面面积，z为气体压缩性因子，r为气体常数，t为气体温度，c为声速，对于天然气如果温度不变，c2恒定，为气体流速平均值；
[0019]
(13)获得天然气管道差分数学模型
[0020][0021][0022]
其中，i，j，k为管网的三个节点，且k<i<j，气体从较小节点流入大节点，下标ij为管道首末节点号，l
ij
为管道i
‑
j的长度，ρ
i,t
和ρ
i,t+1
分别为节点i处气体在t和t+1时刻的密度，ρ
j,t
和ρ
j,t+1
分别为节点j处气体在t和t+1时刻的密度，a
ij
为管道i
‑
j横截面积，d
ij
为管道i
‑
j的直径，δt为时间步长，和分别为t和t+1时刻从节点i流至节点j的流量，和分别为t和t+1时刻从节点j流至节点i的流量，定义正方向为流出，源节点的流量值为负数。
[0023]
步骤(2)中包括以下步骤：
[0024]
(21)汇流节点流量平衡的边界条件
[0025][0026]
其中，表示t时刻从节点i流至节点k的流量，k∈i表示节点k通过管道i
‑
k连接到节点i；i∈n
sink
表示i是汇流节点，n
sink
表示汇流节点数，为节点i在t时刻的流量，和表示从节点i流出和流入节点i的流量之和。
[0027]
(22)源节点气体密度的边界条件
[0028]
ρ
g,t ＝ ρ
sg,0
, g∈n
source
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(8)
[0029]
其中，ρ
sg,0
为初始时刻源节点g的气体密度；n
source
表示源节点数，g∈n
source
表示g为源节点。
[0030]
步骤(3)中包括以下步骤：
[0031]
(31)获得系统方程x
t+1
＝fx
t
+u
t+1
(21)
[0032]
其中，x
t
和x
t+1
分别为t时刻和t+1时刻状态向量，x
t
＝[x
r,t
,x
m,t
]
t
，，i<j，x
r,t
为气体密度的状态向量，x
m,t
为流量的状态向量，ρ
1,t
,ρ
2,t
,分别为t时刻节点1、2和n
n
的密度，x
m,t
的维数是2n
l
，n
n
和n
l
表示节点数量和管道数量，系统方程矩阵控制量控制量n
in
和n
out
分别表示注入节点数和流出节点数；a
11
，a
12
，a
21
，a
22
，b
11
，b
22
为子矩阵，0为零矩阵，下标为矩阵维数；u
r,t+1
，u
m,t+1
分别为密度和流量控制量；(l,i)表示第l行第i列；i
∝
l表示节点i是管道l的一端，为实数集，为实数集，为实数集，为实数集，为实数集，α1，α2，α
l
，为矩阵a
12
中元素，β
l
为矩阵a
21
中元素，γ
l
为矩阵a
22
中元素，j
∝
l表示节点j是管道l的一端，(s,g)表示第s行第g列，(j,2l)表示第j行第2l列,(i,2l
‑
1)表示第i行第2l
‑
1列；
[0033]
(32)获得测量方程
[0034]
z
t+1 ＝ h
t+1
x
t+1
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(22)
[0035][0036]
其中，z
t+1
为t+1时刻测量向量、h
t+1
为t+1时刻测量矩阵；为t+1时刻测量矩阵；i为单位矩阵；z
r1,t+1
,z
r2,t+1
,是t+1时刻压强的测量值，z
m1,t+1
,z
m2,t+1
,
…
,是t+1时刻流量的测量值；
[0037]
(33)获得天然气管网动态状态估计的模型
[0038]
[0039]
其中，v
t
是系统误差向量、w
k+1
是测量误差向量，u
k+1
为控制向量，h为测量方程矩阵，v
t+1
方差矩阵为q
t+1
，w
t+1
方差矩阵为r
t+1
。
[0040]
步骤(4)中包括以下步骤：
[0041]
(41)根据t时刻状态量估计值及其协方差矩阵p
t|t
，获得状态量预测值和其预测协方差矩阵p
t+1|t
[0042][0043]
p
t+1|t ＝ fp
t|t
f
t + q
t
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(26)
[0044]
其中，为t+1时刻状态量预测值、为t时刻状态量估计值；p
t+1|t
和p
t|t
分别为预报误差方差和估计误差方差，q
t
为t时刻测量误差方程阵；
[0045]
(42)修正状态量预测值更新协方差矩阵p
t+1|t+1
[0046]
k
k+1 ＝ p
t+1|t
h
t
(p
e,t+1 )
‑1ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(27)
[0047][0048]
p
t+1|t+1 ＝ p
t+1|t
‑
k
k+1
hp
t+1|t
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(29)
[0049]
其中，k
k+1
为增益矩阵；p
e,t+1
为新息矩阵；为t+1时刻状态量估计值。
[0050]
步骤(42)中p
e,t+1
＝hp
t+1|t
h+μ
t+1
r
t+1
，时变尺度因子，时变尺度因子为新息向量，m
w
为窗口长度。
[0051]
步骤(42)中p
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
是μ
t+1
的第i个对角线元素。
[0052]
有益效果：与现有技术相比，本发明的显著优点：(1)根据新息对测量方差矩阵进行调整，使测量能够准确地修正状态量预测值，提高卡尔曼滤波器对坏数据和有色噪声的鲁棒性，获得大型天然气管网在含有坏数据和有色噪声条件下的准确状态。(2)将汇流节点处的零注入流量作为虚拟测量，提高了测量冗余度。(3)用平均气体流速近似二阶项，以管道两端节点压强和流量为状态量，使状态量不会随着网络规模的增加而急剧增加，提高了效率。
附图说明
[0053]
图1为本发明天然气管道的数学模型；
[0054]
图2为本发明测试系统示意图；
[0055]
图3(a)
‑
图3(b)为本发明坏数据条件下压强和流量估算结果；
[0056]
图4(a)
‑
图4(b)为本发明坏数据条件下压强和流量时变尺度因子；
[0057]
图5(a)
‑
图5(b)为本发明有色噪声条件下压强和流量的估算结果；
[0058]
图6(a)
‑
图6(b)为本发明有色噪声条件下压强和流量的时变尺度因子。
具体实施方式
[0059]
下面结合附图对本发明的技术方案作进一步说明。
[0060]
本发明所述的天然气管网的鲁棒动态估计方法，包括以下步骤：
[0061]
(1)建立天然气管道动态数学模型；
[0062]
(2)建立天然气管网边界条件；
[0063]
(3)建立天然气管网动态状态估计的模型；
[0064]
(4)采用鲁棒卡尔曼滤波算法获得天然气管网的状态量估计值。
[0065]
步骤(1)中包括以下步骤：
[0066]
(11)获得天然气管道偏微分数学模型
[0067][0068][0069]
其中：ρ为管道中气体密度，p管道中气体压强，v
g
管道中气体速度，t为时间，ζ为空间坐标，f为摩擦因子；g为重力常数；θ为管道倾角；d为管道直径。
[0070]
(12)获得天然气输气管道的简化偏微分数学模型为：
[0071]
在气体压强小于10bar，温度低于20℃的条件下，压强，流量与气体密度之间的函数：
[0072]
p ＝ c2ρ
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(2)
[0073][0074]
式中：c2＝zrt，和a分别为管道的流量和截面面积；z，r和t分别是气体压缩系数、气体常数和气体温度。c为声速，对于天然气，如果温度不变，c2是恒定的。这里用气体流速平均值的绝对值代替|u|。
[0075]
由于天然气流速远小于声速，且假设管道是水平的，则在上述假设下，天然气管道的简化偏微分数学模型变为：
[0076][0077][0078]
其中，为气体流速平均值。
[0079]
(13)获得天然气管道差分数学模型为：
[0080]
偏微分方程离散化的有限微分过程为：
[0081][0082][0083][0084]
式中：x表示广义变量；δt和δx分别是时间和空间步长；下标s和t分别为节点号和时间；x
s+1,t+1
和x
s+1,t
分别为t+1时刻和t时刻s+1处的变量值；x
s,t+1
和x
s,t
分别为t+1时刻和t时刻s处的变量值。方程(4)可以转换为天然气输气管道差分数学模型：
[0085][0086][0087]
其中，i，j，k为管网的三个节点，且k<i<j，气体从较小节点流入大节点。下标ij为管道首末节点号；l
ij
为管道i
‑
j的长度，ρ
i,t
和ρ
i,t+1
分别为节点i处气体在t和t+1时刻的密度，ρ
j,t
和ρ
j,t+1
分别为节点j处气体在t和t+1时刻的密度，a
ij
为管道i
‑
j横截面积，d
ij
表示管道i
‑
j的直径，δt为时间步长，和分别表示t和t+1时刻从节点i流至节点j的流量，和分别表示t和t+1时刻从节点j流至节点i的流量，定义正方向为流出，因此源节点的流量值为负数。
[0088]
步骤(2)包括以下步骤：
[0089]
(21)汇流节点流量平衡的边界条件表示
[0090]
边界条件为气源节点压强恒定以及汇流节点流量和为0。汇流节点处流量平衡的边界条件表示为：
[0091][0092]
式中：表示t时刻从节点i流至节点k的流量，k∈i表示节点k通过管道i
‑
k连接到节点i；i∈n
sink
表示i为源节点，n
sink
表示汇流节点数，为节点i在t时刻的流量。和分别表示流出节点i和流入节点i的流量之和。
[0093]
(22)源节点气体密度的边界条件
[0094]
对于源节点，假设气体容量足够大，使得气体密度能够在动态过程中保持恒定，形成了以下边界条件：
[0095]
ρ
g,t ＝ ρ
sg,0
, g∈n
source
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(8)
[0096]
式中：ρ
sg,0
为初始时刻源节点g的气体密度；n
source
表示源节点数，g∈n
source
表示g为源节点。
[0097]
步骤(3)中包括以下步骤：
[0098]
(31)获得系统方程x
t+1
＝fx
t
+u
t+1
[0099]
动态状态估计方程以节点处的气体密度和管道两端的流量为状态量。对于具有n
n
节点和n
l
管道的天然气管网，状态量个数为n
n
+2n
l
。x
t
和x
t+1
分别为t时刻和t+1时刻状态向量，x
t
＝[x
r,t
,x
m,t
]
t
，x
r,t
是t时刻气体密度的状态向量；x
m,t
是t时刻流量的状态向量，ρ
1,t
,ρ
2,t
,分别为t时刻节点1，2和n
n
的密度，x
m,t
的维数是2n
l
，n
n
和n
l
表示节点数量和管道数量。根据上述状态量定义，式(6)可以改写为如下矩阵形式：
[0100]
[0101][0102][0103][0104][0105][0106][0107][0108]
其中：为实数集，x
r,t+1
是t+1时刻气体密度的状态向量；x
m,t+1
是t+1时刻流量的状态向量，a
11
，a
12
，a
21
，a
22
，为子矩阵；α1，α2，α
l
，为矩阵a
12
中元素，β
l
为矩阵a
21
中元素，γ
l
为矩阵a
22
中元素，(l，i)表示第l行和第i列；i
∝
l，j
∝
l表示节点i和节点j是管道l的一端。对于天然气管网，状态量的总数为n
n
+2n
l
。
[0109]
方程数量小于状态量个数，需要边界条件作为补充约束。式(7)和(8)可以写为如下矩阵形式：
[0110][0111][0112][0113][0114]
其中，b
11
，b
22
为子矩阵，0为零矩阵，下标为矩阵维数；u
r,t+1
，u
m,t+1
分别为密度和流量控制量；n
in
和n
out
分别表示注入节点数和流出节点数；(s,g)表示第s行第g列，(j,2l)表示第j行第2l列,(i,2l
‑
1)表示第i行第2l
‑
1列；。
[0115]
式(9)、(17)可以写成以下形式：
[0116]
[0117]
式中：
[0118]
式(20)两边同时左乘则式(20)变为：
[0119]
x
t+1 ＝ fx
t + u
t+1
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(21)
[0120]
式中：系统方程矩阵控制向量
[0121]
式(21)为天然气管网的系统方程，用于表示动态过程。
[0122]
(32)获得测量方程z
t+1
＝h
t+1
x
t+1
[0123]
来自监测控制和数据采集系统的压强和流量测量信息为动态状态估计提供了冗余约束。假设所有节点都配备了流量计和气压计，将这些信息作为测量向量，则测量方程为：
[0124]
z
t+1 ＝ h
t+1
x
t+1
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(22)
[0125][0126]
式中：z
t+1
和h
t+1
分别为t+1时刻测量向量和测量矩阵；分别为t+1时刻测量向量和测量矩阵；i为单位矩阵；z
r1,t+1
,z
r2,t+1
,是t+1时刻压强的测量值，z
m1,t+1
,z
m2,t+1
,
…
,是t+1时刻流量的测量值；0为零矩阵，下标为零矩阵维数；n
n
和n
l
分别为节点数和管道数。
[0127]
(33)获得天然气管网动态状态估计模型
[0128][0129]
式中：v
t
和w
k+1
别是系统和测量误差向量，u
k+1
为控制向量。v
t+1
和w
k+1
的方差矩阵分别为q
t+1
和r
t+1
。
[0130]
步骤(4)中包括以下步骤：
[0131]
所述状态量预测值和其预测协方差矩阵p
t+1|t
为：
[0132][0133]
p
t+1|t ＝ fp
t|t
f
t + q
t
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(26)
[0134]
式中，和分别为状态量预测值和估计值；p
t+1|t
和p
t|t
分别为预报误差方差和估计误差方差。
[0135]
所述对状态量预测值进行修正，更新协方差矩阵p
t+1|t+1
为：
[0136]
k
k+1 ＝ p
t+1|t
h
t
(hp
t+1|t
h
t + r
t+1
)
‑1ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(27)
[0137][0138]
p
t+1|t+1 ＝ p
t+1|t
‑
k
k+1
hp
t+1|t
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(29)
[0139]
式中：k
k+1
为增益矩阵；p
e,t+1
＝hp
t+1|t
h
t
+r
t+1
，为新息协方差矩阵，它表示测量值与
其预测值之间的误差，为t+1时刻估计值。
[0140]
t+1时刻时变尺度因子μ
t+1
满足如下条件：
[0141]
p
e,t+1
＝hp
t+1|t
h+μ
t+1
r
t+1
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(30)
[0142]
采用滑动窗口法估计新息协方差矩阵：
[0143][0144]
其中：为新息向量；m
w
为窗口长度。由(31)得到μ
t+1
[0145][0146]
得到的μ
t+1
可能是非对角占优的，针对该问题，定义对角矩阵μ
t
′
+1
：
[0147][0148]
其中：μ
i
′
＝max(1,μ
t+1,ii
)，i＝1,2,
…
,2n
n
,μ
t+1,ii
为μ
t+1
第i个对角线元素。式(27)变为：
[0149]
k
k+1 ＝ p
t+1|t
h
t
(hp
t+1|t
h
t + μ
t
′
+1
r
t+1
)
‑1ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(34)
[0150]
将所述动态状态估计方法应用于图2所示的30节点天然气管网。用公式(6)～(8)计算节点压强和流量的动态过程。以仿真结果为真值，在真值上加随机误差得到测量值。采用传统的卡尔曼滤波和本发明所述的鲁棒卡尔曼滤波算法对天然气管网在正常、不良数据和有偏噪声三种测量条件下进行状态估计。
[0151]
为了进一步说明本发明，下面再以具体的实例对上述方案进行仿真测试，测试系统参数见表1。在源节点1和2压强恒定，分别为27.8bar和28.5bar。仿真时间为24小时，间隔为15分钟。压强和流量测量的标准偏差分别为0.01bar和2％。
[0152]
表1天然气管道数据
[0153][0154]
通过如下滤波系数对估计结果进行评价：
[0155][0156]
其中，ε为滤波系数；为估计值；为测量量真值。
[0157]
两种dse方法的滤波系数如表2所示。节点1和节点2是源节点，其压强恒定，不计算滤波器系数。由表2可见，滤波系数均小于1，说明基于卡尔曼滤波和本发明所述的鲁棒卡尔曼滤波算法都是有效的。另外，本发明所述的鲁棒卡尔曼滤波的大部分滤波系数都比卡尔曼滤波小，意味着鲁棒卡尔曼滤波的性能更好。
[0158]
表2滤波系数
[0159][0160]
基于正常的测量条件，人为地将不良数据添加到测量向量中。节点30的压强测量值在时间5、5.25和5.5小时分别设置为12、10.7和13.8bar；在时间13.25、13.5和13.75小时分别设置为13、15.5和23bar；节点11的流量测量值在时间7.5、7.75、8和8.25小时分别设置为3、2.1、3和2.2kg/s；在时间15.75、16和16.25小时分别重置为3、2.1和1.7kg/s。
[0161]
卡尔曼滤波和本发明所述的鲁棒卡尔曼滤波算法的估计结果如图3所示，时变尺度因子如图4所示。可见，当不良数据出现时，卡尔曼滤波曲线偏离真值。这是由于在不良数据条件下，测量误差方差矩阵与实际情况不符。然而，鲁棒卡尔曼滤波的曲线和真值曲线始终很好地吻合。当不良数据出现时，尺度因子急剧增加，测量值的校正效果降低，估计值不受不良数据的影响。
[0162]
针对有偏噪声，分别在10～19.75小时和5～12.5小时的压强测量和流量测量中分别加入0.2bar和0.1kg/s的恒定偏差。估算结果如图5和图6所示。可见在有偏噪声条件下，卡尔曼滤波的估计曲线与真值曲线有明显的偏差。当测量均值偏离零时，尺度因子明显增大，基于鲁棒卡尔曼滤波的估计曲线始终与真值相吻合。