王永昆，鮑顏紅，2，方勇杰，鄭 偉，徐 偉

（1.國網(wǎng)電力科學(xué)研究院/南京南瑞集團公司，江蘇 南京 210003；2.華北電力大學(xué) 電氣與電子工程學(xué)院，北京 102206；3.中國電力科學(xué)研究院，北京 100192）

0 引言

區(qū)域電網(wǎng)互聯(lián)，特別是特高壓交直流電網(wǎng)的建設(shè)以及大規(guī)模風(fēng)電的接入使得電網(wǎng)的低頻振蕩問題非常突出，進行電力系統(tǒng)小干擾穩(wěn)定在線分析對預(yù)防互聯(lián)電網(wǎng)弱阻尼引起的低頻振蕩和連鎖故障具有重要的現(xiàn)實意義[1-2]。

自20世紀(jì)80年代以來，國內(nèi)外學(xué)者對大型電力系統(tǒng)特征值計算的方法進行了大量的研究[3-5]，目前小干擾穩(wěn)定特征值離線計算方法已經(jīng)很成熟。但是基于離線數(shù)據(jù)的分析計算并不能反映系統(tǒng)時變的特性，為滿足在線小干擾分析[6-7]的要求，提出了并行計算和分布式計算的方法。文獻[8]介紹了BI迭代法、BI迭代法/逆迭代混合算法的并行化處理，該算法是求解一個特征值的并行實現(xiàn)方法，對原有串行代碼改動較大，實際的加速比并不高。文獻[9]提出了按區(qū)域分網(wǎng)的分布式并行算法，該算法將每個區(qū)域設(shè)為一個子網(wǎng)，子網(wǎng)僅計算本區(qū)域數(shù)據(jù)，通過主控機的協(xié)調(diào)和交互來完成全網(wǎng)的小干擾穩(wěn)定分析，但是該算法也需要對原有串行算法進行大量改動，并且網(wǎng)絡(luò)劃分越不均，則各節(jié)點機相互等待時間越長，計算效率也不高。文獻[10]將AESOPS算法與帶邊界分區(qū)的互聯(lián)電網(wǎng)切分方法相結(jié)合，提出了分布式的特征值算法，但該算法在迭代過程中需要多次交換數(shù)據(jù)，網(wǎng)絡(luò)通信的延遲是分布式計算主要的性能瓶頸。文獻[11]基于多核CPU并行計算技術(shù)提出了多進程的特征值并行搜索算法，該算法將特征值搜索點根據(jù)CPU核的數(shù)量和設(shè)定的任務(wù)分配策略分配給各進程，各CPU核的計算與串行計算完全相同，不但保持了串行程序的完整性，而且保證了計算效率。

本文在文獻[11]的理論基礎(chǔ)上，提出了一種基于位移求逆隱式重啟動Arnoldi算法IRAM（Implicitly Restarted Arnoldi Method）的特征值在線并行計算策略，利用前一時刻斷面分析結(jié)果對特征值位移點（即搜索點）進行了優(yōu)化選取，并且介紹了一種新的任務(wù)并行分配策略，可進一步提高計算速度。

1 位移求逆IRAM

IRAM在收斂性能和計算速度上較優(yōu)越[12-13]，但一次位移求逆IRAM只能計算給定位移點附近一定數(shù)量的特征值，為求解搜索區(qū)域內(nèi)所有關(guān)鍵特征值，需要執(zhí)行多次位移求逆IRAM，以使搜索圓盤盡量覆蓋整個區(qū)域[14]（復(fù)平面上位移點附近的一個圓內(nèi)的特征值已被完全搜索）。目前位移點的選擇方法有中國電力科學(xué)研究院的PSD-SSAP①中國電力科學(xué)研究院系統(tǒng)所.PSD-SSAP電力系統(tǒng)小干擾穩(wěn)定性分析程序用戶手冊.2009.采用的靜態(tài)設(shè)置方法和國網(wǎng)電力科學(xué)研究院的FASTEST②國網(wǎng)電力科學(xué)研究院穩(wěn)定所.FASTEST小擾動穩(wěn)定分析技術(shù)手冊.2007.及上海交通大學(xué)的小干擾穩(wěn)定分析包SSAP[15]采用的動態(tài)設(shè)置方法。前者是對整個頻段等間隔地設(shè)置位移點，它的缺點是需憑經(jīng)驗或反復(fù)嘗試來選擇位移點的數(shù)量，若位移點數(shù)量過少，將導(dǎo)致嚴(yán)重漏根現(xiàn)象，反之則會由于特征值的重復(fù)搜索而浪費時間。后者是在頻段兩側(cè)取位移點，若任意2個圓盤不相交，則在中間點再取1個位移點，直至所有圓盤均相交，但該方法不易在并行計算中實現(xiàn)。

2 特征值在線并行計算策略

定義預(yù)估圓盤為復(fù)平面上以位移點為圓心、半徑r=ω/k的圓，如圖1所示。其中，ω為位移點的虛部，k為給定阻尼比范圍對應(yīng)的恒阻尼線斜率。

其中，ξ為給定阻尼比范圍邊界。

圖1 預(yù)估圓盤示意圖Fig.1 Schematic diagram of predict circle

系統(tǒng)結(jié)構(gòu)和運行參數(shù)的變化對特征值的阻尼影響較大，對特征值的頻率影響較小[16]。假定系統(tǒng)運行方式?jīng)]有發(fā)生大的改變，則可根據(jù)前一時刻斷面的特征值計算結(jié)果，得到指定位移點的預(yù)估圓盤內(nèi)的特征值數(shù)量，進而為當(dāng)前時刻斷面計算時位移點及特征值數(shù)量的選取提供依據(jù)。

根據(jù)電力系統(tǒng)的特點，通常在小干擾穩(wěn)定分析中區(qū)內(nèi)振蕩模式較多，而區(qū)間振蕩模式較少[10]，因此中高頻段位移點的預(yù)估圓盤半徑與實際圓盤半徑相差不大。對于小型電網(wǎng)，特征值分布稀疏，預(yù)估圓盤與實際圓盤誤差較大，可能會導(dǎo)致重復(fù)搜索或遺漏關(guān)鍵根，效率降低。

2.1 主位移點的選擇策略

位移點的選擇和數(shù)量對并行算法的性能起著重要作用。為減少特征值的重復(fù)搜索，選擇的位移點應(yīng)滿足預(yù)估圓盤相交面積較少；為盡可能多地搜索到特征值，選擇的位移點應(yīng)滿足預(yù)估圓盤相交交點距復(fù)平面虛軸較遠。已知第1個位移點位置及預(yù)估圓盤半徑，可將第2個位移點位置的選擇，轉(zhuǎn)化為求解兩圓盤相交面積極小且交點阻尼比絕對值極大的多目標(biāo)優(yōu)化問題。

2.1.1 主位移點選擇的數(shù)學(xué)模型

如圖2所示，搜索區(qū)域內(nèi)兩相鄰位移點的搜索圓盤相交，其圓心（即位移點）分別為B和C，半徑分別為r1和r2，圓心距d=BC，交點A到兩圓心所在直線的距離h=AD，兩圓相交部分（即陰影部分）面積為S。

圖2 主位移點選擇示意圖Fig.2 Schematic diagram of main searching point selection

由圖2可建立如下等式：

在復(fù)平面中又有：

由式（1）—（3）消去 d 和 r2，可定義：

其中，S（h）為兩圓相交部分面積，是關(guān)于h的一元表達式，由式（1）—（3）消去 d 和 r2得來。

計算目標(biāo)是使 f1（h）盡量小，而 f2（h）盡量大。

權(quán)系數(shù)均取 0.5，且對 f1（h）、f2（h）進行無量綱化處理，根據(jù)線性加權(quán)和法建立評價函數(shù)如下：

求解 max f（h）最優(yōu)值 h，進而求出 d，即可得到第2個位移點。

2.1.2 主位移點的選擇策略

根據(jù)第2.1.1節(jié)位移點選擇方法計算指定頻率及阻尼比范圍內(nèi)的各位移點和預(yù)估圓盤半徑，然后基于前一時刻斷面分析結(jié)果計算各預(yù)估圓盤內(nèi)特征值數(shù)量。

當(dāng)特征值數(shù)量較少時預(yù)估圓盤與實際圓盤會存在較大誤差，實際圓盤極有可能出現(xiàn)嚴(yán)重重疊，因此為了減少重復(fù)計算量，需要對位移點進行修正，即合并預(yù)估圓盤內(nèi)特征值數(shù)量較少的位移點，取中點替代被合并的位移點，取預(yù)估特征值數(shù)量之和作為修正位移點的特征值數(shù)量。對于高頻段（0.7～2.0 Hz），特征值分布密集且搜索圓盤較大，建議合并2個位移點；而對于低頻段（0.1～0.7 Hz），特征值分布稀疏且搜索圓盤較小，可以考慮多合并1個位移點。由于預(yù)估圓盤與實際圓盤總會存在誤差，為防止漏根，建議采取保守計算，即對于主位移點實際計算的特征值數(shù)量n稍多于預(yù)估特征值數(shù)量，例如取n=+2。

2.2 輔助位移點的選擇策略

由于小干擾穩(wěn)定較關(guān)注的頻段范圍內(nèi)中高段附近振蕩模式較多，在線計算時預(yù)估圓盤半徑與實際圓盤半徑相差不大，可根據(jù)預(yù)估圓盤交點阻尼比與給定阻尼比來選擇輔助位移點。如圖3所示，輔助位移點取線段AF靠近F的黃金分割點E處作為輔助位移點，預(yù)估圓盤半徑為AE，對于輔助位移點的實際計算特征數(shù)取+3。

對于低頻段，由于特征值分布較少且每次計算時都進行了保守計算，搜索圓盤可以完全覆蓋搜索區(qū)域，因此不必增設(shè)輔助位移點。

圖3 輔助位移點選擇示意圖Fig.3 Schematic diagram of assistant searching point selection

2.3 任務(wù)分配策略

在當(dāng)前的在線動態(tài)安全分析系統(tǒng)中，通常利用分布式并行計算技術(shù)，采用同構(gòu)的計算節(jié)點組成計算集群，在集群計算平臺的管理節(jié)點上形成多個計算方案，分配到計算節(jié)點上進行并行計算，然后在管理節(jié)點上匯總計算結(jié)果。

文獻[9]所提不均分配策略由于受計算任務(wù)隨機分配的影響，CPU利用率也具有很大的隨機性，本文將對其進行改進。

假設(shè)位移點數(shù)量為N，CPU核的數(shù)量為M，當(dāng)N?M時，參與并行計算的CPU核應(yīng)盡量多；當(dāng)N與M相差不大或N＜M時，考慮其他在線分析任務(wù)的需求，參與并行計算的CPU核的數(shù)量以小于N為宜。

IRAM的計算時間與位移點的位置和計算特征值的數(shù)量有關(guān)?？烧J為每個位移點計算量的大小主要取決于每次計算特征值的數(shù)量，數(shù)量越多計算量越大；對于特征值數(shù)量相同的位移點，計算量大小則取決于位移點所處位置，所處位置頻率越低計算量越大。基于以上原則可估計出各位移點的計算量相對大小。

本文按預(yù)估計算量的大小分配任務(wù)，具體為：基于預(yù)估計算量將各任務(wù)壓入棧中，計算量大的位移點處于棧頂，計算量小的位移點處于棧底，首先根據(jù)CPU核的數(shù)量分配計算量大的位移點，剩余的位移點作為預(yù)留任務(wù)；當(dāng)某個CPU核的計算任務(wù)完成后，即到預(yù)留任務(wù)中按出棧順序自動提取1個位移點進行計算，直至所有任務(wù)計算完成。該分配策略可近似實現(xiàn)任務(wù)計算時間均分，最大限度地減少CPU核的閑置時間。

2.4 漏根判斷

由特征值計算結(jié)果很容易得出各位移點的實際圓盤半徑。若兩相鄰主位移點的圓盤不相交，則存在嚴(yán)重漏根現(xiàn)象，取兩相鄰主位移點的中點為新的位移點；若相交但交點不在臨近輔助位移點圓盤內(nèi)（高頻段）或交點阻尼比小于給定阻尼比（低頻段），則交點附近存在漏根，取圓盤交點為新的位移點。

2.5 計算步驟

特征值在線并行計算的整體計算流程如圖4所示。首次特征分布可采用現(xiàn)有的靜態(tài)設(shè)置位移點選擇方法和任務(wù)分配策略（平均分配或不平均分配策略）來得到或者采用離線分析結(jié)果。

圖4 特征值在線并行計算整體步驟Fig.4 Overall flowchart of online parallel eigenvalue calculation

3 仿真分析

測試算例采用南方電網(wǎng)2012年5月7日的在線數(shù)據(jù)：前一斷面時刻t1和當(dāng)前斷面時刻t2分別為23:49:34和23:54:35；搜索頻率范圍為0.1～2.0 Hz；阻尼比為 -0.1～0.1。

對t1時刻斷面進行特征值分析，結(jié)果如圖5所示，圖中橫、縱軸分別為阻尼和振蕩角頻率，ε＝ξ。

圖5 t1時刻特征值分布圖Fig.5 Eigenvalue distribution at t1

根據(jù)第2節(jié)的位移點選擇策略選擇的位移點如表1所示，表中一對實部和虛部（均保留2位有效數(shù)字）組成1個位移點，n為對應(yīng)位移點處欲計算的特征值數(shù)量。

根據(jù)表1中數(shù)據(jù)對當(dāng)前（t2時刻）斷面進行位移求逆IRAM特征值搜索計算，結(jié)果如圖6所示，為便于作圖，將虛軸和各圓盤整體進行了縱向壓縮。

表1 位移點選擇情況Tab.1 Selection of searching points

圖6 t2時刻特征值分析結(jié)果Fig.6 Results of eigenvalue analysis at t2

從圖6可以看出，圓盤幾乎覆蓋整個小干擾穩(wěn)定所關(guān)注的頻段及阻尼比范圍，而且減少了特征值的重復(fù)搜索。僅在較高阻尼比且較高頻段內(nèi)會遺漏個別特征值（阻尼比分別為0.0957和0.0983），在實際中這些特征值是非關(guān)鍵的，而且如果搜索這些非關(guān)鍵特征值會耗費更多時間。

表1中位移點的序號是按預(yù)估計算量由大到小的順序排列的，位移點的實際搜索時間如圖7所示。從圖中可以看出，預(yù)估計算量大的位移點計算耗時較大，預(yù)估計算量小的位移點計算耗時較小。由于特征值分布的不均導(dǎo)致某些位移點的預(yù)估計算量誤差較大，但對本文提出的任務(wù)分配策略影響不大。

圖7 各位移點的計算耗時Fig.7 Time consumption of different searching points

在相同計算環(huán)境下，通過計算機模擬雙CPU核下3種任務(wù)并行分配策略。忽略數(shù)據(jù)處理及數(shù)據(jù)輸出時間，各策略的計算時間取決于計算最慢的CPU核，如表2所示。從表2可以看出，本文所提分配策略耗時最短。任務(wù)分配的具體情況對平均分配策略和不均分配策略的計算時間影響很大，但這2種策略均缺乏有效的任務(wù)分配依據(jù)，只能隨機分配任務(wù)，而本文所提分配策略則能較好地解決任務(wù)分配不合理的問題。

表2 不同策略計算時間對比Tab.2 Comparison of computing time among different strategies

同時基于本文所提分配策略也對四川電網(wǎng)和西北電網(wǎng)進行了在線仿真分析，結(jié)果均驗證了本文提出策略的正確性和有效性。

4 結(jié)論

本文基于位移求逆IRAM提出了一種電力系統(tǒng)小干擾特征值在線并行計算策略，該策略保持現(xiàn)有位移求逆IRAM特征值計算程序不變，充分運用在線分析計算結(jié)果，高效、準(zhǔn)確地完成小干擾穩(wěn)定所關(guān)注的計算區(qū)域內(nèi)關(guān)鍵特征值的計算，程序?qū)崿F(xiàn)簡單可靠，可為其他特征值算法的在線并行實現(xiàn)提供參考。