多时间尺度电热综合能源系统状态估计研究

2021-10-28慈文斌顾海飞朱劲松

热力发电 2021年9期

慈文斌，顾海飞，朱劲松

（1.国网山东省电力公司，山东济南 250000；2.东南大学电气工程学院，江苏南京 210096；3.国网江苏省电力有限公司，江苏南京 210096）

日益增长的能源危机和环境问题引发了构建高效、清洁能源结构的新需求[1-2]。作为能源变革的一种主要形式，电热综合能源系统通过能源梯级利用和协调优化将电力系统和热力系统联系在一起，从而促进可再生能源消纳和运行灵活性[3]，并在工程实践中得到广泛应用[4-5]。因此，电热综合能源系统的统一建模和精确分析，对促进综合能源系统发展十分必要。

状态估计作为一种评估系统运行及预测系统发展的数值方法，为系统在线优化和控制提供理论指导。文献[6]首先利用最小二乘法构建了热力系统状态估计模型，但由于仅考虑用户侧数据，其估计精度较低。文献[7]通过构建基于全量测配置的电热综合能源系统状态模型，有效提高了估计精度，但该方法抗扰动能力较弱。文献[8]引入双线性抗差方法用于解决电热综合能源系统中的非线性问题。文献[9]考虑了多主体特性，实现了综合能源系统的分布式状态估计。文献[10-11]将电热联合状态估计与独立估计对比，并指出联合估计有利于提高量测冗余从而改善估计精度。

目前对于电热综合能源系统的技术应用已经开展了一定研究，如协调优化[12]、新能源消纳[13]等，但这些研究均忽略了系统的动态特性。考虑到不同的运行时间尺度，在一定控制指令下，电力系统能迅速达到稳定状态，而热力系统仍处于动态过程[14]。若采用稳态模型对热力系统进行建模，则忽略了热力状态分布的时空演变特性，将产生较大的估计误差[15]。多时间尺度特性使得传统的针对单一时间断面的状态估计方法不再适用于电热综合能源系统，由于热力系统的状态量在时序上相互影响，计及多时间尺度的状态估计问题需研究一段时间内的状态分布，在满足各状态量在不同控制时刻下的交互特性的前提下，最小化一段时间内量测量与状态量之间的误差，从而提高估计精度。

为此，本文构建了一种考虑多时间尺度的电热综合能源系统状态估计方法，通过刻画热力系统状态量的时序分布特性，对电热综合能源系统的状态量进行准确估计，并利用拉格朗日乘子进行非线性优化问题的求解。最后，通过算例分析验证了本文所提模型的有效性，并与稳态模型所得结果进行了对比。

1 电热综合能源系统动态建模

1.1 电力系统模型

电力系统主要关注节点和支路处的有功和无功分布，本文采用交流模型对其进行建模。节点注入功率平衡和支路传输功率平衡分别为[14]：

式中：Pi和Qi分别为节点i的注入有功和无功功率，Ui和θi分别为节点i的电压幅值和相角，Gij和Bij分别为节点i与节点j之间的电导和电纳，θij为节点i与节点j之间的相角差，Pl和Ql分别为线路传输的有功功率和无功功率。

1.2 热力系统模型

热力系统中，热源产生的热能通过工质流动在供/回水网中传递，因此，热力系统模型可分为水力部分和热力部分。本文主要关注运行在质调节的热网，在此场景下，热负荷需求通过调节在热源处供水温度得到满足，并维持网络中管道流量不变。由于水力状态稳定且操作简单，该调节方式在我国得到大量应用[9,15]。

热力部分主要描述系统内的供/回水温度和功率分布。根据换热器模型，节点注入功率[13]为

式中：cp为水的比热容，φ为节点注入功率，q为节点注入流量，Ts和Tr分别为节点供水和回水温度。

混合后的节点温度可表示为

式中：Tout为节点处的混合温度，Tin为节点注入管道的末端温度，mout为节点流出管道的流量，min为节点流入管道的流量。

管道温度分布的动态过程一般通过常系数偏微分方程描述，即管道温度分布与流速v、管道热阻r以及环境温度Ta相关[16]：

应用微分方程进行求解极大地增加了分析的复杂性，本文采用有限差分法对描述管道温度分布动态特性的偏微分方程进行处理。首先，在区间0≤x≤L，0≤t≤th内将式(5)进行离散化，所得离散时空步长τ和h可分别表示为

式中：L和Γ分别表示管道长度和研究时间区间，N和M分别表示差分的时空步数。

在(xi,tk)处将式(5)的偏微分方程差分化，式(5)中各项偏微分可表示为：

式中O(h)和O(τ)分别为关于空间步长和时间步长的高阶余项。

将式(7)和式(8)代入式(5)，则式(5)表示的管道温度分布方程可转换为：

1.3 耦合环节模型

电力系统和热力系统通过耦合设备实现电能和热能的双向传递。常见设备包括热电联产(CHP)机组和电锅炉(electric boiler，EB)等，电热综合能源系统耦合关系如图1所示。

图1 电热综合能源系统耦合关系示意Fig.1 Schematic diagram of the coupling relationship in the integrated electricity and heat system

由图1 可见，电锅炉通过消耗电能为热力系统供热，可将其视为热源和等效电负荷；热电联产机组同时产生电能和热能，可将其视为等效热源和电源。其能量转换关系为[15]

式中：ϕEB、ϕCHP分别为EB 或CHP 机组的输出热功率，PEB和PCHP分别为EB 和CHP 机组的输出电功率，ηEB、ηCHP分别为EB 和CHP 机组的效率系数。

从耦合特性上看，电力和热力系统主要通过CHP 机组和EB 耦合，耦合机组的电力和热力输出分别根据电力系统和热力系统的运行时间尺度进行运行调整。因此，热力系统的动态特性会通过耦合机组最终传递至电力系统，从而使得电热综合能源系统在整体上呈现动态特性[9-10,14-15,17]。

2 多时间尺度状态估计

2.1 状态量与量测量

电力系统中的量测量ze,t一般包括节点电压幅值Ut,me、节点注入功率Pt,me和Qt,me以及支路传输功率P1,t,me和Q1,t,me；状态量xe,t包括节点电压幅值Ut和相角θt。热力系统中量测量一般布置于热源、用户以及网络中的部分管道，在质调节下，系统的水力状态是确定的[7]，即可认为管道流量mme和节点压力pme全可量测。在此基础上，热力系统中的量测量zh,t还包括供水温度Ts,t,me和回水温度Tr,t,me；状态量xh,t包括节点供水温度Ts,t和节点回水温度Tr,t以及节点热功率φt。电热综合能源系统中量测量和状态量分布如式(12)所示，式中下标me 代表量测量。

2.2 多时间尺度特性

考虑系统运行的动态特性，在一定的控制指令下，电力系统能迅速达到稳定状态而热力系统仍然处于动态过程。因此，电力系统的调度周期通常较短，为分钟级；而热力系统调度周期相对较长，通常为几十分钟甚至数小时[9-10]。图2 为电热综合能源系统耦合关系动态特性，图2 中各网格为电力和热力系统的调度周期。由图2 可见，计及电力系统和热力系统的运行时间尺度特性差异，在相同的调度周期内，电力系统调度周期较短、分布较密集；而热力系统调度周期较长，分布较为稀疏。

图2 电热综合能源系统耦合关系动态特性Fig.2 Dynamic characteristics of the coupling relationship in the integrated electricity and heat system

2.3 状态估计模型

状态估计模型的本质是：在一定约束条件下，最小化状态量估计值和量测值的误差。其目标函数可通过最小二乘法构建，表示为[7]：

式中：Je和Jh分别为研究区间内的电力系统和热力系统的状态估计误差；Re和Rh分别为电力系统和热力系统量测误差的协方差矩阵，由历史数据集确定；re和rh分别为电力系统和热力系统中的各项状态估计误差；ξP,t、ξQ,t和ξU,t分别为t时刻节点有功功率、无功功率和电压幅值的估计误差；ξPl,t和ξQl,t分别为t时刻支路有功功率和无功功率的估计误差；ξcp,t和ξϕ,t分别为t时刻耦合机组功率和节点热功率的估计误差；ξT,t和ξTout,t分别为t时刻管道温度方程和节点温度混合方程的估计误差；ξTs/r,t为t时刻节点供/回水温度的估计误差；Δte和Δth分别为电力和热力系统的调控时间间隔；Γe和Γh分别为电力和热力系统的调控时间标志，通常根据实际系统的运行特征作为常数给定[14-17]，其选择不影响模型本质特性。

倘若所研究时间区间Γ为10 min，电力和热力系统调控间隔Δte和Δth分别为1 min 和5 min，则根据式(13)，Γe和Γh分别取10、2 min。在给定的调控时刻集下，各状态量的估计误差可表示为：

式(14)和式(15)分别表示电力系统中由伪量测量构建的节点、支路功率估计误差，其中Pi、Qi、Pl,i和Ql,i分别根据式(1)和式(2)构建。式(16)表示直接量测的节点电压幅值的估计误差[7]。式(17)反映了电力系统和热力系统在不同时间尺度下，量测量和状态量的耦合，由耦合设备运行特性决定。式(18)反映了热力系统中的管道温度的动态特性，即相邻时间点管道温度的相互关系，管道温度分布的动态特性影响着耦合机组的输出热功率，进而通过式(17)依次影响着机组的输出电功率与电力系统中的状态量分布，其中，TI,k+1根据式(9)构建。式(19)和式(20)表示热力系统中由伪量测量构建的节点热功率、混合温度的估计误差。式(21)表示直接量测的节点供/回水温度估计误差[7-10]。

2.4 求解算法

由于式(14)和式(15)的非线性特征，电热综合能源系统状态估计模型本质上是一个难以直接求解的非线性优化问题。因此，引入拉格朗日乘子对非线性约束进行处理。处理后的目标函数可表示为

式中，F(xe,xh)对应于式(14)—式(21)中的各个等式约束，ω为对应于各个约束的乘子系数向量。

处理后的模型为一种典型的非线性优化问题，而本文所提模型仅在优化模型中考虑了线性的多时间尺度耦合约束和本地量测约束，并不影响状态估计模型和求解算法的自身特性。因此，可通过在MATLAB 下调用求解器IPOPT 进行直接求解，其可解性已得到充分验证[7,9-10]。计算机配置环境为Intel Core i7、8GB RAM。

3 算例分析

本文通过2 个算例对所提模型进行分析：1）热力系统动态模型精度验证算例，用于分析所构建差分化热力系统动态模型的计算精度，并且与实测数据、不计及动态特性的热网模型进行对比验证；2）电热综合能源系统状态估计验证算例用于分析所提多时间尺度状态估计模型的准确性，并在不同的量测误差下与不计及动态特性的热网模型进行对比验证，估计精度通过平均相对误差衡量。

3.1 算例1

算例1 采用石家庄某供热工程进行分析[13]，热源为联供机组，管道全长9.25 km，管道直径为1.4 m，室外温度为–10 ℃，热阻系数为0.67 W/(m·℃)。热源通过对供热温度进行调整来满足热负荷需求。此外，为验证模型适用性，供水网中的管道流量是时变的，即构建差分式动态模型对于不同水力状态分布的热力系统同样适用。

动态模型的差分时空步长影响建模精度，本文通过常用的仿真试探法进行选择，因而在模型精度验证和状态估计问题中作为常数给定[14-15,17]。在热源处添加阶跃扰动，研究不同差分时空步长下的管道末端温度的仿真结果，不同差分步长下的仿真精度如图3所示。

图3 不同差分步长下的仿真精度Fig.3 The simulation accuracy with various differential sizes

由图3 可见：当差分的时间步长为6 min 时，不论差分的空间步长如何选择，仿真结果均出现了不同程度的抖动；而当差分的时间步长为3 min 时，仿真结果的稳定性较高。因此，本算例中最终选择动态模型的差分时间步长和空间步长分别为3 min和370 m。

未计及动态特性所得结果和本文所提动态模型的换热站供水温度计算结果对比如图4所示。由图4可见，动态模型准确地反映了供水温度在一段时间内的变化趋势，其平均相对误差为0.27%。当不考虑动态特性时，系统中的状态分布主要取决于实时的功率平衡且忽略了延时特性。因此，所得的温度分布与动态模型相比提前了若干时刻。此外，不计动态特性时，模型所得结果完全由热源处的供水温度决定，其计算结果误差较大，平均相对误差为1.03%。

图4 未计及动态特性模型和动态模型的供水温度结果对比Fig.4 The supply water temperatures of the models with and without dynamic features

3.2 算例2

算例2 采用文献[18]中电热综合能源系统进行分析，并对其中数据进行修正。算例2 系统结构如图5所示。

图5 算例2 系统结构Fig.5 The structure of IEHS in case 2

图5 中，电力系统节点1 和节点2 为发电机组，节点6 为电热耦合点，通过CHP 机组进行耦合。电力系统的调控间隔为1 min，热力系统的调控间隔为5 min。与算例1 类似，动态模型中的差分时空步长通过提前测试以保证热力系统动态模型具有较高的建模精度，本算例采用时间步长为1 min，每根管道分为25 段可获得较高的仿真精度，并作为常数在状态估计中给定。

算例中，以能流计算结果为真实值，在此基础上叠加50 dB 的高斯白噪声作为量测值进行分析，研究不同量测噪声下的状态估计结果从而验证方法的通用性。研究时间区间为10:00—11:00，共1 h。

利用式(14)模型对研究区间进行状态估计，所得电力系统和热力系统估计结果分别如图6—图9所示。图6 和图8 分别为全时段内电力系统节点5电压幅值和热力系统节点3 的供水温度时序分布；图7 和图9 分别为10:30 时刻电力系统和热力系统中整体的节点电压幅值和供水温度估计结果分布。

图6 和图8 以及图7 和图9 分别反映了电力系统中的电压幅值和热力系统中供水温度在时间和空间上的分布规律，即刻画了电热综合能源系统内状态量分布的时空相关性。

图6 各时刻电力系统节点5 电压幅值估计值Fig.6 The estimated U5 during the whole period

图7 10:30 时刻电力系统电压幅值估计值Fig.7 The estimated voltage amplitude at 10:30

图8 各时刻热力系统节点3 供水温度估计值Fig.8 The estimated Ts3 during the whole period

图9 10:30 时刻热力系统供水温度估计值Fig.9 The estimated supply water temperature at 10:30

由图6—图9 可见，本文所提方法进行状态估计所得结果具有较高精度，且从整体上看，估计值的相对误差小于量测值的相对误差，从而验证了在给定误差的量测数据下，本文多时间尺度状态估计模型的有效性。

此外，通过电力系统和热力系统的结果对比可以看出：电力系统中所得到的状态估计误差明显小于量测误差；而在热力系统中，在部分时刻下的估计误差大于量测误差。这是因为在多时间尺度状态估计模型中，对于电力系统的分析本质上是对多时间断面的研究，且各个时间断面相互独立；而热力系统中的状态分布相互影响，在较早时刻产生的估计误差会通过管道温度动态方程逐步向后续时刻传递，即形成了误差传递链，从而影响了部分时刻的状态估计精度。这是由热力系统的物理特性决定的，但从整体上看，本文所提方法仍取得了较好的估计效果。

由图8 可知，所建立的动态模型能有效反映热力系统内温度变化的动态过程，即各节点的供水温度随热源处的温度调节而变化，但由于前序时刻状态差异，所体现的变化趋势仍有部分偏差。将本文所提多时间尺度状态估计模型与传统的基于稳态模型所获取的状态估计结果进行对比，所得估计结果见表1。

表1 多时间尺度/稳态模型状态估计平均相对误差 %Tab.1 The mean relative error of state estimation results of the multi-timescale and dynamic model

由表1 可见：考虑电热综合能源系统的多时间尺度特性后，状态估计的精度得到显著提高。其改善主要在热力系统部分，电力系统估计精度的改善相对较小。这是因为电力系统和热力系统根据热力系统调控尺度进行交互，热力系统模型的变化仅在部分时刻对电力系统产生影响。此外，在算例2 中，热力系统对电力系统的影响完全由联供机组反映，而联供机组在电力系统中可视为等效PV 节点，其功率变化对电力系统运行状态影响较小，因而热力系统模型精度的提高对于电力系统部分状态估计精度的影响相对较小。对于热力系统，由于动态模型利用了热力系统不同时间、位置上量测数据的相关性，使得状态估计所得结果精度得到明显提高。而在质调节热网中，由于节点注入流量是固定的，节点功率与温度呈线性关系。因此，节点功率的估计误差与温度的估计误差整体上近似。

对不同量测误差对动态状态估计结果的计算精度和计算时间的影响进行进一步研究。取量测误差的噪声分别为40、50、60 dB，不同量测误差下，基于稳态模型状态估计方法和动态状态估计方法所得的计算误差和计算时间对比见表2。

表2 不同量测误差下的状态估计误差和计算时间结果对比Tab.2 The mean relative error and computation time with different measurement errors

由表2 可见，不同的量测误差对于本文所提状态估计方法计算时间差异的影响较小。这是因为状态估计问题的计算时间主要由计算规模决定，而优化问题中所涉及的优化变量数目和约束条件数不随量测误差的增减而变化；而多时间尺度状态估计所需的计算时间明显多于稳态估计方法，其原因在于差分化的热力系统动态模型中引入了较多的分段变量，使得优化问题规模显著增大。在本案例中，多时间尺度状态估计方法所需的计算时间约为稳态状态估计方法的3～4 倍。

对比不同量测噪声下的状态估计结果，多时间尺度状态估计方法和稳态状态估计方法所得结果均随着量测误差的增大而增大，且稳态状态估计结果精度受量测误差的增大下降更为明显。这是因为所研究系统规模较大，系统的延时特性十分明显，而稳态状态估计方法忽略了各个时间尺度下的耦合关系，使得时序的状态估计过程中各个时刻点的估计结果较差；在多时间尺度方法中，除了网络运行约束外额外在优化问题中添加了不同估计时刻之间的数值联系，该处理方法等效在优化模型中添加了约束条件，使得多时间尺度状态估计方法所得结果精度相对较高。

另一方面，由于量测项在优化目标中所对应的权重随着量测误差的增大而减小，因此，量测项在状态估计问题中所产生的影响随着量测误差的增大而减小，使得状态估计结果的误差增长趋势相对于量测误差的变化趋势较为缓慢。