风力发电是利用新能源的主要方式之一，然而，目前电网还是存在一定的弃风现象，究其原因是对于风电技术的研究还不够成熟。实际已经发生一些事故，如2015年新疆哈密风电场发生次同步振荡，导致相邻火电机组保护动作，表明类似的问题亟待进一步研究解决。

研究次同步振荡的方法，主要有特征值分析法、复转矩系数法、阻抗法、时域仿真法等^[1,2]。本文重点研究特征值法，该方法是通过建立系统的小扰动线性化模型，求解特征根、特征向量来分析系统动态响应的方法^[3]。将一个n维动态系统在其运行点线性化后，可得到线性化模型ΔX˙＝AΔX，A为状态矩阵。通过对A求特征值可获得n个特征根(λ＝σ＋jω)，实部的正负可判断系统的稳定性。负实部表示衰减模态，正实部表示发散模态，每一对共轭负特征根对应一种振荡模态。

文献[4,5]应用特征值分析法分析电力系统的稳定性问题；文献[6]针对特征值法存在的“维数灾”问题进行了降维研究，然而特征值分析的详细过程叙述甚少，尤其是对初值计算的说明。本文以双馈感应电机(doubly fed induction generator, DFIG)单机无穷大系统为研究对象，将其在Simulink中搭建小扰动模型，用特征值分析法进行相关研究。

本文所研究系统为工频60 Hz双馈风机单机无穷大系统^[7]，如图1所示。

图1Fig.1

图1中IG为100 MW(50×2 MW)的等值感应电机，定子直接与升压变低压侧690 V母线相连，进而升压经过161 kV串补线路连接到无穷大电网。与发电机端口相连的部分称为网侧换流器(grid side converter, GSC)，与发电机转子相连的称为转子侧换流器(rotor side converter, RSC)。如图中虚线框所示，发电机轴系采用两质量块模型。系统的主要参数参考文献[8]。

当风以速度V_w穿过横扫面积为S的桨叶时，风力机吸收的功率P_t可表示为

(1)Pt=ρSCp(β,λ)Vw3/2

式中：ρ为空气密度；C_p为风能利用系数，C_p是关于桨距角β和叶尖速比λ的函数。

叶轮转矩T_t由吸收的功率和风机的转速ω_t计算得到：

(2)Tt=Pt/ωt

风速、转速和转矩的关系^[8]如表1所示，转矩由对应风速下最大功率所得。空气动力学模型采用Simulink中的Lookup Tables模块描述。

现有的风力发电机轴系模型有单质量块、两质量块和三质量块模型^[9,10]。两质量块模型广泛用于包含风机的电力系统稳定性分析中，因此本文采用两质量块模型，其3阶动态微分方程如下：

(3)X˙tg=AtgXtg+BtgUtg

其中：

(4)XtgT=[ωtωrTtg]

(5)UtgT=[TtTe0]

(6)Atg=[−Dt−Dtg2HtDtg2Ht−12HtDtg2Hg−Dt−Dtg2Hg12HtKtgωb−Ktgωb0]

(7)Btg=diag[12Ht−12Ht1]

轴系方程中，ω_t、ω_r分别为叶轮、发电机转子转速，T_tg为两质量块内部转矩，T_t、T_e分别为叶轮转矩和发电机电磁转矩，H_t、H_g为叶轮和发电机的惯性常数(s)，D_t、D_g为对应的阻尼系数，K_tg、D_tg为两质量块间的弹性系数和阻尼系数。

dq同步旋转坐标系下的6阶异步感应电机方程^[11]如下：

(8)X˙G=AGXG+BGUG

其中：

(9)XGT=[iqsidsi0siqridri0r]

(10)UGT=[uqsudsu0suqrudru0r]

式中：i_qs、i_ds、u_qs和u_ds分别为感应电机定子的dq轴电流、电压；i_qr、i_dr、u_qr和u_dr分别为转子的dq轴电流、电压；其他为零序电流、电压分量。

(11)A1=ωb[Xss00XM000Xss00XM000Xls000XM00Xrr000XM00Xrr000000Xlr]−1

(12){ω1=ωeωbω2=ωe−ωrωb

(13)A2=[Rsω1Xss00ω1XM0−ω1XssRs0−ω1XM0000Rs0000ω2XM0Rrω2Xrr0−ω2XM00−ω2XrrRr000000Rr]

(14){AG=−A1A2BG=A1

式中：X_ls、X_lr和X_M分别为定、转子漏抗和互感电抗；R_s、X_ss分别为定子电阻、电抗；R_r、X_rr分别为转子电阻、电抗；ω_b为基准角速度；ω_e为旋转坐标角速度。

转子和定子侧采用文献[8]的控制策略及参数，其控制回路如图2、图3所示。直流部分电容是连接定转子换流器的关键部分，本文采用下式所示的1阶微分方程描述其动态过程：

(15)CUDCdUDCdt=−(Pr+Pg)

式中：C为换流器电容值；U_DC为电容电压；P_r、P_g为转子侧和网侧有功功率。

图2Fig.2

图3Fig.3

对线路进行串联补偿，能增大电能传输距离与容量，但随着串补度K的增加，也带来了系统谐振的风险，下文对次同步的分析也与串补密切相关。包含4个状态变量的串补网络动态微分方程如下：

(16)X˙T=ATXT+BTUT

其中：

(17)XTT=[iqiducqucd]

(18)UTT=[uqs−EqXLuds−EdXL00]

(19)AT=ωb[−RLXL−ωe−1XL0ωe−RLXL0−1XLXC00−ωe0XCωe0]

(20)BT=ωbI

式中：i_q、i_d分别为线路dq轴电流；u_cq、u_cd分别为串补电容两端电压；E_q、E_d为无穷大母线电压。

本文在Simulink平台下搭建图形化动态模型，主要用到元件库中的Matlab Function block(下文用fcn代替)。在原始动态微分方程不变的情况下，将其以一种可视化的形式呈现。

fcn模块如图4所示，利用它可将方程写到Simulink中来使用，其流程如图5所示。

图4Fig.4

图5Fig.5

在画布中拉入一个fcn模块，根据式(8)—(14)输入后，命名为“A_DFIG fcn”；同理，完成输入矩阵与输入变量的自定义，命名为“B_ DFIG fcn”。根据式(8)，加入求和、积分及变量分离模块后，感应电机的Simulink模型如图6所示。

图6Fig.6

对于所研究系统的轴系模块、串补网络、RSC和GSC等的搭建，参照DFIG模型搭建的过程，结合式(3)、(15)、(16)即可完成。

在fcn模块中，动态微分方程以矩阵的形式输入，这对于方程的检查与修改都是极其便捷的。再者，模型中信号之间的传递采用图7所示的From与Goto模块^[12]，避免了较多信号连接线之间的交错。这都将使轴系、电机和传输线路等模块内部以及模块间的内容和结构更加清晰直观。

图7Fig.7

在风速为7 m/s，串联补偿度为30%的条件下，对系统进行特征值分析。这部分主要用到3个Matlab命令语句：fsolve、linmod与eig，分别对应状态变量的初值求解、模型线性化与求取状态矩阵特征值的功能。

fsolve是Matlab中求解非线性方程组的语句，其基本原理是最小二乘法。这里用来求解状态变量的初值，待求解的微分方程组由式(8)、(16)以及如下各有功、无功之间的关系式构建：

(21){PS0+Pr0+Ptref=0PL0+PS0−Pg0=0QL0+QS0−Qg0=0

式中：P_S0、Q_S0为定子初始有功、无功功率；P_r0、Q_r0、P_g0和Q_g0分别为转子侧、网侧换流器初始有功、无功功率；P_L0、Q_L0为传输线路初始有功、无功功率。各有功、无功功率流向如图1所示。

在Matlab编写．m文件，将上述求解初值的微分方程写入，利用fsolve命令即可得到系统变量的初值。需要注意的是，除表2中所列出的初值外，其他初值，如叶轮转速查表1可得。

linmod是提取系统稳态运行点附近、连续时间线性状态空间模型的语句。在将2.1节所建模型线性化的同时，可得到系统A、B、C和D矩阵。

针对状态矩阵A，利用Matlab中的eig语句便可得系统的特征值。表3为当前工况下系统各模态及其对应的特征值。

求解参与因子离不开特征值所对应的左右特征向量。第i个左特征向量表征了状态变量对第i个模态的贡献大小，体现的是“可观性”；第i个右特征向量表征的是状态变量中第i个模态的活跃程度，体现的是“可控性”。参与因子可通过下式得到：

(22)Pki=|Yki||Zik|∑​k=1n|Yki||Zik|

式中：P_ki为第k个状态变量与第i个模态的参与因子，反映了状态变量对模态的强可观强可控，是一个综合性指标；Z_i和Y_i分别为第i个模态的左右特征向量。

值得注意的是，通常所说的特征向量一般指右特征向量。而针对左特征向量的求解有两种方法：一种是先求状态矩阵A的转置的右特征向量，再将其转置即为左特征向量；另一种方法是直接将右特征向量求逆再转置即为左特征向量。由于代数精度问题，Matlab中二者求得的左特征向量会有细微差别，但是不影响问题的分析。本文中采用第1种方法求取左特征向量。

表4为部分系统状态变量对应的参与因子，为便于后续模态分析，其中较大值已加粗表示。

当风速为7 m/s时，由表1可知叶轮转速为0.75 pu，换成有名值即为45 Hz(60×0.75＝45 Hz)。在该风速下线路串补为30%时，模态λ_7，8对应的特征值为－6.753±j93.821，对应频率为93.821 rad/s，即14.93 Hz，刚好与45 Hz互补(60－14.93＝45.07 Hz)，其他风速下亦如此。再而，由表4可知此模态与定转子dq轴电流呈强相关关系。因此，模态λ_7，8与机械动态(风速、叶轮转速)和电气动态(定、转子电流)相关，故将其识别为机电模态。另外，通过在表5中纵向观察比较，可发现此模态受串补度影响较小。

由表4可知λ_12，13与叶轮转速ω_t以及两质量块间转矩T_tg强相关，因此将其识别为轴系模态，轴系振荡频率很低，约0.95 Hz。

系统的自然振荡频率为

(23)fn=fsXC/XL

式中：f_s为系统工频；X_C为系统容抗；X_L为系统感抗。次同步振荡频率为f_s－f_n，超同步振荡频率为f_s＋f_n。

因此，串补为30%时的自然频率为32.86 Hz，次同步频率为27.14 Hz(60－32.86＝27.14 Hz)，超同步频率为92.86 Hz。与之对应的是模态λ_5，6，频率为28.25 Hz(177.5/2π)，模态λ_3，4频率为91.53 Hz。故将λ_5，6识别为次同步模态，λ_3，4识别为超同步模态。

表6列出了不同风速与串补度下的次/超同步模态，可看出：超同步模态实部均为负；而次同步模态在7 m/s，串补度由45%变为50%时，实部变为正数，振荡发散，变为不稳定。也因此，少有人研究超同步振荡，而把重点放在次同步振荡上。

图8为不同风速、串补度下次同步模态的频率和阻尼比关系图，K由25%开始，间隔为5%依次递增至90%。可观察到K增加或风速减小，都会导致阻尼比减小，使系统存在次同步振荡风险。

图8Fig.8

图9反映改变转子侧换流器控制参数对次同步模态的影响，分别使K_p和K_i参数中的一个不变，增加另一个的值。可发现，当控制参数增加到一定数值后，阻尼比呈现负值，系统次同步模态不稳定。

图9Fig.9

事实上，文献[13,14,15]都已证实风速、串补度及控制参数等对双馈风机次同步的影响，本文对此的分析旨在验证一般性结论，证明了所搭建模型的合理性。

本文列写了双馈风机单机无穷大系统的动态微分方程组，包括轴系、发电机、串补网络和控制环节等在内的22个状态变量。基于Simulink平台搭建了用于双馈风机特征值分析的fcn模型，详细介绍了包括模型搭建，及运行初值、特征值与参与因子计算在内的过程。通过模态分析，识别出了机电模态、轴系模态和次/超同步模态，并通过改变风速、串补度和控制参数，验证了关于次同步振荡的一般性结论。所做工作对于今后的研究有基础性意义。