CN106786514A - 一种电力系统低频振荡模式在线辨识方法 - Google Patents

Abstract

本发明公开了一种电力系统低频振荡模式在线辨识方法，属于电力系统低频振荡领域。所述方法包括：首先读取一段电力系统扰动后多台发电机的角速度量测信号；然后由所述自由振荡信号构造Hankel矩阵，然后对上述Hankel矩阵计算协方差序列构成T矩阵；然后对T矩阵进行奇异值分解，并构造旋转算子；最后通过计算旋转算子的特征值，从而得到系统的振荡频率、阻尼比。本发明采用协方差矩阵代替原始信号的Hankel矩阵，使其抗噪性增强，计算速度更快，计算结果更稳定，能精确地辨识低频振荡模式的频率、阻尼比，具有更高的工程实用价值。

Description

一种电力系统低频振荡模式在线辨识方法

技术领域

本发明涉及电力系统低频振荡技术领域，具体为一种电力系统低频振荡模式在线辨识方法。

背景技术

随着电力系统互联程度提高、负荷的增长、发电机快速励磁系统的大量使用，低频振荡问题益发突出，已成为威胁电力系统稳定和限制交流互联电网传输能力的瓶颈。电力系统低频振荡分析方法可分为基于系统模型的方法和基于量测信号进行辨识的方法。基于模型的方法具有计算量较大、不能及时跟踪系统参数和结构的变化的缺点，如WECC系统标准模型在北美大停电事故发生前并未发现不稳定的低频振荡模式。由于实测数据真实体现了系统当前的运行状态，因此基于量测的低频振荡分析方法弥补基于模型的分析方法的不足，具有广阔的应用前景。

近些年来，基于量测数据辨识低频振荡的方法大量涌现。在基于时域信号的方法中，Prony方法被广泛应用于基于自由振荡信号的低频振荡辨识，但Prony方法对噪声敏感，模型的阶数对结果影响也较大。TLS-ESPRIT方法也常用于电力系统低频振荡模式辨识中，但由于TLS-ESPRIT方法需要对矩阵进行两次奇异值分解，因此计算速度较慢。另外，还有一些典型的基于频域信号的方法，主要包括基于傅里叶变换的方法，基于小波变换的方法和基于希尔伯特黄变换(HHT)的方法。傅里叶变换只能给出信号的频率信息，基于Morlet小波提取小波时频分布脊线，辨识不同时段信号的低频振荡参数；基于希尔伯特黄变换(HHT)的方法，其在电力系统低频振荡模式识别中的应用也比较广泛，但是，其采用EMD没有坚实的理论基础，所得到的模态函数仍需进一步研究。也有卡尔曼滤波的方法实现电力系统低频振荡辨识，但该方法需要构建系统传递函数以求得状态空间矩阵，依赖于系统模型的有效程度。

电力系统的量测、传输环节常常会引入干扰噪声，主要来源于系统负荷的随机波动，通常将这些噪声看成高斯白噪声。高斯白噪声的存在，会影响低频振荡的辨识精度，而目前的方法大多通过奇异值分解进行去噪，这种方法对噪声的处理能力有限。所以，现有的方法存在抗噪性弱、计算速度慢等问题。

发明内容

针对上述问题本发明的目的在于提供一种抗噪性强，计算稳定，计算速率快，能精确地辨识低频振荡模式的频率、阻尼比的电力系统低频振荡模式在线辨识方法，技术方案如下：

一种电力系统低频振荡模式在线辨识方法，包括以下步骤：

步骤1：读取一段电力系统扰动后的发电机角速度或联络线功率的自由振荡信号y：

y＝[y(0),y(1),…y(i)…,y(r)]

其中，y(i)表示第i个测量数据，i＝0,1,…,r，r为所取信号长度；

步骤2：利用上述自由振荡信号构造扩展Hankel矩阵Y：

其中，j＝[r/2]，此处[]表示向下取整，且2k+j-2＝r；

步骤3：对上述Hankel矩阵Y计算协方差序列构成Τ矩阵；

步骤4：对上述Τ矩阵进行奇异值分解，根据奇异值大小确定有效阶数p，生成信号子空间V_s和噪声子空间V_n；

步骤5：构造矩阵[V1V2]，并对其进行奇异值分解：

其中，V₁为去掉V_s最后一行得到的矩阵；V₂为去掉V_s第一行得到的矩阵；

步骤6：把分解为4个p×p的矩阵：

步骤7：计算旋转算子的特征值λ_i(i＝1,2,……,p)，并由此计算振荡频率f_i和阻尼比ζ_i：

进一步的，所述步骤3中对上述Hankel矩阵Y计算协方差序列构成T矩阵的具体方法为：

3.更进一步的，所述步骤4具体包括：

a)对上述T矩阵按下式进行奇异值分解：

式中，Σ为对角矩阵，对角线元素为T的奇异值ξ_i；

b)确定信号子空间阶数p：对角阵Σ中的元素，找出满足下式的最小的整数i，取信号子空间的阶数p＝i；

c)按下式生成信号子空间V_s和噪声子空间V_n：

式中，Σ_S为矩阵X幅值最大的p个奇异值组成的对角阵，Σ_n为矩阵X余下的奇异值组成的对角阵；酉矩阵V按奇异值的大小划分为信号子空间V_s和噪声子空间V_n，V_s的列向量是对应矩阵X的幅值最大的p个奇异值的奇异向量，酉矩阵U按奇异值大小划分为U_s和U_n。

本发明的有益效果是：本发明采用协方差矩阵代替原始信号的Hankel矩阵，使其抗噪性增强，计算速度更快，计算结果更稳定，能精确地辨识低频振荡模式的频率、阻尼比，具有更高的工程实用价值。

附图说明

图1为本发明实施例中一种电力系统低频振荡模式在线辨识方法的流程图。

图2为应用实例中16机系统结构图。

图3为故障后16台发电机角速度振荡曲线。

图4为不同量测噪声水平下本发明方法的辨识结果。

具体实施方式

本发明一种电力系统低频振荡模式在线辨识方法，只需要利用所量测到的系统的时域响应数据便能实现模态分析功能，包括振荡的频率、阻尼比，该方法辨识精度高，抗噪性强。

下面结合附图和具体实施方式对本发明作进一步详细的说明。在实施例中，提供了一种电力系统低频振荡模式在线辨识方法，参考图1，所述方法包括以下步骤：

步骤1：读取一段电力系统扰动后的发电机角速度或联络线功率自由振荡信号y：

y＝[y(0),y(1),…y(i)…,y(r)]

其中，y(i)表示第i个测量数据，i＝0,1,…,r，r为所取信号长度。

步骤2：利用自由振荡信号构造扩展Hankel矩阵Y

利用步骤1中所述自由振荡信号x中的数据构造扩展Hankel矩阵Y：

其中，j＝[r/2]，此处[]表示向下取整，2k+j-2＝r。

步骤3：对上述Hankel矩阵计算协方差序列构成Τ矩阵：

步骤4：对上述Τ矩阵进行奇异值分解，根据奇异值大小确定有效阶数p，生成信号子空间V_s和噪声子空间V_n：

a)对上述T矩阵按下式进行奇异值分解：

式中，Σ为对角矩阵，对角线元素为T的奇异值ξ_i。

b)确定信号子空间阶数p：对角阵Σ中的元素，找出满足下式的最小的整数i，取信号子空间的阶数p＝i。

在方法的应用中，K_c的值可取为0.01。

c)按下式生成信号子空间V_s和噪声子空间V_n

式中，Σ_S为矩阵X幅值最大的p个奇异值组成的对角阵，Σ_n为矩阵X余下的奇异值组成的对角阵。按奇异值的大小划分为信号子空间V_s和噪声子空间V_n，V_s的列向量是对应矩阵X的幅值最大的p个奇异值的奇异向量，相应地，酉矩阵U按奇异值大小可划分为U_s,U_n。

步骤5：构造矩阵[V1V2]，并对其进行奇异值分解

其中，V₁为去掉V_s最后一行得到的矩阵；V₂为去掉V_s第一行得到的矩阵。

步骤6：把分解为4个p×p的矩阵：

步骤7：计算旋转算子的特征值λ_i(i＝1,2,……,p)，并由此计算振荡频率f_i和阻尼比ζ_i：

应用实例：选用16机68节点仿真系统作为算例对本实施例中的方案进行验证，所述16机仿真系统为研究区域间低频振荡问题的经典系统，其结构图详见图2。

16机系统划分为5个区域：区域1包含发电机G1～G9，区域2包含发电机G10～G13，发电机G14、G15、G16分别在区域3、区域4、区域5中。在MATLAB提供的PST(Power Systemtoolbox)中搭建该16机系统模型，并求解系统线性化后系统的状态矩阵的特征值，知系统中存在4个区域间振荡模式，其低频振荡模式的频率和阻尼比如表1所示。

表1 16机系统低频振荡模式真实值

模式	频率/Hz	阻尼比/％
			1	0.3763	11.43
2	0.5214	0.85
			3	0.6497	13.87
4	0.7928	3.56

通过对系统线性化后的状态矩阵的特征值分析可知，模式1主要表现为区域1-2中的发电机相对于区域3-5中的发电机振荡，模式2主要表现为区域1-4中的发电机相对于区域5中的发电机振荡，模式3主要表现为区域1中的发电机相对于区域2中的发电机振荡，模式4主要表现为区域3和区域5中的发电机相对于区域4中的发电机振荡。

本算例扰动设置如下：0.1s时系统1-27输电线路发生3相短路故障(图2中加粗部分)。故障后16台发电机角速度的振荡曲线如图3所示。

表2所示为无噪声情况下改进TLS-ESPRIT对频率和阻尼比的辨识结果。

表2无噪声时的辨识结果

从表2以看出，本发明方法计算的4种低频振荡模式的振荡频率和阻尼比与真实值都很接近，4种模式下的频率和阻尼比的误差都小于1％，表明改进TLS-ESPRIT能非常准确地辨识四个低频振荡模式的频率和阻尼比。

在实测PMU的数据中，往往含有量测噪声，因此，本发明通过向得到的仿真数据中叠加不同分贝高斯白噪声的方式来验证本发明方法的抗噪性能。为了排除偶然因素的影响，在不同噪声水平下均采用蒙特卡洛思路，进行100次试验并记录每次的辨识结果。

图4给出了在不同噪声水平下，100次蒙特卡洛仿真中改进TLS-ESPRIT对4个低频振荡模式的辨识结果和真实值的对比结果。可以看出，在存在量测噪声时，改进TLS-ESPRIT仍能准确地辨识4种振荡模式的频率和阻尼比。

表3给出了在SNR＝20dB时在100次蒙特卡洛仿真中本发明方法辨识结果的统计数据。从表3中以看出，当信噪比SNR＝20dB时，本发明方法计算的频率和阻尼比的均值误差和标准差都很小，对四种低频振荡模式的阻尼比辨识效果都较为准确。

表3信噪比SNR＝30dB时的辨识结果

以最为关注的、对系统影响最大的若阻尼模式2为研究对象，采用本发明方法与传统TLS-ESPRIT方法在相同量测噪声水平下进行辨识，表4和表5分别给出了两种方法对频率和阻尼比的辨识结果。

表4：两种方法在不同噪声水平下对模式2频率辨识结果

表5：两种方法在不同噪声水平下对模式2阻尼比辨识结果

由表4和表5可以看出，在不同噪声水平下，基于本发明方法计算的阻尼比标准差和均值误差明显小于传统TLS-ESPRIT计算结果，随着信号信噪比降低，基于本发明方法计算的阻尼比标准差和均值误差基本不变，而传统TLS-ESPRIT计算的阻尼比标准差和均值误差变化明显，表明本发明方法具有较强的抗噪性。以相同的方法分析模式1、模式3和模式4，该结论仍然成立。

表6本发明与TLS-ESPRIT耗时对比

方法	本发明方法	传统TLS-ESPRIT
			一次辨识平均耗时(s)	8.316	23.532

表6给出了在信噪比为20dB时100次蒙特卡洛试验中本发明方法和传统TLS-ESPRIT方法的计算速度对比，从记录的结果以看到，相比于传统TLS-ESPRIT方法，本发明方法在计算效率方面具有更大的优势。

Claims (3)

1.一种电力系统低频振荡模式在线辨识方法，其特征在于，包括以下步骤：

步骤1：读取一段电力系统扰动后的发电机角速度或联络线功率的自由振荡信号y：

y＝[y(0),y(1),…y(i)…,y(r)]

其中，y(i)表示第i个测量数据，i＝0,1,…,r，r为所取信号长度；

步骤2：利用上述自由振荡信号构造扩展Hankel矩阵Y：

其中，j＝[r/2]，此处[]表示向下取整，且2k+j-2＝r；

步骤3：对上述Hankel矩阵Y计算协方差序列构成Τ矩阵；

步骤4：对上述Τ矩阵进行奇异值分解，根据奇异值大小确定有效阶数p，生成信号子空间V_s和噪声子空间V_n；

步骤5：构造矩阵[V1V2]，并对其进行奇异值分解：

$\[V_1,V_2\]\stackrel{SVD}{=}\stackrel{\‾}{U}\·\stackrel{\‾}{\Σ}\·{\stackrel{\‾}{V}}^H$

其中，V₁为去掉V_s最后一行得到的矩阵；V₂为去掉V_s第一行得到的矩阵；

步骤6：把分解为4个p×p的矩阵：

$\stackrel{\‾}{V}=\left[\begin{array}{cc}{\stackrel{\‾}{V}}_{11}& {\stackrel{\‾}{V}}_{12}\\ {\stackrel{\‾}{V}}_{21}& {\stackrel{\‾}{V}}_{22}\end{array}\right]$

步骤7：计算旋转算子的特征值λ_i(i＝1,2,……,p)，并由此计算振荡频率f_i和阻尼比ζ_i：

$\left\{\begin{array}{c}{f}_i=\frac{\left|\arctan \left(\mathrm{Im}\right({\mathrm{\λ}}_i)/\mathrm{Re}({\mathrm{\λ}}_i\left)\right)\right|}{2\mathrm{\π}\mathrm{\Δ}t}\\ {\mathrm{\ξ}}_i=-\frac{\mathrm{Re}\left(\ln \right({\mathrm{\λ}}_i\left)\right)}{\sqrt{{\left(\mathrm{Re}\right(\ln \left({\mathrm{\λ}}_i\right)\left)\right)}^2/{\left(\mathrm{Im}\right(\ln \left({\mathrm{\λ}}_i\right)\left)\right)}^2}}\end{array}\right..$

2.根据权利要求1所述的电力系统低频振荡模式在线辨识方法，其特征在于，所述步骤3中对上述Hankel矩阵Y计算协方差序列构成T矩阵的具体方法为：

3.根据权利要求1所述的电力系统低频振荡模式在线辨识方法，其特征在于，所述步骤4具体包括：

a)对上述T矩阵按下式进行奇异值分解：

$T\stackrel{SVD}{=}{\mathrm{U\ΣV}}^H$

式中，Σ为对角矩阵，对角线元素为T的奇异值ξ_i；

b)确定信号子空间阶数p：

对角阵Σ中的元素，找出满足下式的最小的整数i，取信号子空间的阶数p＝i；

$\left\{\begin{array}{c}\frac{{\mathrm{\ξ}}_i}{{\mathrm{\ξ}}_1}>K_c\\ \frac{{\mathrm{\ξ}}_i}{{\mathrm{\ξ}}_1}\≤K_c\end{array}\right.;$

c)按下式生成信号子空间V_s和噪声子空间V_n：

$X={\mathrm{U\ΣV}}^H=\[U_s,U_n\]\left[\begin{array}{cc}{\mathrm{\Σ}}_s& 0\\ 0& {\mathrm{\Σ}}_n\end{array}\right]\left[\begin{array}{c}{V}_s^H\\ V_n^H\end{array}\right]$

式中，Σ_S为矩阵X幅值最大的p个奇异值组成的对角阵，Σ_n为矩阵X余下的奇异值组成的对角阵；酉矩阵V按奇异值的大小划分为信号子空间V_s和噪声子空间V_n，V_s的列向量是对应矩阵X的幅值最大的p个奇异值的奇异向量，酉矩阵U按奇异值大小划分为U_s和U_n。