基于Givens正交相似变换的低频振荡分析方法与流程

基于givens正交相似变换的低频振荡分析方法
技术领域
1.本发明涉及电力系统自动化技术领域，尤其涉及一种基于givens正交相似变换的低频振荡分析方法。


背景技术：

2.随着电网互联区域日益扩大、电力系统结构和模型越来越复杂，电网的稳定和低频振荡问题日益严重。低频振荡是威胁电力系统稳定运行和限制互联电网传输能力的首要问题，严重时甚至导致电网的瓦解，因此对低频振荡的分析与控制已成为电力系统稳定研究的重要课题之一。电力系统低频振荡频率范围一般为0.2～2.5hz，prony算法能直接提取振荡信号的的幅值、相位、频率和衰减因子，是当前低频振荡模式识别应用最广泛的方法。该算法属于模态参数辨识的时域方法，数学实质为用一组指数项的多阶线性组合来拟合等间隔采样信号，从中直接计算出信号的幅值、相位、频率、阻尼因子等信息。当振荡模式为多阶且采样率增加时，识别振荡幅值和出相的计算量呈指数增加。
3.目前电力系统低频振荡prony算法中，在求解自回归(autoregressive model，ar)模型的正则方程求得差分方程系数后，是一个求实系数特征多项式求根的问题，其对应为一个求hessenberg矩阵的特征值，通常采用qr算法使用householder正交相似变换求解。prony算法中所形成的hessenberg矩阵是一个高度稀疏的矩阵，householder变换通常会将稀疏矩阵变为满阵，增加求解计算量，且由于prony算法采样中采样点数通常远大于系统的实际阶数，在系统模型不清楚的情况下，为了使样本矩阵尽量包含信号的全面特征，样本矩阵的维数取其最大值，一般选择初始阶数为采样点数的一半，因此矩阵维数较大。基于householder变换的qr算法求解高维数高稀疏的hessenberg矩阵特征值，计算耗时大大增加，难以满足电网实时分析需求，因此选择合适的正交变换矩阵，对提高电力系统低频振荡分析的效率起着举足轻重的作用。


技术实现要素：

4.本发明提出了一种基于givens正交相似变换的低频振荡分析方法，其目的是：提高电力系统低频振荡prony分析的计算效率。
5.本发明技术方案如下：
6.一种基于givens正交相似变换的低频振荡分析方法，包括如下步骤：
7.s1：读取振荡信号，确定输入信号采样点数以及采样间隔；
8.s2：求解低频振荡prony算法自回归模型的差分方程系数，根据差分方程系数形成高次方程对应的hessenberg矩阵，记为矩阵h；
9.s3：通过基于givens相似变换的二重qr算法求解矩阵h的特征值；
10.s4：根据矩阵h的特征值求解低频振荡各个主模式的衰减因子和频率；
11.s5：根据特征值形成雅克比矩阵，运用最小二乘法求解低频振荡各个主模式的幅值和相角。
12.作为本方法的进一步改进，所述步骤s3具体包括如下步骤：
13.s31：将矩阵h作为当前矩阵a；
14.s32：找出当前矩阵a的次对角线倒数第一个零元素的位置，得到当前矩阵的最后一个对角块
‑
不可约准三角矩阵c，所述矩阵c为以所述零元素所在的行和列为分割线，行分割线下侧以及列分割线右侧的所有元素形成的当前矩阵的右下角子矩阵；
15.s33：判断所述对角块
‑
不可约准三角矩阵c的阶数，如果矩阵c的阶数为1或2，则求出矩阵c的特征值，收缩当前矩阵a，并返回步骤s32，否则执行步骤s34；
16.s34：对当前矩阵a作基于givens正交相似变换的二重qr变换，得到准三角矩阵b；
17.s35：将准三角矩阵b作为新的当前矩阵a，转至步骤s32，进行循环迭代运算，直到当前矩阵a的阶数收缩为零，即求出矩阵h的全部特征值为止。
18.作为本方法的进一步改进，步骤s34所述对当前矩阵a作基于givens正交相似变换的二重qr变换得到准三角矩阵b的步骤为：
19.s341：确定一个正交变换矩阵g0，对矩阵a作相似变换：所述正交变换矩阵g0用于化第一列为上三角形，为矩阵a到准三角矩阵b的二重qr变换矩阵：
[0020][0021][0022]
b＝q
t
aq
[0023]
其中q是正交矩阵，原点位移μ1和μ2是a的右下角2
×
2阶子矩阵的特征值；
[0024]
s342：继续利用givens变换构成多个新的相似变换矩阵
[0025]
g1，g2，
…
，g
i
，
…
，g
n
‑2[0026]
通过n
‑
2次givens变换
[0027][0028]
化矩阵b1为准三角矩阵b；
[0029]
各相似变换矩阵之间的关系为：q
t
＝g
n
‑2…
g1g0。
[0030]
作为本方法的进一步改进，步骤s341所述正交变换矩阵g0的确定方法为：
[0031]
第一步，矩阵a的右下角二阶子矩阵的特征值为μ1和μ2，令δ＝μ1 μ2，ρ＝μ1μ2，则有
[0032][0033]
矩阵的第一列只有前三个元素非零，设该列为x0＝(p0，q0，r0，0，...，0)
t
，根据以下公式求出x0的非零元：
[0034][0035]
第二步，利用givens变换矩阵g
01
，通过g
01
x0消去元素q0，将x0变为x
′0；
[0036]
所述变换矩阵g
01
为
[0037][0038]
其中
[0039][0040]
计算结果为x
′0＝(p1，0，r0，0，...，0)
t
；
[0041]
第三步，利用givens变换矩阵g
02
，通过g
02
x
′0消去元素r0；所述变换矩阵g
02
为
[0042][0043]
其中
[0044][0045]
所述变换矩阵g
01
和变换矩阵g
02
之间的关系为：g0＝g
02
g
01
。
[0046]
作为本方法的进一步改进，步骤s4所述低频振荡各个主模式的衰减因子和频率的计算方法为：
[0047][0048]
其中z
i
为矩阵h的特征值对应的高次方程的特征根。
[0049]
相对于现有技术，本发明具有以下有益效果：相对于householder变换，本方法在prony分析计算时应用基于givens正交相似变换的二重qr变换，变换过程中增加的非零注入元素较少，大大降低了计算量，提高了电力系统低频振荡prony分析的计算效率，进一步提高了电力系统低频振荡分析的实用性。
附图说明
[0050]
图1为本发明的流程图。
具体实施方式
[0051]
下面结合附图详细说明本发明的技术方案：
[0052]
如图1，一种基于givens正交相似变换的低频振荡分析方法，包括如下步骤：
[0053]
s1：读取振荡信号，确定输入信号采样点数以及采样间隔。
[0054]
s2：求解低频振荡prony算法自回归模型的差分方程系数，根据差分方程系数形成对应的高次方程hessenberg矩阵，记为矩阵h。
[0055]
具体求解步骤如下：
[0056]
prony算法针对等间隔采样数据的一种信号辨识方法，信号的模型可以表现为一组衰减的正弦分量，认为采样信号x(0)，x(1)，
…
，x(n
‑
1)的估计值可以表示为
[0057][0058]
其中，b
i
，z
i
假定为复数，即：
[0059][0060]
式中，a
i
为振幅；α
i
为衰减因子；f
i
为频率；δt为采样间隔。
[0061]
为了求出拟合参数，最小化误差平方和
[0062][0063]
即可求出，即求解一个复杂的最小二乘问题。
[0064]
即设常系数线性差分方程为：
[0065][0066]
如果把z
i
看作线性常系数差分方程的根，那么z
i
满足特征方程
[0067][0068]
则对x(n)的拟合是一个常系数线性差分方程的齐次解，设估计误差为e(n)即
[0069][0070]
将(5)代入到(3)得到
[0071][0072]
其中
[0073][0074]
求解式(5)就可以得到a
i
、θ
i
、α
i
、f
i
的线性估计。由公式(5)x(n)可以看作是噪声e
(n)激励的一个自回归(autoregressive model，ar)模型产生的输出信号，求解ar模型的正则方程可以求得差分方程的系数a
i
。针对(6)参数的最小二乘估计是使得
[0075][0076]
最小。
[0077]
根据公式(6)有(x(n)用x
n
表示)
[0078][0079]
将公式(6)看作变量为a的代价函数
[0080][0081]
为使其最小化，令则有
[0082][0083]
对应的最小误差能量为
[0084][0085]
定义
[0086][0087]
根据(10)、(11)得到prony算法中求系数a的法方程形式：
[0088][0089]
得到
[0090][0091]
简记为：r
·
a＝r
[0092]
其中矩阵r的各元素计算如下：
[0093][0094]
其中i，j＝1，2，
…
，p
[0095]
而
[0096][0097]
其中i＝1，2，
…
，p.
[0098]
由上述求得差分方程系数a
i
后，代入下式通过多项式求根得到z
i
[0099][0100]
也就是
[0101][0102]
由线性代数知识可知，可以看成是某个实矩阵的特征多项式，因此，求解高次方程的特征根z
i
(i＝1，2，...，p)的问题就变成了求矩阵的全部特征值的问题。可以验证，上述特征多项式所对应的矩阵为上hessenberg矩阵
[0103][0104]
s3：由于矩阵h是一个上hessenberg矩阵，因此可以直接用二重qr算法求出全部特征值。
[0105]
实矩阵二重qr算法的一般步骤，矩阵a到b的二重qr变换，计算过程归结为
[0106][0107][0108]
b＝q
t
aq
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(20)
[0109]
其中a和b都是实矩阵，q是正交矩阵，原点位移μ1和μ2是a的右下角2
×
2阶子矩阵的特征值。
[0110]
显然，变换的关键是如何求正交矩阵q。因为qr分解式的实质是化实矩阵为上三角形：且q
t
可取n
‑
1个正交矩阵的乘积
[0111]
q
t
＝g
n
‑2…
g1g0[0112][0113]
矩阵g0，g1，
…
，g
n
‑2的作用是，用它们左乘时，依次化它的前n
‑
1列为上三角形。因此，正交相似变换b＝q
t
aq可分解成n
‑
1次正交相似变换
[0114][0115]
这就是说，qr分解和相似变换b＝q
t
aq可以同时进行，而并不需要求出矩阵q。
[0116]
本实施例通过基于givens相似变换的二重qr算法求解矩阵h的特征值，具体步骤为：
[0117]
s31：将矩阵h作为当前矩阵a。
[0118]
s32：找出当前矩阵a的次对角线倒数第一个零元素的位置，得到当前矩阵的最后一个对角块
‑
不可约准三角矩阵c，所述矩阵c为以所述零元素所在的行和列为分割线，行分割线下侧以及列分割线右侧的所有元素形成的当前矩阵的右下角子矩阵。
[0119]
s33：判断所述对角块
‑
不可约准三角矩阵c的阶数，如果矩阵c的阶数为1或2，则求出矩阵c的特征值，收缩当前矩阵a(收缩后的矩阵为与步骤s32所述右下角子矩阵相对应的左上角子矩阵)，并返回步骤s32，否则执行步骤s34。
[0120]
s34：对当前矩阵a作基于givens正交相似变换的二重qr变换，得到准三角矩阵b，具体步骤如下：
[0121]
s341：确定一个正交变换矩阵g0，对矩阵a作相似变换：，对矩阵a作相似变换：所述正交变换矩阵g0用于化第一列为上三角形。
[0122]
所述正交变换矩阵g0的确定方法为：
[0123]
第一步，矩阵a的右下角二阶子矩阵的特征值为μ1和μ2，令δ＝μ1 μ2，ρ＝μ1μ2，则有
[0124][0125]
因为a是准三角阵，所以矩阵的第一列只有前三个元素非零，设该列为x0＝(p0，q0，r0，0，...，0)t，根据以下公式求出x0的非零元：
次givens变换次givens变换化矩阵b1为准三角矩阵b；各相似变换矩阵之间的关系为：q
t
＝g
n
‑2…
g1g0。
[0144]
上述初等正交矩阵g
i
(i＝1，2，...，n
‑
2)有与g0类似的简单算式。
[0145]
对于矩阵b1，若令b1的第一列为x1，则矩阵g1应使线性变换g1x1消去x1的元素q1和r1，可见g1的算式以及分块对角形式与g0类同。按照各相似变换矩阵依次变换下去，一般地，b
i
与b1一样，次对角线以下仅有三个非零元素(第i列两个，第i 1列一个)
[0146][0147]
若令其第i列为x
i
＝(
×
，
×
，...，p
i
，q
i
，r
i
，0，...，0)
t
，其参数p
i
、q
i
、r
i
(i＝0，1，...，n
‑
2)从相应矩阵b
i
的第i列取得：
[0148][0149]
当i＝n
‑
2时，有
[0150][0151]
利用givens变换矩阵g
i
，其变换目的是通过对x
i
经线性变换g
i
x
i
后消去其中的元素q
i
和r
i
，其变换步骤分为两步：
[0152]
第一步，利用givens变换矩阵g
i1
，通过g
i1
x
i
消去元素q0，将x
i
变为x
′
i
；
[0153]
所述变换矩阵g
i1
为
[0154][0155]
其中
[0156][0157]
计算的结果是x
i1
＝(p
i1
，0，r
i
，0，...，0)
t
。
[0158]
第二步，利用givens变换矩阵g
i2
，通过g
i2
x
′
i
消去元素r
i
；
[0159]
所述变换矩阵g
i2
为
[0160][0161]
其中
[0162][0163]
所述变换矩阵g
i1
和变换矩阵g
i2
之间的关系为：g
i
＝g
i2
g
i1
。
[0164]
s35：将准三角矩阵b作为新的当前矩阵a，转至步骤s32，进行循环迭代运算，直到当前矩阵a的阶数收缩为零，即求出矩阵h的全部特征值为止。
[0165]
s4：根据矩阵h的特征值求解低频振荡各个主模式的衰减因子和频率，计算公式如下：
[0166][0167]
其中z
i
为矩阵h的特征值对应的高次方程的特征根。
[0168]
s5：根据特征根形成雅克比矩阵，运用最小二乘法求解低频振荡各个主模式的幅值和相角。
[0169]
综上所述，基于givens正交相似变换的电力系统低频振荡prony算法中，对采样值得拟合是一个常系数线性差分方程的齐次解，在差分方程系数已知的情况下，就可以根据差分方程的特征方程求根，也就是求解由差分方程系数构成的hessenberg矩阵的特征值。实际采样中，采样点数通常远大于系统的实际阶数，在对系统模型不清楚的情况下，为了使样本矩阵尽量包含信号的全面特征，样本矩阵的维数取其最大值，一般选择初始阶数为采样点数的一半，即n/2，因此hessenberg矩阵维数较大。
[0170]
本方法在二重qr算法中重复利用相似变换gag
t
，通常左乘是通过givens正交变换实现次对角线下非零元素的消元操作，当矩阵a的对角线以上元素是满阵或者是填充密度较大的矩阵时，利用householder变换计算量相对小，但是prony算法中差分方程系数构成的hessenberg矩阵具有高度稀疏的特点，householder变换后会将原始相应的行列变成满阵，而givens变换一方面在变换的过程中增加的非零注入元素相对较少，每次计算都少一个元素，当计算次数爆炸增长后，相应减少的计算量非常可观。另一方面，编程简单，不用计算householder变换的3
×
3矩阵。
[0171]
通过某电网低频振荡信号实际验证，采样频率取20hz，采样点数取200，阶数取100，形成hessenberg矩阵为100
×
100，采用基于householder变换的二重qr算法求解特征值乘加运算为6189681次，而采用本方法求解特征值乘加运算次数为4933812，运算效率提高21％，减少了浮点数乘加运算次数，提高了计算效率，进一步提高了电力系统低频振荡分析的实用性。