基于约束卡尔曼滤波的区域多能源系统鲁棒状态估计

廖英祺，荆江平，叶婷

(1.国网江苏省电力有限公司南京供电分公司，江苏 南京 210019；2.国网江苏省电力有限公司，江苏 南京 210024)

区域多能源系统(regional multi-energy system，RMES)是集电能、天然气等异质能源的产生、传输、分配与消费于一体的新型能源系统，能够有效提升用户侧终端综合能效和风光发电消纳水平，是清洁低碳、安全高效的能源体系的重要组成部分[1-3]。为了保障RMES的安全可靠运行、提高用户的用能质量，以及适应主动管理的需求，调度管理人员需实时感知RMES的运行状态。作为态势感知的基础，状态估计(state estimation，SE)能够对量测数据进行分析和计算，为RMES提供更加精确的状态信息。

当下配电网估计领域的学术研究已较为广泛[4-6]，众多学者在此基础上展开了一系列面向电-气多能源系统(electric-gas multi-energy system，EGMES)的研究：文献[7]考虑复杂天然气管网的稳态建模及坏数据辨识，提出了稳态EGMES-SE方法；针对EGMES-SE中存在的配气网初值和求解效率问题，文献[8-9]提出了抗差及双线性抗差EGMES-SE方法；文献[10]计及实时量测误差的不确定性，提出一种区间EGMES-SE方法；文献[11]考虑EGMES的多管理主体特征和通信隐私问题，提出基于交替方向乘子法的分布式估计求解方法。上述文献中配气网均根据Weymouth方程进行稳态建模，忽略了管道内天然气的流量和压强随时间和空间发生变化的特性，会产生较大的估计误差。而配气管道相比于输气管道储气能力较弱，在分钟级至小时级时间尺度下，配气管道的动态特征更为明显，因此亟需构建精确的天然气动态模型[12]。

天然气动态模型特性遵循流体力学定律，通常由质量守恒与动量守恒偏微分方程组描述，难以用解析法直接求解，目前较为主流的处理方法是通过有限差分，将连续的偏微分方程转化为若干组离散的代数方程来降低求解难度[13-14]。文献[15]考虑电、气、热子系统间的耦合约束，提出了满足多场景需求的顺序估计算法；文献[16]将扩展卡尔曼滤波(Kalman filter，KF)算法引入EGMES-SE中，并通过线性插值法生成伪量测来解决电、气多时间尺度问题；文献[17]考虑多源量测设备采样周期的延迟问题，提出了计及时序运行特性的EGMES-SE算法；文献[18]提出异步分布式EGMES-SE策略，根据配气网静态安全域自主调整估计周期。以上研究大大促进了动态EGMES-SE的发展，但仍存在一定的不足：①现有的动态估计方法均不具备抗差特性；②配气网有限差分模型会引入大量的时、空微元，显著增加了状态变量的数量，在有限的量测配置下，配气网量测冗余度过低，导致估计精度下降；③研究对象多为跨区级系统，以大型输电、气系统作为骨干网架，缺乏对接入间歇性可再生能源的RMES场景的考虑。

针对上述问题，本研究提出基于约束卡尔曼的鲁棒RMES-SE方法。首先，采用Lax-Wendroff方法对配气网动态过程进行数值求解，推导出配气网状态空间模型；其次，构建配电网状态空间模型，通过耦合约束搭建电、气子系统之间能量转换的桥梁，为RMES-SE提供模型支撑；最后，以KF算法为基础，添加量测噪声自适应算法，增强算法的抗差性能，并计及配气网时序状态约束及电-气边界耦合约束修正估计结果，大幅提升配气网量测冗余度和估计精度。仿真算例根据改进的IEEE 33节点配电网与10节点配气网耦合而成，考虑分布式电源接入的场景，由此验证所提方法的有效性。

1 配气网状态空间建模

相较于配电网而言，配气网具有更长的暂态过程，主要是因为天然气的传输速度相对较慢，导致配气网的响应速度较慢。天然气的流量受到管道首末端压力差的驱动，并且受到管道长度、内径、粗糙度、传输路径高度和边界条件等多方面因素影响；因此，正确地描述配气网管道内气体流量和压强随时间和空间的变化关系，是配气网状态空间建模的关键[12]。

1.1 配气网管道差分建模

在管道水平条件下，不考虑重力加速度和倾角的影响，管道内部的暂态变化过程，可用沿天然气管道轴线的一维等温传输动态方程来描述，即根据质量及动量守恒定律得到的一组偏微分方程，具体为[16]：

(1)

(2)

式(1)、(2)中：ta、x分别为时间、管道距离；p、q分别为天然气的气体压强、质量流量；va为气体声速；S、D和λ分别为管道的横截面积、内径和摩擦系数；w为管道平均流速。

应用偏微分方程形式的配气网管道方程进行状态估计，会大大增加计算复杂度；因此，本研究对配气网管道的暂态方程采用Lax-Wendroff有限元差分格式建模，差分后的配气网管道方程能够表征相邻时刻间流量和压强的变化关系，具体可将式(1)、(2)转化为[16]：

pi+1，t+1+pi，t+1-pi+1，t-pi，t+

(3)

(4)

式(3)、(4)中：i为管道差分节点(包含真实节点和虚拟节点)编号；t为时刻(时间断面)，用下标t表示时刻t的相关变量，下同；Δx、Δt分别为差分空间步长、时间步长。如图1所示，配气网管道差分建模引入了虚拟节点、支路(将第j条真实支路差分成M段)，同时将连续的偏微分方程转化为若干组离散的代数方程，以降低求解难度，便于精确获取短时间尺度下配气网波动的状态参数。

图1 配气网管道有限元差分建模Fig.1 Finite difference model of gas network pipeline

进一步将式(3)、(4)转化为统一的矩阵形式，即

(5)

(6)

式(5)、(6)中A1—A5均为常系数矩阵的元素。

1.2 配气网状态方程推导

假设配气网中的真实节点数为N，真实支路数量为b，虚拟支路数量为v(即所有真实支路差分段数之和)，则时刻t的配气网状态列向量

xg，t=(pt，qt)T=(pn，t，pB，t，qt)T∈RN+2v.

(7)

其中：

(8)

(9)

式(7)—(9)中：pt=(pn，t，pB，t)，pn，t、pB，t分别为真实节点压强行向量、虚拟节点压强行向量；qt为虚拟支路质量流量行向量；pBj，t、qBj，t分别为第j条真实支路的虚拟节点压强行向量、虚拟支路质量流量行向量；qBj，0，t、qBj，M，t可等价为第j条真实支路首、末端质量流量。

每段虚拟支路的配气网暂态特性皆可由式(5)描述，因此xg，t+1与xg，t间的关系式可表示为

J1xg，t+1=J2xg，t，

(10)

式中J1，J2∈R2v×(N+2v)为常系数矩阵，具体的矩阵元素可根据式(5)、(6)计算得到。

进一步考虑配气网节点平衡关系，即任意时刻t下流入、流出任意节点的质量流量之和等于节点注入质量流量，根据拓扑约束关系可构建方程

(11)

联合式(10)、(11)，可得

(12)

式中零矩阵02∈RN×(N+2v)，03∈R2v。

由式(12)可转化得到配气网的暂态状态方程

xg，t+1=Fgxg，t+Gg，t=

(13)

式中：Fg∈R(N+2v)×(N+2v)、Gg，t∈RN+2v分别为配气网的状态转移矩阵、控制变量矩阵。计算过程中：Fg为常系数矩阵，与原始管道参数和配气网动态特性有关；Gg，t与下一时刻各真实节点注入质量流量qn，t+1有关，由于配气网短期负荷波动相对平缓，可通过Holt两参数指数平滑法预测获得。

1.3 配气网状态空间建模

配气网量测向量zg，t+1主要包括真实节点压强pn，t+1，真实支路首、末端质量流量q0，t+1、qM，t+1，以及真实节点注入质量流量qn，t+1，可构建量测方程

(14)

式中：zg，t+1∈R2N+2b；Hg∈R(2N+2b)×(N+2v)、J4∈R(N+2b)×(N+2v)为常系数矩阵，其中J4各行向量中与量测量相对应的状态量位置元素为1，其余元素为0。

由式(13)、(14)可构建配气网状态空间模型

(15)

式中wg，t、vg，t+1分别为配气网过程向量、量测噪声向量，一般服从互不相关的零均值高斯分布。

2 RMES-SE建模

在配气网状态空间建模的基础上，本节进一步构建配电网状态空间模型，并以燃气轮机(gas turbine，GT)和电转气(power to gas，P2G)为接口，建立电-气双向耦合环节，为RMES-SE提供模型支撑。本研究采用的RMES-SE方法以KF算法为基础，通过自适应更新量测噪声统计特性，保证估计的鲁棒性，并计及配气网时序状态约束及电-气边界耦合约束修正估计结果，在提高量测冗余度的同时，确保电-气边界条件的一致性，实现对RMES运行状态的实时感知。

2.1 配电网状态空间建模

本研究建立的配电网状态空间模型为准稳态模型[19]，主要面向配电网负荷波动较小、稳态运行的场景。配电网的状态方程通过时间序列预测方法建立，即采用Holt两参数指数平滑法；而量测方程则与常规配电网稳态估计相同。配电网状态量选取节点α的电压幅值Uα和电压相角θα，量测量包括节点电压幅值Uα，节点注入有功功率Pα、无功功率Qα，以及节点α、β之间支路有功功率Pαβ、无功功率Qαβ。构建配电网状态空间模型为

(16)

式中：xe，t+1、ze，t+1分别为配电网状态向量、量测向量；Fe，t、Ge，t为配电网的状态转移矩阵、控制变量矩阵，可通过Holt两参数指数平滑法求解获得；h(·)为配电网量测方程；we，t、ve，t+1分别为配电网过程向量、量测噪声向量，误差分布原理与上述配气网噪声相同。

2.2 耦合元件建模

在RMES中，天然气能与电能一般通过GT和P2G实现双向传递。GT通过燃烧天然气产生电能，具有启停灵活、能量梯级利用的特点；而P2G通过电解水制备氢气，并以氢气为原料制造甲烷后，注入到现有配气网进行规模化存储和传输，实现节能减排[20-21]。本研究中GT和P2G输入/生成气流量和输出/消耗电功率的函数关系为[20-21]

(17)

式中：PGT、PP2G分别为GT、P2G输出/消耗的电功率；qGT、qP2G分别为GT、P2G燃烧/生成的气流量；GT、P2G分别为GT、P2G的能量转换效率，非热电联产形式下GT效率可达40%左右，P2G完整化学反应综合能量转换效率约为45%～60%；LHV为天然气低热值[21]，标准状况下其值为35.40～39.12 MJ/m3，采用其平均值37.26 MJ/m3，标准状况下天然气密度0.717 4 kg/m3，经转换后LHV=51.94 MJ/kg；ηGT、ηP2G分别为等价后的GT、P2G能量转换系数。

2.3 鲁棒RMES-SE

受现场工况、通信干扰等因素影响，RMES量测量中不可避免地出现非高斯噪声和不良数据，若量测噪声统计特性与真实噪声情况不匹配，将导致常规KF算法无法准确跟踪系统状态。本文所提方法在KF算法的基础上，运用Huber的鲁棒M估计理论[22]，动态修正量测噪声协方差阵，保证量测噪声统计特性符合实际工况。在分别获得电、气子系统本地估计值后，考虑电-气边界耦合约束及配气网时序状态约束，采用带线性等式约束的KF[23]进行融合估计，修正得到全局状态一致解。本文所提方法主要分为以下4个步骤。

a)预测步。

(18)

式中：下标k∈{e，g}，表示电、气子系统之一，当k=“g”时，Fg，t=Fg，t+1=Fg；xk，f，t+1、Pk，f，t+1分别为时刻t+1配电/气网状态预测向量、预测协方差矩阵；xk，p，t、Pk，p，t分别为时刻t配电/气网状态估计向量、估计协方差矩阵；Qk，t为过程噪声协方差矩阵。

b)量测噪声更新[22]。

(19)

(20)

式(19)、(20)中：Rk，t+1、Rk，c，t+1分别为修正前、后的量测噪声方差矩阵；φk，t+1为量测噪声尺度因子，φk，s，t+1为φk，t+1的第s个对角元素；c为残差阈值，通常取值为1.3～2.0；rk，s，t+1为归一化量测预测残差向量rk，t+1的第s个元素。rk，t+1计算公式为：

rk，t+1=diag(zk，t+1-zk，f，t+1)[diag(Pk，f，z，t+1)]-1/2，

(21)

(22)

式(21)、(22)中：zk，f，t+1为量测预测向量，其中ze，f，t+1=h(xe，f，t+1)，zg，f，t+1=Hgxg，f，t+1；Pk，f，z，t+1为量测预测协方差矩阵；He，t+1=∂h(xe，f，t+1)/∂xe，f，t+1为配电网雅可比矩阵。

c)滤波步。

(23)

式中Kk，t+1为卡尔曼增益矩阵。

d)约束KF状态修正。

在配电、气网独立进行鲁棒估计的基础上，将配气网时序状态约束和电-气边界耦合约束添加到RMES-SE中，在协调中心侧采用带线性等式约束的KF对全局状态进行融合估计[23]。根据估计投影法，融合估计可转化为求解下列问题：

(24)

(25)

式(24)、(25)中：xc，t+1、xp，t+1和Pp，t+1分别为电-气融合修正状态向量、估计状态向量和估计协方差矩阵；hc(·)为配电网耦合量测方程；η为RMES能量转换系数矩阵；Hgc为配气网耦合量测系数矩阵；04∈R(2m-1)×(N+2v)，05∈R(N+2v)×(2m-1)，其中m为配电网节点数量。式(24)的解为

xc，t+1=xp，t+1-

(26)

(27)

(28)

式(26)—(28)中：Ht+1、D(·)分别为电-气融合系数矩阵、等式约束向量；06∈R2v×(2k-1)；Hec，t+1为配电网耦合量测对应的雅可比矩阵。

3 算例分析

为验证本文所提方法的性能效果，采用改进的IEEE 33节点配电网和10节点配气网[24]构建RMES算例作为本研究的测试系统，共包含32条配电线路、11条配气管道，并通过1台GT和1台P2G实现电气耦合，系统结构如图2所示，其中：配电网中2台光伏发电装置(PV)经电压源换流器并网，分别接于节点18、31，均采用有功功率-无功功率控制策略[24]，即控制交流侧注入总有功、无功功率为定值；2台风电机组(WG)采用直接并网方式，分别接于节点9、21[25]。

图2 RMES测试系统Fig.2 Test system of RMES

仿真实验中：配电网的基准电压为12.66 kV，基准容量为1 MV·A；配气网的基准气压为1 bar，基准流量为1 kg/s，配气网管道摩擦系数设为0.004，气体声速设为400 m/s；配气网差分空间步长、时间步长分别设置为500 m、15 min。在潮流真值的基础上添加服从高斯分布的随机噪声获得量测值，相较潮流真值，各量测值最大误差分别为：1.5%(电压幅值实时量测)、3%(支路功率、负荷功率实时量测)、30%(根据历史数据、实时天气预测的负荷功率和直接并网分布式电源出力伪量测)、3%(节点压强实时量测)、4.5%(支路首末端质量流量、节点注入质量流量实时量测)。量测配置信息见表1。仿真环境的建立采用Intel Core i7-10700 CPU和16 GiB RAM的计算机，通过MATLAB 2016b进行求解。

表1 量测配置信息Tab.1 Measurement configuration

3.1 滤波效果测试

根据RMES负荷曲线和分布式电源出力曲线，以15 min为采样间隔时长，进行连续96个时刻的测试，在无坏数据干扰情况下使用KF(无量测噪声更新和状态修正步骤)和本文所提的鲁棒约束卡尔曼滤波(robust constrained Kalman filter，RCKF)算法进行滤波效果测试。采用状态量的最大估计误差Emax和平均估计误差Eavg衡量估计精度：

(29)

(30)

式(29)、(30)中：T为时刻总数；n为状态变量数；xt为状态向量真实值，xh，t为其第h个元素；xc，h，t为状态向量融合修正值的第h个元素。表2展示了不同方法的状态量估计精度对比结果。

在某一时刻分别采用KF与RCKF测试，使用蒙特卡洛模拟法进行1 000次抽样，计算电-气状态量的平均估计值，算法精度对比结果如图3所示。结合表2和图3可知：相较KF，RCKF对配气网估计精度有明显的提升作用；而配电网估计精度虽有所下降，但基本处于同一数量级，对整体滤波效果影响不大。

图3 RMES各节点状态估计精度对比Fig.3 Comparison of state estimation accuracy for each node in RMES

为进一步比较KF、RCKF对RMES局部滤波效果的影响，选取算例中与GT相关联的电气耦合量测，即配电网节点7注入有功功率和配气网节点7注入质量流量，并采用相对估计误差来衡量量测估计精度，各时刻测试结果如图4所示。

图4 RMES耦合量测估计精度对比Fig.4 Comparison of coupled measurement estimation accuracy in RMES

配气网估计精度提升明显的原因主要是：差分模型下配气网状态量大幅增加，量测数目有限条件下，量测冗余度较低；而本文所提的RCKF通过添加配气网时序状态约束和电-气边界耦合约束变相提升量测冗余，能够在KF基础上修正得到更为精确的滤波结果。量测冗余度对比见表3。

配电网估计精度下降的原因是：一方面配电网量测冗余改变不足以明显提升估计精度；另一方面RCKF添加了量测噪声更新算法，当部分正常量测数据时序波动较大时，为保障状态估计的鲁棒性，会主动修正数据统计特性，降低其对估计结果的影响，导致估计误差略微增大。

3.2 抗差性能测试

为测试RCKF在坏数据干扰下的抗差性能，假设在第32—第40个采样时刻，RMES中配电网节点5-6支路首端有功功率和配气网节点2-4支路首端质量流量量测通信丢失(即量测值为0)，在此条件下进行动态估计。图5展示了RMES量测通信丢失下的估计结果。

图5 RMES量测丢失下的估计结果Fig.5 Estimation results under measurement loss in RMES

常规KF算法缺乏对量测噪声特性的自适应处理能力，当量测量统计特性偏离先验统计特性时，丢失量测对应的估计值大幅偏离真实值；而本文所提方法在此期间仍能较好跟踪系统状态变化。

除了较为严重的量测数据丢失情况外，当量测受到通信干扰或仪表故障影响，其误差也会远大于正常的量测噪声。为此在各时刻下，设坏数据由叠加10倍标准差的高斯白噪声随机生成，占子系统量测的2%～10%，针对不同比例的坏数据场景，分别进行1 000次蒙特卡洛仿真实验，并使用状态量的Eavg和Emax来评估本文所提方法抗差性能的普适性。表4展示了相应的仿真结果。

表4 不同坏数据比例下的估计结果Tab.4 Estimation results under different bad data ratios

由表4可知，随着坏数据占比的增加，RMES状态量的Eavg和Emax也会相应上升，但仍能够维持较为稳定的数据水平。相较于正常量测下的估计结果，数据差异不大，这主要是因为本文所提方法采用基于Huber的鲁棒M估计理论，能够根据量测新息计算时变多维的量测噪声尺度因子。当发现可疑的量测数据时，本文所提方法会修正噪声统计特性，进而调整滤波增益，降低可疑数据对状态估计的影响，即使在面对坏数据干扰时也能够保持鲁棒性能。

3.3 计算效率测试

本文所提的RMES-SE方法分为在线估计和离线计算2个部分：在线估计指接收到实时量测信息后进行配电网和配气网鲁棒估计，并通过约束修正达成全局一致解和精度提升；离线计算指RMES-SE可根据历史状态信息或负荷预测信息对下一时刻的状态量进行预测，预测步骤无需占用系统实时算力资源。表5列出了本文所提RCKF方法各步骤计算耗时。

表5 RMES-SE计算效率统计数据Tab.5 Statistical data of calculation efficiency for RMES-SE

由表5可知，本文所提方法整体在线计算时间仅需约3×10-3s，能够满足实时跟踪RMES动态变化的需求。相较常规KF方法，RCKF方法增加了量测噪声更新和约束修正步骤，能增强算法鲁棒性，提高电气量测冗余度，达成边界条件约束一致，且不会增加大时间成本，具有良好的工程应用前景。

4 结论

伴随着以天然气分布式能源与可再生能源为主的新型供能技术的飞速发展，RMES需求侧的能源多样化促使子系统间耦合性增强，连锁故障引发的风险系数也在不断提高；因此，需要实时感知RMES的运行状态，预知未来变化趋势，保障RMES安全稳定运行。鉴于此，本文提出了基于约束KF的鲁棒RMES-SE方法，得到以下结论：

a)本文所提方法采用KF算法进行各子系统独立估计，并基于Huber的鲁棒M估计理论实时修正量测噪声信息，原理简单，易于实现，有效保障了估计的鲁棒性。

b)本文所提方法利用耦合元件边界约束和配气网时序状态约束，协同修正各子系统估计值，获得全局一致解，可提升量测冗余度，具有较好的估计精度与计算效率。

本文采用管道平均流速对配气网偏微分方程进行线性化处理，适用于各管道流速在设定流速附近波动的场景；但面对负荷突变等场景，可能会产生较大估计误差。后续将在估计过程中添加兼顾计算效率与结果精度的流速自适应校正步骤。此外，电、气子系统存在的量测传输非同步性及耦合元件非线性建模等因素对RMES-SE的影响，也值得加入到今后的探索研究中。