秦浩宇，蘇勛文，楊福寶，葉圣瞳

(1.黑龍江科技大學(xué)電氣與控制工程學(xué)院，黑龍江 哈爾濱 150022；2.國網(wǎng)浙江省電力有限公司溫州供電公司，浙江 溫州 325000)

0 引言

風(fēng)力發(fā)電是利用新能源的主要方式之一，然而，目前電網(wǎng)還是存在一定的棄風(fēng)現(xiàn)象，究其原因是對于風(fēng)電技術(shù)的研究還不夠成熟。實(shí)際已經(jīng)發(fā)生一些事故，如2015年新疆哈密風(fēng)電場發(fā)生次同步振蕩，導(dǎo)致相鄰火電機(jī)組保護(hù)動作，表明類似的問題亟待進(jìn)一步研究解決。

文獻(xiàn)[4-5]應(yīng)用特征值分析法分析電力系統(tǒng)的穩(wěn)定性問題；文獻(xiàn)[6]針對特征值法存在的“維數(shù)災(zāi)”問題進(jìn)行了降維研究，然而特征值分析的詳細(xì)過程敘述甚少，尤其是對初值計算的說明。本文以雙饋感應(yīng)電機(jī)(doubly fed induction generator, DFIG)單機(jī)無窮大系統(tǒng)為研究對象，將其在Simulink中搭建小擾動模型，用特征值分析法進(jìn)行相關(guān)研究。

1 DFIG模型的建立

本文所研究系統(tǒng)為工頻60 Hz雙饋風(fēng)機(jī)單機(jī)無窮大系統(tǒng)[7]，如圖1所示。

圖1 雙饋風(fēng)機(jī)單機(jī)無窮大系統(tǒng)結(jié)構(gòu)圖Fig.1 DFIG single-machine infinite system

圖1中IG為100 MW(50×2 MW)的等值感應(yīng)電機(jī)，定子直接與升壓變低壓側(cè)690 V母線相連，進(jìn)而升壓經(jīng)過161 kV串補(bǔ)線路連接到無窮大電網(wǎng)。與發(fā)電機(jī)端口相連的部分稱為網(wǎng)側(cè)換流器(grid side converter, GSC)，與發(fā)電機(jī)轉(zhuǎn)子相連的稱為轉(zhuǎn)子側(cè)換流器(rotor side converter, RSC)。如圖中虛線框所示，發(fā)電機(jī)軸系采用兩質(zhì)量塊模型。系統(tǒng)的主要參數(shù)參考文獻(xiàn)[8]。

1.1 風(fēng)機(jī)空氣動力學(xué)模型

當(dāng)風(fēng)以速度Vw穿過橫掃面積為S的槳葉時，風(fēng)力機(jī)吸收的功率Pt可表示為

(1)

式中:ρ為空氣密度；Cp為風(fēng)能利用系數(shù)，Cp是關(guān)于槳距角β和葉尖速比λ的函數(shù)。

葉輪轉(zhuǎn)矩Tt由吸收的功率和風(fēng)機(jī)的轉(zhuǎn)速ωt計算得到：

Tt=Pt/ωt

(2)

風(fēng)速、轉(zhuǎn)速和轉(zhuǎn)矩的關(guān)系[8]如表1所示，轉(zhuǎn)矩由對應(yīng)風(fēng)速下最大功率所得?？諝鈩恿W(xué)模型采用Simulink中的Lookup Tables模塊描述。

表1 不同風(fēng)速下葉輪轉(zhuǎn)速與轉(zhuǎn)矩Table 1 Rotor shaft speed and torque at different wind speeds

1.2 軸系模型

現(xiàn)有的風(fēng)力發(fā)電機(jī)軸系模型有單質(zhì)量塊、兩質(zhì)量塊和三質(zhì)量塊模型[9-10]。兩質(zhì)量塊模型廣泛用于包含風(fēng)機(jī)的電力系統(tǒng)穩(wěn)定性分析中，因此本文采用兩質(zhì)量塊模型，其3階動態(tài)微分方程如下：

(3)

其中：

軸系方程中，ωt、ωr分別為葉輪、發(fā)電機(jī)轉(zhuǎn)子轉(zhuǎn)速，Ttg為兩質(zhì)量塊內(nèi)部轉(zhuǎn)矩，Tt、Te分別為葉輪轉(zhuǎn)矩和發(fā)電機(jī)電磁轉(zhuǎn)矩，Ht、Hg為葉輪和發(fā)電機(jī)的慣性常數(shù)(s)，Dt、Dg為對應(yīng)的阻尼系數(shù)，Ktg、Dtg為兩質(zhì)量塊間的彈性系數(shù)和阻尼系數(shù)。

1.3 感應(yīng)電機(jī)模型

dq同步旋轉(zhuǎn)坐標(biāo)系下的6階異步感應(yīng)電機(jī)方程[11]如下：

(8)

其中：

式中：iqs、ids、uqs和uds分別為感應(yīng)電機(jī)定子的dq軸電流、電壓；iqr、idr、uqr和udr分別為轉(zhuǎn)子的dq軸電流、電壓；其他為零序電流、電壓分量。

(11)

式中：Xls、Xlr和XM分別為定、轉(zhuǎn)子漏抗和互感電抗；Rs、Xss分別為定子電阻、電抗；Rr、Xrr分別為轉(zhuǎn)子電阻、電抗；ωb為基準(zhǔn)角速度；ωe為旋轉(zhuǎn)坐標(biāo)角速度。

1.4 變流器部分模型

轉(zhuǎn)子和定子側(cè)采用文獻(xiàn)[8]的控制策略及參數(shù)，其控制回路如圖2、3所示。直流部分電容是連接定轉(zhuǎn)子換流器的關(guān)鍵部分，本文采用下式所示的1階微分方程描述其動態(tài)過程：

圖2 轉(zhuǎn)子側(cè)控制策略Fig.2 RSC control strategy

圖3 定子側(cè)控制策略Fig.3 GSC control strategy

(15)

式中：C為換流器電容值；UDC為電容電壓；Pr、Pg為轉(zhuǎn)子側(cè)和網(wǎng)側(cè)有功功率。

1.5 串補(bǔ)網(wǎng)絡(luò)模型

對線路進(jìn)行串聯(lián)補(bǔ)償，能增大電能傳輸距離與容量，但隨著串補(bǔ)度K的增加，也帶來了系統(tǒng)諧振的風(fēng)險，下文對次同步的分析也與串補(bǔ)密切相關(guān)。包含4個狀態(tài)變量的串補(bǔ)網(wǎng)絡(luò)動態(tài)微分方程如下：

(16)

其中：

式中：iq、id分別為線路dq軸電流；ucq、ucd分別為串補(bǔ)電容兩端電壓；Eq、Ed為無窮大母線電壓。

2 Simulink下特征值分析

2.1 模型搭建

本文在Simulink平臺下搭建圖形化動態(tài)模型，主要用到元件庫中的Matlab Function block(下文用fcn代替)。在原始動態(tài)微分方程不變的情況下，將其以一種可視化的形式呈現(xiàn)。

fcn模塊如圖4所示，利用它可將方程寫到Simulink中來使用，其流程如圖5所示。

圖4 Matlab中的fcn模塊Fig.4 Fcn block in Matlab

圖5 fcn模塊的使用流程圖Fig.5 Flow chart of fcn block

在畫布中拉入一個fcn模塊，根據(jù)式(8)—(14)輸入后，命名為“A_DFIG fcn”；同理，完成輸入矩陣與輸入變量的自定義，命名為“B_ DFIG fcn”。根據(jù)式(8)，加入求和、積分及變量分離模塊后，感應(yīng)電機(jī)的Simulink模型如圖6所示。

圖6 Simulink感應(yīng)電機(jī)部分結(jié)構(gòu)圖Fig.6 Simulink induction generator part structure

對于所研究系統(tǒng)的軸系模塊、串補(bǔ)網(wǎng)絡(luò)、RSC和GSC等的搭建，參照DFIG模型搭建的過程，結(jié)合式(3)、(15)、(16)即可完成。

在fcn模塊中，動態(tài)微分方程以矩陣的形式輸入，這對于方程的檢查與修改都是極其便捷的。再者，模型中信號之間的傳遞采用圖7所示的From與Goto模塊[12]，避免了較多信號連接線之間的交錯。這都將使軸系、電機(jī)和傳輸線路等模塊內(nèi)部以及模塊間的內(nèi)容和結(jié)構(gòu)更加清晰直觀。

圖7 Matlab中From與Goto模塊Fig.7 From and Goto block in Matlab

2.2 特征值分析

在風(fēng)速為7 m/s，串聯(lián)補(bǔ)償度為30%的條件下，對系統(tǒng)進(jìn)行特征值分析。這部分主要用到3個Matlab命令語句：fsolve、linmod與eig，分別對應(yīng)狀態(tài)變量的初值求解、模型線性化與求取狀態(tài)矩陣特征值的功能。

fsolve是Matlab中求解非線性方程組的語句，其基本原理是最小二乘法。這里用來求解狀態(tài)變量的初值，待求解的微分方程組由式(8)、(16)以及如下各有功、無功之間的關(guān)系式構(gòu)建：

(21)

式中：PS0、QS0為定子初始有功、無功功率；Pr0、Qr0、Pg0和Qg0分別為轉(zhuǎn)子側(cè)、網(wǎng)側(cè)換流器初始有功、無功功率；PL0、QL0為傳輸線路初始有功、無功功率。各有功、無功功率流向如圖1所示。

在Matlab編寫.m文件，將上述求解初值的微分方程寫入，利用fsolve命令即可得到系統(tǒng)變量的初值。需要注意的是，除表2中所列出的初值外，其他初值，如葉輪轉(zhuǎn)速查表1可得。

表2 風(fēng)速7 m/s、串補(bǔ)30%系統(tǒng)狀態(tài)變量初值Table 2 State variable initial value at Ww=7 m/s, K=30%

linmod是提取系統(tǒng)穩(wěn)態(tài)運(yùn)行點(diǎn)附近、連續(xù)時間線性狀態(tài)空間模型的語句。在將2.1節(jié)所建模型線性化的同時，可得到系統(tǒng)A、B、C和D矩陣。

針對狀態(tài)矩陣A，利用Matlab中的eig語句便可得系統(tǒng)的特征值。表3為當(dāng)前工況下系統(tǒng)各模態(tài)及其對應(yīng)的特征值。

2.3 參與因子

求解參與因子離不開特征值所對應(yīng)的左右特征向量。第i個左特征向量表征了狀態(tài)變量對第i個模態(tài)的貢獻(xiàn)大小，體現(xiàn)的是“可觀性”；第i個右特征向量表征的是狀態(tài)變量中第i個模態(tài)的活躍程度，體現(xiàn)的是“可控性”。參與因子可通過下式得到：

表3 風(fēng)速7 m/s、串補(bǔ)30%下的系統(tǒng)特征值Table 3 System eigenvalues at Ww=7 m/s, K=30%

(22)

式中：Pki為第k個狀態(tài)變量與第i個模態(tài)的參與因子，反映了狀態(tài)變量對模態(tài)的強(qiáng)可觀強(qiáng)可控，是一個綜合性指標(biāo)；Zi和Yi分別為第i個模態(tài)的左右特征向量。

值得注意的是，通常所說的特征向量一般指右特征向量。而針對左特征向量的求解有兩種方法:一種是先求狀態(tài)矩陣A的轉(zhuǎn)置的右特征向量，再將其轉(zhuǎn)置即為左特征向量；另一種方法是直接將右特征向量求逆再轉(zhuǎn)置即為左特征向量。由于代數(shù)精度問題，Matlab中二者求得的左特征向量會有細(xì)微差別，但是不影響問題的分析。本文中采用第1種方法求取左特征向量。

表4為部分系統(tǒng)狀態(tài)變量對應(yīng)的參與因子，為便于后續(xù)模態(tài)分析，其中較大值已加粗表示。

表4 部分系統(tǒng)狀態(tài)變量對應(yīng)的參與因子Table 4 Partial participation factors of state variables

3 模態(tài)分析

3.1 機(jī)電模態(tài)和軸系模態(tài)

當(dāng)風(fēng)速為7 m/s時，由表1可知葉輪轉(zhuǎn)速為0.75 pu，換成有名值即為45 Hz(60×0.75=45 Hz)。在該風(fēng)速下線路串補(bǔ)為30%時，模態(tài)λ7,8對應(yīng)的特征值為-6.753±j93.821，對應(yīng)頻率為93.821 rad/s，即14.93 Hz，剛好與45 Hz互補(bǔ)(60-14.93=45.07 Hz)，其他風(fēng)速下亦如此。再而，由表4可知此模態(tài)與定轉(zhuǎn)子dq軸電流呈強(qiáng)相關(guān)關(guān)系。因此，模態(tài)λ7,8與機(jī)械動態(tài)(風(fēng)速、葉輪轉(zhuǎn)速)和電氣動態(tài)(定、轉(zhuǎn)子電流)相關(guān)，故將其識別為機(jī)電模態(tài)。另外，通過在表5中縱向觀察比較，可發(fā)現(xiàn)此模態(tài)受串補(bǔ)度影響較小。

表5 不同風(fēng)速、串補(bǔ)度下模態(tài)λ7,8的特征值Table 5 Eigenvalues of mode λ7,8 at different Ww and K

由表4可知λ12,13與葉輪轉(zhuǎn)速ωt以及兩質(zhì)量塊間轉(zhuǎn)矩Ttg強(qiáng)相關(guān)，因此將其識別為軸系模態(tài)，軸系振蕩頻率很低，約0.95 Hz。

3.2 次同步和超同步模態(tài)

系統(tǒng)的自然振蕩頻率為

(23)

式中：fs為系統(tǒng)工頻；XC為系統(tǒng)容抗；XL為系統(tǒng)感抗。次同步振蕩頻率為fs-fn，超同步振蕩頻率為fs+fn。

因此，串補(bǔ)為30%時的自然頻率為32.86 Hz，次同步頻率為27.14 Hz(60-32.86=27.14 Hz)，超同步頻率為92.86 Hz。與之對應(yīng)的是模態(tài)λ5,6，頻率為28.25 Hz(177.5/2π)，模態(tài)λ3,4頻率為91.53 Hz。故將λ5,6識別為次同步模態(tài)，λ3,4識別為超同步模態(tài)。

表6列出了不同風(fēng)速與串補(bǔ)度下的次/超同步模態(tài)，可看出：超同步模態(tài)實(shí)部均為負(fù)；而次同步模態(tài)在7 m/s，串補(bǔ)度由45%變?yōu)?0%時，實(shí)部變?yōu)檎龜?shù)，振蕩發(fā)散，變?yōu)椴环€(wěn)定。也因此，少有人研究超同步振蕩，而把重點(diǎn)放在次同步振蕩上。

圖8為不同風(fēng)速、串補(bǔ)度下次同步模態(tài)的頻率和阻尼比關(guān)系圖，K由25%開始，間隔為5%依次遞增至90%?？捎^察到K增加或風(fēng)速減小，都會導(dǎo)致阻尼比減小，使系統(tǒng)存在次同步振蕩風(fēng)險。

表6 不同風(fēng)速及串補(bǔ)度下次/超同步模態(tài)特征值Table 6 Eigenvalues of SSR/SupSR modeat different Ww and K

圖8 次同步模態(tài)隨風(fēng)速及串補(bǔ)度的變化關(guān)系圖Fig.8 SSR mode with variation of Ww and K

圖9 次同步模態(tài)隨控制參數(shù)變化關(guān)系圖Fig.9 SSR mode with variation of control parameters

圖9反映改變轉(zhuǎn)子側(cè)換流器控制參數(shù)對次同步模態(tài)的影響，分別使Kp和Ki參數(shù)中的一個不變，增加另一個的值?？砂l(fā)現(xiàn)，當(dāng)控制參數(shù)增加到一定數(shù)值后，阻尼比呈現(xiàn)負(fù)值，系統(tǒng)次同步模態(tài)不穩(wěn)定。

事實(shí)上，文獻(xiàn)[13-15]都已證實(shí)風(fēng)速、串補(bǔ)度及控制參數(shù)等對雙饋風(fēng)機(jī)次同步的影響，本文對此的分析旨在驗(yàn)證一般性結(jié)論，證明了所搭建模型的合理性。

4 結(jié)語

本文列寫了雙饋風(fēng)機(jī)單機(jī)無窮大系統(tǒng)的動態(tài)微分方程組，包括軸系、發(fā)電機(jī)、串補(bǔ)網(wǎng)絡(luò)和控制環(huán)節(jié)等在內(nèi)的22個狀態(tài)變量?；赟imulink平臺搭建了用于雙饋風(fēng)機(jī)特征值分析的fcn模型，詳細(xì)介紹了包括模型搭建，及運(yùn)行初值、特征值與參與因子計算在內(nèi)的過程。通過模態(tài)分析，識別出了機(jī)電模態(tài)、軸系模態(tài)和次/超同步模態(tài)，并通過改變風(fēng)速、串補(bǔ)度和控制參數(shù)，驗(yàn)證了關(guān)于次同步振蕩的一般性結(jié)論。所做工作對于今后的研究有基礎(chǔ)性意義。