CN103336909A - 一种风电接入电网的低频振荡辨识方法 - Google Patents

Abstract

本发明公开了风电并网监测技术领域中的一种风电接入电网的低频振荡辨识方法。包括：对不含风电机组的电力系统全系统线性化状态矩阵进行修正，得到修正后的全系统线性化状态矩阵；求取修正后的全系统线性化状态矩阵的特征值和特征向量；根据所述特征值和特征向量，计算特征值的频率、衰减阻尼比和机电回路相关比；根据特征值的频率、衰减阻尼比和机电回路相关比，对特征值进行分析，确定风电接入电网的低频振荡模式。本发明解决了现有的低频振荡辨识方法未考虑大规模风电接入与无法自动辨识系统主导振荡模式的问题。

Description

一种风电接入电网的低频振荡辨识方法

技术领域

本发明属于风电并网监测技术领域，尤其涉及一种风电接入电网的低频振荡辨识方法。

背景技术

随着风力发电技术的快速发展和国家在政策上对可再生能源发电的重视，我国风力发电建设已进入了一个快速发展的时期。同时，我国大容量风电场的规划、建设和运行使得风力发电在电力系统中的比重逐步增加，风电场与电力系统之间的相互影响也随之加大。

风电场一般在电网末端接入，远离负荷中心。由于风速的随机性与不确定性，大规模风电接入将对整个电网的稳定性有较大影响。其中，低频振荡，或称功率振荡，是较为普遍的一个问题。分析低频振荡常用特征值分析法，该方法将低频振荡归于小干扰稳定研究范畴。首先，在工作点附近将系统线性化，形成系统状态方程矩阵，如文献P.Kundur.Power System Stability and Control.McGraw-Hill,Inc.,1993中（详见第703页，12.1.3节），就给出了电力系统的全系统线性化状态矩阵A=A_D+B_D(Y_N+Y_D)^-1C_D的生成过程，其中的A_D、B_D、C_D和Y_D是参数全系统线性化状态矩阵A的4个参数矩阵，Y_N是节点导纳矩阵。其次，求取矩阵A的特征值和特征向量，并分析系统低频振荡的模式、振型、参与因子和灵敏度，从而考察原始系统的小干扰稳定性。特征值分析法可以准确反映出振荡模式的频率和阻尼以及非振荡模式的衰减率，因此分析人员可以及时发现系统中可能的全部弱阻尼模式或不稳定模式。

但是，特征值分析法依赖准确的全系统状态方程，并且在现有的方法中认为，风电不参与低频振荡，其只影响系统的低频振荡特性，因此需要考虑大规模风电接入以后对全系统状态方程的修正。而且，特征值分析法的计算结果包含大量无物理意义的根，筛选并识别系统中的弱阻尼区间振荡模式仍需要大量人工参与，这将不利于系统调度人员对系统中的低频振荡模式进行快速与有效的辨识。因此，开发一种考虑大规模风电接入的低频振荡辨识方法显得尤为重要。

发明内容

本发明的目的在于，提供一种风电接入电网的低频振荡辨识方法，用于解决电网现有的低频振荡辨识方法未考虑大规模风电接入的问题以及无法自动辨识电网主导振荡模式的问题。

为了实现上述目的，本发明提出的技术方案是，一种风电接入电网的低频振荡辨识方法，其特征是所述方法包括：

步骤1：对不含风电机组的电力系统全系统线性化状态矩阵A进行修正，得到包含风电机组的电力系统全系统线性化状态矩阵A′，并记为修正后的全系统线性化状态矩阵A′；

步骤2：求取修正后的全系统线性化状态矩阵A′的特征值λ_S和特征向量v_S；

步骤3：根据所述特征值λ_S和特征向量v_S，计算特征值的频率f、衰减阻尼比ξ和机电回路相关比ρ；

步骤4：根据特征值的频率f、衰减阻尼比ξ和机电回路相关比ρ，对特征值进行分析，确定风电接入电网的低频振荡模式。

所述风电机组为双馈异步风力发电机组。

所述修正后的全系统线性化状态矩阵A′=A′_D+B′_D(Y′_N+Y′_D)^-1C′_D；

其中， $A_D^{\′}=\left[\begin{array}{cc}{A}_D& \\ & A_{\mathrm{wf}}\end{array}\right];$

$A_{\mathrm{wf}}=\left[\begin{array}{ccccccc}-\frac{{\mathrm{\ω}}_s}{T_r^{\′}}& 0& -\frac{{\mathrm{\ω}}_s}{X_r^{\′}}& 0& 0& 0& 0\\ K_1{U}_t\frac{X_m}{X_s}(\frac{{\mathrm{\ω}}_s}{T_r^{\′}}-\frac{1}{T_1})& 0& K_1{U}_t\frac{X_m}{X_s}\frac{{\mathrm{\ω}}_s}{X_r^{\′}}& 0& 0& 0& 0\\ K_2(\frac{1}{T_2}-\frac{{\mathrm{\ω}}_s}{T_r^{\′}})& -\frac{K_2}{T_2}& -K_2\frac{{\mathrm{\ω}}_s}{X_r^{\′}}& 0& 0& 0& 0\\ 0& 0& 0& -\frac{{\mathrm{\ω}}_s}{T_r^{\′}}& 0& -\frac{{\mathrm{\ω}}_s}{X_r^{\′}}& 0\\ 0& 0& 0& K_1{U}_t\frac{X_m}{X_s}(\frac{{\mathrm{\ω}}_s}{T_r^{\′}}-\frac{1}{T_1})& 0& K_1{U}_t\frac{X_m}{X_s}\frac{{\mathrm{\ω}}_s}{X_r^{\′}}& 0\\ 0& 0& 0& K_2(\frac{1}{T_2}-\frac{{\mathrm{\ω}}_s}{T_r^{\′}})& -\frac{K_2}{T_2}& -K_2\frac{{\mathrm{\ω}}_s}{X_r^{\′}}& 0\\ 0& 0& 0& -U_t\frac{X_m}{2H{X}_s}& 0& 0& 0\end{array}\right];$

$B_D^{\′}=\left[\begin{array}{cc}{B}_D& \\ & B_{\mathrm{wf}}\end{array}\right];$

B_wf=[0]_7×2；

Y′_N是修正后的全系统节点导納矩阵；

$C_D^{\′}=\left[\begin{array}{cc}{C}_D& \\ & C_{\mathrm{wf}}\end{array}\right];$

$C_{\mathrm{wf}}={\left[\begin{array}{cc}\sin \mathrm{\δ}& -\cos \mathrm{\δ}\\ \cos \mathrm{\δ}& \sin \mathrm{\δ}\end{array}\right]}^{-1}\left[\begin{array}{ccccccc}-\frac{X_s}{X_m}& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& -\frac{X_s}{X_m}& 0& 0& 0\end{array}\right];$

$Y_D^{\′}=\left[\begin{array}{cc}{Y}_D& \\ & Y_{\mathrm{wf}}\end{array}\right];$

Y_wf是风电机组对应的节点导纳矩阵；

A_D、B_D、C_D和Y_D是修正前的系统线性化状态矩阵A的参数矩阵；

矩阵A_wf中的ω_s是同步角速度；

矩阵A_wf和矩阵C_wf中的X_s是定子电抗；

矩阵A_wf和矩阵C_wf中的X_m是定子与转子间的互抗；

矩阵A_wf中的X_r是转子电抗；

矩阵A_wf中的K₁是外环比例积分微分控制器的第一参数；

矩阵A_wf中的U_t是发电机机端电压；

矩阵A_wf中的T'_r=X'_r/R_r，R_r是转子电阻；

矩阵A_wf中的T₁是外环比例积分微分控制器的第二参数；

矩阵A_wf中的K₂是内环比例积分微分控制器的第一常数；

矩阵A_wf中的T₂是内环比例积分微分控制器的第二常数；

矩阵A_wf中的H是惯性常数；

矩阵C_wf中的δ是发电机功角。

所述求取修正后的全系统线性化状态矩阵A′的特征值λ_S和特征向量v_S包括：

步骤101：给定2-范数单位迭代初始相量

其中，

是n维复向量空间，n是特征向量的维数；

步骤102：令k=0，根据公式v^(k+1)=z^(k+1)/||z^(k+1)||₂计算v^(k+1)；

其中，z^(k+1)根据(A′-s^(k)I)z^(k+1)=v^(k)计算获得，

s^(k)是给定的初始位移点，为复常数；

A′是修正后的全系统线性化状态矩阵；

I是单位矩阵；

z^(k+1)是计算中间值；

[·]^H是共轭转置运算；

||·||₂是2-范数运算；

步骤103：判断是否达到收敛条件，如果达到收敛条件，则执行步骤104；否则，令k=k+1，返回步骤102；

步骤104：将v^(k+1)作为修正后的全系统线性化状态矩阵A′的特征向量，记为v_S，将s^(k)记为s；

步骤105：根据公式λ^(k+1)=[v^(k+1)]^HA′v^(k+1)计算特征值的迭代值λ^(k+1)，再根据公式

计算得到修正后的全系统线性化状态矩阵A′的特征值λ_S。

所述收敛条件为

其中，λ^(k)为特征值的迭代值，λ^(k)=[v^(k)]^HA′v^(k)，ε为设定误差。

所述计算特征值的频率f采用公式

其中，ω是特征值λ_S虚部的值。

所述计算特征值的衰减阻尼比ξ采用公式

其中，α是特征值λ_S实部的值，ω是特征值λ_S虚部的值。

所述计算特征值的机电回路相关比ρ采用公式

其中，x_k是第k个状态变量，1≤k≤n，n是状态变量的个数；

Δω是角速度变化量；

Δδ是功角变化量；

p_ki为相关因子，1≤i≤n，n为修正后的全系统线性化状态矩阵A′的阶数；

u_S是与特征值λ_S相对应的左特征向量；

v_S是与特征值λ_S相对应的右特征向量；

u_kS是u_S的第k个分量；

v_kS是v_S的第k个分量。

所述步骤4具体包括：

步骤201：选取频率为0.1Hz<f<2.5Hz、衰减阻尼比ξ<3%与机电回路相关比ρ>1的特征值；

步骤202：在选取的特征值下进行机组相位分群判定，即如果存在选取的特征值对应的特征向量集合的子集U，使得对任意i,j∈U，有

则判定在选取的该特征值下机组相位分群；如果机组相位分群且机组相位分群的个数大于等于2，则执行步骤203；否则，结束；

其中，θ_i是第i个特征向量的相位，

α_i和ω_i分别为选取的特征值对应的第i个特征向量的实部与虚部；

θ_j是第j个特征向量的相位，

α_j和ω_j分别为选取的特征值对应的第j个特征向量的实部与虚部；

M是用于表示特征向量相位接近程度的设定值；

步骤203：将每个选取的特征值作为一个风电接入电网的低频振荡模式。

本发明提供的方法，在大规模风电接入电网的背景下，建立了考虑大规模风电接入电网的系统数学模型，采用逆迭代转Rayleigh商迭代法进行小干扰稳定计算，并对振荡模式进行振型图自动分析，能够快速准确辨识的电网中的低频振荡模式，解决现有的低频振荡辨识方法未考虑大规模风电接入与无法自动辨识系统主导振荡模式的问题，为实际电网调度进行低频振荡预警与控制提供了参考。

附图说明

图1是风电接入电网的低频振荡辨识方法流程图；

图2是风力发电机组转子侧有功控制系统框图；

图3是风力发电机组转子侧无功控制系统框图；

图4是甘肃河西地区电网接线图；

图5是低频振荡振型自动分析图；其中（a）是一组特征值下的机组相位分群示意图，（b）是另一组特征值下的机组不呈现相位分群示意图。

具体实施方式

下面结合附图，对优选实施例作详细说明。应该强调的是，下述说明仅仅是示例性的，而不是为了限制本发明的范围及其应用。

实施例1

步骤1：首先建立大规模风电接入的电力系统数学模型。即对不含风电机组的电力系统全系统线性化状态矩阵A进行修正，得到包含风电机组的电力系统全系统线性化状态矩阵A′，并记为修正后的全系统线性化状态矩阵A′。

在电力系统中，每一个动态元件包括发电机及其控制系统的方程可表示为：

${\stackrel{\&CenterDot;}{x}}_i=f_d(x_i,U)$

i_i=g_d(x_i,U)

上式中，x_i为动态元件的状态变量，

为动态元件状态变量x_i的对时间t的导数，U为动态元件的端电压，i_i为注入网络的电流。每一个动态元件的线性化状态方程式可表示化为：

$\mathrm{\&Delta;}{\stackrel{\&CenterDot;}{x}}_i=A_i\mathrm{\&Delta;}{x}_i+B_i\mathrm{\&Delta;U}$

Δi_i=C_DΔx_i-Y_DΔU

其中B_i、Y_i对应于动态元件端电压或用来控制的远端电压的元素为非零以外，其他都是零。电流及电压都具有实部及虚部两个分量。

将所有动态元件的上述方程结合在一起，就形成了全系统状态方程式：

$\mathrm{\&Delta;}{\stackrel{\&CenterDot;}{x}}_i=A_D\mathrm{\&Delta;}{x}_i+B_D\mathrm{\&Delta;U}$

Δi=C_DΔx-Y_DΔU

其中，Δx就是全系统的状态向量，A_D及C_D是由对应于各个动态元件的A_i及C_i形成的分块的对角矩阵。Δi是动态元件对网络的注入电流向量，除了对接有动态元件的母线，其他都是零。网络是以导纳矩阵来表示的，因此Δi=Y_NΔU。

代入上式可得：

ΔU=(Y_N+Y_D)^-1C_DΔx

将ΔU代入全系统状态方程式可得：

$\mathrm{\&Delta;}\stackrel{\&CenterDot;}{x}=[A_D+B_D{(Y_N+Y_D)}^{-1}{C}_D]\mathrm{\&Delta;x}=\mathrm{A\&Delta;x}$

这样就得到了电力系统全系统线性化状态矩阵：

A=A_D+B_D(Y_N+Y_D)^-1C_D

由于风电机组也属于动态元件，因此大规模风电接入电力系统以后，需要针对风电机组的特点对全系统数学模型进行拓展。

因为电力系统中主要接入的风机类型为双馈异步风机，因此下面以双馈异步风力发电机组为例进行说明。确定电网当前实际的运行方式，获得风电机组出力及省际间联络线断面潮流。

双馈异步风机属于异步发电机，首先介绍异步发电机的模型，异步发电机组实用三阶模型可以用下式表示：

$\left\{\begin{array}{c}{T}_{\mathrm{j\&Sigma;}}\frac{{\mathrm{ds}}_{\mathrm{\&Sigma;}}}{\mathrm{dt}}=M_e-M_m\\ T_{d0\mathrm{\&Sigma;}}^{\&prime;}\frac{{\mathrm{dE}}_{\mathrm{\&omega;x}}^{\&prime;}}{2}=-E_{\mathrm{\&omega;x}}^{\&prime;}+(X_{\mathrm{\&Sigma;}}-X_{\mathrm{\&Sigma;}}^{\&prime;})I_{\mathrm{sy}}+T_{d0\mathrm{\&Sigma;}}^{\&prime;}{\mathrm{\&omega;}}_0{s}_{\mathrm{\&Sigma;}}{E}_{\mathrm{\&omega;y}}^{\&prime;}\\ T_{d0\mathrm{\&Sigma;}}^{\&prime;}\frac{d{E}_{\mathrm{\&omega;y}}^{\&prime;}}{2}=-E_{\mathrm{\&omega;y}}^{\&prime;}-(X_{\mathrm{\&Sigma;}}-X_{\mathrm{\&Sigma;}}^{\&prime;})I_{\mathrm{sx}}-T_{d0\mathrm{\&Sigma;}}^{\&prime;}{\mathrm{\&omega;}}_0{s}_{\mathrm{\&Sigma;}}{E}_{\mathrm{\&omega;x}}^{\&prime;}\end{array}\right.$

其中，Σ表示求和，T_jΣ是发电机转子惯性时间常数，s_Σ是转差率，M_m是机械转矩，M_e是电磁转矩。T′_d0Σ是定子绕组开路时转子绕组时间常数且

x₂为转子漏抗标幺值，x_m为励磁电抗标幺值，其中f₀为系统频率基值，E′_ωx是暂态电势实部，X_Σ是同步电抗且X_Σ=x₁+x_m，X′_∑是等效暂态电抗(p.u.)且

x₁为定子漏抗标幺值，x_m为励磁电抗标幺值，I_sy定子电流的虚部，ω₀为同步角速度且ω₀=2πf₀，E′_ωy暂态电势虚部，I_sx定子电流实部。

选取E′_ωx、E′_ωy、s_Σ为状态变量，将上式写成矩阵形式，从而可以直接对全系统线性化状态矩阵进行扩展。

双馈型异步风力发电机的建模主要包括双馈型风力发电机组动态建模和换流器建模。在同步旋转的dq坐标系下，双馈发电机的电压方程和磁链方程分别如下两式所示：

$\left\{\begin{array}{c}{U}_{\mathrm{ds}}=-R_s{I}_{\mathrm{ds}}-{\mathrm{\&psi;}}_{\mathrm{qs}}+\frac{1}{{\mathrm{\&omega;}}_s}\frac{{\mathrm{d\&psi;}}_{\mathrm{ds}}}{\mathrm{dt}}\\ U_{\mathrm{qs}}=-R_s{I}_{\mathrm{qs}}+{\mathrm{\&psi;}}_{\mathrm{ds}}+\frac{1}{{\mathrm{\&omega;}}_s}\frac{{\mathrm{d\&psi;}}_{\mathrm{qs}}}{\mathrm{dt}}\\ U_{\mathrm{dr}}=-R_r{I}_{\mathrm{dr}}-{\mathrm{s\&psi;}}_{\mathrm{qr}}+\frac{1}{{\mathrm{\&omega;}}_s}\frac{{\mathrm{d\&psi;}}_{\mathrm{dr}}}{\mathrm{dt}}\\ U_{\mathrm{qr}}=-R_r{I}_{\mathrm{qr}}+{\mathrm{s\&psi;}}_{\mathrm{dr}}+\frac{1}{{\mathrm{\&omega;}}_s}\frac{{\mathrm{d\&psi;}}_{\mathrm{qr}}}{\mathrm{dt}}\end{array}\right.$

$\left\{\begin{array}{c}{\mathrm{\&psi;}}_{\mathrm{ds}}=-X_s{I}_{\mathrm{ds}}-X_m{I}_{\mathrm{dr}}\\ {\mathrm{\&psi;}}_{\mathrm{qs}}=-X_s{I}_{\mathrm{qs}}-I_m{I}_{\mathrm{qr}}\\ {\mathrm{\&psi;}}_{\mathrm{dr}}=-X_r{I}_{\mathrm{dr}}-X_m{I}_{\mathrm{ds}}\\ {\mathrm{\&psi;}}_{\mathrm{qr}}=-X_r{I}_{\mathrm{qr}}-X_m{I}_{\mathrm{qs}}\end{array}\right.$

其中，U_ds是定子电压d轴分量，U_qs是定子电压q轴分量，U_dr是转子电压d轴分量，U_qr是转子电压q轴分量，R_s是定子电阻，R_r是转子电阻，I_ds是定子电流d轴分量，I_qs是定子电流q轴分量，I_dr是转子电流d轴分量，I_qr是转子电流q轴分量，ψ_qs是定子磁链q轴分量，ψ_ds是定子磁链d轴分量，ψ_qr为是转子磁链q轴分量，ψ_dr是转子磁链d轴分量，ω_s是同步角速度，s是转差率。X_s是定子电抗且X_s=X_sσ+X_m，X_sσ是定子漏抗，X_m是定子与转子间的互抗，X_r是转子电抗且X_r=X_rσ+X_m，X_rσ是转子漏抗，X_m是定子与转子间的互抗。

双馈异步风力发电机组的无功控制部分和有功控制部分可以分别用下两式所示的方程描述。

参考图2和图3，上述方程中，ω_s是同步角速度，X_s是定子电抗，X_m是定子与转子间的互抗，X_r是转子电抗，K₁是外环比例积分微分控制器的第一参数，U_t是发电机机端电压，T'_r=X'_r/R_r，R_r是转子电阻，T1是外环比例积分微分控制器的第二参数，K₂是内环比例积分微分控制器的第一常数，T₂是内环比例积分微分控制器的第二常数，I_{qr_ref}是转子电流q轴分量参考值，P_sref是有功功率参考值，

是将I_qr和I_dr解耦引入的量，其含义为转子电压q轴分量。

通过双馈异步风力发电机组的无功控制部分和有功控制部分的方程，可以得到双馈异步风力发电机组小扰动分析的7阶状态矩阵。

Λ_wf=A_wf+B_wf(Y_wf+D_wf)^-1C_wf

上述矩阵中，

$A_{\mathrm{wf}}=\left[\begin{array}{ccccccc}-\frac{{\mathrm{\&omega;}}_s}{T_r^{\&prime;}}& 0& -\frac{{\mathrm{\&omega;}}_s}{X_r^{\&prime;}}& 0& 0& 0& 0\\ K_1{U}_t\frac{X_m}{X_s}(\frac{{\mathrm{\&omega;}}_s}{T_r^{\&prime;}}-\frac{1}{T_1})& 0& K_1{U}_t\frac{X_m}{X_s}\frac{{\mathrm{\&omega;}}_s}{X_r^{\&prime;}}& 0& 0& 0& 0\\ K_2(\frac{1}{T_2}-\frac{{\mathrm{\&omega;}}_s}{T_r^{\&prime;}})& -\frac{K_2}{T_2}& -K_2\frac{{\mathrm{\&omega;}}_s}{X_r^{\&prime;}}& 0& 0& 0& 0\\ 0& 0& 0& -\frac{{\mathrm{\&omega;}}_s}{T_r^{\&prime;}}& 0& -\frac{{\mathrm{\&omega;}}_s}{X_r^{\&prime;}}& 0\\ 0& 0& 0& K_1{U}_t\frac{X_m}{X_s}(\frac{{\mathrm{\&omega;}}_s}{T_r^{\&prime;}}-\frac{1}{T_1})& 0& K_1{U}_t\frac{X_m}{X_s}\frac{{\mathrm{\&omega;}}_s}{X_r^{\&prime;}}& 0\\ 0& 0& 0& K_2(\frac{1}{T_2}-\frac{{\mathrm{\&omega;}}_s}{T_r^{\&prime;}})& -\frac{K_2}{T_2}& -K_2\frac{{\mathrm{\&omega;}}_s}{X_r^{\&prime;}}& 0\\ 0& 0& 0& -U_t\frac{X_m}{2H{X}_s}& 0& 0& 0\end{array}\right],$

B_wf=[0]_7×2，

$C_{\mathrm{wf}}={\left[\begin{array}{cc}\sin \mathrm{\&delta;}& -\cos \mathrm{\&delta;}\\ \cos \mathrm{\&delta;}& \sin \mathrm{\&delta;}\end{array}\right]}^{-1}\left[\begin{array}{ccccccc}-\frac{X_s}{X_m}& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& -\frac{X_s}{X_m}& 0& 0& 0\end{array}\right],$

D_wf=[0]_2×2。

上述矩阵A_wf中的H是惯性常数，矩阵C_wf中的δ是发电机功角。

形成双馈异步风力发电机组后，对全系统线性化状态矩阵A进行修正，即可得到考虑大规模风电接入的电力系统数学模型。以接入一台双馈异步风机为例，双馈风机的状态变量为：

因此双馈风机的线性化状态方程式为：

$\mathrm{\&Delta;}{\stackrel{\&CenterDot;}{x}}_i=A_{\mathrm{wf}}\mathrm{\&Delta;}{x}_i+B_{\mathrm{wf}}\mathrm{\&Delta;}{U}_i$

Δi_i=C_wfΔx_i-Y_wfΔU_i

由于在系统中,该双馈型风机与系统中其他动态元件无耦合（即无受控关系），则大规模风电接入系统以后，全系统状态方程式为

$\mathrm{\&Delta;}{\stackrel{\&CenterDot;}{x}}^{\&prime;}=A_D^{\&prime;}\mathrm{\&Delta;}{x}^{\&prime;}+B_D^{\&prime;}\mathrm{\&Delta;}{U}^{\&prime;}$

Δi′=C′_DΔx′-Y′_DΔU′

其中， $\mathrm{\&Delta;}{\stackrel{\&CenterDot;}{x}}^{\&prime;}=\left[\begin{array}{c}\mathrm{\&Delta;}\stackrel{\&CenterDot;}{x}\\ \mathrm{\&Delta;}{\stackrel{\&CenterDot;}{x}}_i\end{array}\right],A_D^{\&prime;}=\left[\begin{array}{cc}{A}_D& \\ & A_{\mathrm{wf}}\end{array}\right],\mathrm{\&Delta;}{x}^{\&prime;}=\left[\begin{array}{c}\mathrm{\&Delta;x}\\ \mathrm{\&Delta;}{x}_i\end{array}\right],B_D^{\&prime;}=\left[\begin{array}{cc}{B}_D& \\ & B_{\mathrm{wf}}\end{array}\right],$

$\mathrm{\&Delta;}{U}^{\&prime;}=\left[\begin{array}{c}\mathrm{\&Delta;U}\\ \mathrm{\&Delta;}{U}_i\end{array}\right],\mathrm{\&Delta;}{i}^{\&prime;}=\left[\begin{array}{c}\mathrm{\&Delta;i}\\ \mathrm{\&Delta;}{i}_i\end{array}\right],C_D^{\&prime;}=\left[\begin{array}{cc}{C}_D& \\ & C_{\mathrm{wf}}\end{array}\right],Y_D^{\&prime;}=\left[\begin{array}{cc}{Y}_D& \\ & Y_{\mathrm{wf}}\end{array}\right]$

同样地，由Δi'_i=Y'_nΔU'，代入上式并消去ΔU'可得

$\mathrm{\&Delta;}{\stackrel{\&CenterDot;}{x}}^{\&prime;}={[A}_D^{\&prime;}+B_D^{\&prime;}{(Y_N^{\&prime;}+Y_D^{\&prime;})}^{-1}{C}_D^{\&prime;}]\mathrm{\&Delta;}{x}^{\&prime;}=A^{\&prime;}\mathrm{\&Delta;}{x}^{\&prime;}$

从而得到A'=A'_D+B'_D(Y'_N+Y'_D)^-1C'_D，即修正后的全系统线性化状态矩阵A′。

步骤2：求取修正后的全系统线性化状态矩阵A′的特征值λ_S和特征向量v_S。

根据得到的电力系统数学模型，采用逆迭代转Rayleigh商迭代法进行小干扰稳定计算，从而获得矩阵A′的特征值λ_S和特征向量v_S。

逆迭代转Rayleigh商法是指：对于由状态方程描述的线性系统，其小干扰稳定性由状态矩阵的所有特征值决定。如果所有的特征值实部都为负，则系统在该运行点是稳定的，否则，只要有一个特征值具有正的实部或者有实部为零的重根，则系统在该运行点是不稳定的。因此，分析系统在某运行点的小干扰稳定性问题，可以归结为求解状态矩阵A′的全部特征值的问题。

在电力系统中，系统阶数很高，直接计算状态矩阵的特征值会受到系统规模的限制，通常采用稀疏特征值算法进行计算。逆迭代转Rayleigh商迭代法是稀疏特征值算法的一种。为利用稀疏特征值算法，需要将矩阵A′作如下变换：

S=(A′-sI)^-1

式中，复常数s为给定的初始位移点，通过如上谱变换，将状态矩阵A′关键特征值的求解问题转变为复常数矩阵S模值最大的特征值的计算。

设λ_A和v_A分别为矩阵A′的特征值和特征向量，λ_S和v_S分别为矩阵S的特征值和特征向量，且有：

λ_S=1/(λ_A-s),λ_A=1/λ_S+s

v_s=v_A

在矩阵S的特征值求解过程中，当k→∞时，

${\mathrm{\&lambda;}}^{\left(k\right)}\&RightArrow;1-({\mathrm{\&lambda;}}_p-s),s+\frac{1}{{\mathrm{\&lambda;}}^{\left(k\right)}}\&RightArrow;{\mathrm{\&lambda;}}_P$

v^(k)→v_p

式中，λ_P为矩阵A′的所有特征值中最靠近s的特征值，v_P为相应的特征向量。即将矩阵A′关键特征值映射为矩阵S的最大特征值。

逆迭代转Rayleigh商迭代法是逆迭代法和Rayleigh商迭代法相结合的算法。逆迭代法是计算非奇异矩阵A′模值最大的特征值和相应特征向量的一种非常有效的迭代方法。将逆迭代法用于矩阵A′^-1，可以用来计算非奇异矩阵A′的模数最小的特征值及其相应的特征向量，此时得到的方法称为逆迭代法（又称反幂法）。

对于复矩阵S=(A′-sI)^-1，逆迭代法的迭代过程如下：

1）给定2-范数单位迭代初始向量

2）对k=1,2,...

(A′-sI)z^(k)=v^(k-1)

v^(k)=z^(k)/||z^(k)||₂

λ^(k)=[v^(k)]^HA′v^(k)

则当k→∞时，

${\mathrm{\&lambda;}}^{\left(k\right)}\&RightArrow;\frac{1}{{\mathrm{\&lambda;}}_p-s},s+\frac{1}{{\mathrm{\&lambda;}}^{\left(k\right)}}\&RightArrow;{\mathrm{\&lambda;}}_p,v^{\left(k\right)}\&RightArrow;v_p$

式中，λ_p为矩阵A′的所有特征值中最靠近s的特征值，v_p为相应的特征向量。

逆迭代法具有线性收敛特性。为了提高算法的收敛速度，可用具有立方收敛特性的Rayleigh商来代替该式λ^(k)=[v^(k)]^HA′v^(k)来计算λ^(k)，从而得到逆迭代转Rayleigh商迭代算法。

Rayleigh商定义为：

$r\left(v\right)=\frac{v^H{A}^{\&prime;}v}{v^{H}v}$

其为目标函数的最优解。

$\underset{\mathrm{\&lambda;}}{\min }{{\left|\right|A}^{\&prime;}v-\mathrm{\&lambda;v}\left|\right|}_2^2=\underset{\mathrm{\&lambda;}}{\min }{(A^{\&prime;}v-\mathrm{\&lambda;v})}^H(A^{\&prime;}v-\mathrm{\&lambda;v})$

$=\underset{\mathrm{\&lambda;}}{\min }=[{\left(A^{\&prime;}v\right)}^H\left(A^{\&prime;}v\right)-\frac{{\left(v^H{A}^{\&prime;}v\right)}^H\left(v^H{A}^{\&prime;}v\right)}{v^{H}v}$

$+{(\mathrm{\&lambda;}-\frac{v^H{A}^{\&prime;}v}{v^{H}v})}^H(\mathrm{\&lambda;}-\frac{v^H{A}^{\&prime;}v}{v^{H}v})]$

于是，可得逆迭代转Rayleigh商迭代算法的迭代过程如下：

步骤101：给定2-范数单位迭代初始相量

其中，

是n维复向量空间，n是特征向量的维数，即矩阵A′的阶数。

步骤102：令k=0，根据公式v^(k+1)=z^(k+1)/||z^(k+1)||₂计算v^(k+1)。

上式中，z^(k+1)根据(A′-s^(k)I)z^(k+1)=v^(k)计算获得，s^(k)是给定的初始位移点，为复常数。A′是修正后的全系统线性化状态矩阵，I是单位矩阵，z^(k+1)是计算中间值，[·]^H是共轭转置运算，||·||₂是2-范数运算。修正后的全系统线性化状态矩阵A′中的各个状态变量的值可以通过EMS系统获得。

步骤103：判断是否达到收敛条件，如果达到收敛条件，则执行步骤104；否则，令k=k+1，返回步骤102。

收敛条件的判断准则为

其中，λ^(k)为特征值的迭代值，λ^(k)=[v^(k)]^HA′v^(k)，ε为设定误差。

步骤104：将v^(k+1)作为修正后的全系统线性化状态矩阵A′的特征向量，记为v_S，将s^(k)记为s；

步骤105：根据公式λ^(k+1)=[v^(k+1)]^HA′v^(k+1)计算特征值的迭代值λ^(k+1)，再根据公式

计算得到修正后的全系统线性化状态矩阵A′的特征值λ_S。

步骤3：根据所述特征值λ_S和特征向量v_S，计算特征值的频率f、衰减阻尼比ξ和机电回路相关比ρ。

计算特征值的频率f采用公式

其中，ω是特征值λ_S虚部的值。

计算特征值的衰减阻尼比ξ采用公式

其中，α是特征值λ_S实部的值，ω是特征值λ_S虚部的值。

计算特征值的机电回路相关比ρ采用公式其中，x_k是第k个状态变量，1≤k≤n，n是状态变量的个数。Δω是角速度变化量，Δδ是功角变化量。p_ki为相关因子，

1≤i≤n，n是修正后的全系统线性化状态矩阵A′的阶数。此处的n与系统状态变量的个数相同。因为修正后的全系统状态方程为

其中包含了所有状态变量，而状态变量有n个，因此修正后的全系统线性化状态矩阵A′的阶数为n。u_S是与特征值λ_S相对应的左特征向量，v_S是与特征值λ_S相对应的右特征向量，u_kS是u_S的第k个分量，v_kS是v_S的第k个分量。

步骤4：根据特征值的频率f、衰减阻尼比ξ和机电回路相关比ρ，对特征值进行分析，确定风电接入电网的低频振荡模式，具体为：

步骤201：选取频率为0.1Hz<f<2.5Hz、衰减阻尼比ξ<3%与机电回路相关比ρ>1的特征值。

步骤202：在选取的特征值下进行机组相位分群判定，即如果存在选取的特征值对应的特征向量集合的子集U，使得对任意i,j∈U，有

则判定在选取的该特征值下机组相位分群；如果机组相位分群且机组相位分群的个数大于等于2，则执行步骤203；否则，结束。

其中，θ_i是第i个特征向量的相位，

α_i和ω_i分别为选取的特征值对应的第i个特征向量的实部与虚部。θ_j是第j个特征向量的相位，

α_j和ω_j分别为选取的特征值对应的第j个特征向量的实部与虚部。M是用于表示特征向量相位接近程度的设定值；

步骤203：将每个选取的特征值作为一个风电接入电网的低频振荡模式。

实施例2

下面以2012年甘肃河西地区电网作为本发明的一个实例，对本发明的发明内容做进一步说明。

甘肃河西地区电网接线图如图4所示。2012年甘肃电网风电总装机容量为6915.1兆瓦，火电总装机容量为4500兆瓦。甘肃2012年10月风电出力计划为2800MW，即风电出力的最大值。

获取EMS数据，确定甘肃电网运行方式为2012年冬季小负荷，河西地区750千伏电网全停检修方式。潮流断面潮流情况为新疆外送0MW，南电北送3807MW，西电东送2607MW，风机出力为800MW。

计算使用软件为电力系统分析综合程序(PSASP)6.282版本，风机模型采用PSASP自带的双馈直驱风力发电机通用模型，由此建立大规模风电接入的数学模型。

计算采用PSASP的小干扰稳定计算程序，使用方法为逆迭代转Rayleigh商迭代法，其特征值搜索范围为0.1Hz-2.5Hz，阻尼比搜索范围为0%-5%。小干扰稳定计算结果共得到182个特征值。筛选频率为0.1Hz<f<2.5Hz、衰减阻尼比ξ<3%与机电回路相关比ρ>1的特征值，筛选得到59个特征值，进行低频振荡振型图自动分析，绘制振型图。

以图5为例说明低频振荡振型图自动分析。如图5（a）图，机组呈现相位分群，因此确定图5（a）所对应的特征值为一个振荡模式，放入低频振荡模式集。如图5（b）图，振型图没有明显的相位分群，表明该特征值仅是系统状态方程的数学解，并不具有实际的物理意义，图5（b）所对应的特征值不是一个振荡模式。

经过低频振荡振型图自动分析，确定当前电网的低频振荡模式如下：

模式序号	模式	频率/Hz	阻尼比/%
				1	-0.02811+j2.894975	0.460750	0.9709
2	-0.19659+j8.350194	1.328975	2.3537

以上结果表明，在2012年冬季小负荷，河西地区750千伏电网全停检修方式下，风电接入800MW时，河西地区可能存在频率为0.46Hz与1.32Hz的弱阻尼低频振荡模式。

上述实例分析表明：本方法对大规模风电接入后的对大电网低频振荡模式进行了分析，仿真辨识出甘肃河西地区的低频振荡模式。本方法可以考虑大规模风电接入对低频振荡的影响，并能对低频振荡振型图进行自动分析，能够准确的辨识电网中的振荡模式，为电力系统的低频振荡特性分析提供了指导。

以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应该以权利要求的保护范围为准。

Claims (9)

1.一种风电接入电网的低频振荡辨识方法，其特征是所述方法包括：

步骤1：对不含风电机组的电力系统全系统线性化状态矩阵A进行修正，得到包含风电机组的电力系统全系统线性化状态矩阵A′，并记为修正后的全系统线性化状态矩阵A′；

步骤2：求取修正后的全系统线性化状态矩阵A′的特征值λ_S和特征向量v_S；

步骤3：根据所述特征值λ_S和特征向量v_S，计算特征值的频率f、衰减阻尼比ξ和机电回路相关比ρ；

步骤4：根据特征值的频率f、衰减阻尼比ξ和机电回路相关比ρ，对特征值进行分析，确定风电接入电网的低频振荡模式。

2.根据权利要求1所述的低频振荡辨识方法，其特征是所述风电机组为双馈异步风力发电机组。

3.根据权利要求2所述的低频振荡辨识方法，其特征是所述修正后的全系统线性化状态矩阵A′=A′_D+B′_D(Y′_N+Y′_D)^-1C′_D；

其中， $A_D^{\&prime;}=\left[\begin{array}{cc}{A}_D& \\ & A_{\mathrm{wf}}\end{array}\right];$

$A_{\mathrm{wf}}=\left[\begin{array}{ccccccc}-\frac{{\mathrm{\&omega;}}_s}{T_r^{\&prime;}}& 0& -\frac{{\mathrm{\&omega;}}_s}{X_r^{\&prime;}}& 0& 0& 0& 0\\ K_1{U}_t\frac{X_m}{X_s}(\frac{{\mathrm{\&omega;}}_s}{T_r^{\&prime;}}-\frac{1}{T_1})& 0& K_1{U}_t\frac{X_m}{X_s}\frac{{\mathrm{\&omega;}}_s}{X_r^{\&prime;}}& 0& 0& 0& 0\\ K_2(\frac{1}{T_2}-\frac{{\mathrm{\&omega;}}_s}{T_r^{\&prime;}})& -\frac{K_2}{T_2}& -K_2\frac{{\mathrm{\&omega;}}_s}{X_r^{\&prime;}}& 0& 0& 0& 0\\ 0& 0& 0& -\frac{{\mathrm{\&omega;}}_s}{T_r^{\&prime;}}& 0& -\frac{{\mathrm{\&omega;}}_s}{X_r^{\&prime;}}& 0\\ 0& 0& 0& K_1{U}_t\frac{X_m}{X_s}(\frac{{\mathrm{\&omega;}}_s}{T_r^{\&prime;}}-\frac{1}{T_1})& 0& K_1{U}_t\frac{X_m}{X_s}\frac{{\mathrm{\&omega;}}_s}{X_r^{\&prime;}}& 0\\ 0& 0& 0& K_2(\frac{1}{T_2}-\frac{{\mathrm{\&omega;}}_s}{T_r^{\&prime;}})& -\frac{K_2}{T_2}& -K_2\frac{{\mathrm{\&omega;}}_s}{X_r^{\&prime;}}& 0\\ 0& 0& 0& -U_t\frac{X_m}{2H{X}_s}& 0& 0& 0\end{array}\right];$

$B_D^{\&prime;}=\left[\begin{array}{cc}{B}_D& \\ & B_{\mathrm{wf}}\end{array}\right];$

B_wf=[0]7×2；

Y′_N是修正后的全系统节点导納矩阵；

$C_D^{\&prime;}=\left[\begin{array}{cc}{C}_D& \\ & C_{\mathrm{wf}}\end{array}\right];$

$C_{\mathrm{wf}}={\left[\begin{array}{cc}\sin \mathrm{\&delta;}& -\cos \mathrm{\&delta;}\\ \cos \mathrm{\&delta;}& \sin \mathrm{\&delta;}\end{array}\right]}^{-1}\left[\begin{array}{ccccccc}-\frac{X_s}{X_m}& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& -\frac{X_s}{X_m}& 0& 0& 0\end{array}\right];$

$Y_D^{\&prime;}=\left[\begin{array}{cc}{Y}_D& \\ & Y_{\mathrm{wf}}\end{array}\right];$

Y_wf是风电机组对应的节点导纳矩阵；

A_D、B_D、C_D和Y_D是修正前的系统线性化状态矩阵A的参数矩阵；

矩阵A_wf中的ω_s是同步角速度；

矩阵A_wf和矩阵C_wf中的X_s是定子电抗；

矩阵A_wf和矩阵C_wf中的X_m是定子与转子间的互抗；

矩阵A_wf中的X_r是转子电抗；

矩阵A_wf中的K₁是外环比例积分微分控制器的第一参数；

矩阵A_wf中的U_t是发电机机端电压；

矩阵A_wf中的T′_r=X'_r/R_r，R_r是转子电阻；

矩阵A_wf中的T₁是外环比例积分微分控制器的第二参数；

矩阵A_wf中的K₂是内环比例积分微分控制器的第一常数；

矩阵A_wf中的T₂是内环比例积分微分控制器的第二常数；

矩阵A_wf中的H是惯性常数；

矩阵C_wf中的δ是发电机功角。

4.根据权利要求2所述的低频振荡辨识方法，其特征是所述求取修正后的全系统线性化状态矩阵A′的特征值λ_S和特征向量v_S包括：

步骤101：给定2-范数单位迭代初始相量

其中，

是n维复向量空间，n是特征向量的维数；

步骤102：令k=0，根据公式v^(k+1)=z^(k+1)/||z^(k+1)||₂计算v^(k+1)；

其中，z^(k+1)根据(A′-s^(k)I)z^(k+1)=v^(k)计算获得，

s^(k)是给定的初始位移点，为复常数；

A′是修正后的全系统线性化状态矩阵；

I是单位矩阵；

z^(k+1)是计算中间值；

[·]^H是共轭转置运算；

||·||₂是2-范数运算；

步骤103：判断是否达到收敛条件，如果达到收敛条件，则执行步骤104；否则，令k=k+1，返回步骤102；

步骤104：将v^(k+1)作为修正后的全系统线性化状态矩阵A′的特征向量，记为v_S，将s^(k)记为s；

步骤105：根据公式λ^(k+1)=[v^(k+1)]^HA′v^(k+1)计算特征值的迭代值λ^(k+1)，再根据公式计算得到修正后的全系统线性化状态矩阵A′的特征值λ_S。

5.根据权利要求4所述的低频振荡辨识方法，其特征是所述收敛条件为其中，λ^(k)为特征值的迭代值，λ^(k)=[v^(k)]^HA′v^(k)，ε为设定误差。

6.根据权利要求4所述的低频振荡辨识方法，其特征是所述计算特征值的频率f采用公式其中，ω是特征值λ_S虚部的值。

7.根据权利要求4所述的低频振荡辨识方法，其特征是所述计算特征值的衰减阻尼比ξ采用公式

其中，α是特征值λ_S实部的值，ω是特征值λ_S虚部的值。

8.根据权利要求4所述的低频振荡辨识方法，其特征是所述计算特征值的机电回路相关比ρ采用公式

其中，x_k是第k个状态变量，1≤k≤n，n是状态变量的个数；

Δω是角速度变化量；

Δδ是功角变化量；

p_ki为相关因子，

1≤i≤n，n为修正后的全系统线性化状态矩阵A′的阶数；

u_S是与特征值λ_S相对应的左特征向量；

v_S是与特征值λ_S相对应的右特征向量；

u_kS是u_S的第k个分量；

v_kS是v_S的第k个分量。

9.根据权利要求4所述的低频振荡辨识方法，其特征是所述步骤4具体包括：

步骤201：选取频率为0.1Hz<f<2.5Hz、衰减阻尼比ξ<3%与机电回路相关比ρ>1的特征值；

步骤202：在选取的特征值下进行机组相位分群判定，即如果存在选取的特征值对应的特征向量集合的子集U，使得对任意i,j∈U，有

则判定在选取的该特征值下机组相位分群；如果机组相位分群且机组相位分群的个数大于等于2，则执行步骤203；否则，结束；

其中，θ_i是第i个特征向量的相位，

α_i和ω_i分别为选取的特征值对应的第i个特征向量的实部与虚部；

θ_j是第j个特征向量的相位，

α_j和ω_j分别为选取的特征值对应的第j个特征向量的实部与虚部；

M是用于表示特征向量相位接近程度的设定值；

步骤203：将每个选取的特征值作为一个风电接入电网的低频振荡模式。

Images