CN113010844B - 一种基于子空间动态模式分解的参与因子计算方法 - Google Patents

Abstract

一种基于子空间动态模式分解的参与因子计算方法，属于电力系统稳定与控制技术领域。该方法采用系统正常运行状态下的随机响应数据作为输入，在引入基于正交三角分解的子空间技术基础上，获取正交投影矩阵的奇异向量；再利用动态模式分解法对获得的奇异向量进行计算，得到系统低维近似状态矩阵；再将特征分析法与模式能量矩阵相结合，计算出发电机有功功率的参与因子。本发明能够从系统随机响应数据中实时辨识参与因子。通过与发电有功调度的阻尼比变化结果对比可见，本方法能够有效的提取发电机有功功率的参与因子，进而为实际发电机有功的快速准确调度提供指导依据。本发明避免了建模分析过程，能实时跟踪系统的动态变化，具有重要的实用价值。

Description

一种基于子空间动态模式分解的参与因子计算方法

技术领域

本发明属于电力系统稳定与控制技术领域，特别是涉及到一种基于子空间动态模式分解的参与因子计算方法。

背景技术

近年，随着电力系统的逐步扩展以及高比例可再生能源并网，使得系统的机电振荡现象频发，这对系统的安全稳定运行带来了巨大挑战。在电力系统正常运行过程中，负荷随机波动以及新能源有功输出的不确定性带来的随机小幅扰动无法避免，这使得量测到的系统变量外在表征为杂乱的随机信号。这类信号蕴含丰富的机电振荡动态信息，能够反映系统的实际工况，基于随机信号的小干扰稳定性研究，可以在系统出现主导振荡模式前，分析并挑选出其中的弱阻尼模式，实现电力系统正常运行状态下的弱阻尼模式预警，为运行人员采取适当的阻尼调制措施留有充足的反应时间，这对保障系统在受到扰动冲击时能快速抑制振荡，同时提高电力系统机电小干扰稳定性具有重要研究意义。

基于发电机有功调度的开环调制策略主要应用于区域间振荡模式的阻尼调制，通过减小送端系统中发电机有功输出的同时提高受端系统中发电机有功输出来降低流经联络线的有功功率，以此提高区域间模式阻尼比。开环有功调制过程中的关键环节是参调发电机的选择，利用有功功率参与因子来确定发电机在振荡模式中的参与程度，进而准确快速的选择参调发电机已引起有关学者的重视。

传统参与因子的计算多数以模型为基础，通过小干扰稳定分析获得与状态变量相关的参与因子，但对于大规模系统，小干扰分析过程中存在“维数灾”问题，在线应用较为困难，且对模型的准确性要求较高。

因此现有技术当中亟需要一种新型的技术方案来解决这一问题。

发明内容

本发明所要解决的技术问题是：提供一种基于子空间动态模式分解的参与因子计算方法用于解决传统参与因子的计算多数以模型为基础，在线应用较为困难，且对模型的准确性要求较高等技术问题。

一种基于子空间动态模式分解的参与因子计算方法，包括以下步骤，并且以下步骤顺次进行，

步骤一、将电力系统正常运行状态下的随机响应数据作为输入，输入的随机响应数据利用离散的线性电力系统方程表示，并将离散的线性系统数据序列构造数据矩阵：

X₁＝[x₁,...,x_m]＝[x₁,...,A^m-1x₁]

X₂＝[x₂,...,x_m+1]＝[Ax₁,...,A^mx₁]

式中：X₁为(n×m)阶矩阵，表示为由时刻k＝1至k＝m采集到的状态变量序列构成的数据矩阵；X₂为(n×m)阶矩阵，表示为由时刻k＝2至k＝m+1采集到的状态变量序列构成的数据矩阵；

步骤二、利用步骤一中的数据矩阵获得数据矩阵组(X₁，X₂，X₃，X₄)并构造新的数据矩阵X_p和X_f：

X_p＝[X₁ ^T,X₂ ^T]^T，X_f＝[X₃ ^T,X₄ ^T]^T

式中：X₃为(n×m)阶矩阵，表示由时刻k＝3至k＝m+2采集到的状态变量序列构成的数据矩阵；X₄为(n×m)阶矩阵，表示由时刻k＝4至k＝m+3采集到的状态变量序列构成的数据矩阵；X_p为X₁和X₂组合而成的数据矩阵；X_f为X₃和X₄组合而成的数据矩阵；

步骤三、利用X_p和X_f形成汉克尔矩阵H，对汉克尔矩阵H进行正交三角分解：

式中：L为下三角矩阵；Q^T为正交矩阵；L₁₁，L₂₁和L₂₂均为L的2n阶子矩阵；

和

均为Q^T的(2n×m)阶子矩阵；

将X_f的行向量投影在X_p的行空间上，形成正交投影矩阵O：

式中：(·)^-1表示对矩阵求逆；

对正交投影矩阵O进行奇异值分解：

O＝U_rS_rV_r ^T

式中：r为正交投影矩阵O的秩；U_r为正交投影矩阵O的左奇异向量；V_r为正交投影矩阵O的右奇异向量；S_r为正交投影矩阵O的奇异值矩阵；

步骤四、定义U_r1为U_r的前n行向量组合而成的矩阵，U_r2为U_r的后n行向量组合而成的矩阵，利用U_r1和U_r2替换步骤一中的X₁和X₂作为输入数据矩阵.

对U_r1进行奇异值分解：

U_r1＝U_lS_lV_l ^T

式中：l为U_r1的秩；U_l为U_r1的左奇异向量；V_l为U_r1的右奇异向量；S_l为U_r1的奇异值矩阵；

构造系统低维近似状态矩阵

步骤五、提取振荡模式参与因子

对步骤四中获得的系统低维近似状态矩阵

进行特征值分解：

式中：Λ＝diag[λ₁，...，λ_i，...],diag[·]表示为对角矩阵，λ_i表示为系统模式i的特征值；Φ＝[φ₁,...,φ_i,...]，φ_i为λ_i对应的右特征向量；

结合系统低维近似状态矩阵

公式和

的特征值分解公式，可得：

U_r2≈U_lΦΛΦ^-1S_lV_l ^T

获得空间结构矩阵E：

E＝Φ^-1S_lV_l ^T＝[e₁,...,e_i,...]^T

式中：e_i为矩阵E的第i项；

同时，得到模式i对应的能量矩阵a_i：

a_i＝||e_i||

式中：||·||表示矩阵的欧几里得范数；

则第i个模式中第j个状态变量的参与因子p_ji表示为：

p_ji＝|ψ_jiλ_ia_i|

式中：ψ_ji为ψ_i的第j项，其中ψ_i＝U_lφ_i，振荡模式参与因子提取完毕。

述步骤一中的离散的线性电力系统方程为：

x_k+1＝Ax_k

式中：x_k为采集到的离散线性系统在时刻k的状态变量序列，其中包含n个采样点；x_k+1为采集到的离散线性系统在时刻k+1的状态变量序列，其中包含n个采样点；A为系统状态矩阵。

通过上述设计方案，本发明可以带来如下有益效果：

本发明提供一种随机数据驱动下基于子空间动态模式分解的参与因子计算方法，该方法在引入基于LQ分解的子空间技术对输入随机数据进行预处理的基础上，利用动态模式分解方法得到系统低维近似状态矩阵，再通过计算模式所含能量获得发电机有功功率的参与因子。本发明具有能够从系统随机响应数据中实时辨识参与因子。通过与发电有功调度的阻尼比变化结果对比可见，所提方法能够有效的提取发电机有功功率的参与因子，进而为实际发电机有功的快速准确调度提供指导依据。基于系统随机响应数据的参与因子提取方法避免了建模分析过程，且能实时跟踪系统的动态变化，具有重要的实用价值。

附图说明

以下结合附图和具体实施方式对本发明作进一步的说明：

图1为本发明一种基于子空间动态模式分解的参与因子计算方法的实施例中IEEE16机系统模型图；

图2为本发明一种基于子空间动态模式分解的参与因子计算方法的实施例中系统有功功率响应数据波形图；

图3为本发明一种基于子空间动态模式分解的参与因子计算方法的子空间动态模式分解提取参与因子的求解流程框图。

具体实施方式

本发明提供一种基于子空间动态模式分解的参与因子计算方法，该方法在引入基于LQ分解的子空间技术对输入随机数据进行预处理的基础上，利用动态模式分解方法得到系统低维近似状态矩阵，再通过计算模式所含能量获得发电机有功功率的参与因子，其流程框图见图3。其具体步骤如下：

步骤1：传统动态模式分解算法：

存在离散的线性电力系统方程：

x_k+1＝Ax_k (1)

式中：x_k和x_k+1分别为采集到的离散线性系统在时刻k和k+1的状态变量序列，其中包含n个采样点；A为系统状态矩阵。

利用离散的线性系统数据序列构造数据矩阵：

X₁＝[x₁,...,x_m]＝[x₁,...,A^m-1x₁] (2)

X₂＝[x₂,...,x_m+1]＝[Ax₁,...,A^mx₁] (3)

式中：X₁为(n×m)阶矩阵，表示为由时刻k＝1至k＝m采集到的状态变量序列构成的数据矩阵；X₂为(n×m)阶矩阵，表示为由时刻k＝2至k＝m+1采集到的状态变量序列构成的数据矩阵；。

由式(2)(3)可以得到关系：

X₂＝AX₁ (4)

对数据矩阵X₁进行奇异值分解：

X₁＝USV^T (5)

式中：U和V分别为X₁的左奇异向量和右奇异向量；S为X₁的奇异值矩阵，(·)^T表示为矩阵的转置。

构造系统状态矩阵A的低维近似状态矩阵

通过上述方式提取的低维近似状态矩阵

能够从时域外在表征具有明显振荡过程的大扰动暂态振荡响应信号中获得有效的机电振荡模态信息。然而，环境激励下随机响应过程中并未发生实质性功率振荡，随机信号中包含的机电特征信息不足。因此，上述方法无法从随机响应信号中获得可靠的机电振荡特征参数。

步骤2：基于正交三角分解的子空间动态模式分解算法。

类似于式(2)和(3)，利用数据矩阵组(X₁，X₂，X₃，X₄)构造新的数据矩阵X_p，X_f：

X_p＝[X₁ ^T,X₂ ^T]^T，X_f＝[X₃ ^T,X₄ ^T]^T (7)

式中：X₃为(n×m)阶矩阵，表示由时刻k＝3至k＝m+2采集到的状态变量序列构成的数据矩阵；X₄为(n×m)阶矩阵，表示由时刻k＝4至k＝m+3采集到的状态变量序列构成的数据矩阵；X_p为X₁和X₂组合而成的数据矩阵；X_f为X₃和X₄组合而成的数据矩阵。

利用X_p和X_f形成汉克尔矩阵H，同时对H进行正交三角分解：

式中：L为下三角矩阵；Q^T为正交矩阵；L₁₁，L₂₁和L₂₂为L的2n阶子矩阵；

和

为Q^T的(2n×m)阶子矩阵。

将X_f的行向量投影在X_p的行空间上，形成正交投影矩阵O：

式中：(·)^-1表示对矩阵求逆。

对正交投影矩阵O进行奇异值分解：

O＝U_rS_rV_r ^T (10)

式中：r为正交投影矩阵O的秩；U_r和V_r分别为O的左奇异向量和右奇异向量；S_r为O的奇异值矩阵。

定义U_r1和U_r2分别为U_r的前n行向量和后n行向量所组合而成的矩阵。同时在传统动态模式分解算法的框架下利用U_r1和U_r2替换X₁和X₂作为输入数据矩阵，利用传统动态模式分解法进行辨识。

对U_r1进行奇异值分解：

U_r1＝U_lS_lV_l ^T (11)

式中：l为U_r1的秩；U_l和V_l分别为U_r1的左奇异向量和右奇异向量；S_l为U_r1的奇异值矩阵。

构造系统低维近似状态矩阵

步骤3:提取振荡模式参与因子。

对所获得的系统低维近似状态矩阵

进行特征值分解：

式中：Λ＝diag[λ₁，...，λ_i,...],diag[·]表示为对角矩阵，λ_i表示为系统模式i的特征值；Φ＝[φ₁,...,φ_i,...]，φ_i为λ_i对应的右特征向量。

结合式(12)和(13)，可得：

U_r2≈U_lΦΛΦ^-1S_lV_l ^T (14)

可以获得空间结构矩阵E：

E＝Φ^-1S_lV_l ^T＝[e₁,...,e_i,...]^T (15)

式中：e_i为矩阵E的第i项。

同时，可以得到模式i对应的能量矩阵a_i：

a_i＝||e_i|| (16)

式中：||·||表示矩阵的欧几里得范数。

则第i个模式中第j个状态变量的参与因子p_ji可以表示为：

p_ji＝|ψ_jiλ_ia_i| (17)

式中：ψ_ji为ψ_i的第j项，其中ψ_i＝U_lφ_i，振荡模式参与因子提取完毕。

本方法能够从随机响应信号中有效提取振荡频率、阻尼比和模态振型等特征信息。但更为重要的是，能够有效的提取发电机有功功率的参与因子，从而为后续基于发电机有功调制的互联电网阻尼提升策略中参调发电机组的选择提供指导依据。本发明避免了复杂的系统建模分析过程，仅依赖电力系统正常运行过程中获得的随机响应数据，即可实现电力系统小干扰稳定的在线评估和实时跟踪，具有重要的工程实用价值。

下面通过实施算例，结合数据分析，对本发明的技术方案作进一步具体说明。

实施例：

本发明所提方法普遍适用于电力系统机电振荡模式下发电机有功的参与因子计算，限于篇幅，对IEEE16机电力系统随机响应数据进行计算分析以验证所提计算方法的有效性，并与发电有功调度的阻尼比变化结果进行对比。具体情况如下：

本发明在IEEE16模型上验证所提取出参与因子的准确性，该系统根据发电机之间的振荡关系可分为5个区域，具体模型如图1所示。为了模拟系统中存在的随机响应数据，将基准值3％的随机扰动添加到负荷中，所模拟的发电机有功功率数据如图2所示。利用本文提出的算法对获得的系统随机响应数据进行辨识可获得4种区域间振荡模式，同时可得到这4种区域间模式的参与因子，如表1所示。

表1区域间模式的参与因子

由于区域间振荡模式阻尼比与区域间联络线的功率密切相关，因此可以通过调度参与区域间振荡的各台发电机有功功率来改变区域间联络线功率，区域间联络线功率的改变会导致区域间振荡模式阻尼比的增加或者减少，阻尼比变化的程度可以体现所调度的发电机在该模式下的参与程度。通过比较三种案例来验证所提取参与因子的准确性，各发电机有功调度1p.u为案例1，2p.u为案例2，3p.u为案例3。

表2模式3中区域A的发电机调度结果

表3模式3中区域B的发电机调度结果

以模式3为例，该模式表现为区域A(G1-G9)和区域B(G10-G13)的发电机组相互振荡，且未调度前的基本阻尼为2.80％。在案例1中，送端区域(区域A)发电机有功功率减少1p.u，受端区域(区域B)发电机有功功率增加1p.u，使得流经联络线的功率减少2p.u。其中固定G13有功功率增加1p.u，同时下调区域A其余发电机有功功率，阻尼比变化的统计结果如表2所示。由表2可以看出，区域间联络线功率变化量相同时，调节(G13，G6)发电机对对阻尼比的影响最大，当调节(G13，G5)和(G13，G7)发电机对时，阻尼比的变化仅次于调节(G13，G6)。根据阻尼比的变化顺序，可以得模式3中区域A发电机有功功率参与因子排序，即G6，G5，G7，G4，G9，G3，G2，G8，G1。同样，要验证区域B中发电机有功功率参与因子的准确性，将送端区域A中的发电机G6减小1p.u，受端区域B中的发电机G10-G13分别增加1p.u，阻尼变化量如表3所示。由表中可以看出阻尼比的变化顺序为(G13，G6)(G12，G6)(G11，G6)(G10，G6)，可以得到区域间模式3中区域B发电机有功功率参与因子排序，为G13，G12，G11，G10。阻尼比变化的程度可以体现所调度的发电机在该模式下的参与程度，在案例1中，上述发电机有功调度结果与子空间动态模式分解得到的参与因子结果相同，这也验证了子空间动态模式分解所提取的发电机有功功率参与因子的准确性。分析案例2和案例3得到的结论与案例1相同。

由于发电调度的增加导致区域间联络线的有功功率增加(案例2和案例3与案例1相比)，统计结果表明阻尼比变化与参与因子呈正相关，且阻尼比变化量排序与案例1所得结论相同。以上分析验证了仿真系统中与发电机有功功率相关的参与因子的准确性。该结论可以在系统弱阻尼模式下直接调度对区域间振荡模式影响最大的发电机，从而指导调度人员通过增加或减少发电机有功功率增加系统阻尼比，这对于维持电力系统的稳定性具有重要意义。

Claims (2)

1.一种基于子空间动态模式分解的参与因子计算方法，其特征是：包括以下步骤，并且以下步骤顺次进行，

步骤一、将电力系统正常运行状态下的随机响应数据作为输入，输入的随机响应数据利用离散的线性电力系统方程表示，并将离散的线性系统数据序列构造数据矩阵：

X₁＝[x₁,…,x_m]＝[x₁,…,A^m-1x₁]

X₂＝[x₂,…,x_m+1]＝[Ax₁,…,A^mx₁]

式中：X₁为(n×m)阶矩阵，表示为由时刻k＝1至k＝m采集到的状态变量序列构成的数据矩阵；X₂为(n×m)阶矩阵，表示为由时刻k＝2至k＝m+1采集到的状态变量序列构成的数据矩阵；

步骤二、利用步骤一中的数据矩阵获得数据矩阵组(X₁，X₂，X₃，X₄)并构造新的数据矩阵X_p和X_f：

X_p＝[X₁ ^T,X₂ ^T]^T，X_f＝[X₃ ^T,X₄ ^T]^T

式中：X₃为(n×m)阶矩阵，表示由时刻k＝3至k＝m+2采集到的状态变量序列构成的数据矩阵；X₄为(n×m)阶矩阵，表示由时刻k＝4至k＝m+3采集到的状态变量序列构成的数据矩阵；X_p为X₁和X₂组合而成的数据矩阵；X_f为X₃和X₄组合而成的数据矩阵；

步骤三、利用X_p和X_f形成汉克尔矩阵H，对汉克尔矩阵H进行正交三角分解：

式中：L为下三角矩阵；Q^T为正交矩阵；L₁₁，L₂₁和L₂₂均为L的2n阶子矩阵；

和

均为Q^T的(2n×m)阶子矩阵；

将X_f的行向量投影在X_p的行空间上，形成正交投影矩阵O：

式中：(·)^-1表示对矩阵求逆；

对正交投影矩阵O进行奇异值分解：

O＝U_rS_rV_r ^T

式中：r为正交投影矩阵O的秩；U_r为正交投影矩阵O的左奇异向量；V_r为正交投影矩阵O的右奇异向量；S_r为正交投影矩阵O的奇异值矩阵；

步骤四、定义U_r1为U_r的前n行向量组合而成的矩阵，U_r2为U_r的后n行向量组合而成的矩阵，利用U_r1和U_r2替换步骤一中的X₁和X₂作为输入数据矩阵；

对U_r1进行奇异值分解：

U_r1＝U_lS_lV_l ^T

式中：l为U_r1的秩；U_l为U_r1的左奇异向量；V_l为U_r1的右奇异向量；S_l为U_r1的奇异值矩阵；

构造系统低维近似状态矩阵

步骤五、提取振荡模式参与因子

对步骤四中获得的系统低维近似状态矩阵

进行特征值分解：

式中：Λ＝diag[λ₁,...,λ_i,...],diag[·]表示为对角矩阵，λ_i表示为系统模式i的特征值；Φ＝[φ₁,...,φ_i,...]，φ_i为λ_i对应的右特征向量；

结合系统低维近似状态矩阵

公式和

的特征值分解公式，可得：

U_r2≈U_lΦΛΦ^-1S_lV_l ^T

获得空间结构矩阵E：

E＝Φ^-1S_lV_l ^T＝[e₁,…,e_i,…]^T

式中：e_i为矩阵E的第i项；

同时，得到模式i对应的能量矩阵a_i：

a_i＝||e_i||

式中：||·||表示矩阵的欧几里得范数；

则第i个模式中第j个状态变量的参与因子p_ji表示为：

p_ji＝|ψ_jiλ_ia_i|

式中：ψ_ji为ψ_i的第j项，其中ψ_i＝U_lφ_i，振荡模式参与因子提取完毕。

2.根据权利要求1所述的一种基于子空间动态模式分解的参与因子计算方法，其特征是：所述步骤一中的离散的线性电力系统方程为

x_k+1＝Ax_k

式中：x_k为采集到的离散线性系统在时刻k的状态变量序列，其中包含n个采样点；x_k+1为采集到的离散线性系统在时刻k+1的状态变量序列，其中包含n个采样点；A为系统状态矩阵。

Images