CN105740209B - 一种Givens迭代的Prony低频振荡分析方法 - Google Patents

一种Givens迭代的Prony低频振荡分析方法 Download PDF

Info

Publication number
CN105740209B
CN105740209B CN201610065625.5A CN201610065625A CN105740209B CN 105740209 B CN105740209 B CN 105740209B CN 201610065625 A CN201610065625 A CN 201610065625A CN 105740209 B CN105740209 B CN 105740209B
Authority
CN
China
Prior art keywords
formula
givens
equation
coefficient
signal model
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Expired - Fee Related
Application number
CN201610065625.5A
Other languages
English (en)
Other versions
CN105740209A (zh
Inventor
姚玉斌
李崇
吴志良
王丹
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Dalian Maritime University
Original Assignee
Dalian Maritime University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Dalian Maritime University filed Critical Dalian Maritime University
Priority to CN201610065625.5A priority Critical patent/CN105740209B/zh
Publication of CN105740209A publication Critical patent/CN105740209A/zh
Application granted granted Critical
Publication of CN105740209B publication Critical patent/CN105740209B/zh
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/15—Correlation function computation including computation of convolution operations

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Computational Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Computing Systems (AREA)
  • Complex Calculations (AREA)
  • Measuring Frequencies, Analyzing Spectra (AREA)

Abstract

本发明公开了一种Givens迭代的Prony低频振荡分析方法，采用Givens迭代法来减小信号模型的阶数，提高解算特征方程系数求解方程组的数值稳定性，包括以下步骤：读入滤波后的测量数据；用奇异值分解法计算信号模型实际阶数；用Givens迭代法计算特征方程系数；计算频率和衰减因子；计算幅值和初相。本发明通过Givens迭代的Prony法分析低频振荡，扩展Prony分析时的信号模型的阶数可取信号模型的实际阶数，有效避免了特征方程系数求解方程组的系数矩阵的奇异，提高了数值计算的稳定性。本发明在所取信号模型阶数为信号模型实际阶数的条件下，本发明比传统Prony法计算出的低频振荡信息的结果精度高。

Description

一种Givens迭代的Prony低频振荡分析方法
技术领域
本发明涉及一种电力系统的Prony低频振荡的分析方法，特别是一种Givens迭代的Prony法低频振荡分析方法。
背景技术
电力系统运行时，由于扰动会发生频率为0.1Hz～2.5Hz的低频振荡，严重影响电力系统的稳定性，危及电网及其相关设备的安全运行。因此必须严密监测可能出现的低频振荡现象，及时处理，防止造成破坏系统稳定性的后果。Prony分析是目前分析低频振荡最常用最有效的方法之一，通过对信号进行Prony分析，可以直接得到反映低频振荡的幅值、相位、频率及衰减因子等参数。
Prony算法用一系列的任意频率、衰减因子、幅值和初相的指数函数的线性组合来拟合一个函数，不用再通过频域响应来求解，其计算量大为减少。也就是说该数学模型可以由一组衰减的正弦分量组成，是一种使用线性方程组来求解非线性问题的分析方法。
假设输入信号x(n)有N个采样点，即输入信号为x(0),…,x(N-1)。对该输入信号建立信号模型如下：
式中，是x(n)的近似值，p为信号模型的阶数，p值与信号分量数有关，bi和zi的表达式为：
式中，Ai为幅值，θi为初相，αi为衰减因子，fi为频率，Δt为采样时间间隔。
为了求式(1)中zi的值，由zi构造如下的特征多项式
则zi是特征多项式(4)构成的如下特征方程的根
式中，a0＝1。
根据式(1)和式(5)，可以推导出满足递推的常系数线性差分方程，该差分方程为：
式(6)中是实际测量数据x(n)的近似值，它们之间存在误差e(n)，即
式(6)代入式(7)，得
对式(8)进行最小二乘估计，使得误差平方和最小，则导致一组非线性方程，求解困难。为了实现线性估计，令
则式(8)可写成
式(10)写成矩阵形式为
式(11)为特征方程系数求解方程组，它的方程个数为N–p，未知数个数为p，一般情况下方程个数多于未知数个数，可以用最小二乘法计算估计值。求解式(11)的线性最小二乘法称为扩展Prony法。
采用最小二乘法，使最小，得到最小二乘法的法方程
式中，样本函数r(i,j)为(下式中上标(*)表示共轭)
最小误差能量为
用高斯消去法解式(12)，就可以求出a＝[a1 a2 … ap]T。a的精确求解是Prony法的关键，求出的a代入特征方程(5)，求出z＝[z1 z2 … zp]T后，就可以计算出信号分量的频率和衰减因子。式(5)是一元高次方程，一般采用QR分解法求解。
信号分量的频率和衰减因子计算公式如下：
式中，arctan、ln、Im、Re分别为反正切函数、自然对数函数、取复数虚部函数、取复数实部函数。
计算出信号模型的zi代入式(1)，得到关于未知数b的线性方程如下：
式中
b＝[b1 b2 … bp]T (18)
式(17)的Z是N×p维的范德蒙矩阵。由于zi各不相同，范德蒙矩阵Z的各列线性独立，为满列秩的。求解式(16)计算出b后，就可以计算信号分量的幅值和初相了。式(16)也是方程个数多于未知数个数的方程组，仍然用最小二乘法计算估计值。
幅值和初相计算公式为
由于事先不知道低频振荡中信号分量的个数，即式(11)中的p不确定，因此需要把信号模型的阶数取得大些(即取为pe>p)，再通过奇异值分解法计算出方程组系数矩阵的有效秩r，即可得到信号模型的实际阶数p＝r。这样把式(11)改写为式(21)
式(21)的系数矩阵为扩展阶矩阵
扩展阶矩阵Xe为(N-pe)×pe阶矩阵，用奇异值分解法可以求出该扩展阶矩阵的有效秩r，即可得到信号模型的实际阶数p＝r。
如图1所示，现有电力系统的Prony法低频振荡的分析方法，主要包括以下步骤：
A、读入滤波后的测量数据
测量数据中含有高频信号及噪声，需先对测量数据滤波去掉高频信号及噪声，才能进行Prony法低频振荡分析。
B、奇异值分解法信号模型实际阶数计算
根据测量数据形成扩展阶矩阵Xe，所取信号模型阶数pe大于信号模型实际阶数p(pe取满足的整数，符号为向下取整符号)，然后采用奇异值分解法对扩展阶矩阵Xe进行计算，求出扩展阶矩阵的所有奇异值，此奇异值为非负的，并按下列顺序排列：
σ11≥σ22≥…≥σhh≥0 (23)
式中
h＝min(N-pe,pe) (24)
求出所有奇异值并排顺后，按下式对奇异值归一化
σkk0＝σkk/σ11,1≤k≤h (25)
选择一个较小的正数作为阈值，使σkk0大于此阈值的最大整数k取为扩展阶矩阵Xe的有效秩r，则信号模型的实际阶数p＝r。
C、特征方程系数计算
计算时，一般需要把信号模型的阶数m取得比信号模型的实际阶数p大一些，否则计算结果的精度很差。取信号模型的阶数m大于信号模型实际阶数p(m取满足的整数，符号为向下取整符号)，重新列方程组(11)得到新的特征方程系数求解方程组如下：
此方程组的方程数(N–m)大于未知数个数m，可以采用最小二乘法对未知数进行估计，计算出对应特征方程的系数。由于方程组(26)是线性方程组，形成相应法方程后采用高斯消去法求解。
D、频率和衰减因子计算
信号模型的阶数为m时的特征方程为：
式中，a0＝1。
式(27)是一元高次方程，其解为实数或共轭复数，一般采用QR分解法求解，得到z后，就可以计算出信号分量的频率和衰减因子。
E、幅值和初相计算
信号模型的阶数为m时，求解信号模型的系数bi的公式(17)和(18)相应变为：
b＝[b1 b2 … bm]T (29)
把求出的z代入式(28),并根据式(16)计算出b后，就可以计算信号分量的幅值和初相。式(16)的方程数也大于未知数个数，仍然可以用最小二乘法对未知数进行估计。计算出信号分量按幅值由大到小排序，取前p个信号分量作为最终的信号分量结果。
计算出的信号分量为如下复指数函数，
一般情况下，复指数函数成对出现，即复指数函数和它的共轭值成对出现。函数xi′(t)的共轭值为
信号分量xi′(t)和信号分量xi″(t)叠加为正弦分量如下：
由于计算特征方程系数矩阵时，未知数个数m大于信号模型实际阶数p，计算出信号分量比实际的信号分量多。这些多出的分量幅值很小，可以根据幅值剔除这些多余的分量。
目前采用Prony模型，用最小二乘法计算低频振荡时，所取的信号模型阶数m要比信号模型的实际阶数p大得较多，否则计算精度很差。但取信号模型阶数m大于信号模型的实际阶数p时，特征方程系数求解方程组的系数矩阵已经奇异，数值条件变得很差，造成求解困难，信号模型的阶数比信号模型的实际阶数大得越多，数值条件就越差，特征方程系数求解方程组的求解就更困难。
申请人在申请号为201610008988.5专利申请中披露了残差迭代的Prony低频振荡分析方法来减小扩展Prony分析时的信号模型的阶数，提高解算特征方程系数求解方程组的数值稳定性，但仍然有进一步提高的余地。
发明内容
为解决现有技术存在的上述问题，本发明要提出一种Givens迭代的Prony法低频振荡分析方法，以便有效改善特征方程系数求解方程组的系数矩阵的数值条件，提高特征方程系数求解方程组数值稳定性。
为了实现上述目的，本发明提出了Givens迭代的Prony低频振荡分析方法来减小扩展Prony分析时的信号模型的阶数，使所取信号模型的阶数等于信号模型的实际阶数，提高解算特征方程系数求解方程组的数值稳定性。
本发明的技术方案如下：一种Givens迭代的Prony低频振荡分析方法，采用Givens迭代的方法来减小信号模型的阶数，提高解算特征方程系数求解方程组的数值稳定性。
特征方程系数求解方程组(11)是方程数多于未知数的一组方程，由于右端向量未知，用最小二乘法求解，就是把式(11)写成式(33)所示的不相容方程组，而令误差平方和最小。
把式(33)改写成下式：
上述方程式(34)的方程个数为N-p，未知数个数为p，一般情况下方程个数多于未知数个数，可以用最小二乘法计算估计值。
式(34)简写为
Xa＝c (35)
式中，
a＝[a1 a2 … ap]T (37)
c＝-[x(p) x(p+1) … x(N-1)]T (38)
对式(35)进行最小二乘估计，除了形成法方程后用高斯消去法求解的法方程法外，也可以用Givens变换的方法直接求解式(35)进行最小二乘估计。
方程组的求解存在数值稳定性问题，如果方程组的系数矩阵奇异，方程组则无解；如果方程组的系数矩阵接近奇异，方程组则属于病态，求解非常困难。通常用条件数来衡量方程组是否病态，条件数越大，方程组的病态越严重。用法方程法进行最小二乘估计会增大方程的条件数，使其变为原条件数的平方，不利于病态方程组的求解；Givens变换为正交变换，通过正交变换求解方程组时不会增加系数矩阵的条件数，能够比较准确反映方程组的奇异程度，有利于提高计算精度，但Givens变换也对接近奇异的系数矩阵非常敏感，可能导致方程组无法求解。
Givens变换的原理如下。
对式(35)两边进行Givens变换如下：
GXa＝Gc (39)
式中，X为式(36)所示的(N-p)×p阶矩阵，G为(N-p)×(N-p)阶矩阵，由一系列Givens变换矩阵相乘得到。
G＝Gp,N-p…Gp,p+1…Gi,N-p…Gi,j…Gi,i+1…G1,N-p…G1,2 (40)
式中，Gij为对单位阵E的元素eii、eij、eji、ejj修改后得到的矩阵
式中，cosθ和sinθ分别被变换矩阵X的元素的函数，为：
式中，上标(n)表示前次Givens变换得到的新元素。
由式(39)和式(40)可知，式(40)中各矩阵依次左乘X和c，对应的式(42)中cosθ和sinθ中X的元素都是进行前次Givens变换得到的新元素。
在实际采用扩展Prony方法分析低频振荡信号时，如果所取信号模型阶数m为信号模型的实际阶数p时，计算结果的精度很低。如果所取信号模型阶数m大于信号模型的实际阶数p，特征方程系数求解方程组的系数矩阵则接近奇异(由于数值计算的舍入误差，未真正奇异)，数值条件变得很差。Givens变换方法能够比较准确反映其所求解方程组系数矩阵的奇异程度，因而Givens变换方法对方程组的系数矩阵的奇异程度非常敏感，很难求解系数矩阵接近奇异的方程组，这造成Givens变换解算特征方程系数求解方程组比传统的法方程解算该方程组更困难。虽然法方程进行最小二乘估计，对方程组的系数矩阵的奇异程度不太敏感，但因为数值条件太差，也常常出现分析失败的情况。为此可以采用下面介绍的Givens迭代法进行低频振荡分析，此方法能够在所取信号模型阶数m为信号模型的实际阶数p时，得到非常精确的分析结果。
Givens迭代法低频振荡分析的原理如下：
由于a是用最小二乘法估计出的值，代入方程(35)的左端得到的值c′
Xa＝c′ (43)
显然不等于方程右端的值c，二者的差即为残差Δc，即
Δc＝c-c′ (44)
可以根据如下关于Δc的方程求残差所对应的a的偏差Δa
XΔa＝Δc (45)
用Givens变换方法求解式(45)，得到变量Δa估计值。
用得到的Δa对a的原值进行修正，得到a的新值作为估计值，效果可能更好，修正公式如下：
a(k+1)＝a(k)+Δa(k) (46)
式中，a(k+1)为第k+1次迭代得到的a值。
计算信号模型的系数bi的公式(16)中的可以表示为：
按照对式(11)的处理方法，式(16)也写成式(48)所示的不相容方程组，而令误差平方和最小。
Zb＝x (48)
式中，右端向量为
x＝[x(0) x(1) … x(N-1)]T (49)
本发明方案包括以下步骤：
A、读入滤波后的测量数据
测量数据中含有高频信号及噪声，需先对测量数据滤波去掉高频信号及噪声，才能进行Prony法低频振荡分析。
B、用奇异值分解法计算信号模型实际阶数
根据测量数据形成扩展阶矩阵Xe，所取信号模型阶数pe大于信号模型实际阶数p(pe取满足的整数，符号为向下取整符号)，然后采用奇异值分解算法计算扩展阶矩阵Xe的有效秩r，即可得到信号模型的实际阶数p＝r。
C、用Givens迭代法计算特征方程系数
根据得到的信号模型p，列方程组(11)，方程组(11)的方程数(N-p)大于未知数个数p，不能直接求解，可以用Givens变换对未知数进行最小二乘估计，计算出特征方程的系数a＝[a1 a2 … ap]T。
为了提高计算精度，采用迭代法求解a。其步骤如下：
C1、设迭代计数k＝0，收敛精度δ＝0.00001，最大迭代次数km＝30；
C2、根据测量数据形成方程的系数矩阵X和右端向量c；
C3、用Givens变换求解式Xa＝c，得到a的初值a(0)；
C4、把a(k)代入式(c′)(k)＝Xa(k)，得到(c′)(k)；
C5、按照式Δc(k)＝c-(c′)(k)，计算残差Δc(k)；
C6、用Givens变换求解式XΔa(k)＝Δc(k)，得到Δa(k)；
C8、判断Δamax是否小于δ，如果Δamax小于δ，转至步骤C11；
C9、令a(k+1)＝a(k)+Δa(k)；
C10、令k＝k+1，如果k<km，转至步骤C4；
C11、结束；
D、频率和衰减因子计算
特征方程的系数a求出后，代入特征方程(5)，求出z后，就可以计算出信号分量的频率和衰减因子。特征方程(5)是一元高次方程，其解为实数或共轭复数，一般采用QR分解法求解。
E、幅值和初相计算
把求出的z代入式(17),并根据式(48)计算出b后，就可以计算信号分量的幅值和初相。式(48)的方程数也大于未知数个数，不能直接求解，可以采用Givens变换对未知数进行最小二乘估计。
与现有技术相比，本发明具有以下有益效果：
1、本发明通过Givens迭代的Prony法分析低频振荡，扩展Prony分析时的信号模型的阶数可以取信号模型的实际阶数，有效避免了特征方程系数求解方程组的系数矩阵的奇异，提高了数值计算的稳定性。
2、本发明在所取信号模型阶数为信号模型实际阶数的条件下，本发明比传统Prony法计算出的低频振荡信息的结果精度高。
附图说明
本发明共有附图3张。其中：
图1是Prony法低频振荡分析的流程图。
图2是本发明Givens迭代的Prony法低频振荡分析的流程图。
图3是本发明Givens迭代的Prony法计算特征方程系数的迭代流程图。
具体实施方式
下面结合附图对本发明进行进一步地说明。采用图1所示的传统Prony法低频振荡分析的流程图和图2-3所示的本发明方法的流程图，对下列算例信号进行了分析。
x(t)＝10e-0.0311tcos(2π×0.2473t+20°)+2e-0.2652tcos(2π×0.42t+13°)+
10e-0.0027tcos(2π×0.423t+110°)+e-0.2936tcos(2π×1.0349t+60°)+ (50)
1.6e-45.8878tcos(2π×2.4t+10°)
上述信号包含5个正弦分量，等效为10个复指数分量，采用式(1)对上述信号建模时，信号模型的实际阶数p为10。算例信号的参数见表1。
表1算例信号的参数
对上述信号每间隔0.1s计算一个值，共计算200个值，用来模拟测量值。则此测量数据的采样间隔为0.1s，时长为20s，测量值个数为200。计算时，信号模型的阶数m取为p。本发明方法及其他对比方法计算结果见表2～表5。
本发明方法的计算结果见表2，本发明方法迭代3次。
表2信号模型阶数m取p时本发明方法的计算结果
分量序号 频率f(Hz) 衰减因子α 幅值A 初相θ(°)
1 0.24730 -0.03110 10.00000 20.00000
2 0.42000 -0.26520 2.00000 13.00000
3 0.42300 -0.00270 10.00000 110.00000
4 1.03490 -0.29360 1.00000 60.00000
5 2.40000 -45.87880 1.60000 10.00016
由表2参照表1可见，当所取信号模型阶数m恰好为实际阶数p时，本发明方法所得慢速衰减正弦分量误差很小，计算精度达到小数点后5位，快速衰减的正弦分量5误差也很小。
不迭代的Givens法的计算结果见表3。
表3信号模型阶数m取p时不迭代Givens方法的计算结果
分量序号 频率f(Hz) 衰减因子α 幅值A 初相θ(°)
1 0.24732 -0.03118 10.01283 19.93054
2 0.42210 -0.21853 2.23935 16.85023
3 0.42274 -0.00209 9.90945 111.74451
4 1.03491 -0.29367 0.99378 59.99634
5 0.00000 -53.76538 0.15105 0.00000
6 0.00000 -38.79178 3.08698 0.00000
由表3参照表1可见，当所取信号模型阶数m恰好为实际阶数p时，不迭代Givens法所得慢速衰减正弦分量误差较大，快速衰减信号分量5没有辨识出来，还多了原低频振荡信号根本不存在的两个具有不同衰减因子的直流分量5和直流分量6。
传统方法的计算结果见表4。
表4信号模型阶数m取p时传统方法的计算结果
分量序号 频率f(Hz) 衰减因子α 幅值A 初相θ(°)
1 0.24657 -0.03029 9.86541 22.52293
2 0.32772 -0.96931 1.64788 -2.34213
3 0.42390 -0.00418 10.13640 104.91954
4 1.03486 -0.29176 0.99660 62.12383
5 5.00000 -78.17693 1.38551 0.00000
注：表中分量5是复指数分量。
由表4参照表1可见，当所取信号模型阶数m恰好为实际阶数p时，传统方法所得慢速衰减正弦分量误差很大，比不迭代的Givens法的误差还要大；快速衰减信号分量5未能够辨识出来；计算结果中多出了一个复指数分量。
专利201610008988.5方法的计算结果见表5，由于信号模型阶数m取实际阶数p时，专利201610008988.5方法计算的一个衰减因子为62.57118，很大，按式(17)计算Z时，部分元素超过IEEE双精度浮点数的上溢值1.797×10308，无法形成Z矩阵，也就无法计算各信号分量幅值和初相，导致低频振荡分析失败。表5是m＝11时201610008988.5方法的计算结果，迭代37次，计算时间比本发明方法要长得多。
表5信号模型阶数m取11时专利201610008988.5方法的计算结果
分量序号 频率f(Hz) 衰减因子α 幅值A 初相θ(°)
1 0.24730 -0.03110 10.00000 20.00000
2 0.42000 -0.26520 2.00000 13.00001
3 0.42300 -0.00270 10.00000 110.00000
4 1.03490 -0.29360 1.00000 60.00000
5 2.39999 -45.87880 1.60000 10.00029
由表5参照表1可见，当所取信号模型阶数m恰好为实际阶数p时，专利201610008988.5方法不收敛；当所取信号模型阶数m大于实际阶数p时，专利201610008988.5方法所得的慢速衰减正弦分量误差很小，快速衰减的正弦分量5误差也很小。
信号模型阶数m分别取为10和11时，不迭代Givens法计算的特征方程系数a的结果见表6和表7。
表6信号模型阶数m取p时不迭代Givens方法计算出的特征方程系数值
序号 1 2 3 4 5
a -7.3586453 24.1891262 -46.4426937 57.0386443 -45.9678269
序号 6 7 8 9 10
a 23.8199795 -7.3247505 1.0694481 -0.0233223 0.0000857
表7信号模型阶数m取11时不迭代Givens方法计算出的特征方程系数值
序号 1 2 3 4
a -6.3216960×1011 4.6314930×1012 -1.5146110×1013 2.8897514×1013
序号 5 6 7 9
a -3.5203338×1013 2.8052201×1013 -1.4284045×1013 4.2517231×1012
序号 9 10 11
a -5.6839466×1011 1.1554009×109 -5.8127158×107
对比表6和表7可见，当信号模型阶数m取实际阶数p时，由于特征方程系数求解方程组的系数矩阵不奇异，不迭代Givens法计算的特征方程系数a的值比较正常；当所取信号模型阶数m大于实际阶数p时，由于特征方程系数求解方程组的系数矩阵接近奇异(由于数值计算的舍入误差，未真正奇异)，不迭代Givens法计算的特征方程系数a的值很大。另外，信号模型阶数m取11时，计算出一个信号分量的衰减因子为271.72424，很大，按式(28)计算Z时，部分元素超过IEEE双精度浮点数的上溢值1.797×10308，无法形成Z矩阵，也就无法计算各信号分量幅值和初相，导致低频振荡分析失败。说明Givens变换解算方程组的精度高，更能反映方程组的奇异程度。因此采用Givens迭代法分析低频振荡，可以使信号模型阶数m取实际阶数p，既能避免特征方程系数求解方程组的系数矩阵奇异，又能达到理想的求解精度。
本发明方法可以采用任何一种编程语言和编程环境实现，如C语言、C++、FORTRAN、Delphi等。开发环境可以采用Visual C++、Borland C++Builder、Visual FORTRAN等。

Claims (1)

1.一种Givens迭代的Prony低频振荡分析方法，其特征在于：包括以下步骤：
A、读入滤波后的测量数据
读入测量数据，并对测量数据滤波去掉高频信号及噪声；
B、用奇异值分解法计算信号模型实际阶数
根据测量数据形成扩展阶矩阵Xe，所取信号模型阶数pe大于信号模型的实际阶数p，pe取值为满足的整数，然后采用奇异值分解法计算扩展阶矩阵Xe的有效秩r，即得到信号模型的实际阶数p＝r；
扩展阶矩阵Xe为：
式中，pe为所取信号模型的阶数，x(i)表示第i个测量数据，i＝0、1、2、…、N-2，N为测量数据总个数；
C、用Givens迭代法计算特征方程系数
根据计算出的信号模型的实际阶数p，列出特征方程系数求解方程组，由于特征方程系数求解方程组的方程数大于未知数个数，不能直接求解，用Givens变换对未知数进行最小二乘估计，计算出特征方程的系数
a＝[a1 a2 … ap]T (2)
特征方程系数求解方程组为：
式中：
e(n)为估计值与实际测量数据x(n)的误差，定义如下：
信号模型的阶数为p时的特征方程为：
式中，a0＝1；
为了提高计算精度，采用迭代法求解a；其步骤如下：
C1、设迭代计数k＝0，收敛精度δ＝0.00001，最大迭代次数km＝30；
C2、根据测量数据形成方程组的系数矩阵X和右端向量c；
用最小二乘法对特征方程系数求解方程组(3)求解，就是把式(3)写成下式所示的不相容方程组，并令误差平方和最小，
X和c分别为式(5)经过变化后得到的方程组的简写形式的系数矩阵和右端向量：
Xa＝c (6)
式中，
a＝[a1 a2 … ap]T (8)
c＝-[x(p) x(p+1) … x(N-1)]T (9)
C3、用Givens变换求解式(6)，得到a的初值a(0)；
C4、把a(k)代入式(c′)(k)＝Xa(k)，得到(c′)(k)；
C5、按照式Δc(k)＝c-(c′)(k)，计算残差Δc(k)；
C6、用Givens变换求解式XΔa(k)＝Δc(k)，得到Δa(k)；
C7、求
C8、判断Δamax是否小于δ，如果Δamax小于δ，转至步骤D；
C9、令a(k+1)＝a(k)+Δa(k)；
C10、令k＝k+1，如果k<km，转至步骤C4；
D、计算频率和衰减因子
特征方程的系数a求出后，代入特征方程(4)，得到z后，并计算出信号分量的频率和衰减因子；特征方程(4)是一元高次方程，其解为实数或共轭复数，采用QR分解法求解；
信号分量的频率和衰减因子计算公式如下：
式中，arctan、ln、Im、Re分别为反正切函数、自然对数函数、取复数虚部函数、取复数实部函数，Δt为采样时间间隔；
E、计算幅值和初相
把求出的z代入以下公式(11)，并根据该式计算出b后，再计算信号分量的幅值和初相；式(11)的方程数也大于未知数个数，不能直接求解，采用Givens变换对未知数进行最小二乘估计；
计算b的公式为：
Zb＝x (11)
式中
b＝[b1 b2 … bp]T (13)
x＝[x(0) x(1) … x(N-1)]T (14)
幅值和初相计算公式为
步骤C中所述的Givens变换方法如下：
对式(6)两边进行Givens变换如下：
GXa＝Gc (16)
式中，X为式(7)所示的(N-p)×p阶矩阵，G为(N-p)×(N-p)阶矩阵，由一系列Givens变换矩阵相乘得到；
G＝Gp,N-p…Gp,p+1…Gi,N-p…Gi,j…Gi,i+1…G1,N-p…G1,2 (17)
式中，Gij为对单位阵E的元素eii、eij、eji、ejj修改后得到的矩阵
式中，cosθ和sinθ分别为 被变换矩阵X的元素的函数，为：
式中，上标(n)表示前次Givens变换得到的新元素；
由式(16)和式(17)可知，式(17)中各矩阵依次左乘X和c，对应的式(19)中cosθ和sinθ中X的元素都是进行前次Givens变换得到的新元素。
CN201610065625.5A 2016-01-28 2016-01-28 一种Givens迭代的Prony低频振荡分析方法 Expired - Fee Related CN105740209B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201610065625.5A CN105740209B (zh) 2016-01-28 2016-01-28 一种Givens迭代的Prony低频振荡分析方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201610065625.5A CN105740209B (zh) 2016-01-28 2016-01-28 一种Givens迭代的Prony低频振荡分析方法

Publications (2)

Publication Number Publication Date
CN105740209A CN105740209A (zh) 2016-07-06
CN105740209B true CN105740209B (zh) 2018-06-29

Family

ID=56247105

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201610065625.5A Expired - Fee Related CN105740209B (zh) 2016-01-28 2016-01-28 一种Givens迭代的Prony低频振荡分析方法

Country Status (1)

Country Link
CN (1) CN105740209B (zh)

Families Citing this family (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108152651A (zh) * 2017-12-27 2018-06-12 重庆水利电力职业技术学院 基于gmapm和som-lvq-ann的输电线路故障综合识别方法
CN111293706A (zh) * 2018-12-06 2020-06-16 中国移动通信集团山东有限公司 一种电力系统低频振荡参数辨识方法和装置
CN111046327B (zh) * 2019-12-18 2021-10-19 河海大学 适用于低频振荡与次同步振荡辨识的Prony分析方法
CN113517686B (zh) * 2021-05-06 2022-09-20 东方电子股份有限公司 基于Givens正交相似变换的低频振荡分析方法

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104217112A (zh) * 2014-09-02 2014-12-17 南京工程学院 一种基于多类型信号的电力系统低频振荡分析方法
CN104333020A (zh) * 2014-10-17 2015-02-04 广西电网有限责任公司 一种电力系统实时低频振荡分析及最优校正控制方法
CN104504257A (zh) * 2014-12-12 2015-04-08 国家电网公司 一种基于双重并行计算的在线Prony分析方法
CN104865497A (zh) * 2015-04-30 2015-08-26 国电南瑞科技股份有限公司 基于扩展Prony算法的低频振荡就地化在线辨识方法
CN104953583A (zh) * 2015-07-01 2015-09-30 河海大学 基于变点探测和Prony方法相结合的电力系统低频振荡在线监测方法

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104217112A (zh) * 2014-09-02 2014-12-17 南京工程学院 一种基于多类型信号的电力系统低频振荡分析方法
CN104333020A (zh) * 2014-10-17 2015-02-04 广西电网有限责任公司 一种电力系统实时低频振荡分析及最优校正控制方法
CN104504257A (zh) * 2014-12-12 2015-04-08 国家电网公司 一种基于双重并行计算的在线Prony分析方法
CN104865497A (zh) * 2015-04-30 2015-08-26 国电南瑞科技股份有限公司 基于扩展Prony算法的低频振荡就地化在线辨识方法
CN104953583A (zh) * 2015-07-01 2015-09-30 河海大学 基于变点探测和Prony方法相结合的电力系统低频振荡在线监测方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
基于Prony算法的电力系统低频振荡分析;董航 等;《高电压技术》;20060630;第32卷(第6期);第97-100页 *
基于经验模态分解滤波的低频振荡Prony分析;侯王宾 等;《物理学报》;20100531;第59卷(第5期);第3531-3537页 *

Also Published As

Publication number Publication date
CN105740209A (zh) 2016-07-06

Similar Documents

Publication Publication Date Title
CN106845010B (zh) 基于改进SVD降噪和Prony的低频振荡主导模式辨识方法
CN105740209B (zh) 一种Givens迭代的Prony低频振荡分析方法
Grant et al. Comparison of matrix pencil and prony methods for power system modal analysis of noisy signals
Su et al. Power harmonic and interharmonic detection method in renewable power based on Nuttall double‐window all‐phase FFT algorithm
CN105548739B (zh) 一种避雷器运行状态信号处理方法
CN105158723A (zh) 一种数字化电能计量系统的误差评估系统及方法
CN102222911A (zh) 基于ar模型和卡尔曼滤波的电力系统间谐波估计方法
Chang et al. Measurement techniques for stationary and time-varying harmonics
CN104849545B (zh) 一种电力系统频率测量方法及测量装置
Borkowski et al. Frequency estimation in interpolated discrete fourier transform with generalized maximum sidelobe decay windows for the control of power
Petrović et al. Computational effective modified Newton–Raphson algorithm for power harmonics parameters estimation
CN103207931A (zh) 基于mm和arma算法的次同步振荡模态辨识方法
CN105528496B (zh) 一种残差迭代的Prony低频振荡分析方法
CN117233687B (zh) 一种基于历史数据的cvt初始误差评估方法、介质及终端
CN108459992A (zh) 一种基于Prony算法的配电网同步相量测量方法
CN104407197B (zh) 一种基于三角函数迭代的信号相量测量的方法
Nanda et al. Field programmable gate array implementation of fuzzy variable step size adaptive linear element for adaptive frequency estimation<? show [AQ ID= Q1]?>
Zhou et al. Multi‐innovation stochastic gradient method for harmonic modelling of power signals
CN116826735A (zh) 一种针对新能源场站宽频振荡辨识方法、装置
CN107732940B (zh) 一种基于adpss的电力系统稳定器参数优化试验方法
Yang et al. Hybrid extended‐cubature Kalman filters for non‐linear continuous‐time fractional‐order systems involving uncorrelated and correlated noises using fractional‐order average derivative
US20140174946A1 (en) System and method for estimating impedance in electrochemical impedance spectroscopy
Juhlin et al. Fast gridless estimation of damped modes
Zhuang et al. Four harmonic analysis and energy metering algorithms based on a new cosine window function
CN105468880B (zh) 一种低频振荡快速衰减信号分量参数的提取方法

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20180629

Termination date: 20190128

CF01 Termination of patent right due to non-payment of annual fee