小干擾穩(wěn)定分析中特定模式降階法的快速實現(xiàn)

劉 暢，崔 惟，王克文

（鄭州大學(xué) 電氣工程學(xué)院，河南 鄭州 450001）

0 引言

在電力系統(tǒng)小干擾穩(wěn)定分析中，特征值分析法[1]是得到廣泛應(yīng)用的有效方法之一。通過求解系統(tǒng)狀態(tài)矩陣的特征值、特征向量，不僅可直接判斷系統(tǒng)的穩(wěn)定性，還能得到更多的系統(tǒng)信息，為穩(wěn)定器設(shè)置提供指導(dǎo)[1-4]。

求解全維特征值的QR算法具有魯棒性強(qiáng)、收斂速度快等優(yōu)點，但矩陣階數(shù)較高時，結(jié)果可能存在較大誤差。因此，多種部分特征值方法[5-10]相繼提出。子空間迭代法[5-6]計算模數(shù)較大的一組特征值及特征向量，包括同時迭代法和 Arnoldi法[6]；而降階選擇模式法[7-10]則是將系統(tǒng)矩陣降階處理，使得降階系統(tǒng)的特征值包含原系統(tǒng)的主導(dǎo)特征值，以較小計算量得到所需的特征子集，如AESOPS法[7]、選擇模態(tài)分析法[8-9]等。

概率特征值分析法可以考慮更多的系統(tǒng)運行方式，并得到特征值的統(tǒng)計特性。為了和概率特征值分析法[11]配合使用，文獻(xiàn)[12]描述了特定模式降階法。其大體計算思路與選擇模態(tài)法[9]類似，由于采用了改進(jìn)Rayleigh商逆迭代的完整表達(dá)，可以直接得到較準(zhǔn)確的特征值和特征向量，而迭代初值的確定仍基于全系統(tǒng)狀態(tài)矩陣的QR計算。文獻(xiàn)[13]進(jìn)而考核了降階系統(tǒng)下特征值靈敏度的計算誤差。文獻(xiàn)[12]主要在算法方面實現(xiàn)了降階系統(tǒng)特征值、特征向量的準(zhǔn)確計算，計算效率方面考慮不多。

本文進(jìn)一步完善特定模式降階分析法，提高該方法計算振蕩模式特征值的效率，通過研究降階過程，從共享保留機(jī)組集、迭代算式調(diào)整以及并行計算等方面改進(jìn)，并用70機(jī)和105機(jī)算例進(jìn)行驗證。

1 特定模式降階法

小干擾穩(wěn)定分析中重點關(guān)注的是機(jī)電振蕩模式或阻尼不足的模式。省級電網(wǎng)狀態(tài)矩陣的計算規(guī)模是幾千階，用QR算法計算包括機(jī)電模式在內(nèi)的全維特征值、特征向量，作為準(zhǔn)確計算的初值。

計算第k個特定模式λk時，將狀態(tài)空間方程中的變量分為兩部分：

其中，X1為與模式λk強(qiáng)相關(guān)機(jī)組的狀態(tài)列向量，X2為待消去的狀態(tài)列向量，A1、A2、A3、A4為狀態(tài)矩陣的分塊子矩陣。

由式（1）消去 X2，得：

其中，I為單位矩陣，p為微分算子。

通過 B（p）將式（2）簡記為：

文獻(xiàn)[9]證明了當(dāng)p=λk時降階系統(tǒng)的2個重要性質(zhì)：λk同時為原系統(tǒng)和降階系統(tǒng)的特征值；若降階前、后λk的右特征向量分別為Uk0和Ukr，則Ukr等于Uk0中對應(yīng)X1的元素。因此，降階系統(tǒng)不會畸變原系統(tǒng)的振蕩模式特征值及特征向量相應(yīng)元素[9]。因此，當(dāng) p=λk時，有：

顯然 B（λk）是 λk的函數(shù)，所以由 B（λk）只能準(zhǔn)確計算特定模式λk，其余特征值存在不同程度的偏差。

振蕩模式的特征值 λk是復(fù)數(shù)，B（λk）也是復(fù)矩陣。為簡化表達(dá)，記降階矩陣如式（5）所示，并略去特征值的下標(biāo)k。此外，下文中V和U均對應(yīng)于降階后矩陣的左、右特征向量。

文獻(xiàn)[12]采用三重迭代格式，主要迭代算式如式（6）—（8）所示。

其中，l、m和n分別為外層、次外層和內(nèi)層的迭代次數(shù)，λ（l）為特定模式特征值， τmn為歸一化系數(shù)。

外層迭代中，先按式（6）計算降階矩陣，然后用QR 法求解實部 B1（l），得到對應(yīng) λ（l）的特征值 λ（0）和左、右特征向量 V（n-1）、U（n-1）。

而次外層迭代中，將虛部 jB2（l）作為B的復(fù)數(shù)形式的變化量，按式（7）計算特征值修正量。其中，初始值 V（0）的采用保證了算式的完整表達(dá)[12]。

內(nèi)層迭代，按式（8）計算對應(yīng)的右特征向量 U（n）。

2 降階中保留機(jī)組集的共享

計算降階矩陣的外層迭代式（6）中計及特征值的影響，且式（7）采用完整表達(dá)，因此特定模式降階法得到的特征值具有較高的計算精度。但對于省網(wǎng)規(guī)模電力系統(tǒng)，需關(guān)注的機(jī)電模式或弱阻尼模式較多，計算量仍很可觀。為減少計算量，進(jìn)行如下嘗試。

嘗試1：共享降階矩陣。初值接近的多個特征值采用式（6）形成的同一降階矩陣；導(dǎo)致不同初值往往收斂于同一結(jié)果，不可行。

嘗試2：修正降階矩陣。先按一個特征值形成降階矩陣；分析初值接近的其他模式時，按特征值偏差量修正降階矩陣，由于泰勒級數(shù)固有的截斷誤差，計算誤差偏大。

嘗試3：共享保留機(jī)組。降階計算時需要確定各模式的保留機(jī)組，對應(yīng)于式（1）中的X1。分析各模式的保留機(jī)組情況，對保留機(jī)組相同或相近的模式，采用同一保留機(jī)組，即共享式（6）中的 A1、A2、A3、A4。 部分減少降階矩陣的重復(fù)計算量。

與λk強(qiáng)相關(guān)的機(jī)組構(gòu)成λk的保留機(jī)組。強(qiáng)相關(guān)機(jī)組的確定可依據(jù)參與因子或特征值對PSS增益的靈敏度[12]。出于計算速度的考慮，本文采用參與因子評判各機(jī)組的參與度。

3 降階算法的改進(jìn)實現(xiàn)

3.1 計算效率的改善

降階過程中的主要計算量由三部分組成：式（6）的復(fù)矩陣求逆，式（8）的線性方程組求解，各算式中的矩陣乘法和加法運算。

a.實數(shù)化表達(dá)。

為避免復(fù)數(shù)運算中類型自動轉(zhuǎn)換產(chǎn)生的運算量和相應(yīng)的累計誤差，將算式轉(zhuǎn)換為實數(shù)表達(dá)。

b.采用稀疏技術(shù)。

小干擾穩(wěn)定分析中系統(tǒng)矩陣A的稀疏度與系統(tǒng)結(jié)構(gòu)參數(shù)、運行參數(shù)和控制方式等因素有關(guān)。盡管A的稀疏度不如節(jié)點導(dǎo)納矩陣，但隨著系統(tǒng)規(guī)模的擴(kuò)大，稀疏度有所提高。

除特征向量、逆陣（A4-λ（l）I）-1和部分原始參數(shù)矩陣外，對矩陣進(jìn)行稀疏存儲。

在稀疏存儲基礎(chǔ)上進(jìn)行矩陣乘法和加法運算，并利用因子表技術(shù)[1]求解式（8）的線性方程組。

c.利用OpenMP功能方便實現(xiàn)并行計算。

OpenMP是為共享存儲的多處理機(jī)編寫并行程序而開發(fā)的應(yīng)用編程接口[15-16]，由編譯器自動完成串行程序的并行化處理，可與Fortran、C或C＋＋結(jié)合工作，只要在源代碼中插入編譯指導(dǎo)語句，并進(jìn)行少量修改，就可以方便實現(xiàn)并行計算[15-16]。

對于省級電力系統(tǒng)，分析的機(jī)電模式較多，利用OpenMP的并行計算模式便于同時使用多核CPU。

3.2 用順序高斯消去法進(jìn)行復(fù)矩陣求逆

式（6）中的（A4-λ（l）I）-1表示對復(fù)矩陣求逆，將其按實、虛部分開，并設(shè)其逆矩陣為C＋jD，則有：

寫成實矩陣形式為：

其中，σ、ω 分別為特征值 λ（l）的實部、虛部數(shù)值。

對式（10）中的矩陣C和D，可通過反復(fù)求解線性方程組得到。然而，式（10）的系數(shù)矩陣雖然稀疏，卻不能保證對角優(yōu)勢，因此文獻(xiàn)[12]、[13]中按主元消去法解線性方程組。

特征值的虛部對應(yīng)振蕩頻率，當(dāng)考慮0.1～2.5 Hz范圍的機(jī)電振蕩模式[1]時，ω 取值 0.628～15.7 rad／s，虛部的絕對值不至于太小。所以對式（10）重新排列，將特征值虛部置于對角位置，如式（11）所示。

這樣排列之后，可以采用順序高斯消去法形成因子表，提高計算效率。

3.3 順序消去和主元搜索相結(jié)合的內(nèi)迭代計算

式（8）為內(nèi)迭代算式，其作用是通過迭代得到特征值 λ（m）的右特征向量 U（n）。

由于在內(nèi)迭代時，矩陣 B-λ（m）I恒定，利用因子表技術(shù)提高求解速度。仿照式（11），將式（8）按實、虛部展開，對應(yīng)特征值虛部的子矩陣置于對角元位置，有：

與式（11）相比，式（12）中的對角元素不是僅由特征值虛部構(gòu)成，不能保證對角優(yōu)勢。

將順序消去和主元消去結(jié)合，在某對角位置出現(xiàn)元素絕對值較小的情況下，搜索該列剩余行的主元，將待消去行與主元所在行交換。如此形成的因子表，既避免了計算失敗，又不需要每次消去時都搜索主元。編程實現(xiàn)時考慮到數(shù)據(jù)類型的精度，將消去法的主元門檻值設(shè)為10-5，并利用一維數(shù)組記錄行交換信息。

3.4 矩陣相乘的稀疏處理

針對次外層迭代算式（7），編程時按三元組（矩陣非零元素值、對應(yīng)行號及列號）壓縮存儲矩陣虛部jB2，然后和特征向量相乘；式（6）中，（A4-λ（l）I）-1是稠密矩陣，將 A2、A3稀疏存儲后再和（A4-λ（l）I）-1進(jìn)行計算。由于僅有非零元素參與運算，提高了速度和精度。

3.5 并行計算的采用

隨著電力系統(tǒng)規(guī)模的擴(kuò)大，機(jī)電振蕩模式的數(shù)目相應(yīng)增加。本文借助OpenMP實現(xiàn)多模式特征值的并行計算，對計算程序稍加修改即可實現(xiàn)。

下列示例為Fortran程序中采用OpenMP三線程并行計算的一種格式，其他線程數(shù)類似。

OpenMP功能使并行計算的實現(xiàn)變得簡單，尤其是在多核微機(jī)上同時使用多個CPU。當(dāng)BLOCK1、BLOCK2和BLOCK3計算程序相同但特征值初值不同時，能同時進(jìn)行3個降階模式計算。對于在一個降階模式計算中使用并行計算技術(shù)，需要進(jìn)一步的算式處理，本文尚未考慮。

4 算例分析

算例數(shù)據(jù)來源于某省網(wǎng)系統(tǒng)等值，算例1中發(fā)電機(jī)數(shù)為70臺；算例2為105臺，接近省網(wǎng)規(guī)模。發(fā)電機(jī)采用五階模型，勵磁調(diào)節(jié)系統(tǒng)、原動機(jī)調(diào)速系統(tǒng)采用原有模型。系統(tǒng)矩陣A的有關(guān)信息如表1所示。

表1給出了狀態(tài)矩陣階數(shù)、矩陣非零元個數(shù)等信息，而算例中稀疏度較高的原因，不僅與運行方式有關(guān)，還在于程序中將絕對值小于10-12的元素當(dāng)作零元素處理。

表1 系統(tǒng)狀態(tài)矩陣信息Tab.1 Information of system state matrix

本文計算程序在Visual Studio 2008平臺的Intel Fortran 11環(huán)境下編寫，所有變量都采用雙精度數(shù)據(jù)類型；在個別模塊（如QR計算）采用記錄有效數(shù)字位的方法進(jìn)行乘除運算，減小舍入誤差。程序在主頻為2.13 GHz的雙核計算機(jī)上編譯運行。

用QR法算出狀態(tài)矩陣的初始特征值、特征向量，用特征方程校核，最大偏差小于10-7，因此本文按準(zhǔn)確值處理。所以算例中只需按最外層迭代算式（6）形成降階矩陣一次，然后反復(fù)進(jìn)行次外層迭代和內(nèi)層迭代。

4.1 共享保留機(jī)組的分析

研究機(jī)電模式時，先根據(jù)閾值（本文取各模式最大參與因子絕對值的10%）確定各模式中的強(qiáng)相關(guān)機(jī)組；然后將強(qiáng)相關(guān)機(jī)組相同或較為相近的模式歸為一類；對每一類包含的若干模式，同時保留與之強(qiáng)相關(guān)的所有機(jī)組，得到該類各模式降階分析時可共享的保留機(jī)組。

在算例1中，依據(jù)強(qiáng)相關(guān)機(jī)組相同或相近模式歸為一類的原則，將69個機(jī)電模式分成7類，如表2所示。以表中第1行數(shù)據(jù)為例，當(dāng)計算第1類的15個機(jī)電模式時，只需采用可共享的17臺強(qiáng)相關(guān)機(jī)組，而不必每次重新確定需保留的機(jī)組。

表2 算例1中保留機(jī)組信息Tab.2 Information of retained generators in case 1

由上述結(jié)果可知，只需要較少的保留機(jī)組，就可以實現(xiàn)對所有機(jī)電模式的強(qiáng)相關(guān)機(jī)組劃分。在減少降階計算量的同時，使各模式的分析更加簡潔。

4.2 稀疏技術(shù)的使用效果

形成降階矩陣的外迭代算式（6）中，計算量是矩陣求逆和相乘。為驗證改進(jìn)效果，用2種方法實現(xiàn)。

方法1：按本文所述，用稀疏技術(shù)計算逆陣，用三元組法處理矩陣運算。

方法2：以主元消去法計算逆陣，滿陣存儲矩陣并進(jìn)行相應(yīng)運算。

表3列出算例1中弱阻尼模式的計算數(shù)據(jù)。比較計算結(jié)果，2種方法所得到的降階陣數(shù)值一致；而稀疏技術(shù)的采用，使方法1能較好地節(jié)省時間，算例1中求逆和乘法運算的平均時間為5 s（表3中第3、4列之和）。

表3 算例1中弱阻尼模式降階陣計算數(shù)據(jù)Tab.3 Calculation data of order reduction matrix under low damping mode in case 1

對于降階陣實部的QR計算、次外層迭代式（7）、內(nèi)層迭代式（8），由于矩陣階數(shù)低，耗時較短。算例1中平均用時1.5 s。

因此，對于包含外層迭代、降階陣實部QR計算、次外層迭代、內(nèi)層迭代的3重運算，按方法1的稀疏技術(shù)進(jìn)行計算時，算例1的平均耗時約為6.5 s。

4.3 并行計算效果分析

方法3：在采用稀疏技術(shù)的基礎(chǔ)上，按3.5節(jié)的格式編程，實現(xiàn)并行計算。由于算例分析中采用雙核計算機(jī)，僅可以實現(xiàn)雙線程計算。

表4給出了在算例1中方法1和方法3的計算耗時。算例1計算69個模式，在總耗時中計及了確定共享保留機(jī)組的時間，但不含程序運行中數(shù)據(jù)輸入、輸出時間。

表4 算例1中方法1和3的用時比較Tab.4 Comparison of computation time between method 1 and method 3 in case 1

采用并行計算后，降階法的計算時間有明顯改善。

4.4 算例2

算例2含105臺發(fā)電機(jī)，狀態(tài)矩陣為1177階。

類似地，將104個機(jī)電模式共分成9類，各類中保留機(jī)組的數(shù)目在9～21臺之間，即計算單個機(jī)電模式時的發(fā)電機(jī)數(shù)不超過21臺。

僅采用稀疏技術(shù)時，即方法1，平均1個振蕩模式的計算時間是25.1 s，計算104個機(jī)電模式的總耗時為2 612.6s。

進(jìn)一步結(jié)合雙線程并行計算，即方法3，平均1個模式的計算時間是13.1 s，計算104個機(jī)電模式的總耗時是1358.3 s。

為考核計算效果，本文在算例分析中計算了所有機(jī)電模式。對實際系統(tǒng)，僅需計算關(guān)注的振蕩模式，計算用時會成比例降低；如果進(jìn)一步采用更多線程，例如利用4核或8核微機(jī)，計算時間將進(jìn)一步減少。

5 結(jié)論

特定模式降階法采用完整表達(dá)形成降階矩陣，能夠得出精度較高的特征值。本文對該方法進(jìn)行了改進(jìn)。通過研究降階過程，提出了機(jī)電模式共享保留機(jī)組集以減少降階計算量。調(diào)整迭代算式采用順序消去或?qū)㈨樞蛳ズ椭髟ソY(jié)合，提高計算速度；在矩陣運算中稀疏處理，以節(jié)省計算時間。引入OpenMP的并行計算功能，實現(xiàn)多核計算機(jī)上對降階的快速計算。針對省級電力系統(tǒng)的規(guī)模，在70機(jī)和105機(jī)算例中具體討論了機(jī)電模式共享保留機(jī)組集的情況，并詳細(xì)分析了采用稀疏技術(shù)和并行計算的改進(jìn)效果。