CN105140957A - 基于风电场和光伏电站聚合模型的机电振荡模式估算方法 - Google Patents

Abstract

本发明属于新能源电力系统动态稳定性分析技术领域，尤其涉及一种基于大规模风电场和光伏电站聚合模型的电力系统机电振荡模式估算方法，首先建立单台直驱风机或光伏发电单元的传递函数模型，整个风电场和光伏电站的聚合模型可以通过单台发电设备的模型直接倍乘，或相加后降阶的方法得到；对于参数未知的风电场或光伏电站，聚合的传递函数模型通过现场频率响应试验得到，电力系统机电振荡模式由系统的开环特征值加上聚合传递函数值与残差的乘积来估算。本发明能够在保证一定准确度的情况下，大大降低模型阶数和计算复杂度，从而为新能源发电电源接入后的电力系统小干扰振荡稳定性分析提供帮助。

Description

基于风电场和光伏电站聚合模型的机电振荡模式估算方法

技术领域

本发明属于新能源电力系统动态稳定性分析技术领域，尤其涉及一种基于大规模风电场和光伏电站聚合模型的电力系统机电振荡模式估算方法。

背景技术

近年来，随着风电机组技术和光伏发电技术的不断进步、风电场和光伏电站规模的不断扩大、区域电网互联的飞速发展，新能源发电电源的接入对电网稳定性的影响越来越大。在这种发展趋势下，全国各地区电网公司对大规模新能源发电电源接入电网所带来的影响，尤其是小干扰稳定性问题非常重视。因此，对于大规模风电场和光伏电站接入电力系统后进行小干扰稳定性及阻尼特性的研究成为了当今最迫切的课题。

目前，在分析电力系统小干扰稳定性问题时，特征值分析法得到了广泛应用。特征值分析法是有效解决计算电力系统机电振荡模式的定量信息问题最有效的手段。它是借助线性系统的稳定性理论来研究非线性系统的稳定性。该方法严格建立在现代控制理论基础上，把电力系统视为一般控制系统，用标准的线性化状态方程来描述。系统的振荡模式完全由状态矩阵的特征根的对数来决定。特征根的实部和虚部分别对应振荡模式的阻尼和频率，而特征根对应的特征向量反映了振荡模式在整个系统中的动态行为。因此，不论系统有多么复杂，都可以采用特征值分析法对系统中任一振荡模式做详细分析，而与其它振荡模式无关。特征值分析法不仅可以用来刻画系统的振荡模式，而且可以用来确定控制器最佳安装地点。

当大规模风电场或光伏电站接入电力系统后，使用特征值分析法分析电力系统机电振荡模式时，一个最大的障碍就是状态矩阵阶数的增加大大增加了计算特征值的复杂度和计算时间。因此对可降低系统阶数的风电场和光伏电站的聚合模型的研究成为新能源电力系统小干扰稳定性分析中很重要的工作。

风电场的聚合模型是指将风电场等效成为单台或几台机组。显然，当风电场包含几百台风电机组时，采用详细模型会使计算规模过于庞大，计算能力往往不能满足要求。因此，当研究风电场对电网的影响时，可将风电场看作一个整体进行聚合，以减小仿真规模，提高计算效率和收敛性。目前，建立风电场聚合模型的方法可以归纳为3种：

1)单机法(完全聚合法)。将风电场等值成一台风力发电机组。等值风电机组的容量等于所有风电机组容量的代数和。其余参数都由等值计算获得。

2)复合模型法。将风电场电气部分等值为一台发电机，而空气动力学部分、传动链和最大功率跟踪模型则根据运行条件分区分别聚合成几台并联，构成一台复合的风电机组。

3)分群法(多机法)。按一定分群原则，将风电场等值成为多台机组。

上述三种方法中，单机法简化程度最大，但往往被认为精度较低；而复合模型法和分群法是为了体现风场中机组运行状态的差异性对输出特性带来的影响，等效机组过多时，简化程度不够。而且上述方法都用于新能源接入后电力系统的动态分析中，而用在小干扰分析中时不免“大材小用”，因为它们都提供了一些无用于小干扰分析的状态量与参数。

发明内容

为了在新能源电源接入电力系统后，建立用于新能源电力系统小干扰稳定分析的风电场和光伏电站聚合模型，以降低特征值分析中的模型阶数和计算复杂度，本发明提出了一种基于大规模风电场和光伏电站聚合模型的电力系统机电振荡模式估算方法，包括：

步骤1、建立大规模风电场和光伏电站的聚合传递函数模型；

步骤2、将电力系统中风电场和光伏电站的输出电流作为控制量，建立开环的多机电力系统的线性化模型；

步骤3、把风电场或光伏电站的传递函数模型看作是控制器，作为连接开环系统输入与输出的反馈回路；开环系统和反馈回路共同构成了多机电力系统的闭环状态空间模型；

步骤4、将风电场和光伏电站看作恒电流源，即令风电场或光伏电站输出电流为零，忽略掉风电场和光伏电站与多机电力系统间的动态交互影响，计算开环系统的特征值；通过系统的开环特征值加上聚合传递函数值与残差的乘积来估算电力系统闭环的机电振荡模式。

所述步骤1若在各台直驱风力发电机或光伏发电单元的型号、参数和运行工况都已知的情况下，则建立的各台直驱风机和光伏发电单元的传递函数模型为：

$\begin{array}{c}{\mathrm{\ΔI}}_{wdj}^k=b_1^k\left(s\right)\mathrm{\Δ}\left|V_{wj}\right|\\ {\mathrm{\ΔI}}_{wqj}^k=b_2^k\left(s\right)\mathrm{\Δ}\left|V_{wj}\right|+b_3^k\left(s\right){\mathrm{\Δu}}_j^k\end{array}---\left(23\right)$

其中，和分别为第j个风电场或光伏电站内第k台直驱风机或光伏发电单元输出电流的d轴和q轴分量，|V_wj|为第j个风电场或光伏电站内低压母线的电压幅值，为第j个风电场或光伏电站内第k个直驱风机或光伏发电单元的外部变量；对于直驱风机，外部变量为机械转矩，对于光伏发电单元，外部变量为光照强度；和分别为所对应变量间的传递函数。

所述步骤1若在一个风电场或光伏电站内所有直驱风机或光伏发电单元的型号、参数和运行工况都相同的情况下，那么该风电场或光伏电站的聚合传递函数模型由单台风机或光伏发电单元的传递函数倍乘直接得到：

$\begin{array}{c}{G}_{dj}\left(s\right)=\frac{{\mathrm{\ΔI}}_{wdj}}{\mathrm{\Δ}\left|V_{wj}\right|}=L_j{b}_1^k\left(s\right)\\ G_{qj}\left(s\right)=\frac{{\mathrm{\ΔI}}_{wqj}}{\mathrm{\Δ}\left|V_{wj}\right|}=L_j{b}_2^k\left(s\right)\end{array}---\left(26\right)$

若一个风电场或光伏电站内直驱风机或光伏发电单元的型号、参数或运行工况不同，那么整个风电场或光伏电站的聚合传递函数模型由所有风机或光伏发电单元的传递函数相加得到：

$\begin{array}{c}{G}_{dj}\left(s\right)=\frac{{\mathrm{\ΔI}}_{wdj}}{\mathrm{\Δ}\left|V_{wj}\right|}=\underset{k=1}{\overset{L}{\Σ}}{b}_1^k\left(s\right)\\ G_{qj}\left(s\right)=\frac{{\mathrm{\ΔI}}_{wqj}}{\mathrm{\Δ}\left|V_{wj}\right|}=\underset{k=1}{\overset{L}{\Σ}}{b}_2^k\left(s\right)\end{array}---\left(25\right)$

其中，I_wdj和I_wqj分别为第j个风电场或光伏电站内总输出电流的d轴和q轴分量，L_j为第j个风电场或光伏电站内所包含直驱风机或光伏发电单元的数量，G_dj(s)和G_qj(s)分别为第j个风电场或光伏电站内输出量ΔI_wdj和ΔI_wqj对输入量Δ|V_wj|的聚合传递函数；

对于式(25)，传递函数G_dj(s)和G_qj(s)的阶数会很高，通过目前已有的模型降阶方法，将高阶的传递函数模型降阶成低阶的传递函数模型。

所述步骤1中在直驱风机或光伏发电单元的型号、参数或运行工况未知的情况下，通过现场频率响应试验来得到风电场和光伏电站的聚合传递函数模型；通过在低压母线稳态电压幅值上叠加某频率的正弦信号作为输入，检测风电场或光伏电站的输出电流，得到该频率下的输出电流对输入电压的传递函数值；对各个频率下传递函数值的幅值图和相位图，进行不同阶数的拟合，从而近似得到风电场和光伏电站的聚合模型。

所述步骤2中建立开环的多机电力系统的线性化模型为：

$\begin{array}{c}{\mathrm{s\ΔX}}_g=A_0{\mathrm{\ΔX}}_g+\underset{j=1}{\overset{M}{\Σ}}(b_{0dj}{\mathrm{\ΔI}}_{wdj}+b_{0qj}{\mathrm{\ΔI}}_{wqj})\\ \mathrm{\Δ}\left|V_{wj}\right|=c_{0Vj}^T{\mathrm{\ΔX}}_g+\underset{j=1}{\overset{M}{\Σ}}(d_{0dj}{\mathrm{\ΔI}}_{wdj}+d_{0qj}{\mathrm{\ΔI}}_{wqj})\end{array}---\left(44\right)$

其中，X_g为系统状态变量，A₀为开环系统的状态矩阵，b_0dj和b_0qj为控制矩阵，为输出矩阵，d_0dj和d_0qj为反馈系数，M为多机电力系统所包含风电场和光伏电站的数量，I_wdj和I_wqj分别为第j个风电场或光伏电站总输出电流的d轴和q轴分量，|V_wj|为第j个风电场或光伏电站内低压母线的电压幅值。

所述步骤4中电力系统机电振荡模式为：

${\stackrel{\‾}{\mathrm{\λ}}}_{ci}\≈\stackrel{\‾}{{\mathrm{\λ}}_i}+\frac{\underset{j=1}{\overset{M}{\Σ}}\[{\stackrel{\‾}{R}}_{dji}{G}_{dj}\left(\stackrel{\‾}{{\mathrm{\λ}}_i}\right)+{\stackrel{\‾}{R}}_{qji}{G}_{qj}\left(\stackrel{\‾}{{\mathrm{\λ}}_i}\right)\]}{1-\underset{j=1}{\overset{M}{\Σ}}\[d_{0dj}{G}_{dj}\left(\stackrel{\‾}{{\mathrm{\λ}}_i}\right)+d_{0qj}{G}_{qj}\left(\stackrel{\‾}{{\mathrm{\λ}}_i}\right)\]}---\left(54\right)$

残差为：

$\begin{array}{c}{\stackrel{\‾}{R}}_{dji}={\stackrel{\‾}{p}}_i^T{b}_{0dj}{c_{oVj}}^T\stackrel{\‾}{r_i}\\ {\stackrel{\‾}{R}}_{qji}={\stackrel{\‾}{p}}_i^T{b}_{0qj}{c_{oVj}}^T\stackrel{\‾}{r_i}\end{array}---\left(55\right)$

其中，为开环系统的状态矩阵的特征值，和分别为对应的左右特征向量，为考虑了风电场和光伏电站动态后闭环阵的特征值估计值，G_dj(s)和G_qj(s)分别为第j个风电场或光伏电站的输出量ΔI_wdj和ΔI_wqj对输入量Δ|V_wj|的聚合传递函数，b_0dj和b_0qj为控制矩阵，为输出矩阵，d_0dj和d_0qj为反馈系数，M为多机电力系统所包含风电场和光伏电站的数量。

本发明的有益效果在于：用本发明中提出的方法计算大规模风电场和光伏电站接入电力系统后的机电振荡模式时，不需要推导闭环系统的状态矩阵。只需要得出把风电场和光伏电站等效为恒电流源后的开环状态矩阵的特征值，再计算残差，即可近似估算闭环系统的机电振荡模式。这种方法可以大大降低矩阵的阶数，减小模态计算中计算负担和复杂度。

附图说明

图1为直驱风机的结构图；

图2为直驱风机机侧换流器的控制方案；

图3为直驱风机网侧换流器的控制方案；

图4为光伏发电单元的结构图；

图5为光伏发电单元DC/DC换流器的控制方案；

图6为大规模风电场和光伏电站接入后的多机电力系统的结构图；

图7为x-y坐标系与d-q坐标系的旋转关系；

图8为大规模风电场和光伏电站接入后的多机电力系统的线性化模型示意图；

图9为纽约-新英格兰互联电力系统的结构图；

图10为连接在母线48处的光伏电站的频率响应曲线。

具体实施方式

下面结合附图，对实施例作详细说明。

本发明提出了一种基于大规模风电场和光伏电站聚合模型的电力系统机电振荡模式估算方法，包括：

步骤1、建立大规模风电场和光伏电站的聚合传递函数模型；

步骤2、每台直驱风机或光伏发电单元通过变压器连接到低压母线上，然后再通过升压变压器连接到公共连接点处的高压母线上；建立开环的多机电力系统的线性化模型；

步骤3、把不考虑风电场和光伏电站动态的多机电力系统的状态空间模型，看作是开环的多机电力系统模型；风电场或光伏电站的传递函数模型看作是控制器，作为连接开环系统输入与输出的反馈回路；开环系统和反馈回路共同构成了多机电力系统的闭环状态空间模型；

步骤4、令风电场或光伏电站输出电流的d轴和q轴分量均为零，忽略掉风电场和光伏电站与多机电力系统间的动态交互影响，计算开环系统、反馈回路、闭环系统的特征值，并估算电力系统机电振荡模式。

所述步骤1，建立大规模风电场和光伏电站的聚合传递函数模型，具体包括下面两种情况：

(1)在各台直驱风力发电机或光伏发电单元的型号、参数和运行工况都已知的情况下，首先建立各台直驱风机和光伏发电单元的传递函数模型。

直驱风机的结构如附图1所示。

在d-q坐标系下，直驱风机的定子电压方程、定子磁链方程和转子运动方程分别为：

$\begin{array}{c}{V}_{pmd}=-R_a{I}_{pmd}-\frac{1}{{\mathrm{\ω}}_0}\frac{{\mathrm{d\ψ}}_{pmd}}{dt}+{\mathrm{\ω}}_{pm}{\mathrm{\ψ}}_{pmq}\\ V_{pmq}=-R_a{I}_{pmq}-\frac{1}{{\mathrm{\ω}}_0}\frac{{\mathrm{d\ψ}}_{pmq}}{dt}-{\mathrm{\ω}}_{pm}{\mathrm{\ψ}}_{pmd}\end{array}---\left(1\right)$

ψ_pmd＝X_dI_pmd-λ₀

(2)

ψ_pmq＝X_qI_pmq

$\frac{J_{pm}}{{\mathrm{\ω}}_0}\frac{{\mathrm{d\ω}}_{pm}}{dt}=T_{wm}-T_{pme}---\left(3\right)$

其中，ω₀为同步角速度，λ₀是永磁体的恒定磁链，R_a、X_z、ψ_pmz、V_pmz和I_pmz(z＝d,q)分别是定子电阻、定子漏抗、定子磁通、定子电压和定子电流，J_pm为直驱风机的惯性常数，ω_pm为直驱风机的转子转速，T_wm和T_pme分别为机械转矩和电磁转矩，其中：

T_pme＝ψ_pmqI_pmd-ψ_pmdI_pmq(4)

忽略定子电阻，同时假设X_d＝X_q＝X_pm，并把式(2)代入式(1)，可得：

V_pmd＝ω_pmX_pmI_pmq

(5)

V_pmq＝-ω_pm(X_pmI_pmd-λ₀)

直驱风机的机侧换流器的控制策略如附图2所示。忽略其中的电磁暂态快速过程，那么控制策略中的内环控制可以忽略，由附图2的外环控制可得：

I_pmq＝I_pmq ^ref＝K_pm(s)(ω_pm-ω_pm ^ref)

(6)

I_pmd＝I_pmd ^ref＝I_pmd0

在网侧换流器的控制中，将低压母线电压的方向定为d-q坐标系的q轴方向。由附图1可得：

$\begin{array}{c}\frac{X_w}{{\mathrm{\ω}}_0}\frac{{\mathrm{dI}}_{wd}}{dt}=V_{wd}-V_{cd}+X_w{I}_{wq}\\ \frac{X_w}{{\mathrm{\ω}}_0}\frac{{\mathrm{dI}}_{wq}}{dt}=V_{wq}-V_{cq}-X_w{I}_{wd}\end{array}---\left(7\right)$

直驱风机的网侧换流器的控制策略如附图3所示。同样忽略其中的电磁暂态快速过程，那么控制策略中的内环控制可以忽略，由附图3的外环控制可得：

$\begin{array}{c}{I}_{wq}=I_{wq}^{ref}=K_{dc}\left(s\right)({V_{dc}}^{ref}-V_{dc})\\ I_{wd}=I_{wd}^{ref}=K_q\left(s\right)({Q_w}^{ref}-Q_w)\end{array}---\left(8\right)$

电容器的动态方程为：

$\frac{C_{dc}{V}_{dc}}{{\mathrm{\ω}}_0}\frac{{\mathrm{dV}}_{dc}}{dt}=P_{pme}+P_w---\left(9\right)$

其中，

P_pme＝I_pmdV_pmd+I_pmqV_pmq

P_w＝|V_w|I_wq(10)

Q_w＝|V_w|I_wd

对式(1)-(6)进行线性化，可得：

ΔI_pmq＝K_pm(s)Δω_pm

J_pmsΔω_pm＝ΔT_wm-λ₀ΔI_pmq(11)

ΔP_pme＝λ₀ω_pm0ΔI_pmq+λ₀I_pmq0Δω_pm

对式(8)-(10)进行线性化，可得：

ΔI_wq＝-K_dc(s)ΔV_dc

ΔI_wd＝-K_q(s)ΔQ_w

$\frac{{\mathrm{d\ΔV}}_{dc}}{dt}=\frac{{\mathrm{\ω}}_0}{C_{dc}{V}_{dc0}}({\mathrm{\ΔP}}_{pme}+{\mathrm{\ΔP}}_w)---\left(12\right)$

ΔP_w＝I_wq0Δ|V_w|+|V_w0|ΔI_wq

ΔQ_w＝I_wd0Δ|V_w|+|V_w0|ΔI_wd

由式(11)(12)可得：

ΔP_w＝a₁(s)Δ|V_w|+a₂(s)ΔT_wm

ΔQ_w＝a₃(s)Δ|V_w|(13)

ΔV_dc＝a₄(s)Δ|V_w|+a₅(s)ΔT_wm

其中，

$a_1\left(s\right)=\frac{{\mathrm{sC}}_{dc}{V}_{dc0}{I}_{wq0}}{{\mathrm{sC}}_{dc}{V}_{dc0}+{\mathrm{\ω}}_0{K}_{dc}\left(s\right)\left|V_{w0}\right|}$

$a_2\left(s\right)=\frac{-{\mathrm{\ω}}_0{K}_{dc}\left(s\right)\left|V_{w0}\right|}{{\mathrm{sC}}_{dc}{V}_{dc0}+{\mathrm{\ω}}_0{K}_{dc}\left(s\right)\left|V_{w0}\right|}\frac{{\mathrm{\λ}}_0{\mathrm{\ω}}_{pm0}{K}_{pm}\left(s\right)+{\mathrm{\λ}}_0{I}_{pmq0}}{J_{pm}s+{\mathrm{\λ}}_0{K}_{pm}\left(s\right)}$

$a_3\left(s\right)=\frac{I_{wd0}}{1+\left|V_{w0}\right|K_q\left(s\right)}$

$a_4\left(s\right)=\frac{{\mathrm{\ω}}_0}{{\mathrm{sC}}_{dc}{V}_{dc0}}{a}_1$

$a_5\left(s\right)=\frac{{\mathrm{\ω}}_0}{{\mathrm{sC}}_{dc}{V}_{dc0}}(a_2+\frac{{\mathrm{\λ}}_0{\mathrm{\ω}}_{pm0}{K}_{pm}\left(s\right)+{\mathrm{\λ}}_0{I}_{pmq0}}{J_{pm}s+{\mathrm{\λ}}_0{K}_{pm}\left(s\right)})$

将式(13)代入式(12)中，可得：

ΔI_wd＝b₁(s)|ΔV_w|

(14)

ΔI_wq＝b₂(s)|ΔV_w|+b₃(s)ΔT_wm

其中，

b₁(s)＝-K_q(s)a₃(s)

b₂(s)＝-K_dc(s)a₄(s)

b₃(s)＝-K_dc(s)a₅(s)

光伏发电单元的结构如附图4所示。

光伏发电单元的电压-电流(V_pv-I_pv)特性和出力分别为：

$V_{pv}=\frac{N_{s}nkT}{q}ln(\frac{\frac{N_p{I}_{sc}{I}_r}{100}-I_{pv}}{N_p{I}_0}+1)---\left(15\right)$

P_pv＝I_pvV_pv(16)

其中，T为光伏板工作温度，I_r为光照强度，N_s和N_p分别是串并联光伏板的单元数量，n为理想化因子，k为波尔兹曼常数，q为电子的电荷量，I_sc为短路电流，I₀为饱和电流。

附图4中左侧电感的动态方程为：

$\frac{{\mathrm{dI}}_{pv}}{dt}=\frac{{\mathrm{\ω}}_0}{L_{dc}}(V_{pv}-V_{dc1})---\left(17\right)$

光伏发电单元DC/DC换流器的控制策略如附图5所示，由附图5可得：

$\begin{array}{c}{I}_{pv}^{ref}=K_{pv1}\left(s\right)(P_{pv}^{ref}-P_{pv})\\ V_{dc1}=-K_{pv2}\left(s\right)(I_{pv}^{ref}-I_{pv})+V_{pv}\end{array}---\left(18\right)$

对式(15)-(18)进行线性化，可得：

ΔV_pv＝a_IΔI_pv+a_rΔI_r

ΔP_pv＝V_pv0ΔI_pv+I_pv0ΔV_pv

${\mathrm{s\ΔI}}_{pv}=\frac{{\mathrm{\ω}}_0}{L_{dc}}({\mathrm{\ΔV}}_{pv}-{\mathrm{\ΔV}}_{dc1})---\left(19\right)$

${\mathrm{\ΔI}}_{pv}^{ref}=-K_{pv1}\left(s\right){\mathrm{\ΔP}}_{pv}$

ΔV_dc1＝K_pv2(s)ΔI_pv+ΔV_pv

其中，

$a_I=\frac{N_{s}nkT}{q}\frac{-100}{N_p{I}_{sc}{I}_r-100I_{pv}+100N_p{I}_0}$

$a_r=\frac{N_{s}nkT}{q}\frac{N_p{I}_{sc}}{N_p{I}_{sc}{I}_r-100I_{pv}+100N_p{I}_0}$

继而有：

ΔV_dc1＝a_VIΔI_pv+a_VrΔI_r

${\mathrm{\ΔI}}_{pv}=\frac{{\mathrm{\ω}}_0(a_r-a_{Vr})}{{\mathrm{sL}}_{dc}-{\mathrm{\ω}}_0(a_I-a_{VI})}{\mathrm{\ΔI}}_r$

${\mathrm{\ΔV}}_{pv}=\[a_r+\frac{a_I{\mathrm{\ω}}_0(a_r-a_{Vr})}{{\mathrm{sL}}_{dc}-{\mathrm{\ω}}_0(a_I-a_{VI})}\]{\mathrm{\ΔI}}_r---\left(20\right)$

${\mathrm{\ΔP}}_{pv}=\[\frac{V_{pv0}{\mathrm{\ω}}_0(a_r-a_{Vr})}{{\mathrm{sL}}_{dc}-{\mathrm{\ω}}_0(a_I-a_{VI})}+I_{pv0}{a}_r+\frac{I_{pv0}{a}_I{\mathrm{\ω}}_0(a_r-a_{Vr})}{{\mathrm{sL}}_{dc}-{\mathrm{\ω}}_0(a_I-a_{VI})}\]{\mathrm{\ΔI}}_r$

其中，

a_VI＝K_pv2(s)+K_pv1(s)K_pv2(s)V_pv0+(1+K_pv1(s)K_pv2(s)V_pv0)a_I

a_Vr＝(1+K_pv1(s)K_pv2(s)I_pv0)a_r

光伏发电单元的DC/AC换流器的控制策略与直驱风机的网侧换流器的控制策略相同，如附图3所示。其线性化模型的推导见式(12)-(14)，只需要将其中的ΔP_pme和ΔT_wm替换为ΔP_pv和ΔI_r即可：

ΔP_w＝a₁(s)Δ|V_w|+a₂(s)ΔI_r

ΔQ_w＝a₃(s)Δ|V_w|(21)

ΔV_dc＝a₄(s)Δ|V_w|+a₅(s)ΔI_r

其中，a₁(s)、a₃(s)和a₄(s)与式(13)中的相同，其余参数为：

$a_2\left(s\right)=\frac{-{\mathrm{\ω}}_0{K}_{dc}\left(s\right)\left|V_{w0}\right|}{{\mathrm{sC}}_{dc}{V}_{dc0}+{\mathrm{\ω}}_0{K}_{dc}\left(s\right)\left|V_{w0}\right|}(\frac{V_{pv0}{\mathrm{\ω}}_0(a_r-a_{Vr})}{{\mathrm{sL}}_{dc}-{\mathrm{\ω}}_0(a_I-a_{VI})}+I_{pv0}{a}_r+\frac{I_{pv0}{a}_I{\mathrm{\ω}}_0(a_r-a_{Vr})}{{\mathrm{sL}}_{dc}-{\mathrm{\ω}}_0(a_I-a_{VI})})$

$a_5\left(s\right)=\frac{{\mathrm{\ω}}_0}{{\mathrm{sC}}_{dc}{V}_{dc0}}(a_2+\frac{V_{pv0}{\mathrm{\ω}}_0(a_r-a_{Vr})}{{\mathrm{sL}}_{dc}-{\mathrm{\ω}}_0(a_I-a_{VI})}+I_{pv0}{a}_r+\frac{I_{pv0}{a}_I{\mathrm{\ω}}_0(a_r-a_{Vr})}{{\mathrm{sL}}_{dc}-{\mathrm{\ω}}_0(a_I-a_{VI})})$

将式(21)代入到式(12)中，可得：

ΔI_wd＝b₁(s)|ΔV_w|

(22)

ΔI_wq＝b₂(s)|ΔV_w|+b₃(s)ΔI_r

其中，b₁(s)、b₂(s)和b₃(s)与式(14)中相同。

自此，对于如附图6所示的多机电力系统，在第j个风电场或光伏电站内第k个直驱风机或光伏发电单元的输出电流与低压母线电压幅值之间的关系可以统一表示为：

$\begin{array}{c}{\mathrm{\ΔI}}_{wdj}^k=b_1^k\left(s\right)\mathrm{\Δ}\left|V_{wj}\right|\\ {\mathrm{\ΔI}}_{wqj}^k=b_2^k\left(s\right)\mathrm{\Δ}\left|V_{wj}\right|+b_3^k\left(s\right){\mathrm{\Δu}}_j^k\end{array}---\left(23\right)$

其中，对于直驱风机，外部变量Δu为机械转矩ΔT_wm；对于光伏发电单元，外部变量Δu为光照强度ΔI_r。

第j个风电场或光伏电站的总输出电流为：

$\begin{array}{c}{\mathrm{\ΔI}}_{wdj}=\underset{k=1}{\overset{L_j}{\Σ}}{\mathrm{\ΔI}}_{wdj}^k=G_{dj}\left(s\right)\mathrm{\Δ}\left|V_{wj}\right|\\ {\mathrm{\ΔI}}_{wqj}=\underset{k=1}{\overset{L_j}{\Σ}}{\mathrm{\ΔI}}_{wqj}^k=G_{qj}\left(s\right)\mathrm{\Δ}\left|V_{wj}\right|+\underset{k=1}{\overset{L_j}{\Σ}}{b}_3^k\left(s\right){\mathrm{\Δu}}_j^k\end{array}---\left(24\right)$

G_dj(s)和G_qj(s)分别为第j个风电场或光伏电站内，输出量ΔI_wdj和ΔI_wqj对输入量Δ|V_wj|的聚合传递函数：

$\begin{array}{c}{G}_{dj}\left(s\right)=\frac{\mathrm{\Δ}\left|V_{wj}\right|}{{\mathrm{\ΔI}}_{wdj}}=\underset{k=1}{\overset{L_j}{\Σ}}{b}_1^k\left(s\right)\\ G_{qj}\left(s\right)=\frac{\mathrm{\Δ}\left|V_{wj}\right|}{{\mathrm{\ΔI}}_{wqj}}=\underset{k=1}{\overset{L_j}{\Σ}}{b}_2^k\left(s\right)\end{array}---\left(25\right)$

当风电场或光伏电站内所有直驱风机或光伏发电单元的型号、参数和运行工况完全相同时，式(25)所示的聚合传递函数可以写为：

$\begin{array}{c}{G}_{dj}\left(s\right)=\frac{\mathrm{\Δ}\left|V_{wj}\right|}{{\mathrm{\ΔI}}_{wdj}}=L_j{b}_1^k\left(s\right)\\ G_{qj}\left(s\right)=\frac{\mathrm{\Δ}\left|V_{wj}\right|}{{\mathrm{\ΔI}}_{wqj}}=L_j{b}_2^k\left(s\right)\end{array}---\left(26\right)$

显然当风电场或光伏电站内有部分直驱风机或光伏发电单元的型号、参数和运行工况近似相同时，式(25)同样可以用来降低聚合传递函数的阶数。

当风电场或光伏电站内直驱风机或光伏发电单元的型号、参数和运行工况不同时，由式(25)所得的聚合传递函数的阶数可能会非常高。此时利用成熟的模型降阶方法，可以将高阶的传递函数模型降阶成低阶的传递函数模型，这样可以大大降低电力系统小干扰稳定分析的计算压力和复杂度。

(2)在直驱风机或光伏发电单元的型号、参数或运行工况未知的情况下，可以通过现场频率响应试验来得到风电场和光伏电站的聚合传递函数模型。通过在低压母线稳态电压幅值上叠加某频率的正弦信号作为输入，检测整个风电场或光伏电站的输出电流，可以得到该频率下的输出电流对输入电压的传递函数值G_dj(jω_i)＝|G_dj(jω)|∠G_dj(jω)和G_qj(jω_i)＝|G_qj(jω)|∠G_qj(jω)。将一组连续频率下的传递函数值的幅值和相位作图，形成传递函数幅值图和相位图。然后对幅值图和相位图可以进行不同阶数的曲线拟合，从而可以近似得到风电场和光伏电站的聚合模型。模型降阶的程度可以在曲线拟合过程中给予考虑。

所述步骤2具体包括以下步骤：

建立开环的多机电力系统的线性化模型：

如附图6所示的多机电力系统，有N台同步发电机和M个大规模风电场或光伏电站。第j个风电场或光伏电站内有L_j台直驱风机或光伏发电单元。每台直驱风机或光伏发电单元通过变压器连接到低压母线上，然后再通过升压变压器连接到公共连接点(PCC)处的高压母线上。

含有N台同步发电机的多机电力系统中同步发电机的线性化模型为：

sΔX_g＝A_gΔX_g+B_gΔV_g

(27)

ΔI_g＝C_gΔX_g+D_gΔV_g

其中，X_g为同步机的状态变量，V_g和I_g为N台同步发电机的端电压和输入网络的电流，其每个元素值为：

ΔV_gj＝[ΔV_gxjΔV_gyj]^T

(28)

ΔI_gj＝[ΔI_gxjΔI_gyj]^T

其中x和y表示在x-y坐标系下相应物理量的x分量和y分量。

用节点导纳矩阵表示的网络方程为：

$\left[\begin{array}{c}\mathrm{\Δ}{I}_g\\ {\mathrm{\ΔI}}_w\end{array}\right]=\left[\begin{array}{cc}{Y}_{gg}& Y_{gw}\\ Y_{wg}& Y_{ww}\end{array}\right]\left[\begin{array}{c}\mathrm{\Δ}{V}_g\\ {\mathrm{\ΔV}}_w\end{array}\right]---\left(29\right)$

其中，V_w和I_w为M个风电场或光伏电站的端电压和输入网络的电流，其每个元素值为：

ΔV_wj＝[ΔV_wxjΔV_wyj]^T

(30)

ΔI_wj＝[ΔI_wxjΔI_wyj]^T

电压相量对其进行线性化可得：

Δ|V_wj|＝a_1xjΔV_wxj+a_1yjΔV_wyj

(31)

Δθ_wj＝a_2xjΔV_wxj+a_2yjΔV_wyj

将式(31)写成矩阵形式为：

Δ|V_w|＝A_1xΔV_w

(32)

Δθ_w＝A_2xΔV_w

其中Δ|V_w|和Δθ_w的第j个元素为Δ|V_wj|和Δθ_wj。

对于第j个风电场和光伏电站，令的方向为d-q坐标系中q轴的方向，如附图7所示。由附图7可得：

$\left[\begin{array}{c}{\mathrm{\ΔI}}_{wxj}\\ {\mathrm{\ΔI}}_{wyj}\end{array}\right]=\left[\begin{array}{cc}{\mathrm{sin\θ}}_{wj0}& {\mathrm{cos\θ}}_{wj0}\\ -{\mathrm{cos\θ}}_{wj0}& {\mathrm{sin\θ}}_{wj0}\end{array}\right]\left[\begin{array}{c}{\mathrm{\ΔI}}_{wdj}\\ {\mathrm{\ΔI}}_{wqj}\end{array}\right]+\left[\begin{array}{c}{I}_{wdj0}{\mathrm{cos\θ}}_{cj0}-I_{wqj}{\mathrm{sin\θ}}_{wj0}\\ I_{wdj0}{\mathrm{sin\θ}}_{cj0}+I_{wqj0}{\mathrm{cos\θ}}_{wj0}\end{array}\right]{\mathrm{\Δ\θ}}_{wj}---\left(33\right)$

将式(33)写成矩阵形式为：

${\mathrm{\ΔI}}_{wj}=C_{dqj}\left[\begin{array}{c}{\mathrm{\ΔI}}_{wdj}\\ {\mathrm{\ΔI}}_{wqj}\end{array}\right]+C_{v\mathrm{\θ}j}{\mathrm{\ΔV}}_{wj}---\left(34\right)$

因此有：

ΔI_w＝C_dqΔI_wdq+C_vθΔV_w(35)

其中ΔI_wdq的第j个元素值为

ΔI_wdqj＝[ΔI_wdjΔI_wqj]^T(36)

由式(27)和式(29)可得：

ΔI_g＝C_gΔX_g+D_gΔV_g＝Y_ggΔV_g+Y_gwΔV_w

(37)

ΔV_w＝Y_ww ^-1ΔI_w-Y_ww ^-1Y_wgΔV_g

由式(37)可得：

ΔV_g＝(Y_gg-Y_gwY_ww ^-1Y_wg-D_g)^-1(C_gΔX_g-Y_gwY_ww ^-1ΔI_w)

(38)

＝C_ggΔX_g+C_gwΔI_w

将式(38)代入式(37)的第二式中，可得：

ΔV_w＝(Y_ww ^-1-Y_ww ^-1Y_wgC_gw)ΔI_w-Y_ww ^-1Y_wgC_ggΔX_g(39)

将式(39)代入式(35)中可得：

ΔI_w＝B_dqΔI_wdq+B_wgΔX_g(40)

其中，

B_dq＝[I-C_vθ(Y_ww ^-1-Y_ww ^-1Y_wgC_gw)]^-1C_dq

B_wg＝-[I-C_vθ(Y_ww ^-1-Y_ww ^-1Y_wgC_gw)]^-1C_vθY_ww ^-1Y_wgC_gg

将式(39)(40)代入式(32)中，可得：

Δ|V_w|＝C_0VΔX_g+D_0VΔI_wdq

(41)

Δθ_w＝C_0θΔX_g+D_0θΔI_wdq

其中，

C_0V＝A_1x[(Y_ww ^-1-Y_ww ^-1Y_wgC_gw)B_wg-Y_ww ^-1Y_wgC_gg]

D_0V＝A_1x(Y_ww ^-1-Y_ww ^-1Y_wgC_gw)B_dq

C_0θ＝A_2x[(Y_ww ^-1-Y_ww ^-1Y_wgC_gw)B_wg-Y_ww ^-1Y_wgC_gg]

D_0θ＝A_2x(Y_ww ^-1-Y_ww ^-1Y_wgC_gw)B_dq

由式(38)和式(40)可得：

ΔV_g＝(C_gg+C_gwB_wg)ΔX_g+C_gwB_dqΔI_wdq(42)

将式(42)代入式(27)的第一式中，可得：

sΔX_g＝A₀ΔX_g+B₀ΔI_wdq(43)

其中，

A₀＝A_g+B_g(C_gg+C_gwB_wg)

B₀＝B_gC_gwB_dq

根据式(36)中的定义，式(43)和式(41)可以写为：

$\begin{array}{c}{\mathrm{s\ΔX}}_g=A_0{\mathrm{\ΔX}}_g+\underset{j=1}{\overset{M}{\Σ}}(b_{0dj}{\mathrm{\ΔI}}_{wdj}+b_{0qj}{\mathrm{\ΔI}}_{wqj})\\ \mathrm{\Δ}\left|V_{wj}\right|=C_{0Vj}{\mathrm{\ΔX}}_g+\underset{j=1}{\overset{M}{\Σ}}(d_{0dj}{\mathrm{\ΔI}}_{wdj}+d_{0qj}{\mathrm{\ΔI}}_{wqj})\end{array}---\left(44\right)$

使用上文所述的方法建立各个风电场和光伏电站的聚合的传递函数模型：

ΔI_wdj＝G_dj(s)Δ|V_wj|

${\mathrm{\ΔI}}_{wqj}=G_{qj}\left(s\right)\mathrm{\Δ}\left|V_{wj}\right|+\underset{k=1}{\overset{L_j}{\Σ}}{b}_3^k\left(s\right){\mathrm{\Δu}}_j^k---\left(45\right)$

当风电场或光伏电站的参数或运行工况未知时，采用下面的近似来得到聚合的传递函数模型：

$\begin{array}{c}{G}_{dj}\left(\stackrel{\‾}{{\mathrm{\λ}}_i}\right)\≈G_{dj}\left(j{\stackrel{\‾}{\mathrm{\ω}}}_i\right)=\left|G_{dj}\left(j{\stackrel{\‾}{\mathrm{\ω}}}_i\right)\right|\∠G_{dj}\left(j{\stackrel{\‾}{\mathrm{\ω}}}_i\right)\\ G_{qj}\left(\stackrel{\‾}{{\mathrm{\λ}}_i}\right)\≈G_{qj}\left(j{\stackrel{\‾}{\mathrm{\ω}}}_i\right)=\left|G_{qj}\left(j{\stackrel{\‾}{\mathrm{\ω}}}_i\right)\right|\∠G_{qj}\left(j{\stackrel{\‾}{\mathrm{\ω}}}_i\right)\end{array}---\left(46\right)$

步骤3具体包括：

式(44)和式(45)构成了含M个风电场或光伏电站的多机电力系统的线性化模型，如附图8所示。附图8中的上半部分为不考虑风电场和光伏电站动态的多机电力系统的状态空间模型，它可以看作是开环的多机电力系统模型。由式(45)表示的风电场或光伏电站的传递函数模型可以看作是“控制器”，作为连接开环系统输入与输出的反馈回路。开环系统和反馈回路共同构成了考虑M个风电场或光伏电站动态的多机电力系统的闭环状态空间模型。

大规模风电场和光伏电站的接入对多机电力系统的影响体现在两个方面：一是风电场和光伏电站向原多机电力系统提供一定的注入电流，改变了原多机电力系统的潮流分布，这部分影响体现在式(44)中的A₀矩阵中；另一方面是风电场和光伏电站与多机电力系统的动态交互。如附图8所示，风电场和光伏电站向开环系统输入ΔI_wdj和ΔI_wqj，ΔI_wdj和ΔI_wqj正是风电场和光伏电站对Δ|V_wj|的响应，而Δ|V_wj|又是开环系统对ΔI_wdj和ΔI_wqj的响应。显然，开环系统与反馈回路的动态交互过程与外部变量没有任何关系。因此与外部变量相关的动态过程不会影响到电力系统的小干扰功角稳定，在下面的小信号分析中可以直接忽略掉

步骤4具体包括：

令ΔI_wdj＝ΔI_wqj＝0，此时风电场和光伏电站褪化为恒电流源，此时系统的线性化模型褪化为：

sΔX_g＝A₀ΔX_g(47)

A₀矩阵中包含了大规模风电场和光伏电站的接入对原多机电力系统带来的潮流影响，而忽略了风电场和光伏电站与多机电力系统间的动态交互影响。

对于输入为u_j、输出为y_j的线性系统：

sX＝AX+b_ju_j

(48)

y_j＝c_j ^TX

作为反馈回路的控制器的传递函数模型为：

u_j＝h_j(s)y_j(49)

如果i＝1,2,…,N为式(48)中开环矩阵A的第i个特征值，和为其对应的左右特征向量，由控制器造成的特征值的变化量为：

$\mathrm{\Δ}\stackrel{\‾}{{\mathrm{\λ}}_i}\≈{\stackrel{\‾}{p}}_i^T{b}_i{c}_i^T\stackrel{\‾}{r_i}{h}_j\left(\stackrel{\‾}{{\mathrm{\λ}}_i}\right)={\stackrel{\‾}{R}}_i{h}_j\left(\stackrel{\‾}{{\mathrm{\λ}}_i}\right)---\left(50\right)$

其中是特征值i＝1,2,…,N对应的残差。

如果式(48)表示的开环系统变为：

sX＝AX+b_ju_j

(51)

y_j＝c_j ^TX+d_ju_j

那么由控制器造成的特征值的变化量为：

$\mathrm{\Δ}\stackrel{\‾}{{\mathrm{\λ}}_i}\≈\frac{{\stackrel{\‾}{R}}_i}{1-d_j{\stackrel{\‾}{h}}_j\left(\stackrel{\‾}{{\mathrm{\λ}}_i}\right)}{\stackrel{\‾}{h}}_j\left(\stackrel{\‾}{{\mathrm{\λ}}_i}\right)---\left(52\right)$

闭环系统的特征值为：

${\stackrel{\‾}{\mathrm{\λ}}}_{ci}\≈\stackrel{\‾}{{\mathrm{\λ}}_i}+\mathrm{\Δ}\stackrel{\‾}{{\mathrm{\λ}}_i}=\stackrel{\‾}{{\mathrm{\λ}}_i}+\frac{{\stackrel{\‾}{R}}_i}{1-d_j{\stackrel{\‾}{h}}_j\left(\stackrel{\‾}{\mathrm{\λ}}\right)}{\stackrel{\‾}{h}}_j\left(\stackrel{\‾}{{\mathrm{\λ}}_i}\right)---\left(53\right)$

对比式(44)(45)和式(51)(49)，可以看出，式(44)中的ΔI_wdj和ΔI_wqj相当于式(51)中的控制量u_j，Δ|V_wj|相当于式(51)中的输出量y_j，式(45)中的G_dj(s)和G_qj(s)相当于式(49)中控制器的传递函数h_j(s)。根据式(53)，大规模风电场和光伏电站接入后的多机电力系统的机电振荡模式近似为：

${\stackrel{\‾}{\mathrm{\λ}}}_{ci}\≈\stackrel{\‾}{{\mathrm{\λ}}_i}+\frac{\underset{j=1}{\overset{M}{\Σ}}\[{\stackrel{\‾}{R}}_{dji}{G}_{dj}\left(\stackrel{\‾}{{\mathrm{\λ}}_i}\right)+{\stackrel{\‾}{R}}_{qji}{G}_{qj}\left(\stackrel{\‾}{{\mathrm{\λ}}_i}\right)\]}{1-\underset{j=1}{\overset{M}{\Σ}}\[d_{0dj}{G}_{dj}\left(\stackrel{\‾}{{\mathrm{\λ}}_i}\right)+d_{0qj}{G}_{qj}\left(\stackrel{\‾}{{\mathrm{\λ}}_i}\right)\]}---\left(54\right)$

残差为

$\begin{array}{c}{\stackrel{\‾}{R}}_{dji}={\stackrel{\‾}{p}}_i^T{b}_{0dj}{c_{oVj}}^T\stackrel{\‾}{r_i}\\ {\stackrel{\‾}{R}}_{qji}={\stackrel{\‾}{p}}_i^T{b}_{0qj}{c_{oVj}}^T\stackrel{\‾}{r_i}\end{array}---\left(55\right)$

其中，为式(44)中开环阵A₀的特征值，和为其对应的左右特征向量，为考虑了风电场和光伏电站动态后闭环阵的特征值估计值。

如附图9所示为纽约-新英格兰互联电力系统，下面通过该系统来说明本发明提出的基于大规模风电场和光伏电站聚合模型的电力系统机电振荡模式估算方法。

该电力系统中网络参数和同步发电机的参数见文献(G.Rogers,PowerSystemOscillations.Norwell,MA:Kluwer,2000.)。系统中有两个风电场(连接在节点9和26处)和两个光伏电站(连接在节点16和48处)。每个风电场内有50台直驱风力发电机，每个光伏电站内有50台光伏发电单元。

假设节点9处风电场内全部直驱风机的参数都相同，为：

ω_pm＝0.2,X_d＝0.01,X_q＝0.01,X_w＝0.3,λ₀＝1,J_pm＝6s,C_dc＝50,V_dc＝1,

$K_{pm}\left(s\right)=0.1+\frac{5}{s},K_{dc}\left(s\right)=0.4+\frac{1}{s},K_q\left(s\right)=0.4+\frac{5}{s}$

通过式(1)到式(14)的推导，可以得到该风电场内第k台直驱风机的传递函数模型可以表示为：

${\mathrm{\ΔI}}_{wd}^k=b_{1d}^k\left(s\right)\left|{\mathrm{\ΔV}}_w\right|$

${\mathrm{\ΔI}}_{wq}^k=b_{1q}^k\left(s\right)\left|{\mathrm{\ΔV}}_w\right|+b_{1e}^k\left(s\right){\mathrm{\ΔT}}_{wm}^k$

其中，

$b_{1d}^k\left(s\right)=\frac{-0.0102s-0.1269}{1.4143s+5.1788}$

$b_{1q}^k\left(s\right)=\frac{-9.7060s-24.2651}{50s^2+130.1569s+325.3923}$

则整个风电场的传递函数模型为：

$G_{wd1}\left(s\right)=50b_{1d}^k\left(s\right)$

$G_{wq1}\left(s\right)=50b_{1q}^k\left(s\right)$

节点26处风电场内，有20台直驱风机的参数相同，为：

ω_pm＝0.2,X_d＝0.01,X_q＝0.01,X_w＝0.3,λ₀＝1,J_pm＝6s,C_dc＝50,V_dc＝1,

$K_{pm}\left(s\right)=0.1+\frac{5}{s},K_{dc}\left(s\right)=0.4+\frac{1}{s},K_q\left(s\right)=0.4+\frac{5}{s}$

这20台风机中每台风机的传递函数为：

${\mathrm{\ΔI}}_{wd}^k=b_{3d1}^k\left(s\right)\left|{\mathrm{\ΔV}}_w\right|$

${\mathrm{\ΔI}}_{wq}^k=b_{3q1}^k\left(s\right)\left|{\mathrm{\ΔV}}_w\right|+b_{3e1}^k\left(s\right){\mathrm{\ΔT}}_{wm}^k$

其中，

$b_{3d1}^k\left(s\right)=\frac{-0.0063s-0.0787}{1.4177s+5.2211}$

$b_{3q1}^k\left(s\right)=\frac{-6.0171s-15.0427}{50s^2+131.2208s+328.0521}$

另30台直驱风机的参数相同，为：

ω_pm＝0.2,X_d＝0.01,X_q＝0.01,X_w＝0.3,λ₀＝1,J_pm＝6s,C_dc＝50,V_dc＝1,

$K_{pm}\left(s\right)=0.1+\frac{5}{s},K_{dc}\left(s\right)=0.5+\frac{1.1}{s},K_q\left(s\right)=0.5+\frac{5.1}{s}$

这30台风机中每台风机的传递函数为：

${\mathrm{\ΔI}}_{wd}^k=b_{3d2}^k\left(s\right)\left|{\mathrm{\ΔV}}_w\right|$

${\mathrm{\ΔI}}_{wq}^k=b_{3q2}^k\left(s\right)\left|{\mathrm{\ΔV}}_w\right|+b_{3e2}^k\left(s\right){\mathrm{\ΔT}}_{wm}^k$

其中，

$b_{3d2}^k\left(s\right)=\frac{-0.0063s-0.0787}{1.4177s+5.2211}$

$b_{3q2}^k\left(s\right)=\frac{-6.0171s-15.0427}{50s^2+131.2208s+328.0521}$

那么该风电场的聚合传递函数为：

$G_{wd3}\left(s\right)=20b_{3d1}^k\left(s\right)+30b_{3d2}^k\left(s\right)$

$G_{wq3}\left(s\right)=20b_{3q1}^k\left(s\right)+30b_{3q2}^k\left(s\right)$

节点16处光伏电站内，有20台光伏发电单元的参数相同，为：

$T=293.15,I_0=1.2987\×{10}^{-4},n=1.8405,I_p=\frac{N_p{I}_{sc}{I}_r}{100}=1.82,I_r=100\%,$

$N_s=9,N_p=8,L_{dc}=0.1,C_{dc}=50,V_{dc}=1,K_{pv1}\left(s\right)=0.1+\frac{1}{s},K_{pv2}\left(s\right)=0.1+\frac{10}{s},$

$K_{dc}\left(s\right)=0.5+\frac{1.1}{s},K_q\left(s\right)=0.5+\frac{5.1}{s}$

通过式(15)到式(22)的推导，可以得到这20台光伏发电单元中每台光伏发电单元的传递函数为：

${\mathrm{\ΔI}}_{wd}^k=b_{2d1}^k\left(s\right)\left|{\mathrm{\ΔV}}_w\right|$

${\mathrm{\ΔI}}_{wq}^k=b_{2q1}^k\left(s\right)\left|{\mathrm{\ΔV}}_w\right|+b_{2e1}^k\left(s\right){\mathrm{\ΔI}}_r^k$

其中，

$b_{2d1}^k\left(s\right)=\frac{-0.0081s-0.0826}{1.5076s+5.1771}$

$b_{2q1}^k\left(s\right)=\frac{-7.7370s-17.0214}{50s^2+159.4551s+350.8012}$

另30台光伏发电单元的参数相同，为：

$T=293.15,I_0=1.2987\×{10}^{-4},n=1.8405,I_p=\frac{N_p{I}_{sc}{I}_r}{100}=1.82,I_r=100\%,$

$N_s=9,N_p=8,L_{dc}=0.1,C_{dc}=50,V_{dc}=1,K_{pv1}\left(s\right)=0.1+\frac{1}{s},K_{pv2}\left(s\right)=0.1+\frac{10}{s},$

$K_{dc}\left(s\right)=0.4+\frac{1}{s},K_q\left(s\right)=0.4+\frac{5}{s}$

这30台光伏发电单元中每台光伏发电单元的传递函数为：

${\mathrm{\ΔI}}_{wd}^k=b_{2d2}^k\left(s\right)\left|{\mathrm{\ΔV}}_w\right|$

${\mathrm{\ΔI}}_{wq}^k=b_{2q2}^k\left(s\right)\left|{\mathrm{\ΔV}}_w\right|+b_{2e2}^k\left(s\right){\mathrm{\ΔI}}_r^k$

其中，

$b_{2d2}^k\left(s\right)=\frac{-0.0052s-0.0648}{1.4060s+5.0756}$

$b_{2q2}^k\left(s\right)=\frac{-4.9517s-12.3792}{50s^2+127.5641s+318.9102}$

那么该光伏电站的聚合传递函数为：

$G_{wd2}\left(s\right)=20b_{2d1}^k\left(s\right)+30b_{2d2}^k\left(s\right)$

$G_{wq2}\left(s\right)=20b_{2q1}^k\left(s\right)+30b_{2q2}^k\left(s\right)$

节点48处光伏电站内光伏发电单元的参数未知，需要通过现场频率响应试验来得出其聚合模型。要得到其聚合传递函数模型，需要令输入信号为光伏电站内低压母线的电压幅值Δ|V_wj|，输出信号为整个光伏电站输出到网络中的电流ΔI_wd和ΔI_wq。通过在输入的低压母线稳态电压幅值上叠加某频率的正弦信号，可以检测到带有一定噪声的风电场或光伏电站的输出电流，即得到该频率下的输出电流ΔI_wd和ΔI_wq对输入电压Δ|V_wj|的传递函数值。将各个频率下传递函数的幅值和相位作图，进行适当降噪处理后，即得到如附图10所示的幅值图和相位图。最后对附图10所示的曲线进行不同阶数的曲线拟合，从而可以近似得到节点48处连接的光伏电站的聚合传递函数模型：

ΔI_wd4＝G_wd4(s)|ΔV_w|

ΔI_wq4＝G_wq4(s)|ΔV_w|+G_we4(s)ΔI_r

其中，

$G_{wd4}\left(s\right)=\frac{-0.0052s-0.0648}{1.4060s+5.0756}$

$G_{wq4}\left(s\right)=\frac{-4.9517s-12.3792}{50s^2+127.5641s+318.9102}$

下面验证本发明提出的大规模风电场和光伏电站聚合模型的正确性和计算机电振荡模式的方法的正确性。

表1不同方法求得的特征值

(1)当如附图9所示的互联电力系统中所有直驱风机和光伏发电单元全部使用各自的全模型，即直驱风机使用5阶模型，光伏发电单元使用6阶模型时，此时通过求闭环状态矩阵，可以求得特征值如表1的第一列所示。用直驱风机和光伏发电单元的全模型时，计算最为准确，但是状态矩阵阶数和计算复杂度也是最高的，此时由于风电场和光伏电站的接入而增加的状态矩阵阶数为5×100+6×100＝1100。

(2)当如附图9所示的互联电力系统中所有直驱风机和光伏发电单元全部使用本发明中提出的聚合模型时，即节点9处的风电场为3阶模型、节点26处的风电场为6阶模型、节点16处的光伏电站为6阶模型、节点48处的光伏电站为3阶模型，此时由于风电场和光伏电站的接入而增加的状态矩阵阶数为3+6+6+3＝18。此时通过求闭环状态矩阵而得到的系统特征值如表1中的第二列所示。

(3)通过式(27)到式(44)的推导，可以得到开环特征矩阵，开环特征值即可以求得；再利用式(54)和式(55)，可以得到闭环特征矩阵的估算值，如表1的第三列所示。这种方法并没有由于风电场和光伏电站的接入而增加状态矩阵的阶数。

通过对比表1中三列的数值，可以看出，本发明中提出的大规模风电场和光伏电站的聚合模型是正确的，而且本发明中提出的机电振荡模式估算方法也是比较准确的，可以在工程实际中应用。

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

Claims (6)

1.一种基于大规模风电场和光伏电站聚合模型的电力系统机电振荡模式估算方法，其特征在于，包括：

步骤1、建立大规模风电场和光伏电站的聚合传递函数模型；

步骤2、将电力系统中风电场和光伏电站的输出电流作为控制量，建立开环的多机电力系统的线性化模型；

步骤3、把风电场或光伏电站的传递函数模型看作是控制器，作为连接开环系统输入与输出的反馈回路；开环系统和反馈回路共同构成了多机电力系统的闭环状态空间模型；

步骤4、将风电场和光伏电站看作恒电流源，即令风电场或光伏电站输出电流为零，忽略掉风电场和光伏电站与多机电力系统间的动态交互影响，计算开环系统的特征值；通过系统的开环特征值加上聚合传递函数值与残差的乘积来估算电力系统闭环的机电振荡模式。

2.根据权利要求1所述方法，其特征在于，所述步骤1若在各台直驱风力发电机或光伏发电单元的型号、参数和运行工况都已知的情况下，则建立的各台直驱风机和光伏发电单元的传递函数模型为：

$\begin{array}{c}{\mathrm{\ΔI}}_{wdj}^k=b_1^k\left(s\right)\mathrm{\Δ}\left|V_{wj}\right|\\ {\mathrm{\ΔI}}_{wqj}^k=b_2^k\left(s\right)\mathrm{\Δ}\left|V_{wj}\right|+b_3^k\left(s\right){\mathrm{\Δu}}_j^k\end{array}---\left(23\right)$

其中，和分别为第j个风电场或光伏电站内第k台直驱风机或光伏发电单元输出电流的d轴和q轴分量，|V_wj|为第j个风电场或光伏电站内低压母线的电压幅值，为第j个风电场或光伏电站内第k个直驱风机或光伏发电单元的外部变量；对于直驱风机，外部变量为机械转矩，对于光伏发电单元，外部变量为光照强度；和分别为所对应变量间的传递函数。

3.根据权利要求1所述方法，其特征在于，所述步骤1若在一个风电场或光伏电站内所有直驱风机或光伏发电单元的型号、参数和运行工况都相同的情况下，那么该风电场或光伏电站的聚合传递函数模型由单台风机或光伏发电单元的传递函数倍乘直接得到：

$\begin{array}{c}{G}_{dj}\left(s\right)=\frac{{\mathrm{\ΔI}}_{wdj}}{\mathrm{\Δ}\left|V_{wj}\right|}=L_j{b}_1^k\left(s\right)\\ G_{qj}\left(s\right)=\frac{{\mathrm{\ΔI}}_{wqj}}{\mathrm{\Δ}\left|V_{wj}\right|}=L_j{b}_2^k\left(s\right)\end{array}---\left(26\right)$

若一个风电场或光伏电站内直驱风机或光伏发电单元的型号、参数或运行工况不同，那么整个风电场或光伏电站的聚合传递函数模型由所有风机或光伏发电单元的传递函数相加得到：

$\begin{array}{c}{G}_{dj}\left(s\right)=\frac{{\mathrm{\ΔI}}_{wdj}}{\mathrm{\Δ}\left|V_{wj}\right|}=\underset{k=1}{\overset{L}{\Σ}}{b}_1^k\left(s\right)\\ G_{qj}\left(s\right)=\frac{{\mathrm{\ΔI}}_{wqj}}{\mathrm{\Δ}\left|V_{wj}\right|}=\underset{k=1}{\overset{L}{\Σ}}{b}_2^k\left(s\right)\end{array}---\left(25\right)$

其中，I_wdj和I_wqj分别为第j个风电场或光伏电站内总输出电流的d轴和q轴分量，L_j为第j个风电场或光伏电站内所包含直驱风机或光伏发电单元的数量，G_dj(s)和G_qj(s)分别为第j个风电场或光伏电站内输出量ΔI_wdj和ΔI_wqj对输入量Δ|V_wj|的聚合传递函数；

对于式(25)，传递函数G_dj(s)和G_qj(s)的阶数会很高，通过目前已有的模型降阶方法，将高阶的传递函数模型降阶成低阶的传递函数模型。

4.根据权利要求1所述方法，其特征在于，所述步骤1中在直驱风机或光伏发电单元的型号、参数或运行工况未知的情况下，通过现场频率响应试验来得到风电场和光伏电站的聚合传递函数模型；通过在低压母线稳态电压幅值上叠加某频率的正弦信号作为输入，检测风电场或光伏电站的输出电流，得到该频率下的输出电流对输入电压的传递函数值；对各个频率下传递函数值的幅值图和相位图，进行不同阶数的拟合，从而近似得到风电场和光伏电站的聚合模型。

5.根据权利要求1所述方法，其特征在于，所述步骤2中建立开环的多机电力系统的线性化模型为：

$\begin{array}{c}{\mathrm{s\ΔX}}_g=A_0{\mathrm{\ΔX}}_g+\underset{j=1}{\overset{M}{\Σ}}(b_{0dj}{\mathrm{\ΔI}}_{wdj}+b_{0qj}{\mathrm{\ΔI}}_{wqj})\\ \mathrm{\Δ}\left|V_{wj}\right|=c_{0vj}^T{\mathrm{\ΔX}}_g+\underset{j=1}{\overset{M}{\Σ}}(d_{0dj}{\mathrm{\ΔI}}_{wdj}+d_{0qj}{\mathrm{\ΔI}}_{wqj})\end{array}---\left(44\right)$

其中，X_g为系统状态变量，A₀为开环系统的状态矩阵，b_0dj和b_0qj为控制矩阵，为输出矩阵，d_0dj和d_0qj为反馈系数，M为多机电力系统所包含风电场和光伏电站的数量，I_wdj和I_wqj分别为第j个风电场或光伏电站总输出电流的d轴和q轴分量，|V_wj|为第j个风电场或光伏电站内低压母线的电压幅值。

6.根据权利要求1所述方法，其特征在于，所述步骤4中电力系统机电振荡模式为：

${\stackrel{\‾}{\mathrm{\λ}}}_{\mathrm{ci}}\≈{\stackrel{\‾}{\mathrm{\λ}}}_i+\frac{\underset{j=1}{\overset{M}{\mathrm{\Σ}}}[{\stackrel{\‾}{R}}_{\mathrm{dji}}{G}_{\mathrm{dj}}\left({\stackrel{\‾}{\mathrm{\λ}}}_i\right)+{\stackrel{\‾}{R}}_{\mathrm{qji}}{G}_{\mathrm{qj}}\left({\stackrel{\‾}{\mathrm{\λ}}}_i\right)]}{1-\underset{j=1}{\overset{M}{\mathrm{\Σ}}}[d_{0\mathrm{dj}}{G}_{\mathrm{dj}}\left({\stackrel{\‾}{\mathrm{\λ}}}_i\right)+d_{0\mathrm{qj}}{G}_{\mathrm{qj}}\left({\stackrel{\‾}{\mathrm{\λ}}}_i\right)]}---\left(54\right)$

残差为：

$\begin{array}{c}{\stackrel{\‾}{R}}_{\mathrm{dji}}={\stackrel{\‾}{p}}_i^T{b}_{0\mathrm{dj}}{c_{\mathrm{oVj}}}^T{\stackrel{\‾}{r}}_i\\ {\stackrel{\‾}{R}}_{\mathrm{qji}}={\stackrel{\‾}{p}}_i^T{b}_{0\mathrm{qj}}{c_{\mathrm{oVj}}}^T{\stackrel{\‾}{r}}_i\end{array}---\left(55\right)$

其中，为开环系统的状态矩阵的特征值，和分别为对应的左右特征向量，为考虑了风电场和光伏电站动态后闭环阵的特征值估计值，G_dj(s)和G_qj(s)分别为第j个风电场或光伏电站的输出量ΔI_wdj和ΔI_wqj对输入量Δ|V_wj|的聚合传递函数，b_0dj和b_0qj为控制矩阵，为输出矩阵，d_0dj和d_0qj为反馈系数，M为多机电力系统所包含风电场和光伏电站的数量。