馬 靜,彭明法,王 形,楊奇遜

（華北電力大學(xué) 新能源電力系統(tǒng)國(guó)家重點(diǎn)實(shí)驗(yàn)室,北京 102206）

0 引言

在故障切除后，由于電力系統(tǒng)的時(shí)變和非線性因素顯著增強(qiáng)，傳統(tǒng)的低頻振蕩分析方法很難準(zhǔn)確反映事故后特征值的變化情況[1-2]。因此，如何從故障系統(tǒng)中提取時(shí)變非線性的振蕩特性是亟待解決的難題[3-4]。

平衡點(diǎn)特征值法可以描述定長(zhǎng)線性系統(tǒng)的全局動(dòng)態(tài)行為，并且可以定量解釋低頻振蕩的機(jī)理[5-6]，但是隨著系統(tǒng)規(guī)模的日益擴(kuò)大以及結(jié)構(gòu)日趨復(fù)雜，利用平衡點(diǎn)特征值法反映系統(tǒng)的振蕩特性愈加困難。特別是系統(tǒng)在強(qiáng)時(shí)變、強(qiáng)非線性運(yùn)行條件下，利用平衡點(diǎn)特征值法很可能得到完全錯(cuò)誤的結(jié)論。為此，眾多專家和學(xué)者提出了與受擾軌跡密切相關(guān)的軌跡特征值法。有別于平衡點(diǎn)特征值法，軌跡特征值法沿著實(shí)際受擾軌跡，在不同時(shí)間窗口中獲取振蕩頻率和阻尼系數(shù)的時(shí)間序列[7]。軌跡特征值的時(shí)間序列可以反映時(shí)變非線性系統(tǒng)在指定擾動(dòng)下被激發(fā)的振蕩模式及其隨時(shí)間的變化情況[8]。但軌跡特征值法的研究主要集中于確定性擾動(dòng)場(chǎng)景，對(duì)于系統(tǒng)在故障類隨機(jī)因素下（如故障持續(xù)時(shí)間和故障位置），特征值變化情況的研究則相對(duì)較少。且該方法對(duì)于窗口寬度的選取依賴性較強(qiáng)，當(dāng)窗口過(guò)寬時(shí)，將會(huì)掩蓋系統(tǒng)的時(shí)變振蕩特性，而過(guò)窄則將會(huì)大幅度增加計(jì)算量。

本文將平衡點(diǎn)特征值法與軌跡特征值法相結(jié)合，提出了一種適用于分析故障系統(tǒng)低頻振蕩特征值的模型軌跡法。首先，利用共享因子將故障后的系統(tǒng)進(jìn)行等值，再利用原網(wǎng)絡(luò)支路電流和故障分量網(wǎng)絡(luò)支路電流，求取故障后各支路電流以及功率，在此基礎(chǔ)上，建立系統(tǒng)的模型并在各測(cè)點(diǎn)處分段線性化，將分段線性化模型與系統(tǒng)的受擾軌跡相結(jié)合，求解系統(tǒng)特征值的時(shí)間序列。最后，利用2機(jī)系統(tǒng)和16機(jī)系統(tǒng)，分析故障切除后系統(tǒng)特征值和阻尼比的變化情況，算例驗(yàn)證了該方法的正確性和有效性。

1 故障系統(tǒng)的穩(wěn)定性分析

1.1 故障模型的等值

電力系統(tǒng)的節(jié)點(diǎn)電壓方程可描述為[9]：

其中，下標(biāo)N表示發(fā)電機(jī)節(jié)點(diǎn)；下標(biāo)M表示聯(lián)絡(luò)節(jié)點(diǎn)（包含負(fù)荷節(jié)點(diǎn)）；UN為發(fā)電機(jī)節(jié)點(diǎn)電壓；UM為聯(lián)絡(luò)節(jié)點(diǎn)（包含負(fù)荷節(jié)點(diǎn)）電壓；IN為發(fā)電機(jī)節(jié)點(diǎn)注入電流。

消去聯(lián)絡(luò)節(jié)點(diǎn)[10]，可得僅包含發(fā)電機(jī)節(jié)點(diǎn)的節(jié)點(diǎn)電壓方程：

其中，Y′NN=YNN-YNMY-1MMYMN為系統(tǒng)的縮減導(dǎo)納陣。

由式（2）可得各發(fā)電機(jī)注入節(jié)點(diǎn)的電流。網(wǎng)絡(luò)支路電流Ib與節(jié)點(diǎn)注入電流IN的關(guān)系如下：

其中，Y為支路導(dǎo)納陣；A為節(jié)點(diǎn)關(guān)聯(lián)陣；Yn為全網(wǎng)絡(luò)節(jié)點(diǎn)導(dǎo)納陣；C=YATYn-1為共享因子矩陣。

圖1為等值2機(jī)系統(tǒng)電路圖。在t0時(shí)刻，系統(tǒng)聯(lián)絡(luò)線1上某點(diǎn)發(fā)生三相短路，其中短路點(diǎn)距節(jié)點(diǎn)B1的距離占整條線路全長(zhǎng)的比例為α。經(jīng)過(guò)短路持續(xù)時(shí)間tc后將故障線路切除，切除后的網(wǎng)絡(luò)如圖2所示。

圖1 2機(jī)系統(tǒng)線路圖Fig.1 Two-generator system

圖2 系統(tǒng)故障支路切除后的線路圖Fig.2 Two-generator system after faulty line is tripped

設(shè)故障前發(fā)電機(jī)的節(jié)點(diǎn)注入電流為IN1，線路切除后發(fā)電機(jī)的注入電流發(fā)生變化，記為IN2。根據(jù)疊加原理，線路切除后網(wǎng)絡(luò)中的潮流分布可以看成由線路切除前的支路潮流和線路切除所引起的故障分量?jī)刹糠纸M成，即圖2可由圖1和圖3疊加等值而成[11]。需要注意的是，此時(shí)線路切除前的（原）網(wǎng)絡(luò)注入電流為IN2。

圖3 僅含故障分量的等值網(wǎng)絡(luò)Fig.3 Equivalent system only containing fault components

當(dāng)節(jié)點(diǎn)注入電流為IN2時(shí)，由式（3）可知，支路i上的電流為：

其中，I′i為支路 i上的電流；Ci為矩陣 C 的第 i行。

當(dāng)系統(tǒng)支路i發(fā)生故障切除后，系統(tǒng)的故障分量網(wǎng)絡(luò)僅有一個(gè)激勵(lì)源，即I′i，因此其他線路上的潮流僅與該激勵(lì)源有關(guān)。此時(shí)，網(wǎng)絡(luò)的節(jié)點(diǎn)電壓方程可寫(xiě)為：

其中，Y′n為故障等值網(wǎng)絡(luò)的節(jié)點(diǎn)導(dǎo)納陣。

由式（5）可知，該激勵(lì)源注入各支路的電流為：

其中，I″為故障分量網(wǎng)絡(luò)中各支路的電流；Y′為該網(wǎng)絡(luò)的支路導(dǎo)納陣；A′為節(jié)點(diǎn)關(guān)聯(lián)陣；D為共享因子矩陣C=YATYn-1的最后一列。

根據(jù)共享因子矩陣，當(dāng)故障支路切除后，支路k的電流可用下式表示：

其中，Dk為矩陣C第k行的最后一列。

將式（2）代入式（7），可得由發(fā)電機(jī)節(jié)點(diǎn)電壓表示的支路電流：

對(duì)于圖1所示的2機(jī)系統(tǒng)，當(dāng)支路1發(fā)生故障切除后，支路2上的電流發(fā)生變化，由前面推導(dǎo)可得支路2的電流為：

此時(shí)，兩發(fā)電機(jī)之間傳送的功率可由式（10）表示：

將故障切除后的網(wǎng)絡(luò)進(jìn)行等值，通過(guò)原網(wǎng)絡(luò)的節(jié)點(diǎn)電壓和故障分量網(wǎng)絡(luò)的支路電流，求取故障切除后網(wǎng)絡(luò)各支路的電流和功率，以此對(duì)系統(tǒng)的特征值進(jìn)行分析。

1.2 故障系統(tǒng)時(shí)變振蕩特性分析

以2機(jī)系統(tǒng)為例，分析系統(tǒng)在故障切除后的特征值問(wèn)題，系統(tǒng)的結(jié)構(gòu)如圖1所示。引入2機(jī)系統(tǒng)經(jīng)典轉(zhuǎn)子運(yùn)動(dòng)方程，并忽略機(jī)械功率的變化，則系統(tǒng)經(jīng)過(guò)小擾動(dòng)線性化后的方程可用下式表示[12]：

其中，δi（i=1，2）為第 i臺(tái)發(fā)電機(jī)轉(zhuǎn)子旋轉(zhuǎn)的電角度；ω0為同步速；ωi（i=1，2）為第 i臺(tái)發(fā)電機(jī)轉(zhuǎn)子的電角速度；Mi（i=1，2）為第 i臺(tái)發(fā)電機(jī)的慣性時(shí)間常數(shù)；Pei（i=1，2）為第 i臺(tái)發(fā)電機(jī)的電磁功率；KDi（i=1，2）為第i臺(tái)發(fā)電機(jī)的阻尼系數(shù)。

根據(jù)功率平衡，發(fā)電機(jī)輸出的電磁功率可表示為[13-14]：

其中，PL1、PL2分別為負(fù)荷 L1、L2的有功功率；P12、P21為發(fā)電機(jī)G1與G2之間的傳輸功率。

發(fā)電機(jī)輸出的功率經(jīng)過(guò)小擾動(dòng)線性化后的表達(dá)式如下：

將式（13）代入式（11），可得故障線路切除后系統(tǒng)的狀態(tài)方程：

由式（15）可得，此時(shí)系統(tǒng)的狀態(tài)矩陣A為：

其中，I 為 2×2 維單位陣；k11t、k12t、k22t、KD1t、KD2t分別為 k11、k12、k22、KD1、KD2的時(shí)變值。

將系統(tǒng)相應(yīng)時(shí)刻電壓的實(shí)測(cè)值代入式（16），可得該時(shí)刻系統(tǒng)的振蕩特征值，據(jù)此，可得特征值的完整時(shí)間序列。

需要注意的是，當(dāng)系統(tǒng)發(fā)生擾動(dòng)的時(shí)候，發(fā)電機(jī)節(jié)點(diǎn)的電壓幅值和相角不再是常量，由式（16）可知，系統(tǒng)的狀態(tài)矩陣A中部分元素也會(huì)發(fā)生改變，同時(shí)，系統(tǒng)的特征值亦發(fā)生變化。再者，隨著電壓幅值和相角變化幅度的增加，特征值的變化幅度也將相應(yīng)增加。

2 算例分析

2.1 等值2機(jī)系統(tǒng)算例仿真

2機(jī)系統(tǒng)采用經(jīng)典二階模型，系統(tǒng)在0.1 s時(shí)聯(lián)絡(luò)線1上某點(diǎn)處發(fā)生三相短路。故障點(diǎn)距節(jié)點(diǎn)1的距離占線路全長(zhǎng)的比例為α，經(jīng)過(guò)故障持續(xù)時(shí)間tc后，將故障線路切除?？紤]以下幾種故障情況，分析故障系統(tǒng)的時(shí)變振蕩特性。

a.比例α取0.5，故障持續(xù)時(shí)間tc為0.05 s。

b.比例α取0.9，故障持續(xù)時(shí)間tc為0.05 s。

c.比例α取0.5，故障持續(xù)時(shí)間tc為0.1 s。

d.比例α取0.9，故障持續(xù)時(shí)間tc為0.1 s。

通過(guò)仿真得到系統(tǒng)受擾軌跡后，將相應(yīng)的實(shí)測(cè)值代入式（16），計(jì)算系統(tǒng)實(shí)時(shí)特征值的變化情況[15]，如圖4—7所示。

由圖4和圖5可知，當(dāng)短路持續(xù)時(shí)間為0.05 s，故障位置改變時(shí)兩者故障切除后的特征值變化影響較小。其中，精確特征值的實(shí)部變化區(qū)間為-0.1395～-0.1380，虛部變化區(qū)間為3.1～6.7，而利用本文方法計(jì)算得到的實(shí)部區(qū)間為-0.1412～-0.1378，虛部區(qū)間為2.5～6.7。由圖6和圖7可知，當(dāng)短路持續(xù)時(shí)間為0.1 s時(shí)，此時(shí)故障位置改變對(duì)故障切除后特征值的影響較大，特征值在實(shí)數(shù)與復(fù)數(shù)之間交替變化：當(dāng)α為0.5時(shí)，精確特征值實(shí)部區(qū)間為-0.14～2.1，虛部區(qū)間為0～7，而用本文方法計(jì)算得到的實(shí)部區(qū)間為-0.14～2.2，虛部區(qū)間為 0～7；當(dāng) α 為 0.9時(shí)，精確特征值實(shí)部區(qū)間為-0.14～3.8，虛部區(qū)間為0～7，而用本文方法計(jì)算得到的實(shí)部區(qū)間為-0.14～3.8，虛部區(qū)間為0～7，雖然兩者的變化區(qū)間相近，然而，在同一時(shí)刻兩者的特征值差異較大。由此可知，不同故障類隨機(jī)因素對(duì)故障切除后的特征值將產(chǎn)生不同影響。一方面，系統(tǒng)在故障切除后，由本文方法計(jì)算得到故障切除后的特征值不再是常量，并且變化幅度較為明顯。另一方面，系統(tǒng)故障切除后的特征值隨著短路持續(xù)時(shí)間和故障位置的不同，其變化差異較大。

圖4 α=0.5、tc=0.05 s時(shí)仿真結(jié)果Fig.4 Simulative results when α is 0.5 and tcis 0.05 s

圖5 α=0.9、tc=0.05 s時(shí)仿真結(jié)果Fig.5 Simulative results when α is 0.9 and tcis 0.05 s

圖6 α=0.5、tc=0.1 s時(shí)仿真結(jié)果Fig.6 Simulative results when α is 0.5 and tcis 0.1 s

圖7 α=0.9、tc=0.1 s時(shí)仿真結(jié)果Fig.7 Simulative results when α is 0.9 and tcis 0.1 s

由圖4—7還可以看出，利用本文方法計(jì)算得到的特征值能較精確地反映系統(tǒng)實(shí)際特征值的變化情況，而利用平衡點(diǎn)算法計(jì)算得到的特征值為常數(shù)，不能實(shí)時(shí)反映系統(tǒng)特征值。由平衡點(diǎn)算法計(jì)算得到的偏差可由下式表示：

其中，Ri、Mi、ξi分別為利用第 i時(shí)刻的實(shí)測(cè)值計(jì)算得到的特征值實(shí)部、虛部和阻尼比；R0、M0、ξ0分別為利用平衡點(diǎn)算法得到的特征值的實(shí)部、虛部和阻尼比；Rerr、Merr、ξerr分別表示平衡點(diǎn)法與本文方法的偏差；m為實(shí)測(cè)點(diǎn)個(gè)數(shù)。

表1給出了4種情況下，利用平衡點(diǎn)法與本文方法計(jì)算得到的特征值和阻尼比的相對(duì)偏差。由表1可見(jiàn)，短路持續(xù)時(shí)間越長(zhǎng)，利用2種方法計(jì)算出的特征值和阻尼比偏差越大，故障位置α越大，兩者的偏差也越大。該表同樣表明未考慮故障類隨機(jī)因素的平衡點(diǎn)法不能正確反映故障系統(tǒng)的實(shí)時(shí)特征值變化情況。

表1 平衡點(diǎn)算法的均方誤差Tab.1 Mean square error of equilibrium-point method

2.2 16機(jī)系統(tǒng)仿真算例分析

采用IEEE 16機(jī)68節(jié)點(diǎn)[16-18]的新英格蘭—紐約互聯(lián)系統(tǒng)進(jìn)一步驗(yàn)證本文方法的正確性和有效性，系統(tǒng)結(jié)構(gòu)如圖8所示，該系統(tǒng)可分為五大區(qū)域。其中發(fā)電機(jī)采用六階詳細(xì)模型，勵(lì)磁采用IEEE-DC1型勵(lì)磁。系統(tǒng)在0.1 s時(shí)刻，線路28-29某處發(fā)生三相短路故障，故障點(diǎn)距節(jié)點(diǎn)28的距離占線路全長(zhǎng)的比例為α。分析以下4種故障情況下，故障切除后系統(tǒng)的特征值變化情況。

a.比例α取0.5，故障持續(xù)時(shí)間tc為0.02 s。

b.比例α取0.9，故障持續(xù)時(shí)間tc為0.02 s。

c.比例α取0.5，故障持續(xù)時(shí)間tc為0.05 s。

圖8 16機(jī)68節(jié)點(diǎn)電網(wǎng)結(jié)構(gòu)圖Fig.8 Structure of 16-generator 68-bus system

d.比例α取0.9，故障持續(xù)時(shí)間tc為0.05 s。

將系統(tǒng)仿真得到的軌跡代入式（16）中，由此得到特征值的時(shí)間序列。其中區(qū)域4與區(qū)域5之間振蕩模式的時(shí)間序列如圖9—12所示。

圖9 α=0.5、tc=0.02 s時(shí)仿真結(jié)果Fig.9 Simulative results when α is 0.5 and tcis 0.02 s

圖10 α=0.9、tc=0.02 s時(shí)仿真結(jié)果Fig.10 Simulative results when α is 0.9 and tcis 0.02 s

由圖9—12可知，基于本文方法計(jì)算得到的特征值亦不再為常量，而是隨時(shí)間發(fā)生相應(yīng)的變化，且隨故障持續(xù)時(shí)間和故障位置的不同，其變化情況也各異。同時(shí)，本文方法能夠較精確地跟蹤系統(tǒng)的特征值變化。

表2給出了不同短路持續(xù)時(shí)間和故障位置時(shí)，利用平衡點(diǎn)法與本文方法計(jì)算得到的特征值和阻尼比的相對(duì)偏差。由表2可知，當(dāng)短路持續(xù)時(shí)間越長(zhǎng)時(shí)，利用平衡點(diǎn)法與利用本文方法得到的特征值和阻尼比的偏差越大，這同樣也表明平衡點(diǎn)法不能正確反映故障系統(tǒng)的振蕩特性。

圖11 α=0.5、tc=0.05 s時(shí)仿真結(jié)果Fig.11 Simulative results when α is 0.5 and tcis 0.05 s

圖12 α=0.9、tc=0.05 s時(shí)仿真結(jié)果Fig.12 Simulative results when α is 0.9 and tcis 0.05 s

表2 平衡點(diǎn)算法的均方誤差Tab.2 Mean square error of equilibrium-point method

3 結(jié)論

本文提出了一種分析故障切除后系統(tǒng)低頻特征值的方法，即模型軌跡法，該方法首先將故障切除后的系統(tǒng)進(jìn)行等值，再利用原網(wǎng)絡(luò)支路電流和故障分量網(wǎng)絡(luò)支路電流，求取故障后各支路電流以及功率，在此基礎(chǔ)上，建立等值系統(tǒng)的運(yùn)動(dòng)方程，并在各測(cè)點(diǎn)處將等值模型分段線性化，結(jié)合實(shí)際受擾軌跡，求解系統(tǒng)特征值的時(shí)間序列。最后，利用2機(jī)系統(tǒng)和16機(jī)系統(tǒng)，分析故障類隨機(jī)因素對(duì)系統(tǒng)振蕩特性的影響，算例結(jié)果驗(yàn)證了該方法較平衡點(diǎn)特征值法，能夠更準(zhǔn)確地跟蹤系統(tǒng)特征值的變化情況，同時(shí)，也有效避免了利用軌跡特征根法求解電力系統(tǒng)時(shí)變特征值的窗口選取問(wèn)題。