考虑边界等式约束的电气热综合能源系统鲁棒状态估计方法

1.本发明涉及综合能源系统鲁棒状态估计方法的技术领域，更具体地说涉及电气热综合能源系统鲁棒状态估计方法的技术领域。


背景技术：

2.随着能源危机和环境污染日益严重，如何高效、清洁利用各类能源已成为近年来研究的重点。在这种情况下，综合能源系统(integrated energy system，ies)应运而生。
3.综合能源系统ies内可含电、气、热和冷等多种形式能流，集多种能源的生产、输送、分配、转换、储存和消费于一体，可为综合能源的高效、清洁利用提供重要解决方案。一方面，ies中的热电联产、热泵和电锅炉等各种能量转换设备，可使得各能源子系统之间耦合紧密。多种能源的转换与互动，可充分提高能源系统的灵活性和能源使用效率。另一方面，不同能源子系统之间相互协调，可促进可再生能源的消费，减少化石能源消耗，对环境友好。因此，综合能源系统ies已成为能源结构调整的重要方向。
4.随着综合能源系统ies的快速发展，综合能源系统ies能量管理系统迫切需要获取系统的实时运行状态，为后续的在线分析和多能流优化调度提供完整可靠的实时数据支持。因此，有必要为ies引入状态估计，通过充分利用测量冗余来消除测量误差和不良数据干扰，从而提供完整准确的全局状态信息。
5.状态估计广泛用于电力系统，相关理论也已有系统的阐述。相比之下，状态估计在供热系统和燃气系统中的应用相对空白。现有集成电、热、气三个能源子系统的ies状态估计，没有考虑量测在非高斯噪声下的情况，无法在非高斯噪声下获得更高精度的全局运行状态。同时，现有研究对耦合元件的处理也未将其既纳入严格的等式约束又作为大权重的虚拟量测方程，以便既满足严格的边界等式约束条件又提高系统量测冗余度。


技术实现要素：

6.本发明为了解决现有技术之不足而提供一种针对电气热综合能源系统，首先厘清各子系统的量测量与状态量，并在此基础上建立非高斯测量噪声下的量测-状态方程，然后，考虑耦合元件边界条件，将其既作为等式约束又作为大权重的虚拟量测方程，基于极大似然鲁棒状态估计方法，提出了一种考虑边界等式约束的电气热综合能源系统鲁棒状态估计方法。
7.本发明为了解决上述技术问题而采用的技术解决方案如下：
8.考虑边界等式约束的电气热综合能源系统鲁棒状态估计方法，其步骤如下：
9.a、建立电气热综合能源系统量测-状态方程
10.a1、电力子系统量测模型
11.作为综合能源系统的子系统，电力系统中的量测量包含各节点注入功率、电压模值和支路功率，用向量形式可表示为：
12.ze＝[p q v p
b qb]
t
；状态量为各电气节点电压模值和相角，用向量形式可表示
为：xe＝[v θ]
t
，
[0013]
量测量ze与状态量xe的方程可表示如下：
[0014][0015]
式中：pi、qi表示节点i注入有功、无功；vi、vk表示节点i、k电压模值；δ
ik
表示节点i、k间相角差；g
ik
、b
ik
表示节点电导、电纳；ne表示电气总节点数；p
ik
、q
ik
表示节点i、k间支路有功、无功；
[0016]
a2、天然气子系统量测模型
[0017]
在气网稳态模型中，类似电路原理中的基尔霍夫电流定律、基尔霍夫电压定律和欧姆定律，有流量连续方程、压力回路方程和阻抗特性方程，其表达式如下所示：
[0018][0019]
式中：f
ij
为支路i、j间的天然气流量；l
gi
为节点i的天然气注入量；σ
gi
为与气节点i相连的管道集合；θg为任一气网回路的管道集合；pi和pj为气节点i、j的压力；π
ij
为气节点i、j压力的平方差；k
ij
为与管道有关的系数；k为与气体输送压力有关的系数；
[0020]
依据以上气网稳态方程，在气网中进行量测布置时，量测量包含气节点压力平方、管道支路天然气流量和各节点天然气注入量，用向量形式可表示为：zg＝[π f l]
t
；状态量为各气节点压力平方，用向量形式可表示为：xg＝[π]
t
，
[0021]
量测量zg与状态量xg的方程可表示为：
[0022][0023]
式中：s
ij
为天然气流向标识，当pi大于pj时，s
ij
＝+1；当pi小于pj时，s
ij
＝-1；
[0024]
a3、热力子系统量测模型
[0025]
在热力系统中，热力系统稳态模型包含水力模型和热力模型，
[0026]
在水力模型中，节点注水量与管道水流量的方程可表示为：
[0027][0028]
式中：li为节点i注水量；m
ij
为管道i、j的水流量；pipei为节点i相关联的管道集合；hi、hj为管道i、j节点压强；k
ij
为与管道和液体有关的参数，
[0029]
在热力模型中，热负荷方程表达式为：
[0030]
φi＝c
p
li(t
si-t
ri
)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(5)
[0031]
式中：φi为节点i热负荷；t
si
、t
ri
为供应温度和返回温度；c
p
为水比热容，
[0032]
依据以上水力和热力模型，在热力系统中进行量测布置时，量测量包含节点压强、节点注水量、管道流量、热负荷、供应温度和返回温度，用向量形式可表示为：zh＝[h l m φ t
s ti]
t
；状态量为各节点压强、供应温度和返回温度，用向量形式可表示为：xh＝[h t
s tr]
t
，
[0033]
量测量zh与状态量xh的方程可表示为：
[0034][0035]
a4、非高斯噪声下电气热综合能源系统量测模型
[0036]
由上可知，在电气热ies中，量测量为：z＝[z
e z
g zh]
t
，状态量为：x＝[x
e x
g xh]
t
，
[0037]
在非高斯噪声下，量测模型可表示为：
[0038]
z＝h(x)+v
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(7)
[0039]
式中：z、x分别为ies量测量、状态量；h(x)为量测方程；v为非高斯噪声，
[0040]
非高斯噪声v由高斯噪声和未知分布噪声混合而成，表达式如下所示：
[0041]
v＝(1-ε)φ(x)+εψ(x)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(8)
[0042]
式中：φ(x)表示高斯分布噪声；ψ(x)表示未知分布噪声，如大方差的高斯或拉普拉斯分布噪声；ε表示噪声污染调节系数；
[0043]
b、耦合元件模型
[0044]
在ies中，耦合元件包含热电联产机组、热泵和电锅炉，各耦合元件将各系统紧密连接在一起，在本方法中，考虑ies中最常见的耦合元件，热电联产机组(combined heat and power，chp)，
[0045]
chp的能量转换关系可用下式描述：
[0046]
[0047]
式中：φ
chp
为chp的产热功率；p
chp
为chp的产电功率；cm为chp产热产电功率比；f
chp
为天然气消耗量；c
gas
为天然气热值；η为chp转换效率，
[0048]
在ies状态估计模型中，耦合元件的边界条件可用下式表示：
[0049][0050]
式中：下标i表示第i个耦合元件，
[0051]
本方法将式(10)的耦合元件边界条件既作为边界等式约束，又作为大权重的虚拟量测方程，以便既满足严格的边界等式约束条件又提高系统量测冗余度，提升状态估计结果精度；
[0052]
c、综合能源系统鲁棒状态估计模型
[0053]
c1、极大似然鲁棒状态估计模型
[0054]
极大似然鲁棒状态估计模型如下：
[0055][0056]
式中：表示权重，d＝1.5，psi表示投影统计参数；ρ(rsi)为huber函数，
[0057]
ρ(rsi)函数表达式如下：
[0058][0059]
式中：为电力量测的标准化残差；s为鲁棒估计尺度参数；c为取值1至3的系数；
[0060]
c2、ies鲁棒状态估计模型
[0061]
本方法考虑量测的非高斯噪声，同时将耦合元件边界条件既作为等式约束又作为大权重的虚拟量测方程，建立ies联合鲁棒状态估计模型，如下所示：
[0062][0063]
式中：wi表示ies量测权重，r
si
为ies量测的标准化残差；c(x)＝0表示耦合元件等式约束；
[0064]
c3、ies鲁棒状态估计算法
[0065]
基于lagrange乘子法，消去等式约束，可得：
[0066][0067]
式中：λ表示lagrange乘子，
[0068]
最小化上式，得其一阶偏导表达式为：
[0069][0070]
式中：ai表示雅可比矩阵的第i行元素；c(x)表示等式约束c(x)对应的雅可比矩阵，
[0071]
将上式改写为如下的矩阵形式：
[0072][0073]
式中：q＝diag(q(r
si
))为q(r
si
)表示的对角阵，)表示的对角阵，为表示的对角阵，
[0074]
基于迭代重新加权，求解以上一阶偏导表达式，可得如下修正量表达式：
[0075][0076]
状态量x的迭代更新表达式为：
[0077]
x
k+1
＝xk+δxkꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(18)
[0078]
基于以上电气热ies鲁棒状态估计方法，可将电、气和热系统进行联合统一估计，增加系统量测冗余，获得全局一致解，其计算步骤如下：
[0079]
(1)输入基础数据，初始化电、气和热子系统状态量x，设置最大迭代次数kmax、收敛阈值ε，置初始迭代k＝0；
[0080]
(2)计算矩阵h、q和r；
[0081]
(3)计算状态量增量δxk；
[0082]
(4)若max|δxk|＜ε，则迭代收敛，计算结束；反之，转步骤(5)；
[0083]
(5)x
k+1
＝xk+δxk，k＝k+1，转步骤(2)。
[0084]
本发明采用上述技术解决方案所能达到的有益效果是：
[0085]
1、详细厘清了电、气和热各子系统的量测量与状态量，并在此基础上建立了非高斯测量噪声下的量测-状态方程；
[0086]
2、考虑耦合元件边界条件，将其既作为等式约束又作为大权重的虚拟量测方程，可保证所得结果严格满足边界等式约束条件，同时增加系统量测冗余度，获得全局状态一致解，提高了估计结果精度。
附图说明
[0087]
图1为考虑边界等式约束的ies鲁棒状态估计流程图；
[0088]
图2为电气热综合能源系统网络结构图。
具体实施方式
[0089]
考虑边界等式约束的电气热综合能源系统鲁棒状态估计方法，其步骤如下：
[0090]
a、建立电气热综合能源系统量测-状态方程
[0091]
a1、电力子系统量测模型
[0092]
作为综合能源系统的子系统，电力系统中的量测量包含各节点注入功率、电压模值和支路功率，用向量形式可表示为：
[0093]
ze＝[p q v p
b qb]
t
；状态量为各电气节点电压模值和相角，用向量形式可表示为：xe＝[v θ]
t
，
[0094]
量测量ze与状态量xe的方程可表示如下：
[0095][0096]
式中：pi、qi表示节点i注入有功、无功；vi、vk表示节点i、k电压模值；δ
ik
表示节点i、k间相角差；g
ik
、b
ik
表示节点电导、电纳；ne表示电气总节点数；p
ik
、q
ik
表示节点i、k间支路有功、无功；
[0097]
a2、天然气子系统量测模型
[0098]
在气网稳态模型中，类似电路原理中的基尔霍夫电流定律、基尔霍夫电压定律和欧姆定律，有流量连续方程、压力回路方程和阻抗特性方程，其表达式如下所示：
[0099][0100]
式中：f
ij
为支路i、j间的天然气流量；l
gi
为节点i的天然气注入量；σ
gi
为与气节点i相连的管道集合；θg为任一气网回路的管道集合；pi和pj为气节点i、j的压力；π
ij
为气节点i、j压力的平方差；k
ij
为与管道有关的系数；k为与气体输送压力有关的系数；
[0101]
依据以上气网稳态方程，在气网中进行量测布置时，量测量包含气节点压力平方、管道支路天然气流量和各节点天然气注入量，用向量形式可表示为zg＝[π f l]
t
；状态量为各气节点压力平方，用向量形式可表示为xg＝[π]
t
，
[0102]
量测量zg与状态量xg的方程可表示为：
[0103][0104]
式中：s
ij
为天然气流向标识，当pi大于pj时，s
ij
＝+1；当pi小于pj时，s
ij
＝-1；
[0105]
a3、热力子系统量测模型
[0106]
在热力系统中，热力系统稳态模型包含水力模型和热力模型，
[0107]
在水力模型中，节点注水量与管道水流量的方程可表示为：
[0108][0109]
式中：li为节点i注水量；m
ij
为管道i、j的水流量；pipei为节点i相关联的管道集合；hi、hj为管道i、j节点压强；k
ij
为与管道和液体有关的参数，
[0110]
在热力模型中，热负荷方程表达式为：
[0111]
φi＝c
p
li(t
si-t
ri
)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(5)
[0112]
式中：φi为节点i热负荷；t
si
、t
ri
为供应温度和返回温度；c
p
为水比热容，
[0113]
依据以上水力和热力模型，在热力系统中进行量测布置时，量测量包含节点压强、节点注水量、管道流量、热负荷、供应温度和返回温度，用向量形式可表示为zh＝[h l m φ t
s ti]
t
；状态量为各节点压强、供应温度和返回温度，用向量形式可表示为xh＝[h t
s tr]
t
，
[0114]
量测量zh与状态量xh的方程可表示为：
[0115][0116]
a4、非高斯噪声下电气热综合能源系统量测模型
[0117]
由上可知，在电气热ies中，量测量为z＝[z
e z
g zh]
t
，状态量为x＝[x
e x
g xh]
t
，
[0118]
在非高斯噪声下，量测模型可表示为：
[0119]
z＝h(x)+v
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(7)
[0120]
式中：z、x分别为ies量测量、状态量；h(x)为量测方程；v为非高斯噪声，
[0121]
非高斯噪声v由高斯噪声和未知分布噪声混合而成，表达式如下所示：
[0122]
v＝(1-ε)φ(x)+εψ(x)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(8)
[0123]
式中：φ(x)表示高斯分布噪声；ψ(x)表示未知分布噪声，如大方差的高斯或拉普拉斯分布噪声；ε表示噪声污染调节系数；
[0124]
b、耦合元件模型
[0125]
在ies中，耦合元件包含热电联产机组、热泵和电锅炉，各耦合元件将各系统紧密连接在一起，在方法中，考虑ies中最常见的耦合元件，热电联产机组(combined heat and power，chp)，
[0126]
chp的能量转换关系可用下式描述：
[0127][0128]
式中：φ
chp
为chp的产热功率；p
chp
为chp的产电功率；cm为chp产热产电功率比；f
chp
为天然气消耗量；c
gas
为天然气热值；η为chp转换效率，
[0129]
在ies状态估计模型中，耦合元件的边界条件可用下式表示：
[0130][0131]
式中：下标i表示第i个耦合元件，
[0132]
本方法将式(10)的耦合元件边界条件既作为边界等式约束，又作为大权重的虚拟量测方程，以便既满足严格的边界等式约束条件又提高系统量测冗余度，提升状态估计结果精度；
[0133]
c、综合能源系统鲁棒状态估计模型
[0134]
c1、极大似然鲁棒状态估计模型
[0135]
极大似然鲁棒状态估计模型如下：
[0136][0137]
式中：表示权重，d＝1.5，psi表示投影统计参数；ρ(rsi)为huber函数,
[0138]
ρ(rsi)函数表达式如下：
[0139][0140]
式中：为电力量测的标准化残差；s为鲁棒估计尺度参数；c为取值1至3的系数；
[0141]
c2、ies鲁棒状态估计模型
[0142]
本方法考虑量测的非高斯噪声，同时将耦合元件边界条件既作为等式约束又作为大权重的虚拟量测方程，建立ies联合鲁棒状态估计模型，如下所示：
[0143][0144]
式中：wi表示ies量测权重，r
si
为ies量测的标准化残差；c(x)＝0表示耦合元件等式约束；
[0145]
c3、ies鲁棒状态估计算法
[0146]
基于lagrange乘子法，消去等式约束，可得：
[0147][0148]
式中：λ表示lagrange乘子，
[0149]
最小化上式，得其一阶偏导表达式为：
[0150][0151]
式中：ai表示雅可比矩阵的第i行元素；c(x)表示等式约束c(x)对应的雅可比矩阵，
[0152]
将上式改写为如下的矩阵形式：
[0153][0154]
式中：q＝diag(q(r
si
))为q(r
si
)表示的对角阵，)表示的对角阵，为表示的对角阵，
[0155]
基于迭代重新加权，求解以上一阶偏导表达式，可得如下修正量表达式：
[0156][0157]
状态量x的迭代更新表达式为：
[0158]
x
k+1
＝xk+δxkꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(18)
[0159]
基于以上电气热ies鲁棒状态估计方法，可将电、气和热系统进行联合统一估计，增加系统量测冗余，获得全局一致解，其计算步骤如下：
[0160]
(1)输入基础数据，初始化电、气和热子系统状态量x，设置最大迭代次数kmax、收敛阈值ε，置初始迭代k＝0；
[0161]
(2)计算矩阵h、q和r；
[0162]
(3)计算状态量增量δxk；
[0163]
(4)若max|δxk|＜ε，则迭代收敛，计算结束；反之，转步骤(5)；
[0164]
(5)x
k+1
＝xk+δxk，k＝k+1，转步骤(2)。
[0165]
计算流程由图1所示。
[0166]
算例仿真与分析
[0167]
1、基础数据与仿真条件
[0168]
由图2所示的电气热综合能源系统网络结构，对本方法所提考虑边界等式约束的电气热综合能源系统鲁棒状态估计方法的准确性、鲁棒性和有效性进行验证，
[0169]
图2中，热电联产机组chp有两个，分别为chp1和chp2，起到将电、气和热三个系统耦合的作用，node
ei
、node
gi
和node
hi
分别表示电、气和热子系统中的节点，
[0170]
在该系统中，量测为完全测量，基于多能流计算的结果，并将其作为真值，基于蒙特卡罗模拟，将量测误差设为高斯和非高斯噪声，在真值基础上叠加噪声来模拟实际测量，蒙特卡洛模拟次数为2
×
103次。
[0171]
以状态量估计值相对真值的总体偏差erms
x
作为估计结果精度的评价指标，计算表达式如下：
[0172][0173]
式中：erms
x
表示状态量估计值相对真值的总体偏差；n表示网络节点总数；表示第i个状态量估计值；x
i,ture
表示第i个状态量真值；t表示蒙特卡罗模拟次数，
[0174]
此外，考虑耦合元件的边界等式约束条件，以s1和s2表示耦合元件边界的匹配指标，其计算表达式如下：
[0175][0176][0177]
式中：n
chp
表示chp元件个数；ρ1和ρ2表示单次模拟的耦合元件边界匹配指标，s1和s2值越小，表示边界匹配结果越好，越符合实际情况。
[0178]
在高斯噪声和非高斯噪声下，将基于wls的电、气和热各子系统单独状态估计结果与本方法的估计结果进行对比。
[0179]
2、仿真结果与分析
[0180]
(1)高斯噪声下状态估计结果比较
[0181]
在真值基础上，叠加方差为10-3
pu的高斯分布噪声，电、气和热各子系统基于wls的单独状态估计结果与本方法所提ies鲁棒联合状态估计结果对比，如表1、表2和表3所示：
[0182]
表1高斯噪声下电网状态估计结果
[0183][0184]
表2高斯噪声下气网状态估计结果
[0185][0186]
表3高斯噪声下热网状态估计结果
[0187][0188]
由表1、表2和表3可知，在蒙特卡洛模拟次数为2
×
103次情况下，基于wls的单独电网估计，ermsv和erms
θ
分别为2.31
×
10-3
和4.79
×
10-3
，而本方法分别为3.27
×
10-4
和1.05
×
10-4
；基于wls的单独气网估计，erms
π
为5.79
×
10-3
，而本方法为1.33
×
10-3
；基于wls的单独热网估计，ermsm、erms
ts
和erms
tr
分别为5.12
×
10-3
、4.35
×
10-3
和3.17
×
10-3
，而本方法分别为2.57
×
10-3
、2.98
×
10-3
和2.25
×
10-3
。
[0189]
相比于基于wls的单独状态估计，本方法所提的状态量估计值相对真值的总体偏差erms
x
均更小，故而本方法将各系统联合估计，有效增大了各个系统的量测冗余度，获取全局状态一致解，状态估计结果精度更高。
[0190]
(2)非高斯噪声下状态估计结果比较
[0191]
将噪声污染系数设为ε＝0.15，φ(x)方差设为10-3
pu，ψ(x)取为高斯分布，方差为φ(x)的100倍，电、气和热各子系统基于wls的单独状态估计结果与本方法所提ies鲁棒联合状态估计结果对比，如表4、表5和表6所示：
[0192]
表4非高斯噪声下电网状态估计结果
[0193][0194]
表5非高斯噪声下气网状态估计结果
[0195][0196]
表6非高斯噪声下热网状态估计结果
[0197]
[0198]
由表4、表5和表6可知，在蒙特卡洛模拟次数为2
×
103次情况下，基于wls的单独电网估计，ermsv和erms
θ
分别为8.73
×
10-2
和9.54
×
10-2
，而本方法分别为6.21
×
10-4
和7.15
×
10-4
；基于wls的单独气网估计，erms
π
为4.37
×
10-2
，而本方法为3.23
×
10-3
；基于wls的单独热网估计，ermsm、erms
ts
和erms
tr
分别为4.17
×
10-2
、4.25
×
10-2
和2.89
×
10-2
，而本方法分别为3.54
×
10-3
、4.21
×
10-3
和4.52
×
10-3
。
[0199]
相比于高斯噪声，在非高斯噪声下基于wls的单独状态估计结果erms
x
指标均显著增大，而本方法erms
x
指标变化较小，且均显著小于基于wls的单独状态估计结果，说明本方法状态估计结果精度更高，鲁棒性好，不易受未知噪声分布的影响。
[0200]
(3)边界条件匹配情况比较
[0201]
在高斯噪声情况下和非高斯噪声情况下，基于wls的单独状态估计和本方法的耦合元件边界的匹配指标如表7和表8所示：
[0202]
表7高斯噪声下边界匹配指标
[0203][0204][0205]
表8非高斯噪声下边界匹配指标
[0206][0207]
由表7和表8可知，在高斯噪声下，基于wls的单独估计，耦合元件边界匹配指标s1和s2分别为3.72
×
10-2
和4.83
×
10-2
，单次模拟的最大边界匹配指标分别max(ρ1)和max(ρ2)分别为0.1375和0.1528；而本方法耦合元件边界匹配指标s1和s2分别为4.78
×
10-7
和6.92
×
10-7
，单次模拟的最大边界匹配指标分别max(ρ1)和max(ρ2)分别为5.93
×
10-6
和7.85
×
10-6
。在非高斯噪声下，基于wls的单独估计，耦合元件边界匹配指标s1和s2分别为8.77
×
10-2
和9.58
×
10-2
，单次模拟的最大边界匹配指标分别max(ρ1)和max(ρ2)分别为0.2257和0.2597；而本方法耦合元件边界匹配指标s1和s2分别为5.15
×
10-6
和3.92
×
10-6
，单次模拟的最大边界匹配指标分别max(ρ1)和max(ρ2)分别为1.23
×
10-5
和2.47
×
10-5
。