基于多新息辨识的电力系统节点惯量估计方法

2022-08-09李元臣文云峰蒋小亮林晓煌

电力自动化设备 2022年8期

李元臣，文云峰，叶希，蒋小亮，林晓煌

（1. 湖南大学电气与信息工程学院，湖南长沙 410082；2. 国网四川省电力公司，四川成都 610041；3. 国网河南省电力公司经济技术研究院，河南郑州 450000）

0 引言

随着新能源的规模化开发和大容量直流输电工程的兴建，电力系统中大量同步电源逐步被非同步电源所替代［1］。由于新能源和直流均通过电力电子换流器并网，其传输功率与电网频率解耦，致使系统惯量大幅度降低，恶化了系统的频率稳定性［2］。为此，有必要对系统惯量水平进行评估，以便运行人员及时感知和调控系统惯量储备情况，避免惯量水平过低而削弱系统稳定性［3］。

与此同时，自然风光资源分布存在明显的地域差异，使得原本相对均衡的惯量资源分布格局被打破。新能源集中接入区域低惯量特征显著，与同步电源富集区域的高惯量态势形成强烈对比［4-5］。大容量直流在一定程度上阻断了异步互联区域电网之间功率快速支援，促使惯量空间分布差异愈加明显。惯量空间分布特性的凸显将限制系统对新能源的消纳能力。高比例非同步电源接入的惯量薄弱区域在有功冲击下频率变化速率RoCoF（Rate of Change of Frequency）更快，为确保频率偏差在安全范围内，电网运营商通常不得不被动弃风弃光以维持区域电网的惯量水平［6］。因此，对于非同步电源高占比的电力系统，惯量评估的研究工作不仅需着眼于系统整体惯量水平，还需关注节点／区域层面的惯量空间分布情况，以便发现惯量薄弱节点／区域，有针对性地提升新能源／直流接入的适应性。准确评估系统惯量空间分布情况，对研究系统频率动态特性，指导新能源并网与配置惯量补偿措施具有重要意义。

目前，电力系统惯量评估方法的相关研究主要集中于系统层级，分为3 类：①基于开停机状态监测的系统惯量计算；②基于频率测量数据分析的系统惯量评估；③基于机电扰动传播的系统惯量估计。第①类方法通过监测同步发电机出口断路器状态，将所有开机的同步发电机惯量与其容量的乘积进行累加得到系统总体旋转动能或等效惯量；该评估方法需提前给定系统中所有同步发电机惯量的准确值，不能体现异步电动机负荷和虚拟惯量资源的惯量贡献，因而只适用于特定系统的同步惯量评估。第②类方法利用频率／功率量测数据，基于频率曲线多项式拟合和RoCoF 估算系统惯量；该方法难以消除频率信号振荡分量的影响，评估精度有限，由于一般需假定已知有功扰动大小和事故发生时间，通常仅适用于事故后离线分析［7-9］。第③类方法依据机电扰动传播波信号提取系统等值惯量参数，可分析小扰动下的惯量分布特征，但评估精度可能受到扰动位置以及测量装置布点情况的影响［10］。

文献［11-12］将系统分区并建立区域等值频率响应模型，利用系统辨识方法实现区域惯量的估计。通过分区域进行惯量评估，可在一定程度上反映各区域电网的惯量水平，但如何更为精细地量化节点层级的惯量空间分布情况仍是目前需解决的问题。文献［13］基于惯量中心频率指标分析了系统惯量对不同节点的惯量支撑力度；文献［14］建立扰动传播速度与惯量之间的映射关系，进而评估电网惯量分布情况；文献［15］从小扰动下节点频率频谱中提取指标对惯量时空特性进行可视化，实现对电力系统惯量空间分布的快速估计。上述工作在惯量空间分布估计方面进行了有益尝试，但适用于不同扰动规模的惯量估计方法仍需深入研究。

本文提出一种基于多新息辨识的电力系统节点惯量估计方法。通过分析系统惯量资源对节点的惯量支撑作用，构建节点惯量表征模型，统一描述不同扰动规模下各节点有功-频率动态过程，有效克服扰动规模对惯量空间分布估计结果的影响；采用多新息辨识方法，将标量新息扩展为向量新息，加速算法收敛，增强辨识精度，快速、准确求解惯量辨识模型中未知参数，进而得出系统内所有节点的等效惯量，使得调度员可清晰纵览惯量空间分布情况。最后通过仿真计算验证了所提方法的有效性和适应性。

1 电力系统惯量基本概念

电力系统等效惯量可由同步电源转动惯量、新能源／储能虚拟惯量等共同提供。电力系统惯量响应实质上是系统内部能量转移的过程。当系统有功平衡被打破时，各类惯量资源自主响应频率变化，通过主动释放或吸收能量来减少系统有功波动、减缓频率变化速率［16］。利用等值转子运动方程可描述系统遭受有功扰动后的有功-频率动态变化过程，如式（1）所示。

式中：Hsys为系统等效惯量；Δf为系统频率偏差；ΔPm为机械功率增量；ΔPe为电磁功率增量；D为阻尼系数；t为时间。式（1）表明，当系统发生有功不平衡时，以同步发电机为代表的惯量资源将立即响应系统频率变化，向系统提供／吸收电磁功率，以抑制系统频率变化。由于同步发电机一次调频具有动作延迟且调整速率较慢，惯量响应对于维持故障瞬间的频率稳定性具有重要支撑作用。

1.1 电力系统惯量空间分布

大量弱／零惯量非同步电源接入电网，打破了原本相对均衡的惯量分布格局，使得电力系统惯量空间分布特征显化。电力系统惯量空间分布情况可由不同子区域的等效惯量集合表示。在电力系统运行过程中，若能识别出惯量资源对系统中各节点的惯量支撑力度，解析出节点惯量值，则当全网各节点的等值惯量确定后，系统惯量空间分布态势自然明晰。系统各节点惯量的集合ΦH为：

式中：Hi（i=1，2，…，n，n为节点总数）为节点i的惯量值，表示有功扰动发生后系统惯量对该节点频率变化的阻碍作用。ΦH内包含一般节点（不与发电机直接相连的节点）以及同步发电机节点（与同步发电机直接相连的节点）的节点惯量，其可反映整个系统惯量空间分布的具体情况。在提出惯量空间分布估计模型和方法前，下面先对节点惯量水平的表征方法进行阐述。

1.2 节点惯量水平表征

节点惯量是在电力系统能量波动过程中，系统惯量阻碍节点频率变化的固有属性［15］，其大小表征该节点在所有惯量资源共同作用下抵抗有功扰动的能力。具体而言，节点惯量较大表示系统内惯量资源能够为该节点提供更多的惯量支撑作用，在有功扰动下能够更好地抑制该节点频率变化，降低其频率变化速率；反之，节点惯量较小，则说明系统为该节点提供惯量支撑的力度相对有限，该节点在有功扰动下频率变化速率相对较快。图1 为不同系统惯量支撑作用下的节点频率变化情况。由图可知，相同扰动故障下，同一节点在不同系统惯量支撑作用下具有不同的频率变化趋势，随着系统惯量的下降，节点频率变化速率越快，暂态频率极值可能越低，节点所呈现出的抗扰动能力弱化。

图1 不同系统惯量值下的节点频率变化曲线Fig.1 Variation curves of nodal frequency under different system inertia values

节点惯量反映了系统惯量资源对节点的惯量支撑能力，该支撑作用体现于惯量响应时节点的有功-频率变化过程之中。为使一般节点的惯量表征与发电机惯量处于同一量纲框架下，将式（1）推广至一般节点，建立一般节点k的节点惯量、频率偏差和功率波动的关系模型，忽略阻尼后可得：

进一步整理可得Hk为：

式中：Δfk为一般节点k的频率偏差；ΔPk为一般节点k的功率增量，描述电力系统暂态或准稳态过程中流过节点的功率波动。式（4）表明，一般节点k的惯量可由该节点的功率波动与节点频率偏差来表征，是对节点阻抗有功扰动作用的定量描述。对同步发电机节点j而言，在惯量响应阶段，同样有：

式中：Hj为同步发电机j的惯量；Δfj为同步发电机j的频率偏差；ΔPm，j为同步发电机j的机械功率增量；ΔPe，j为同步发电机j的电磁功率增量。在惯量响应阶段由于一次调频尚未动作，ΔPm，j=0，式（5）可简写为：

对与同步发电机j相连的发电机节点j而言，ΔPj=-ΔPe，j，此时式（3）与式（6）具有一致性，即式（3）适用于描述一般节点和发电机节点的有功-频率动态过程，且式（3）所表征的节点惯量值可以使得一般节点的惯量与同步发电机惯量表征具备统一性。

由上述分析可知，节点惯量值蕴含于惯量响应阶段的节点功率波动与节点频率偏差的时间序列数据中。后文的工作重点在于构建评估模型和高效求解方法从该时序数据中提取节点惯量值，进而估计系统惯量空间的分布情况。

2 惯量空间分布的估计模型

为实现惯量空间分布的精细化评估，本文通过估计系统内所有节点惯量来获取惯量的空间分布，以协助运行人员清晰纵览整个系统的惯量空间分布情况，明确不同节点频率在有功扰动情况下受系统惯量资源的支撑力度。由于有功扰动后，一次调频等频率调节控制措施存在死区设置，其动作速率相对同步发电机的惯量响应具有明显的滞后性。本文提出的节点惯量估计方法考虑频率波动过程中惯量响应速度最快的特点来估计节点惯量。对应于有功扰动后系统惯量响应，式（3）所示节点有功-频率变化过程可用图2所示节点惯量响应模型表示。

图2 节点惯量响应模型Fig.2 Model of nodal inertia response

频域下节点惯量响应模型的传递函数表达式为：

节点惯量响应模型中描述节点惯量响应特性的单位脉冲响应可表示为：

由式（8）可知，连续传递函数的零极点增益值（即初始时刻单位脉冲响应值）等于2 倍节点惯量Hi的倒数。则上述过程可归纳为首先辨识出节点有功-频率过程的连续传递函数模型，然后获取其单位脉冲响应，得到初始时刻的响应值，进而估计节点惯量。

受测量噪声、高次谐波等大量随机因素的干扰，单纯基于式（3）采用机理建模结合模型优化的惯量估计方法鲁棒性较差。为避免受单一时刻数据测量精度的影响，采用统计辨识方法构建输出误差滑动平均OEMA（Output Error Moving Average）模型［17］，用以描述节点频率偏差与节点注入功率增量之间的动态关系，从而将节点惯量估计问题转化为OEMA模型中的参数辨识问题。OEMA模型可表示为：

式中：u(t)为OEMA 模型的输入序列，即节点注入功率增量；y(t)为OEMA 模型的输出序列，即节点频率偏差；v(t)为零均值高斯随机白噪声序列；x(t)为内部变量序列，表示系统实际不可观测信息；A(z)、B(z)、D(z)表达式如式（10）所示。

式中：a1、…、ana、b1、…、bnb、d1、…、dnd为OEMA模型中多项式待定系数；na、nb、nd分别为A(z)、B(z)、D(z)的阶次。OEMA 模型为离散模型，其对应的离散传递函数可表示为：

采用双线性变换法［18］（Tustin 变换法）将传递函数模型式（11）从离散形式转化为连续形式。双线性变换法是自然对数函数的一阶估计法，可采用式（12）所示映射函数实现z平面到s平面的映射。

式中：T为采样周期。OEMA模型转换后的连续传递函数一般形式为：

式中：b′n-1、…、b′1、b′0和a′n、…、a′1、a′0为连续传递函数模型中的待定系数。由式（13）可得求取Hi的一般方法，如式（14）所示。

式中：G′i(s)、g′i(t)分别为OEMA 模型转换为连续模型后节点i有功-频率变化过程的连续传递函数、描述节点惯量响应特性的单位脉冲响应。

至此已建立数学模型描述节点有功-频率变化过程，构建了细化至节点层级的惯量空间分布估计模型。为实现对待定参数向量b=[bn-1，…，b2，b1]、a=[an，…，a2，a1]的快速求解，满足运行调度过程中的快速监测需求，下面将提出基于多新息辨识的节点惯量估计方法。

3 基于多新息辨识的节点惯量估计方法

3.1 多新息辨识

由前文分析可知，评估节点惯量首先需要辨识出OEMA 模型中的未知参数，本文基于多新息辨识方法实现该功能。新息指能够改善参数估计精度或状态估计精度的有用信息［17］。相比于传统系统辨识，多新息辨识将原本标量新息拓展为向量新息（多新息）。多新息辨识基于可变的新息长度，可以充分利用同步相量测量单元PMU（Phasor Measurement Unit）当前及历史量测数据，提高数据利用效率，加速算法收敛，增强惯量空间分布的辨识精度，快速、有效地辨识出OEMA 模型中的未知参数，进而得出整个系统的惯量空间分布情况。

在多新息辨识方法中，将节点频率偏差、节点功率波动单一新息拓展为数据窗［t-p+1，t］共p组数据的向量新息，建立包含内部变量x(t)、输出向量y(t)、噪声向量v(t)和信息向量φ(t)的多新息矩阵，如式（15）所示。

式中：X(p，t)为堆积内部变量；Y(p，t)为堆积输出向量；V(p，t)为堆积噪声变量；ψ(p，t)为堆积信息矩阵；x、y、v、φ分别为向量序列x(t)、y(t)、v(t)、φ(t)中的元素。式（15）中信息向量φ(t)的表达式为：

在电力系统实际运行中，PMU 采集的频率数据通常包含随机噪声，难以完全反映节点实际的频率动态。因此采用辅助模型辨识思想，借助系统可测得的节点功率，用辅助模型x̂(t)的输出代替实际频率偏差的动态变化，进一步提高OEMA 模型中未知参数的辨识精度。辅助模型为：

综合式（9）、（15）、（17），建立计及辅助模型辨识、多新息辨识方法的OEMA模型，如式（18）所示。

式中：ψS(p，t)为辅助模型中堆积信息矩阵。设数据长度为t，采用辅助模型多新息增广最小二乘算法［17］，以节点功率增量作为输入量、节点频率偏差作为输出量，在每步迭代计算中，将信息向量φ(t)中不可测项x(t)用辅助模型输出估计值代替，反过来估计值由前一次迭代的参数估计计算，二者进行递阶计算，更新未知参数的估计结果，快速有效跟踪节点惯量的变化。每一次更新后，参数估计结果可表示为：

为使预测误差尽可能小，给定准则函数为：

3.2 惯量空间分布估计流程

基于多新息辨识的惯量空间分布估计整体流程见附录A 图A1。首先收集系统内节点频率与输电线路传输功率等时序数据，对所获数据进行预处理，识别扰动发生时刻。然后解析节点惯量、节点频率偏差和节点功率波动三者之间的关系；构建OEMA模型，解析该模型中未知参数与节点惯量之间的耦合关系。其次采用多新息辨识方法辨识模型中的未知参数，为求解节点惯量、估计惯量空间分布情况做准备。最后计算得到各节点惯量的大小，进而估计整个电力系统惯量空间分布情况。具体步骤如下。

1）步骤1：获取PMU 实时采集的系统惯量空间分布情况所需数据（包括节点频率、输电线路传输有功功率、负荷大小）；对获取的数据进行预处理，利用样条插值方法弥补缺失数据，采用低通滤波器消除数据测量噪声，分别基于系统额定频率、稳态时／扰动初始时刻流过节点的有功潮流对该节点预处理后的频率、功率数据进行归一化处理，形成标准形式的数据集。

2）步骤2：基于小波多分辨率分析原理［20］，利用多贝西小波基（Daubechies Wavelet）方法将所获取的频率信号在一系列相异的空间上进行分解，提取各层小波细节系数确定模极大值点，从而检测出频率信号的奇异点，得到大扰动故障发生的初始时刻近似值，获取大扰动初始时段，进而估计系统惯量空间分布情况；在小扰动场景下，直接利用量测频率-有功功率数据估计系统惯量空间分布情况。

3）步骤3：构建OEMA 模型描述节点频率偏差与节点注入功率增量之间动态关系，将节点惯量估计问题转化为模型中参数辨识问题。通过求取表征节点惯量响应过程的连续传递函数G(s)的零极点增益值来估计节点惯量大小。

4）步骤4：基于步骤1 得到的时间序列数据，采用多新息辨识方法求解步骤3 中OEMA 模型中的未知参数。将OEMA 模型中标量新息拓展为向量新息，如内部向量x(t)、输出向量y(t)、噪声向量v(t)等，建立多新息矩阵；然后建立式（17）所示辅助模型预测系统内部不可测变量；以式（20）为目标函数，利用式（19）快速更新未知参数的辨识结果，辨识出OEMA 模型中的未知参数，得到对应的式（11）所示离散传递函数。

5）步骤5：采用双线性变化法构建z平面与s平面之间的映射关系，将步骤4 中离散传递函数连续化，求解连续传递函数初始时刻激响应值，即可得到2 倍节点惯量的倒数值，进而估计系统各节点惯量的大小。

6）步骤6：在求解出系统所有节点惯量值后，即可明晰细化至节点层级的惯量空间分布情况。通过融合计算求得节点惯量大小与节点地理位置信息，进行可视化展示，得到系统惯量空间分布结果，反映系统不同区域惯量支撑力度大小，辅助运行人员实时感知和预警系统惯量空间分布情况。

4 算例分析

为了验证本文所提出的基于多新息辨识的电力系统节点惯量估计方法的有效性，通过PSD-BPA 仿真软件搭建IEEE 39 节点系统。该系统的拓扑结构如附录A图A2所示，系统总负荷为6150 MW。所用计算机配置为CPU Intel Core i5，内存为8 GB。

4.1 本文方法的有效性验证

为分析本文方法对不同扰动故障规模的适应性，设置不同规模的故障，其中大扰动故障包括：①大负荷脱网，在10 s 时切除节点20 处负荷（总计615 MW，约占总负荷的10%）；②机组跳闸，在10 s时设置同步发电机G5跳闸，系统有功出力改变量为508 MW。考虑到系统各同步发电机有功出力均大于200 MW，准稳态运行时仅设置小负荷波动，在节点4、7、15、20、23、27处设置小负荷变化共计100 MW。

对量测得到的同步发电机有功功率和母线频率信号注入高斯白噪声（其幅值为测量信号的1%），以模拟实际系统中PMU 的测量噪声。将所提方法与基于最小二乘一次完成算法的系统辨识工具箱SIT（System Identification Toolbox）方法［11］以及最小二乘迭代（LSI）算法［19］进行对比分析，结果见附录A表A1—A3，同步发电机惯量以其额定容量为基准。

对比由本文方法得到的各台同步发电机惯量估计值与实际值可知：各工况下所有同步发电机惯量估计误差均在6%以内，大扰动故障下同步发电机惯量最小误差为0.11%，最大误差为5.74%；准稳态运行时同步发电机惯量最小误差为1.30%，最大误差为5.62%；所有工况下同步发电机惯量估计的平均误差为3.06%。这说明本文所提出的基于多新息辨识的节点惯量估计方法在不同规模扰动下能够准确估计电力系统的惯量空间分布情况。而基于SIT、LSI 方法的同步发电机惯量估计结果在大扰动故障下的最大估计误差分别达到17.26%、19.52%，平均误差分别为9.40%和11.44%；准稳态运行时基于SIT、LSI 方法同步发电机惯量估计的平均误差分别达到12.8%、19.12%；且部分工况下同步发电机惯量估计结果精度严重不足，最大误差分别达到26.21%和39.02%。由此可知，相较一般系统辨识算法，基于多新息辨识算法估计惯量节点惯量大小的准确度更高，且对不同大小和不同类型的扰动适应性更好。

4.2 所提方法对不同扰动位置和大小的适应性分析

4.1 节已说明基于多新息辨识方法估计惯量空间分布情况的有效性以及对不同规模扰动的适应性。进一步，通过在测试系统中的不同位置设置相同类型扰动，验证不同扰动位置下本文方法的适应性。故障场景设置如下：场景1，10 s 时切除节点3 处200 MW 负荷；场景2，10 s 时切除节点8 处200 MW 负荷；场景3，10 s 时切除节点15 处200 MW负荷；场景4，在节点4、7、15、20、23、27 处总计添加100 MW的小扰动负荷。

图3 为4 种场景下IEEE 39 节点系统中各节点惯量的估计值。由图可知，针对系统所有节点而言，在不同扰动发生位置下，由本文方法估计得出的节点惯量估计结果的最大值与最小值平均偏差为0.25 s，以同步发电机一般惯量（5 s）为衡量标准，该误差在5%以内。这说明相比于仅从频率信息得到节点惯量的方法，本文方法对扰动发生位置具有很好的适应性。

图3 扰动故障不同位置及规模下节点惯量估计结果Fig.3 Estimation results of nodal inertia for disturbance faults of different locations and sizes

4.3 系统惯量空间分布情况的可视化展示

将IEEE 39 节点系统拓扑结构映射到二维平面，以此代表实际地理接线图信息，在该二维空间平面上计算节点惯量大小并做可视化处理。根据附录A 表A4 所示同步发电机惯量设置2 种场景：场景5中左侧区域5 台同步发电机惯量较大，右侧区域同步发电机惯量较小；场景6 与场景5 相反，右侧区域同步发电机惯量均大于左侧区域。

附录A 图A3（a）、（b）分别展示了场景5、6 下IEEE 39 节点系统惯量空间分布可视化效果。图A3（a）中左侧区域节点惯量高于4 s，中间区域节点惯量介于3～4 s 之间，右侧区域惯量低于3 s，节点惯量从左到右逐渐减弱；惯量水平高于4 s 的节点约占12.8%，惯量水平在3～4 s之间的节点约占53.8%，惯量水平低于3 s的节点约占33.3%。图A3（b）中，节点惯量在空间上呈现出左翼小右翼大的分布态势。均与预设场景下同步发电机惯量空间分布相同，这说明了本文方法对估计系统惯量空间分布情况的有效性。对比分析图A3（a）、（b）可知，系统惯量虽然由发电机惯量所决定，但其在空间上并非均匀分布，大惯量发电机集群联网区域高惯量特征明显，与小惯量发电机集中并网区域的低惯量特征形成鲜明对比，二者之间的区域形成高惯量到低惯量的逐步过渡。

将场景5 中G3替换为同容量的风电场（风电场聚合等效惯量为2.7 s）后的系统惯量空间分布评估结果如附录A 图A4（a）所示，与图A3（a）对比可知，G3向系统提供的惯量下降时，其对周围区域的惯量支撑能力减弱，相应节点的频率抗扰动能力降低，惯量水平低于3 s 的节点比例由33.3%升高为48.7%。在场景5 中增强同步发电机组G5惯量（惯量由2.6 s变为3.6 s）后的系统惯量空间分布评估结果如附录A 图A4（b）所示，与图A3（a）对比可知，同步发电机组G5所具有的惯量增强时，其周围区域的频率抗扰动能力显著提升，惯量水平低于3 s 的节点比例由33.3%降低为23.0%。这说明节点惯量空间分布随系统惯量资源配置的变化而变化，从而在时间上、空间上呈现出不同的态势，当运行周期内所有时刻的节点惯量空间分布形成后，可以综合获取系统惯量时空分布特性。