随机数据驱动的电力系统小干扰稳定在线评估方法

周书宇，蔡国伟，杨德友，王丽馨

（东北电力大学电气工程学院，吉林省吉林市 132012）

0 引言

随着现代电力系统的快速发展，大型区域电网互联使得电力系统小干扰稳定问题日益突出，这给电力系统的安全稳定运行带来了巨大挑战。准确评估小干扰稳定特征参数对系统稳定运行状态下的低频振荡预警、及时采取有效调控措施抑制低频振荡现象具有重要意义［1-4］。

近年来，同步相量测量单元（PMU）逐步成为监测电力系统动态运行的主要工具，将PMU 采集到的量测数据按照数据性质可分为暂态振荡数据和随机响应数据。暂态振荡数据指系统受到大扰动冲击后的自由振荡响应。随机响应数据由负荷随机波动等自然扰动所激发，外在表现为类噪声数据，这种数据在系统中时刻存在，且包含系统小干扰稳定特征参数，利用这类数据辨识特征参数更符合系统振荡预警要求，避免了传统方法只能在大扰动特殊工况下进行特征参数提取，针对这类数据的辨识受到国内外研究人员的广泛关注［5-14］。

基于随机响应数据的辨识方法大致可分为预处理方法和直接处理方法。预处理方法先通过预处理技术将随机响应信号处理为自由振荡信号，再利用针对暂态振荡信号的方法进行参数辨识。文献［8］利用随机量减技术（RDT）从系统随机响应信号中获取自由振荡信号，再通过Prony 方法辨识振荡模式。文献［9］将自然激励技术（NExT）与总体最小二乘-旋转不变技术的信号参数估计（TLSESPRIT）方法相结合，通过滑窗的方式获取模式参数。 文献［10］引入NExT 和多参考点复指数（PRCE）法在获得振荡频率和阻尼比的基础上得到模式模态，同时讨论了信噪比对辨识结果的影响。预处理方法中预处理技术和参数辨识方法是2 个独立算法，两者在数据运算过程中误差的叠加可能导致难以保证辨识精度。直接处理方法则可避免算法之间的累加误差，直接利用随机响应数据进行参数辨识。文献［11］将自回归滑动平均模型（ARMA）法应用于特征参数辨识，但这一辨识过程中存在模型定阶问题。文献［12］提出一种贝叶斯方法在频域中辨识模态参数，大大简化了识别过程。文献［13］利用随机子空间（SSI）算法获得系统状态矩阵，通过特征分解求得模式特征信息。文献［14］则在SSI 算法基础上引入递归过程来降低LQ 分解计算的复杂度，提高了算法的在线辨识能力。

传统的动态模式分解（DMD）法能从全局角度提取振荡特征参数［15］，但由于随机响应数据中小干扰稳定特征参数较弱，DMD 法难以对系统随机响应数据进行有效辨识，且DMD 法只能选择在固定的低维正交空间下构建系统近似状态矩阵，加大了系统低维模型的误差。为了克服传统DMD 法的不足，本文提出了一种子空间最优模式分解（subspace optimal dynamic mode decomposition，Sub-OpMD）算法，该方法将基于正交投影的矩阵线性变换方法以及共轭梯度优化算法引入DMD 法，用于在随机响应数据中辨识系统小干扰稳定特征参数。利用IEEE 16 机68 节点系统和实际电网PMU 数据进行仿真验证。相较于SSI 算法，本文所提方法辨识准确度更高，充分满足了现代大型区域互联系统对安全稳定评估的精度要求。

1 系统随机响应信号模型

在电力系统小干扰稳定分析中，运行点处每台发电机的振荡行为可表示为［12］:

式中:M、D、K分别为惯性时间常数、阻尼系数和同步力矩系数；x(t)为发电机在平衡点附近的功角偏差；F(t)为系统中的自然激励，主要表现为负荷的随机波动。

在自然激励条件下式（1）的解x(t)可表示为:

式中:ψi为模式i的模态；ηi(t)为模式i满足运动解耦方程的标量模态。

式中:ξi为模式i的阻尼比；ωi=2πfi，其中fi为模式i的振荡频率；pi(t) 为模式i的激励，pi(t)=(t)/(ψiTMψi)。

在实际电力系统中由于受到负荷投切、新能源出力不稳定等自然激励的影响，PMU 量测的系统随机响应数据外在表现为混乱的类噪声信号。通过观察式（2）和式（3）可知，在这种自然激励条件下，系统看似杂乱的随机响应数据中实则包含了丰富的系统动态信息，利用可靠的随机数据辨识方法可以有效地从中提取系统小干扰特征信息（振荡频率、阻尼比、模态）。

2 小干扰稳定评估方法

2.1 DMD 法

DMD 法作为一种高维系统降阶算法已广泛应用于流体力学等研究领域［16］，在2015 年首次由Barocio 等学者引入电力系统振荡模式辨识［15］。下面简单介绍一下DMD 法原理。

在研究电力系统小干扰稳定分析时，可以将电力系统用状态空间模型表示为［7］:

式中:A为系统高维映射矩阵。对A进行特征值计算可得到系统特征信息。

由于PMU 采集的信号均为离散量，则可以将式（4）离散化表示为xk+1=Axk，其中k为离散化后的一个采样点，设置连续状态变量的采样时间间隔为Δt，则存在xk=x(kΔt)。利用m个所采集的数据序列构造2 个输入矩阵X1=[x1，x2，…，xm]和X2=[x2，x3，…，xm+1]。根据式（4）的离散化表示可得出X1和X2之间存在一个映射矩阵A，使得:

将高维状态矩阵A投影到固定的低维正交空间上，得到A的低维近似状态矩阵。

式中:Uq为低维正交空间，可以通过对X1进行奇异值分解（singular value decomposition，SVD）得到，q为子矩阵X1的秩。

式中:Uq和Vq分别为X1的左奇异矩阵和右奇异矩阵；Sq为X1的奇异值矩阵。

通过求解Frobenius 范数的最小值优化问题以计算低维近似矩阵A～。

则系统低维近似状态矩阵A～ 为:

DMD 法仅仅适用于从大扰动数据中辨识小干扰特征参数，对外在表现为类噪声的随机响应数据适应性较差，无法有效辨识。此外，DMD 法利用SVD 求取固定维数的低维正交空间，而低维正交空间的维数选择不当会导致A～ 无法获取最佳的系统动态信息。

2.2 Sub-OpMD 算法

本文提出的Sub-OpMD 算法在DMD 法的基础上引入随机子空间算法中基于正交投影的矩阵线性变换［17］以及共轭梯度算法［18-19］，来弥补DMD 法的不足。

构造数据矩阵X3=[x3，x4，…，xm+2]和X4=[x4，x5，…，xm+3]，并形成2 个新的数据矩阵Xp和Xf。

将Xf的行向量投影到Xp的行空间上，构造正交投影矩阵O。

对正交投影矩阵O进行SVD 分析:

式中:U和V分别为O的左奇异矩阵和右奇异矩阵；S为O的奇异值矩阵。

正交投影变换能有效压缩数据，使数据降维，并减少数据降维过程中的信息丢失，SVD 可以在数据降维的同时提取数据所含有的特征，这2 种矩阵线性变换的应用使算法可以有效辨识随机响应数据。本文定义U1和U2分别为U的前n行子矩阵和后n行子矩阵，U1和U2之间同样存在一个包含系统动态信息的高维映射矩阵，满足U2U-11→。

引入正交矩阵L，与式（8）类似，则系统低维近似矩阵构建问题可表示为:

式（13）可进一步改写为:

式中:Z=UT2L(LTU2UT2L)-1LTU2。

则式（13）的双变量（L，Aˉ*）优化问题可以转化为单变量（L）优化问题:

共轭梯度算法首先计算初始梯度G0与初始搜索方向H0。

式中:L0为正交矩阵L的初始矩阵。

在梯度方向G上将L参数化表示为:

Lk+1(t)≐LkQΓQT+PTQTt∈[0，1]（17）式中:Lk为第k次迭代下的正交矩阵；P和Q分别为第k次迭代下的搜索方向矩阵Hk的左奇异矩阵和右奇异矩阵，即Hk=PNQT，其中N为Hk的奇异值矩阵；Γ和T分别为矩阵Nt对角线元素的余弦值和正弦值构成的矩阵。

迭代更新梯度G和搜索方向H:

其中，修正项Δk+1为:

式中:· 表示求内积；τ(Gk)=Gk-[LkQTk+P(I-Γk)]PTGk，其中Γk和Tk分别为矩阵Ntk对角线元素的余弦值和正弦值构成的矩阵，tk为第k次迭代优化t的最小值；τ(Hk)=-(LkQTk+PΓk)NQT。

利用式（18）和式（19）迭代更新得到L*，直至g(Lk+1)和g(Lk)之间的差值在误差允许范围内，即g(Lk+1)-g(Lk)＜ε时，停止迭代。利用迭代优化得到L*，可进一步得到最优低维近似状态矩阵。

利用共轭梯度算法可以实现式（15）单变量优化问题的求解，目的是计算最佳维度的低维正交空间L*，并在此正交空间上找到系统低维近似模型Aˉ*与L*的最优组合［18-19］。所构建的最优系统模型Aˉ*可以更好地描述系统动态特征。

3 系统特征参数辨识

小干扰稳定性分析过程中特征参数可以确定某振荡模式的频率、振荡幅值衰减速度（阻尼比）和状态变量之间的相对振荡关系（模态）。通过对Sub-OpMD 算法得到的最优近似状态矩阵Aˉ*进行特征值分解，可以得到上述振荡特征参数。

式中:Λ=diag(λ1，λ2，…，λi，…，λl)，其中λi为系统模式i的特征值，l为最优近似状态矩阵Aˉ*的维数；Φ=[φ1，φ2，…，φi，…，φl]，其中φi为系统模式i的右特征向量。

计算系统模式i的振荡频率fi和阻尼比ξi［20］。

模式i的模态ψi表示为［18］:

利用Sub-OpMD 算法计算小干扰特征参数的详细过程如图1 所示。

图1 Sub-OpMD 算法辨识小干扰特征参数过程Fig.1 Process of Sub-OpMD algorithm for identifying small signal characteristic parameters

4 算例分析

4.1 IEEE 16 机68 节点系统

本文采用IEEE 16 机68 节点系统作为仿真算例，该系统可划分为5 个区域，详细接线图见附录A图A1。通过将系统在工作点附近线性化并进行特征分解可以得到15 种振荡模式，其中4 种区域间振荡模式是本文研究的重点。为了全面评价Sub-OpMD 算法的性能，将基于模型的小信号稳定分析（SSSA）所得振荡频率、阻尼比结果作为参考值，即模式1 下振荡频率为0.375 2 Hz，阻尼比为11.35%；模式2 下振荡频率为0.504 7 Hz，阻尼比为7.25%；模式3 下振荡频率为0.554 9 Hz，阻尼比为2.80%；模式4 下振荡频率为0.749 1 Hz，阻尼比为6.99%。

为了模拟系统中负荷随机波动所激发的系统随机响应，本文在系统所有负荷中加入基准值为2%的随机扰动，以16 台发电机的角频率、功角、发电机出口功率作为输入信号，在OptiPlex 7080 工作站（Intel Core i7-10700 CPU，2.90 GHz）上验证所提算法的有效性。

利用Sub-OpMD 算法对10 min 内系统随机响应信号进行辨识，采用滑动数据窗方法，分析窗长为2 min，滑动窗长为0.2 s，共得到2 400 个窗口的辨识结果。4 种模式辨识结果波形图见附录A 图A2。由图A2 可以看出，系统振荡频率和阻尼比在均值附近波动，这种波动具有一定随机性，且阻尼比的波动程度明显大于振荡频率。

对10 min 内系统随机响应信号辨识结果进行统计分析，将统计结果与SSI 算法和SSSA 进行比较，比较结果如表1 所示。从表1 可知，Sub-OpMD算法得到的4 种区域间振荡模式的频率和阻尼比均值更接近SSSA 得到的参考结果，同时标准差的统计结果小于SSI 算法得出的结果，说明本文提出的方法在振荡特征参数辨识上具有较高的可靠性，这种现象在强阻尼模式下更为明显。

表1 区域间小干扰特征参数辨识结果对比Table 1 Comparison of identification results for inter-area small signal characteristic parameters

本文验证了Sub-OpMD 算法在不同幅值负荷波动下参数辨识的鲁棒性。由于强阻尼模式表示结果波动较大，验证效果更加显著，故以阻尼比较高的模式1 为例，在原有负荷波动2%的基础上，加入0.5%、1.0%和3.0%的负荷波动下Sub-OpMD 算法以及SSI 算法的辨识结果如表2 所示。从表2 可以看出，Sub-OpMD 算法对这4 种不同幅值负荷波动所激励的随机响应数据均具有鲁棒性，且辨识结果精度优于SSI 算法。此外，根据大量仿真实验，对于负荷波动0.5%以下的信号，本文所提算法的辨识误差明显增大，甚至不能有效辨识小干扰特征参数。

表2 不同负荷波动小干扰特征参数辨识结果对比Table 2 Comparison of identification results for small signal characteristic parameters with different load fluctuations

本文针对上文所采用的10 min 数据还验证了分析窗长和滑动窗长的选择对辨识结果的影响。首先研究分析窗长对辨识结果的影响，固定滑动窗长为0.2 s，并分别对分析窗长为0.5、1.0、2.0、3.0 min的数据进行辨识。由附录A 图A2 可知，阻尼比波动程度明显大于振荡频率，故此次验证以模式1 的阻尼比辨识结果为例，验证所得结果如图2 所示。从图2 可以看出，随着分析窗长长度的增加阻尼比辨识结果均值更接近参考值，且标准差更小。而对于滑动窗长而言，固定分析窗长为2 min，分别调整滑动窗长为0.2、0.3、0.4、0.8 s。经过仿真验证，对于相同一段10 min 数据的参数辨识，不同滑动窗长的选择对辨识结果精度影响不大。

图2 不同分析窗长下阻尼比辨识结果Fig.2 Identification results of damping ratio with different analysis window lengths

在辨识系统振荡频率和阻尼比的同时，还可以得到振荡模式的模态，即发电机组之间的振荡关系，它反映了机群发电机转子之间的非同步摇摆。Sub-OpMD 算法辨识得到了4 个区域间模式下发电机角频率之间的振荡关系，统计结果为每个分析窗口辨识结果的平均值，详细模态图见附录A 图A3。模式1中区域A、B 发电机群（G1 至G13）与区域C、D、E 发电机群（G14 至G16）相互振荡；模式2 中区域C 发电机（G14）与区域E 发电机（G16）相互振荡；模式3 中区域A 发电机群（G1 至G9）与区域B 发电机群（G10至G13）相互振荡；模式4 中区域D 发电机（G15）与区域E 发电机（G16）相互振荡。以上Sub-OpMD 算法模态辨识结果与SSSA 计算得到的发电机振荡分组一致。

4.2 实际系统算例

以某实际系统中PMU 实测数据为例，进一步验证Sub-OpMD 算法在振荡模式辨识的有效性。选取该系统中4 个PMU 采集的母线频率、母线电压幅值、电压相角数据作为输入数据，采集到的10 min内电压幅值波形如图3（a）所示，对这些采集到的数据进行功率谱分析，得到功率谱密度波形如图3（b）所示，可以看到功率谱密度在0.319 4 Hz 处有明显尖峰，表示该振荡模式在这4 个PMU 测量信号中具有很强的可观性［21-22］，本文重点对这一可观性较强的模式进行分析。

图3 实测数据波形和功率密度谱Fig.3 Waveform and power density spectrum of measured data

利用Sub-OpMD 算法对10 min 内实测数据采用滑动数据窗方法进行辨识，分析窗长为2 min，滑动窗长为0.5 s，共得到1 000 个窗口的辨识结果。时域辨识结果以及相应的频次直方图如图4 所示，从图4 可以看出，辨识结果呈正态分布，且多集中于均值附近，这在振荡频率上体现得更为明显。

图4 实测数据振荡频率和阻尼比频次Fig.4 Oscillation frequency and damping ratio frequency of measured data

对实测数据频率和阻尼比的辨识结果进行统计，如表3 所示。从表3 可以看出，所提方法辨识出的模式振荡频率均值为0.332 2 Hz，阻尼比均值为2.90%，相较于SSI 算法更接近功率谱分析得出的真实值（0.319 4 Hz），且标准差更小。该振荡模式在振荡关系上表现为PMU4 相对于PMU1 至PMU3振荡。

表3 实测数据辨识结果对比Table 3 Comparison of identification results of measured data

5 结语

本文在分析系统随机响应信号特性的基础上，将基于正交投影的矩阵线性变换与共轭梯度算法引入DMD 法，提出了随机数据驱动的Sub-OpMD 算法。基于IEEE 16 机68 节点系统和某实际测量数据的辨识结果表明，Sub-OpMD 算法可用于多通道数据，在小幅随机激励的环境下能从全局角度一次性获得全部小干扰特征参数，且辨识结果可以较准确地反映系统运行状态，有利于系统振荡模式的在线监测。通过与SSI 算法进行对比验证，表明Sub-OpMD 算法的辨识结果均值更靠近SSSA 的结果，且标准差更小。

考虑实时评估系统振荡模式参数的要求，缩短输入数据矩阵线性变换时间，提高Sub-OpMD 实时应用能力将是后续研究的重点。

附录见本刊网络版（http：//www.aeps-info.com/aeps/ch/index.aspx），扫英文摘要后二维码可以阅读网络全文。