CN104077480A - 基于Matrix Pencil的电力系统低频振荡模态辨识方法 - Google Patents

Abstract

本发明涉及电力系统低频振荡模态辨识领域，特别是一种基于Matrix Pencil的电力系统低频振荡模态辨识方法。该方法针对实际系统得到的测量数据通常受到现场环境等因素的影响，是带有一定信噪比的信号数据，提出利用旋转不变技术(ESPRIT)改进Matrix Pencil算法，直接以测量数据构成的数据矩阵为基础，将信号空间分解成信号子空间和噪声子空间，准确估计模型阶数，并检测出电力系统低频振荡信号不同振荡模态的振荡频率、衰减因子、振荡幅值和相位等信息，能够有效的提高计算效率和低频振荡辨识能力。该方法适用于电力系统等相关部门，用于电力系统低频振荡模态辨识。

Description

基于Matrix Pencil的电力系统低频振荡模态辨识方法

技术领域

本发明涉及电力系统低频振荡模态辨识领域，特别是一种基于Matrix Pencil的电力系统低频振荡模态辨识方法。

背景技术

电力系统稳定运行是电力系统根本问题。随着我国电网的联网建设，在解决我国能源分布地区与用电需求地区之间不平衡的矛盾同时，一些诸如远距离、弱电气联系等因素导致电力系统的低频振荡时有发生。电力系统低频振荡也越来越成为继自然灾害、电气设备故障之后的严重危害电网安全稳定运行的重大危害。因此有必要深入研究电力系统中低频振荡问题和抑制低频振荡的方法，实现电网的稳定运行。

研究在线的振荡特征辨识算法是实现电力系统低频振荡在线监视以及广域阻尼控制的重要理论基础。电力系统是一个复杂的大规模非线性系统，传统的模态分析方法需要建立详细的数学模型并列写方程、求解大规模矩阵的特征值。通常，由于系统过于复杂或缺乏足够准确的参数，并存在不同程度的维数灾难现象，使特征值计算困难。基于辨识的振荡模态识别方法能够从振荡信号中提取所需要的模态信息，是低频振荡的一种分析方法。现有的辨识方法有Prony算法、基于希尔伯特-黄(Hilbert-Huang，HHT)变换和小波分解等方法。传统的Prony算法利用复指数函数对信号进行强行拟合得到模态信息，会产生大量的虚假模态，同时对噪声敏感；基于HHT变换的模态识别能够提取出信号的瞬时模态信息，但是所需算法可能出现漏辨识现象。使用小波分解处理时变振荡信号有重要意义，但在提取多频率成分信号时存在分辨率不高的问题。

本发明在综合比较各方法之后，利用测量数据构造矩阵，并利用旋转不变技术(ESPRIT)改进Matrix Pencil算法，进行电力系统低频振荡模态辨识，检测出电力系统低频振荡信号不同振荡模态的振荡频率、衰减因子、振荡幅值和相位等信息。

发明内容

本发明的目的在于提供一种基于Matrix Pencil的电力系统低频振荡模态辨识方法，该方法能够快速准确的检测出低频振荡信号参数，从而实现对低频振荡模态的辨识。

为实现上述方案，本发明的技术方案是：一种基于Matrix Pencil的电力系统低频振荡模态辨识方法，利用测量数据构造矩阵，基于旋转不变技术(ESPRIT)改进Matrix Pencil算法，并利用改进之后的Matrix Pencil算法进行电力系统低频振荡模态辨识，检测出电力系统低频振荡信号不同振荡模态的振荡频率、衰减因子、振荡幅值和相位等信息，其具体步骤如下：

步骤1:设理想采样数据为x(n)，n＝0,1,…,N-1，用M阶的指数模型进行估计，如下：

$x\left(n\right)=\underset{k=1}{\overset{M}{\mathrm{\Σ}}}{b}^k{z}_k^{n-1}$

式中，z_k为包含振荡模式k的振荡频率和衰减因子信息的参数，b_k为对应振荡模式k的包含振荡幅值和初始相位信息的参数。

步骤2:根据采样数据x(0),x(1),…,x(N-1)，构造Hankel数据矩阵，如下：

$X=\left(\begin{array}{cccc}x\left(L\right)& x(L-1)& ...& x\left(0\right)\\ x(L+1)& x\left(L\right)& ...& x\left(1\right)\\ ...& ...& ...& ...\\ x(N-1)& x(N-2)& ...& x(N-L-1)\end{array}\right)$

式中，L为矩阵束参数，恰当的选择L可以抑制噪声干扰，通常取L＝N/4～N/3，假设L+1≤N-L。

步骤3:对X进行奇异值分解，得到由奇异值矩阵所组成的对角矩阵，如下：

X＝UDV^T

式中，U为主导左特征值向量矩阵，且为N-L阶正交矩阵，V^T为主导右特征值向量矩阵，且为L+1阶正交矩阵，D为(N-L)×(L+1)阶对角阵，具体表示如下：

式中，d₁,d₂,…,d_L+1为对X进行奇异值分解得到的奇异值，满足d₁≥d₂≥…≥d_L+1，对于理想的M阶信号，有如下等式：

$\left\{\begin{array}{c}{d}_1,d_2,...,d_M\≠0\\ d_{M+1}=...=d_{L+1}=0\end{array}\right.$

而对于非理想型信号，也许d_M+1,…,d_L+1不为零，但是他们的值相对于d₁,d₂,…,d_M比较小。

步骤4:设置阀值令取满足等式最大的i为模型的阶数，即M＝i。

步骤5:重新构造矩阵D′、D′为(N-L)×L阶矩阵，前M行由D的前M个奇异值组成，后N-L-M行为0，这样得到的矩阵D′可以有效的消除噪声的影响，具体表示如下：

同理，为X奇异值分解后的前M个主导右特征向量矩阵V^T的第1行～第L行，为X奇异值分解后的前M个主导右特征向量矩阵V^T的第2行～第L+1行.

步骤6:根据重新构造之后的矩阵D′、重新构造两个样本矩阵 $X_0^{\′}={\mathrm{UD}}^{\′}{V}_0^T,X_1^{\′}={\mathrm{UD}}^{\′}{V}_1^T,$ 表示如下：

$X_0^{\′}={\left(\begin{array}{cccc}x(L-1)& x(L-2)& ...& x\left(0\right)\\ x\left(L\right)& x(L-1)& ...& x\left(1\right)\\ ...& ...& ...& ...\\ x(N-2)& x(N-3)& ...& x(N-L-1)\end{array}\right)}_{(N-L)\×L}$

$X_1^{\′}={\left(\begin{array}{cccc}x\left(L\right)& x(L-1)& ...& x\left(1\right)\\ x(L+1)& x\left(L\right)& ...& x\left(2\right)\\ ...& ...& ...& ...\\ x(N-1)& x(N-2)& ...& x(N-L)\end{array}\right)}_{(N-L)\×L}$

步骤7:定义矩阵B、Z,B为包含M阶信号的所有幅值和相位信息的矩阵，Z为包含M阶信号的所有振荡频率和衰减因子等信息的矩阵，如下:

B＝diag(b₁,b₂,…,b_M)

Z＝diag(z₁,z₂,…,z_M)

根据M阶的指数模型将X₀′、X₁′与B、Z用矩阵的形式联系起来，求解得出(X₀′)^-1X₁′的M个特征值z_k(k＝1,2,….M)和Z矩阵。

步骤8:求得Z矩阵之后，根据数据时间间隔T_s，可以求得相应的衰减因子α_i和振荡频率ω_i，如下：

${\mathrm{\α}}_i+{\mathrm{j\ω}}_i=\frac{\ln Z_i}{T_s},(i=\mathrm{1,2},..,M)$

${\mathrm{\α}}_i=\mathrm{Re}\left(\frac{\ln Z_i}{T_s}\right)$

${\mathrm{\ω}}_i=\mathrm{Im}\left(\frac{\ln Z_i}{T_s}\right)$

由x＝z·b，得：

b＝z^-1·x

式中，x＝(x(0),x(1),…,x(N-1))^T为理想采样数据矩阵，z为由特征值z_k(k＝1,2,….M)组成的N×M阶范德蒙德矩阵，如下：

$z=\left(\begin{array}{cccc}1& 1& ...& 1\\ z_1& z_2& ...& z_M\\ ...& ...& ...& ...\\ z_1^{N-1}& z_2^{N-1}& ...& z_M^{N-1}\end{array}\right)$

进一步求出振荡幅值A_i和相位θ_i，如下：

$\left\{\begin{array}{c}{A}_i=\left|b_i\right|\\ {\mathrm{\θ}}_i=\arctan (\mathrm{Im}\left(b_i\right)/\mathrm{Re}\left(b_i\right))\end{array}\right.,(i=\mathrm{1,2},...M)$

至此，利用改进的Matrix Pencil算法对电力系统低频振荡信号完成模态辨识。

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

1、改进之后的算法速度快、抗噪声能力强。

2、能够准确的辨识出电力系统低频振荡的各个模态，精度得到了很大的提高。

附图说明

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

图2是本发明实施例的奇异值曲线。

图3是本发明实施例时的拟合曲线。

具体实施方式

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

请参见图1，图1是本发明实施例的工作流程图。

步骤1:设理想采样数据为x(n)，n＝0,1,…,N-1，用M阶的指数模型进行估计，如下：

$x\left(n\right)=\underset{k=1}{\overset{M}{\mathrm{\Σ}}}{b}^k{z}_k^{n-1}$

式中，z_k为包含振荡模式k的振荡频率和衰减因子信息的参数，其具体表达式可表示如下：

$z_k=e^{({\mathrm{\α}}_k+{\mathrm{j\ω}}_k)\mathrm{\Δt}}$

式中，α_k为衰减因子，ω_k为振荡频率。b_k为对应振荡模式k的包含振荡幅值和初始相位信息的参数，其具体表达式如下：

b_k＝A_ke^jθk

式中，A_k是振荡幅值，θ_k为初始相位。

步骤2:根据采样数据x(0),x(1),…,x(N-1)，构造Hankel数据矩阵，如下：

$X=\left(\begin{array}{cccc}x\left(L\right)& x(L-1)& ...& x\left(0\right)\\ x(L+1)& x\left(L\right)& ...& x\left(1\right)\\ ...& ...& ...& ...\\ x(N-1)& x(N-2)& ...& x(N-L-1)\end{array}\right)$

式中，L为矩阵束参数，恰当的选择L可以抑制噪声干扰，通常取L＝N/4～N/3，假设L+1≤N-L。

步骤3:对X进行奇异值分解，得到由奇异值矩阵所组成的对角矩阵，如下：

X＝UDV^T，图2是本实施例的奇异值曲线。

式中，U为主导左特征值向量矩阵，且为N-L阶正交矩阵，V^T为主导右特征值向量矩阵，且为L+1阶正交矩阵，D为(N-L)×(L+1)阶对角阵，具体表示如下：

式中，d₁,d₂,…,d_L+1为对X进行奇异值分解得到的奇异值，满足d₁≥d₂≥…≥d_L+1，对于理想的M阶信号，有如下等式：

$\left\{\begin{array}{c}{d}_1,d_2,...,d_M\≠0\\ d_{M+1}=...=d_{L+1}=0\end{array}\right.$

而对于非理想型信号，也许d_M+1,…,d_L+1不为零，但是他们的值相对于d₁,d₂,…,d_M比较小。

步骤4:设置阀值令取满足等式最大的i为模型的阶数，即M＝i。

步骤5:重新构造矩阵D′、D′为(N-L)×L阶矩阵，前M行由D的前M个奇异值组成，后N-L-M行为0，这样得到的矩阵D′可以有效的消除噪声的影响，具体表示如下：

同理，为X奇异值分解后的前M个主导右特征向量矩阵V^T的第1行～第L行，为X奇异值分解后的前M个主导右特征向量矩阵V^T的第2行～第L+1行.

步骤6:根据重新构造之后的矩阵D′、重新构造两个样本矩阵 $X_0^{\′}={\mathrm{UD}}^{\′}{V}_0^T,X_1^{\′}={\mathrm{UD}}^{\′}{V}_1^T,$ 表示如下：

$X_0^{\′}={\left(\begin{array}{cccc}x(L-1)& x(L-2)& ...& x\left(0\right)\\ x\left(L\right)& x(L-1)& ...& x\left(1\right)\\ ...& ...& ...& ...\\ x(N-2)& x(N-3)& ...& x(N-L-1)\end{array}\right)}_{(N-L)\×L}$

$X_1^{\′}={\left(\begin{array}{cccc}x\left(L\right)& x(L-1)& ...& x\left(1\right)\\ x(L+1)& x\left(L\right)& ...& x\left(2\right)\\ ...& ...& ...& ...\\ x(N-1)& x(N-2)& ...& x(N-L)\end{array}\right)}_{(N-L)\×L}$

步骤7:定义矩阵B、Z,B为包含M阶信号的所有幅值和相位信息的矩阵，Z为包含M阶信号的所有振荡频率和衰减因子等信息的矩阵，如下:

B＝diag(b₁,b₂,…,b_M)

Z＝diag(z₁,z₂,…,z_M)

根据M阶的指数模型将X₀′、X₁′与B、Z用矩阵的形式联系起来，有如下关系式：

$X_0^{\′}=\left(\begin{array}{cccc}{z}_1^0& z_2^0& ...& z_M^0\\ z_1^1& z_2^1& ...& z_M^1\\ ...& ...& ...& ...\\ z_1^{N-L-1}& z_2^{N-L-1}& ...& z_M^{N-L-1}\end{array}\right)\·B\·\left(\begin{array}{cccc}{z}_1^{L-1}& z_1^{L-2}& ...& z_1^0\\ z_2^{L-1}& z_2^{L-2}& ...& z_2^0\\ ...& ...& ...& ...\\ z_M^{L-1}& z_M^{L-2}& ...& z_M^0\end{array}\right)$

$X_1^{\′}=\left(\begin{array}{cccc}{z}_1^0& z_2^0& ...& z_M^0\\ z_1^1& z_2^1& ...& z_M^1\\ ...& ...& ...& ...\\ z_1^{N-L-1}& z_2^{N-L-1}& ...& z_M^{N-L-1}\end{array}\right)\·B\·Z\·\left(\begin{array}{cccc}{z}_1^{L-1}& z_1^{L-2}& ...& z_1^0\\ z_2^{L-1}& z_2^{L-2}& ...& z_2^0\\ ...& ...& ...& ...\\ z_M^{L-1}& z_M^{L-2}& ...& z_M^0\end{array}\right)$

令Z_L、Z_R为范德蒙德矩阵，如下：

$Z_L=\left(\begin{array}{cccc}{z}_1^0& z_2^0& ...& z_M^0\\ z_1^1& z_2^1& ...& z_M^1\\ ...& ...& ...& ...\\ z_1^{N-L-1}& z_2^{N-L-1}& ...& z_M^{N-L-1}\end{array}\right)$

$Z_R=\left(\begin{array}{cccc}{z}_1^{L-1}& z_1^{L-2}& ...& z_1^0\\ z_2^{L-1}& z_2^{L-2}& ...& z_2^0\\ ...& ...& ...& ...\\ z_M^{L-1}& z_M^{L-2}& ...& z_M^0\end{array}\right)$

则X₀′、X₁′可以简化为：

X₀′＝Z_L·B·Z_R

X₁′＝Z_L·B·Z·Z_R

设为Z_R的M-L阶广义逆矩阵，即令：

$Z_R^{-1}=(Q_1,Q_2,...,Q_k,...,Q_M)$

其中，Q_k为第k列L维列向量，有如下性质：

Z_R·Q_k＝(0,…,1,…,0)^T

其中，上式中第k个元素为1，其余为零，满足正交关系：

X₁′-z_kX₀′＝Z_LB(Z-z_kI_M)Z_R

Z-z_kI_M＝(z₁-z_k,…,z_k-1-z_k,0,z_k+1-z_k,…,z_M-z_k)

继续推导得：

(X₁′-z_kX₀′)Q_k＝Z_LB(Z-z_kI_M)Z_RQ_k＝[0]

进一步得到：

X₁′Q_k＝z_kX₀′Q_k

上式两边同时乘以(X₀′)^-1，得：

(X₀′)^-1X₁′Q_k＝z_k(X₀′)^-1X₀′Q_k＝z_kQ_k

至此，可得到z_k(k＝1,2,…，M)为(X₀′)^-1X₁′的M个特征值，由于(X₀′)^-1X₁′为L阶方阵，所以还存在L-M个特征值，因此由两个样本矩阵X₀′、X₁′的关系可以求解包含M阶信号的所有振荡频率和衰减因子信息的Z矩阵。

步骤8:求得Z矩阵之后，根据数据时间间隔T_s，可以求得相应的衰减因子α_i和振荡频率ω_i，如下：

${\mathrm{\α}}_i+{\mathrm{j\ω}}_i=\frac{\ln Z_i}{T_s},(i=\mathrm{1,2},..,M)$

${\mathrm{\α}}_i=\mathrm{Re}\left(\frac{\ln Z_i}{T_s}\right)$

${\mathrm{\ω}}_i=\mathrm{Im}\left(\frac{\ln Z_i}{T_s}\right)$

由x＝z·b，得：

b＝z^-1·x

式中，x＝(x(0),x(1),…,x(N-1))^T为理想采样数据矩阵，z为由特征值z_k(k＝1,2,….M)组成的N×M阶范德蒙德矩阵，如下：

$z=\left(\begin{array}{cccc}1& 1& ...& 1\\ z_1& z_2& ...& z_M\\ ...& ...& ...& ...\\ z_1^{N-1}& z_2^{N-1}& ...& z_M^{N-1}\end{array}\right)$

进一步求出振荡幅值A_i和相位θ_i，如下：

$\left\{\begin{array}{c}{A}_i=\left|b_i\right|\\ {\mathrm{\θ}}_i=\arctan (\mathrm{Im}\left(b_i\right)/\mathrm{Re}\left(b_i\right))\end{array}\right.,(i=\mathrm{1,2},...M)$

至此，利用改进的Matrix Pencil算法对电力系统低频振荡信号完成模态辨识。

请参见表1，表1是本发明实施例时的Matrix Pencil计算结果。

振荡模态	幅值	衰减因子	相角/角度	频率/Hz
					1、2	0.0250	-0.323	36.000	1.000
3、4	0.150	-0.200	0.000	0.500
					5、6	0.200	-0.162	10.588	0.250
7、8	0.500	-0.100	60.000	0.100

表1

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

Claims (1)

1.一种基于Matrix Pencil的电力系统低频振荡模态辨识方法，其特征在于：利用测量数据构造矩阵，基于旋转不变技术(ESPRIT)改进Matrix Pencil算法，并利用改进之后的Matrix Pencil算法进行电力系统低频振荡模态辨识，检测出电力系统低频振荡信号不同振荡模态的振荡频率、衰减因子、振荡幅值和相位信息，其具体步骤如下：

步骤1:设理想采样数据为x(n)，n＝0,1,…,N-1，用M阶的指数模型进行估计，如下：

$x\left(n\right)=\underset{k=1}{\overset{M}{\mathrm{\Σ}}}{b}_k{z}_k^{n-1}$

式中，z_k为包含振荡模式k的振荡频率和衰减因子信息的参数，b_k为对应振荡模式k的包含振荡幅值和初始相位信息的参数；

步骤2:根据采样数据x(0),x(1),…,x(N-1)，构造Hankel数据矩阵，如下：

$X=\left(\begin{array}{cccc}x\left(L\right)& x(L-1)& ...& x\left(0\right)\\ x(L+1)& x\left(L\right)& ...& x\left(1\right)\\ ...& ...& ...& ...\\ x(N-1)& x(N-2)& ...& x(N-L-1)\end{array}\right)$

式中，L为矩阵束参数，恰当的选择L可以抑制噪声干扰，通常取L＝N/4～N/3，假设L+1≤N-L；

步骤3:对X进行奇异值分解，得到由奇异值矩阵所组成的对角矩阵，如下：

X＝UDV^T

式中，U为主导左特征值向量矩阵，且为N-L阶正交矩阵，V^T为主导右特征值向量矩阵，且为L+1阶正交矩阵，D为(N-L)×(L+1)阶对角阵，具体表示如下：

式中，d₁,d₂,…,d_L+1为对X进行奇异值分解得到的奇异值，满足d₁≥d₂≥…≥d_L+1，对于理想的M阶信号，有如下等式：

$\left\{\begin{array}{c}{d}_1,d_2,...,d_M\≠0\\ d_{M+1}=...=d_{L+1}=0\end{array}\right.;$

步骤4:设置阀值令取满足等式最大的i为模型的阶数，即M＝i；

步骤5:重新构造矩阵D′、D′为(N-L)×L阶矩阵，前M行由D的前M个奇异值组成，后N-L-M行为0，这样得到的矩阵D′可以有效的消除噪声的影响，具体表示如下：

同理，为X奇异值分解后的前M个主导右特征向量矩阵V^T的第1行～第L行，为X奇异值分解后的前M个主导右特征向量矩阵V^T的第2行～第L+1行；

步骤6:根据重新构造之后的矩阵D′、重新构造两个样本矩阵 $X_0^{\′}={\mathrm{UD}}^{\′}{V}_0^T,X_1^{\′}={\mathrm{UD}}^{\′}{V}_1^T,$ 表示如下：

$X_0^{\′}={\left(\begin{array}{cccc}x(L-1)& x(L-2)& ...& x\left(0\right)\\ x\left(L\right)& x(L-1)& ...& x\left(1\right)\\ ...& ...& ...& ...\\ x(N-2)& x(N-3)& ...& x(N-L-1)\end{array}\right)}_{(N-L)\×L}$

$X_1^{\′}={\left(\begin{array}{cccc}x\left(L\right)& x(L-1)& ...& x\left(1\right)\\ x(L+1)& x\left(L\right)& ...& x\left(2\right)\\ ...& ...& ...& ...\\ x(N-1)& x(N-2)& ...& x(N-L)\end{array}\right)}_{(N-L)\×L};$

步骤7:定义矩阵B、Z,B为包含M阶信号的所有幅值和相位信息的矩阵，Z为包含M阶信号的所有振荡频率和衰减因子信息的矩阵，如下:

B＝diag(b₁,b₂,…,b_M)

Z＝diag(z₁,z₂,…,z_M)

根据M阶的指数模型将X₀′、X₁′与B、Z用矩阵的形式联系起来，求解得出(X₀′)^-1X₁′的M个特征值z_k(k＝1,2,….M)和Z矩阵；

步骤8:求得Z矩阵之后，根据数据时间间隔T_s，求得相应的衰减因子α_i和振荡频率ω_i，如下：

${\mathrm{\α}}_i+{\mathrm{j\ω}}_i=\frac{\ln Z_i}{T_s},(i=\mathrm{1,2},..,M)$

${\mathrm{\α}}_i=\mathrm{Re}\left(\frac{\ln Z_i}{T_s}\right)$

${\mathrm{\ω}}_i=\mathrm{Im}\left(\frac{\ln Z_i}{T_s}\right)$

由x＝z·b，得：

b＝z^-1·x

式中，x＝(x(0),x(1),…,x(N-1))^T为理想采样数据矩阵，z为由特征值z_k(k＝1,2,….M)组成的N×M阶范德蒙德矩阵，如下：

$z=\left(\begin{array}{cccc}1& 1& ...& 1\\ z_1& z_2& ...& z_M\\ ...& ...& ...& ...\\ z_1^{N-1}& z_2^{N-1}& ...& z_M^{N-1}\end{array}\right)$

进一步求出振荡幅值A_i和相位θ_i，如下：

$\left\{\begin{array}{c}{A}_i=\left|b_i\right|\\ {\mathrm{\θ}}_i=\arctan (\mathrm{Im}\left(b_i\right)/\mathrm{Re}\left(b_i\right))\end{array}\right.,(i=\mathrm{1,2},...M).$