CN104993480B - 基于递推随机子空间的电力系统低频振荡在线辨识方法 - Google Patents

Abstract

本发明涉及一种基于递推随机子空间的电力系统低频振荡在线辨识方法。本发明方法针对随机子空间辨识算法中要进行算法复杂度高的SVD分解，从而造成低频振荡模态辨识的实时性和动态辨识效果差这一问题，提出引入遗忘因子对Hankel协方差矩阵进行更新，利用投影近似子空间跟踪方法对子空间进行递推，避免SVD运算，显著减少计算复杂度，在每一次递归计算中，本发明算法的复杂度为远小于SVD计算的复杂度，能够有效提高辨识的实时性，适合低频振荡模态的在线辨识，也可以为多时空尺度下电力系统的在线监测和稳定性分析提供有效支持。

Description

基于递推随机子空间的电力系统低频振荡在线辨识方法

技术领域

本发明涉及电力系统低频振荡分析技术领域，特别是一种基于递推随机子空间的电力系统低频振荡在线辨识方法。

背景技术

随着大区电力系统互联以及大规模风、光发电系统并网运行且在电网装机容量中比重不断提升，电网弱阻尼低频振荡问题日益凸显，严重影响电网的稳定运行和电能质量控制。实时、准确辨识振荡各模态信息，定位振荡区域和机组，进而施加有效低频振荡抑制措施就显得极为关键。当前，广域测量系统(WAMS)的广泛应用，特别相量测量单元(PMU)采样频率和精度的不断提高，为低频振荡在线辨识提供了支持。

近年来，多种基于WAMS量测数据的系统模态辨识方法在电力系统低频振荡在线辨识中得到广泛应用和发展。Prony及其改进算法可准确提取主导振荡模态，获得低频振荡详细参数信息，但算法对噪声敏感，这一缺点使得算法定阶困难，计算过程极易产生虚假特征根或漏根，不利于工程运用。ARMA算法是一种基于线性模型的信号处理方法，具有一定的抗噪性能，但是在低频振荡模态辨识中，无法直接辨识振型，与Prony算法一样存在着模型定阶困难问题。希尔伯特-黄变换(HHT)适用于非线性、非平稳信号处理和辨识，但其理论框架尚需进一步完善，算法本身存在着模态混叠、边界处理等问题，在噪声干扰下经验模态(EMD)筛选结果也不尽人意。随机子空间辨识算法在振型辨识方面具有优势，是系统辨识领域最成功的算法之一。该算法直接以WAMS量测数据构成数据矩阵，利用信号子空间和噪声子空间的正交性，通过奇异值分解将信号空间分解为信号子空间和噪声子空间，具有很高的抗噪性能，是一种高精度的系统辨识方法。但该算法在每次辨识中都需要进行高阶数据矩阵奇异值分解(SVD)，算法费时间、冗余大，在电力系统低频振荡模态辨识应用中，不仅振荡模态动态辨识效果差，而且在线辨识工程实现困难。因此需要寻求能够显著减少算法复杂度、提高辨识实时性的方法，能够实现电力系统低频振荡模态信息的准确在线辨识，也可以为多时空尺度下电力系统的稳定性分析提供了有效支持。

发明内容

本发明的目的在于提供一种基于递推随机子空间的电力系统低频振荡在线辨识方法，具有良好的抗干扰能力和实时性，从而实现对电力系统低频振荡模态信息的在线辨识，为电力系统实时监测和稳定性分析提供有效支持。

为实现上述目的，本发明的技术方案是：一种基于递推随机子空间的电力系统低频振荡在线辨识方法，利用随机子空间方法辨识电力系统低频振荡模态信息，并引入遗忘因子对Hankel协方差矩阵更新，基于投影近似子空间跟踪方法对子空间进行递推，避免SVD运算，降低了算法复杂度、提升辨识实时性，包括如下步骤，

步骤1：由有限长度的振荡电气量实测数据构造Hankel协方差矩阵：

式中，为振荡电气量实测时间序列数据，，，为数据协方差，为数据窗的长度；

步骤2：对矩阵H做一次SVD分解计算：

将奇异值从大到小排列：

令，，和分别为左、右奇异向量，、和、分别为相应信号子空间和噪声子空间的左右奇异向量，是包含全部奇异值的对角矩阵，包含个主奇异值，为张成的信号子空间，为张成的噪声子空间；

步骤3：进行以下计算，其结果作为步骤4随机子空间迭代递归运算的初始值：

其中，符号表示广义逆；为权值矩阵，为协方差矩阵估计；

步骤4：利用时间序列数据，广义可观测矩阵的计算更新可由以下基于投影近似子空间跟踪方法的递归计算实现：

其中，遗忘因子可取值0.95-0.99；

步骤5:依据系统的最小实现理论，系统的广义可观测矩阵可表示为：

其中，C为线性状态空间输出矩阵，为系统矩阵，矩阵由下式求出：

上式中，和分别表示矩阵去除最后一行和第一行后的矩阵；

步骤6:求取矩阵的特征值，设采样时间间隔为，可辨识出各振荡模态的频率衰减因子和阻尼比为：

步骤7:随机子空间模态辨识方法由下式可以直接确定振荡模态振型：

上式中，和分别表示算法模型振型和系统振型，也是矩阵特征值对应的右特征向量；

步骤8:通过最小二乘法求取各模态的幅值和相角；对于个采样数据构造向量，模态幅值和相角辨识通过求解以下线性系统方程：

其中，，最小二乘解为：

由此，各振荡模态幅值和相位为：

步骤9：通过时间序列数据的更新，进行步骤4-8递推计算实现振荡模态信息的在线辨识。

相较于现有技术，本发明具有以下有益效果：

1、利用信号子空间和噪声子空间的正交性，递推随机子空间辨识方法能有效辨识出噪声环境下电力系统低频振荡主导模态：幅值、频率、相位、阻尼比、衰减因子、振型，具有很好的抗噪性能和辨识精度高；

2、引入遗忘因子对Hankel协方差矩阵进行更新，利用投影近似子空间跟踪方法对子空间进行递推，避免SVD运算，显著减少计算复杂度，在每一次递归计算中，本发明算法的复杂度为远小于SVD计算的复杂度，能够有效提高辨识的实时性，适合低频振荡模态的在线辨识，也可以为多时空尺度下电力系统的稳定性分析提供了有效支持。

附图说明

图1为本发明实施例的工作流程图。

图2是本发明方法对振荡模态辨识的效果图。

具体实施方式

下面结合附图，对本发明的技术方案进行具体说明。

本发明的一种基于递推随机子空间的电力系统低频振荡在线辨识方法，具体即电力系统低频振荡模态信息辨识方法方法结合图1进行说明并对电力系统低频振荡理想信号上叠加信噪比为10dB的白噪声，数据采样时间0-20秒，采样频率20Hz。采用本发明方法对此振荡模态的主要信息进行辨识。其中选择5s采样初始化数据构造Hankel矩阵，窗数据为0.1s。具体步骤如下：

步骤1:由采样离散时间序列数据构造协方差Hankel矩阵:

式中，为振荡电气量实测时间序列数据，，，为数据协方差，为数据窗的长度。

步骤2:对矩阵H做一次SVD分解计算：

将奇异值从大到小排列：

令，，和分别为左、右奇异向量，、和、分别为相应信号子空间和噪声子空间的左右奇异向量，是包含全部奇异值的对角矩阵，包含个主奇异值，为张成的信号子空间，为张成的噪声子空间。

步骤3：进行以下计算，其结果作为步骤4随机子空间迭代递归运算的初始值，其中符号“”表示广义逆。

其中， 为权值矩阵，为协方差矩阵估计；

步骤4：利用时间序列数据，广义可观测矩阵的计算更新可由以下基于投影近似子空间跟踪方法的递归计算实现：

其中遗忘因子可取值0.95-0.99。

步骤5:依据系统的最小实现理论，系统的广义可观测矩阵可表示为：

系统矩阵由下式求出：

上式中，和分别表示矩阵去除最后一行和第一行后的矩阵；符号“”表示广义逆。

步骤6: 求取矩阵的特征值，设采样时间间隔为，可辨识出各振荡模态的频率衰减因子和阻尼比为：

步骤7:随机子空间模态辨识方法由式(9-10)可以直接确定振荡模态振型：

上式中，和分别表示算法模型振型和系统振型，也是矩阵特征值对应的右特征向量。

步骤8:通过最小二乘法求取各模态的幅值和相角。对于个采样数据构造向量，模态幅值和相角辨识通过求解以下线性系统方程。

其中，最小二乘解为：

由此，各振荡模态幅值和相位为：

步骤9：通过时间序列数据的更新，通过步骤4-8递推计算实现了振荡模态信息的在线辨识。

为了让本领域人员更容易理解本发明，利用本发明提出的方法对下列振荡理想信号进行辨识，并叠加信噪比为10dB的白噪声，数据采样时间0-20秒，采样频率20Hz。图1给出了本发明所提出方法的流程图，图2给出了振荡理想信号模态信息辨识的仿真图。

以上是本发明的较佳实施例，凡依本发明技术方案所作的改变，所产生的功能作用未超出本发明技术方案的范围时，均属于本发明的保护范围。

Claims (1)

1.一种基于递推随机子空间的电力系统低频振荡在线辨识方法，利用随机子空间方法辨识电力系统低频振荡模态信息，并引入遗忘因子对Hankel协方差矩阵更新，基于投影近似子空间跟踪方法对子空间进行递推，避免SVD运算，降低了算法复杂度、提升辨识实时性，其特征在于：包括如下步骤，

步骤1：由有限长度的振荡电气量实测数据构造Hankel协方差矩阵：

$H=\left[\begin{array}{cccc}{R}_0& R_1& ...& R_I\\ R_1& R_2& ...& R_{I+1}\\ ...& ...& ...& ...\\ R_I& R_{I+1}& ...& R_{2I}\end{array}\right]=\underset{k}{\Σ}{y}_k^{+}{y}_k^{-T}$

式中，y_i为振荡电气量实测时间序列数据， R_i为数据协方差，2×I为数据窗的长度；

步骤2：对矩阵H做一次SVD分解计算：

$H={\mathrm{USV}}^T=\left[\begin{array}{cc}{U}_1& U_2\end{array}\right]\left[\begin{array}{cc}{S}_1& 0\\ 0& S_2\end{array}\right]\left[\begin{array}{c}{V_1}^T\\ {V_2}^T\end{array}\right]$

将奇异值从大到小排列：λ₁≥λ₂≥…≥λ_n>λ_n+1＝…λ_i+I＝σ²

令S₁＝diag(λ₁,λ₁,…,λ_n)，S₂＝diag(λ_n+1,λ_n+2,…λ_I+1)，U和V分别为左、右奇异向量，U₁、V₁和U₂、V₂分别为相应信号子空间和噪声子空间的左右奇异向量，S是包含全部奇异值的对角矩阵，S₁∈R^n×n包含n个主奇异值，span(S₁)为张成的信号子空间，span(S₂)⊥span(S₁)为张成的噪声子空间；

步骤3：进行以下计算，其结果作为步骤4随机子空间迭代递归运算的初始值：

其中，符号表示广义逆；W为权值矩阵，ξ为协方差矩阵估计；

步骤4：利用时间序列数据广义可观测矩阵O的计算更新可由以下基于投影近似子空间跟踪方法的递归计算实现：

$w\left(t\right)={\mathrm{\ξ}}_{t-1}{y}_t^{-}\∈R^{n\×1}$

$h\left(t\right)=W^T(t-1)y_t^{+}\∈R^{n\×1}$

Φ(t)＝[w(t)h(t)]∈R^n×2

${\mathrm{\μ}}^2\mathrm{\Λ}\left(t\right)=\left[\begin{array}{cc}-y_t^{-T}{y}_t^{-}& \mathrm{\μ}\\ \mathrm{\μ}& 0\end{array}\right]\∈R^{2\×2}$

K(t)＝[μ²Λ(t)+Φ^H(t)P(t-1)Φ(t)]^-1×Φ^H(t)P(t-1)∈R^2×n

$v\left(t\right)=\left[\begin{array}{cc}{\hat{H}}_{t-1}{y}_T^{-}& y_t^{+}\end{array}\right]\∈R^{lI\×n}$

$P\left(t\right)=\frac{1}{{\mathrm{\μ}}^2}\[P(t-1)-P(t-1)\mathrm{\Φ}\left(t\right)K\left(t\right)\]\∈R^{n\×n}$

$\mathrm{\ξ}\left(t\right)=\mathrm{\μ}\mathrm{\ξ}(t-1)+h\left(t\right)y_t^{-T}\∈R^{n\×lI}$

$H_t={\mathrm{\μH}}_{t-1}+y_t^{+}{y}_T^{-t}\∈R^{lI\×lI}$

W(t)＝W(t-1)+[v(t)-W(t-1)Φ(t)]×∈R^lI×nK(t)

O_t＝W(t)

其中，遗忘因子μ取值0.95-0.99；

步骤5:依据系统的最小实现理论，系统的广义可观测矩阵O可表示为：

$O=\left[\begin{array}{c}C\\ CA\\ \mathrm{...}\\ {\mathrm{CA}}^n\end{array}\right]=U_1{\left(S_1\right)}^{1/2}$

其中，C为线性状态空间输出矩阵，A为系统矩阵，矩阵A由下式求出：

上式中，O和分别表示矩阵O去除最后一行和第一行后的矩阵；

步骤6:求取矩阵A的特征值λ_i(i＝1,2,…,n)，设采样时间间隔为△t，可辨识出各振荡模态的频率f_i衰减因子α_i和阻尼比ζ_i为：

$\left.\begin{array}{c}{f}_i={\mathrm{arg\λ}}_i/\left(2\mathrm{\π}\mathrm{\Δ}t\right)\\ {\mathrm{\α}}_i=\ln \left|{\mathrm{\λ}}_i\right|/\mathrm{\Δ}t\\ {\mathrm{\ζ}}_i=-{\mathrm{\α}}_i/\sqrt{{{\mathrm{\α}}_i}^2+{\left(2{\mathrm{\πf}}_i\right)}^2}\end{array}\right\},i=1,2,\mathrm{...}n$

步骤7:随机子空间模态辨识方法由下式可以直接确定振荡模态振型：

$\left\{\begin{array}{c}C=O(1:l,:)\\ {\mathrm{\φ}}_i^s=C\×{\mathrm{\φ}}_i^m\end{array}\right.$

上式中，φ_i ^m和φ_i ^s分别表示算法模型振型和系统振型，φ_i ^m也是矩阵A特征值λ_i对应的右特征向量；

步骤8:通过最小二乘法求取各模态的幅值和相角；对于N个采样数据构造向量Y＝[y₀y₁…y_N-1]^T，模态幅值和相角辨识通过求解以下线性系统方程：

Y＝Λc

$\mathrm{\Λ}=\left[\begin{array}{cccc}1& 1& \mathrm{...}& 1\\ {\mathrm{\λ}}_1& {\mathrm{\λ}}_2& \mathrm{...}& {\mathrm{\λ}}_3\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ {{\mathrm{\λ}}_1}^{N-1}& {{\mathrm{\λ}}_2}^{N-1}& \mathrm{...}& {{\mathrm{\λ}}_3}^{N-1}\end{array}\right]$

其中，c＝[c₁ c₂ … c_n]^T，最小二乘解为：c＝(Λ^HΛ)^-1Λ^HY

由此，各振荡模态幅值A_i和相位为：

步骤9：通过时间序列数据的更新，进行步骤4-8递推计算实现振荡模态信息的在线辨识。