江 涵，江全元

(浙江大學(xué)電氣工程學(xué)院，浙江 杭州 310027)

0 引言

電力系統(tǒng)暫態(tài)穩(wěn)定仿真是進(jìn)行電力系統(tǒng)規(guī)劃、安全穩(wěn)定分析和電力科學(xué)研究的重要工具。電網(wǎng)規(guī)模的擴(kuò)大和仿真模型的日益精細(xì)使其計算量不斷增長。與此同時，為進(jìn)一步提高工作效率以及滿足對電網(wǎng)進(jìn)行更快的安全控制，對仿真速度的要求也不斷提升。在此背景下，并行計算技術(shù)被引入暫態(tài)穩(wěn)定仿真領(lǐng)域，并不斷發(fā)展[1-2]。

近幾十年來，隨著并行計算技術(shù)的革新和進(jìn)步，眾多科研工作者提出了多種暫態(tài)穩(wěn)定并行算法，一般可分為空間并行算法，時間并行算法和引自VLSI領(lǐng)域的波形松弛算法三類[3]。算法的數(shù)學(xué)本質(zhì)都是使用并行方法解決一個大規(guī)模微分-代數(shù)方程組初值問題，對于空間并行算法，一個主要思路是基于串行聯(lián)立算法，首先將整個系統(tǒng)形成一個統(tǒng)一的非線性方程組，進(jìn)而通過牛頓法將方程組線性化，再進(jìn)行并行處理，代表性的算法有并行VDHN法，SOR-Newton 法[4]， Maclaurin-Newton 法[5]等。另一條思路是基于交替算法，直接將微分與代數(shù)方程分開處理，兩部分交替計算，其并行研究主要集中在任務(wù)的平衡分配和網(wǎng)絡(luò)矩陣的處理[6-7]。這些算法的共同之處是都需要處理一個線性代數(shù)方程組。

對線性方程組的求解，直接法（如LU分解）占據(jù)了主流地位。但是，直接法由于其序列性，不能完全地并行化。在很多科研領(lǐng)域，迭代法已經(jīng)顯示出其優(yōu)秀的應(yīng)用性質(zhì)。在眾多迭代算法中，共軛梯度（Conjugate Gradient Method，CG）及其改進(jìn)算法[8]具有快速收斂、存儲量小、不必預(yù)先估計參數(shù)等特點(diǎn)，其計算過程主要是矩陣向量相乘。此類方法近些年在電力系統(tǒng)潮流、仿真等計算中得到了很多研究應(yīng)用。文獻(xiàn)[9]在超立方計算機(jī)上實(shí)現(xiàn)了預(yù)處理 CG算法，應(yīng)用于暫態(tài)穩(wěn)定計算。文獻(xiàn)[10]針對分塊式電力系統(tǒng)暫態(tài)穩(wěn)定并行算法，采用喬列斯基分解的預(yù)處理共軛梯度法用于邊界系統(tǒng)方程組的求解。文獻(xiàn)[11]將ILU共軛梯度法與矩陣方程直接求解法相結(jié)合，求解暫態(tài)穩(wěn)定計算中的整個雅可比矩陣，獲得了較好的加速效果。文獻(xiàn)[12]在Beowulf集群平臺上，提出一種采用多項式預(yù)處理的共軛梯度迭代法求解電力系統(tǒng)分析計算中的稀疏線性方程組，收到了一定的效果。但是，共軛梯度算法針對的是實(shí)對稱矩陣，很多電力系統(tǒng)矩陣并不具備這個特點(diǎn)（如包含移相調(diào)節(jié)變壓器的網(wǎng)絡(luò)導(dǎo)納矩陣），會影響算法的收斂性，導(dǎo)致迭代次數(shù)增加。針對求解不對稱矩陣問題，文獻(xiàn)[13]提出的穩(wěn)定雙共軛梯度法（BiConjugate Gradient Stabilized Method，BiCGSTAB）是一種有效的解決方案。

基本的并行計算硬件平臺主要基于PC集群和共享內(nèi)存的多核計算機(jī)等多 CPU體系結(jié)構(gòu)。近年來，圖形處理器（Graphics Processing Unit，GPU）發(fā)展迅猛。自1999年誕生后, 其發(fā)展的速度是CPU更新速度的3倍多[14]。廣泛應(yīng)用于圖像與視頻處理等領(lǐng)域。2010年，文獻(xiàn)[15]首次將GPU引入電力系統(tǒng)計算，在包含GPU的多核計算平臺上（本文簡稱為GPU計算平臺）使用直接法計算電力系統(tǒng)暫態(tài)穩(wěn)定性。我國學(xué)者同時也展開了相應(yīng)研究：文獻(xiàn)[16-17]分析了GPU的特性，并使用其計算電力系統(tǒng)潮流和進(jìn)行諧波分析，提升了計算速度，提高了分析的實(shí)用性。

本文提出了一種應(yīng)用于 GPU計算平臺的電力系統(tǒng)暫態(tài)穩(wěn)定仿真并行算法。首先，在聯(lián)立算法框架下，依據(jù)非誠實(shí)牛頓法(Very Dishonest Newton Method, VDHN)原理，將聯(lián)立矩陣解耦，使動態(tài)元件與網(wǎng)絡(luò)方程可以分別計算。進(jìn)一步，將網(wǎng)絡(luò)矩陣重排為對角加邊（block bordered diagonal form，BBDF）形式，使此部分計算分解為子分區(qū)系統(tǒng)和邊界系統(tǒng)兩部分。子分區(qū)系統(tǒng)與相連的動態(tài)元件交由多核CPU并行處理。對于邊界系統(tǒng)，本文在GPU上使用基于稀疏近似逆預(yù)處理的 BiCGSTAB法求解。針對一個12823節(jié)點(diǎn)的大型電力系統(tǒng)進(jìn)行測試，驗證了算法的準(zhǔn)確性和有效性。結(jié)果顯示算法提高了邊界系統(tǒng)計算速度，可獲得較高加速比。

1 暫態(tài)穩(wěn)定并行算法

1.1 聯(lián)立算法基本原理

電力系統(tǒng)暫態(tài)穩(wěn)定模型的數(shù)學(xué)本質(zhì)為大規(guī)模微分代數(shù)方程初值問題

式中：x為狀態(tài)向量，描述發(fā)電機(jī)組等元件的動態(tài)行為；V為電壓向量，描述節(jié)點(diǎn)電壓變化；I為電流向量，描述節(jié)點(diǎn)注入電流變化。YN為輸電網(wǎng)絡(luò)導(dǎo)納矩陣。公式(1)描述電網(wǎng)各種動態(tài)元件的行為，如發(fā)電機(jī)組。各動態(tài)元件僅由網(wǎng)絡(luò)連接，因此，式(1)由多組相互獨(dú)立的子集組成。

采用聯(lián)立算法時，求解過程如下：

1）應(yīng)用數(shù)值積分算法將微分方程組差分化，獲得聯(lián)立代數(shù)方程組。本文采用隱式梯形積分法：當(dāng)t時步狀態(tài)向量xt，電壓向量Vt已知時，令h為積分步長，由式(1)得t+1時步的狀態(tài)向量x表達(dá)式為

由式(2)得，t+1時步網(wǎng)絡(luò)方程為

定義

2）聯(lián)立計算

式(7)為非線性代數(shù)方程組，1tx+與1tV+為未知量，可采用Newton法求解。第k次迭代計算公式為

YD為輸電網(wǎng)絡(luò)導(dǎo)納YN（固定不變）與動態(tài)元件自導(dǎo)納 YG之和。為第k次迭代時的殘差向量。，為迭代修正量。雅可比矩陣僅在迭代次數(shù) k達(dá)到設(shè)定閾值時才更新。式(8)迭代收斂后，進(jìn)入t+2時步計算。

當(dāng)計及發(fā)電機(jī)凸極效應(yīng)時，如圖1所示。

圖1 發(fā)電機(jī)網(wǎng)絡(luò)接口圖Fig. 1 Interface of the generator with the network

圖中，YGN、YGS分別為 YG中的固定部分和可變部分，即為發(fā)電機(jī)注入電流中與機(jī)端電壓無關(guān)部分。

由圖及VDHN原理[18]，式(8)的計算可簡化為

式中， Y =(YN+YGN)，其在計算過程中為固定值，僅需在網(wǎng)絡(luò)拓?fù)渥兓瘯r重新形成，十分利于算法并行化。

1.2 雙層BBDF并行算法

由電力系統(tǒng)特點(diǎn)，可對上述聯(lián)立算法進(jìn)行雙層BBDF并行分解。

第1層：電力系統(tǒng)中，發(fā)電機(jī)組通過電網(wǎng)連接，這種特點(diǎn)使式(8)的具體結(jié)構(gòu)如式(12)所示。

其中：G代表發(fā)電機(jī)組總數(shù)。g代表第g臺發(fā)電機(jī)；Ag ,Bg,Cg，，分別為A,B,C，，F(xiàn)tk+1的子矩陣。定義 B?g為 B g中非零向量，C?g為 C g中非零向量。

因此，動態(tài)元件（發(fā)電機(jī)組）相關(guān)計算，即式(9)、式(10)可分解為

第2層：由于整個電網(wǎng)為大型稀疏網(wǎng)絡(luò)，可將電網(wǎng)分為q個子分區(qū)系統(tǒng)和一個邊界系統(tǒng)，使式(11)成如下結(jié)構(gòu)[9]：

其中：b代表邊界系統(tǒng)；j代表第j子分區(qū)，j∈[1,q]；I?，Vk+1分別表示第j分區(qū)注入電流向量與電壓向

j t+1,1量。

因此，式(11)可分解為

邊界系統(tǒng)計算

子分區(qū)系統(tǒng)計算

本文中子分區(qū)系統(tǒng)式(16)的計算與相連的發(fā)電機(jī)組分配到相同CPU核心并行計算。同時需要說明的是，為便于處理式(16)，本文使用實(shí)數(shù)導(dǎo)納。

1.3 基于穩(wěn)定雙共軛梯度的邊界網(wǎng)絡(luò)方程計算

應(yīng)用1.2節(jié)并行算法后，式(16)成為唯一需要串行計算的部分，傳統(tǒng)方法使用LU分解法進(jìn)行計算。文獻(xiàn)[10，19]中引入CG算法進(jìn)行并行迭代求解, 簡化式(16)成式(18)形式。

針對不對稱實(shí)矩陣的迭代求解，BiCGSTAB法被證明是一種十分有效的改進(jìn)方案，其基于雙邊Lanczos算法和殘差正交子空間的理論[13], 具有雙正交和平方收斂的特點(diǎn)，運(yùn)算量小，收斂性優(yōu)于CG算法[21]。同時，從算法流程上，其繼承了CG法涉及大量的行(列) 的乘積與求和的特點(diǎn)，并行性很高。因此，非常適合用于邊界系統(tǒng)的并行求解。

對式(17)進(jìn)行求解時，矩陣Y的條件數(shù)也將影響算法的收斂性。而預(yù)處理技術(shù)可以降低矩陣的條件數(shù)。因此，BiCGSTAB算法一般與預(yù)處理算法相結(jié)合，算法一般流程如下。

步驟1)給定系數(shù)矩陣Y,向量I，及預(yù)處理矩陣M（Y?1的近似矩陣）。設(shè)初始值 v0，最大迭代次數(shù)imax，最大容許誤差εmax，計算r0= I ? Y v0，并令 i =1， r=r0。

步驟2)若 i ≤imax且ε≥εmax，則轉(zhuǎn)到3）；否則，得到v的值，結(jié)束迭代；

步驟3) ρi?1=rTri?1

步驟4)若 i = 1，那么 pi= r0；

此計算過程都為矩陣向量間乘加操作，可完全并行化。其中，M的構(gòu)造目前較常用的是ILU預(yù)處理方法，但并行性有限。近年來，稀疏近似逆預(yù)處理技術(shù)[22-23]以其天然的并行性而受到重視，得到迅速的發(fā)展。而且，其具有較少的自定義的參數(shù)，不確定性小。其原理為：若能構(gòu)造兩個相互關(guān)于Y共軛的向量集合 { zl}l=1,,n和 { wl}l=1,,n，組成兩個矩陣則

其中，d=wTY z≠0。那么

l l l

Z,W ,D的完整構(gòu)造過程如文獻(xiàn)[24]所述。我們?yōu)榱说玫揭粋€盡量稀疏的預(yù)條件子， Z ,W,D可被不完全地計算，即扔掉那些絕對值小于給定值的元素，即Z≈Z，W≈W，D≈D，則稀疏近似逆預(yù)條件子

2 GPU計算平臺簡介及算法實(shí)現(xiàn)

2.1 計算平臺架構(gòu)

圖2顯示了GPU計算平臺的架構(gòu)。計算平臺包含兩顆主頻為2.66 Hz的Intel Xeon E5650處理器，每顆CPU包含6個核心(core)，共享24G內(nèi)存。GPU與CPU通過PCIe 2.0總線相連。GPU型號為GeForce GTX460，其包含多個流處理器組（stream multi-processor，SM），每個SM中包含32（或48）個流處理器（stream Processor，SP）。SP是GPU上進(jìn)行數(shù)據(jù)計算的基本單元。

在NVIDIA公司開發(fā)編程環(huán)境——統(tǒng)一計算設(shè)備架構(gòu)（Compute Unified Device Architecture ，CUDA）下，如圖3所示，一個程序分為兩部分：host 端和 device 端。Host 端是指在 CPU 上執(zhí)行的部分， Host端可使用MPI等技術(shù)并行。

圖 2 GPU計算平臺硬件架構(gòu)Fig. 2 Hardware architecture of GPU plugged into a PC motherboard

圖3 CUDA計算模式Fig. 3 Hierarchy of computing structure of CUDA

device端是在GPU芯片上執(zhí)行的部分，稱為內(nèi)核（kernel）函數(shù)。在kernel中，當(dāng)一個指令被觸發(fā)時，將被大量GPU線程并行執(zhí)行，每個線程被分配給一個SP，而且線程的創(chuàng)建、撤銷和切換3種操作僅需要極少的系統(tǒng)開銷[25]。每個線程都具有唯一的編號（index），供編程時使用。為便于組織和使用這些GPU線程，CUDA中還定義了兩層結(jié)構(gòu)：網(wǎng)格（grid）和塊（block）。一個grid包含多個block，每個block管理多個線程（Thread）。

2.2算法實(shí)現(xiàn)

2.2.1 基于GPU的邊界系統(tǒng)計算

本文中邊界系統(tǒng)由1.3節(jié)所述BiCGSTAB算法計算，此方法計算過程包含大量的矩陣-矢量，矢量-矢量計算。以兩個矢量相加為例， vi=vi?1+Δ vi。在CPU上，此類計算的實(shí)現(xiàn)模式是

循環(huán)計數(shù)器index從0增加到n，向量的相應(yīng)元素順序進(jìn)行加運(yùn)算。但在此計算過程中，每個元素的計算實(shí)際是相互獨(dú)立的，沒有先后依存關(guān)系。因此在GPU上，整個向量一次性使用n個GPU線程同時操作，如下所示

其中，index 既表示向量元素的編號，同時定位了GPU線程，所有的數(shù)據(jù)是同步并行操作的。GPU可對成百乃至上萬的元素并行操作，這是CPU所不能比擬的。

2.2.2 基于CPU的分區(qū)系統(tǒng)及動態(tài)元件計算的分配根據(jù)1.2節(jié)算法，使用MPI將式(13)、式(14)、式(16)的計算分配到多個CPU核心中計算。在多核平臺上，為了獲得較高的并行效率，需根據(jù)計算量，進(jìn)行任務(wù)分配。本文依據(jù)文獻(xiàn)[7]所述方法進(jìn)行任務(wù)分配，分配指標(biāo)為

其中：q為任務(wù)數(shù)(=CPU核心數(shù)=分區(qū)總數(shù))，CompCostj代表了任務(wù) j的總計算量，必須說明的是此計算量由網(wǎng)絡(luò)子分區(qū)和相連的發(fā)電機(jī)組計算量共同確定，即式(13)、式(14)、式(17)； C ompCostB代表了邊界系統(tǒng)的計算量。由于本文中使用共享內(nèi)存的多核計算機(jī)，其通訊損耗小于PC集群，因此本文算法相當(dāng)于放寬了對邊界系統(tǒng)的約束，而更加重視子分區(qū)計算量的平均分配。

從并行計算的角度，算法一次迭代流程如圖 4所示。

3 算例分析

3.1 測試方案

本文以聯(lián)立算法為框架，在 GPU計算平臺上實(shí)現(xiàn)了并行暫態(tài)穩(wěn)定算法。

圖4 并行算法流程Fig. 4 Flow chart of an iteration of the parallel algorithm

本文算法本質(zhì)上是一種組合算法，為驗證其有效性，使用2種算法組合進(jìn)行效率比較測試，兩種算法外層框架都使用雙層BBDF并行算法，邊界系統(tǒng)分別采用基于稀疏近似逆預(yù)條件子的 CG法（簡稱為PCG）和基于稀疏近似逆預(yù)條件子的 BiCGSTAB法（簡稱為PBiCGSTAB）計算，如表1所示。

表 1 算法列表Table 1 Lists of algorithms

每時步迭代收斂條件均為 10-4。算法步長為0.01 s，CG與BiCGSTAB法收斂精度為10-5。測試方案為仿真10 s暫態(tài)過程，0.0 s系統(tǒng)中某條線路近母線處發(fā)生三相接地短路，0.083 s故障切除。

3.2 測試結(jié)果及分析

本文使用某12823節(jié)點(diǎn)系統(tǒng)，系統(tǒng)由1 431臺發(fā)電機(jī)組（發(fā)電機(jī)采用6階模型，并考慮勵磁控制及原動機(jī)調(diào)速系統(tǒng)），12365條線路，5 612臺變壓器(包含部分移相調(diào)節(jié)變壓器)組成，串行算法仿真時間serialT 為68.41 s。

首先，選取系統(tǒng)中3個節(jié)點(diǎn)電壓曲線，對本文算法（即算法2）精度進(jìn)行對比分析, 從圖5中可以看出，串行算法與并行算法電壓曲線基本相同，可以保證計算精度。

圖5 算法電壓曲線比較Fig. 5 Comparison of voltage curves between different algorithms

算法1、2并行框架同為雙層BBDF法，因此邊界系統(tǒng)規(guī)模是相同的。邊界系統(tǒng)的劃分結(jié)果及計算時間如表2所示，其中P__C___G__ 表示每時步內(nèi)PCG法___________________求解的平均迭代次數(shù)，PBiCGSTAB 表示每時步內(nèi)PBiCGSTAB法求解的平均迭代次數(shù)。NB為

__________________________PBiCGSTAB 與PCG之比， TB為 PBiCGSTAB法與PCG法對邊界系統(tǒng)計算總時間之比。

表 2 邊界系統(tǒng)計算結(jié)果對比Table 2 Comparison of the simulation results of the boundary system

由表 2可知，分區(qū)數(shù)從2增加到12的過程中，邊界系統(tǒng)由2節(jié)點(diǎn)增長到131節(jié)點(diǎn)。從迭代次數(shù)的角度，PCG法的增加較快，平均次數(shù)由2.5次增加到7.1次，而PBiCGSTAB法迭代次數(shù)僅由2.1次增加到3.1次。從實(shí)際計算時間的角度，當(dāng)邊界系統(tǒng)增大到131節(jié)點(diǎn)時，PBiCGSTAB的時間僅為PCG法的0.62。

仿真整體加速比和總時間如圖6、圖7所示。

從原理上，兩種算法的動態(tài)元件與網(wǎng)絡(luò)子分區(qū)部分的計算策略是相同的，而邊界系統(tǒng)算法不同。因此，圖中兩種算法在并行分區(qū)較少時加速比幾乎相同，但隨著分區(qū)的增多，算法1加速比接近飽和，算法2的加速比逐漸高于算法1。這是由于算法2較好的處理了分區(qū)增多導(dǎo)致的邊界系統(tǒng)增大的情況。

圖6 算法加速比對比Fig. 6 Comparison of speed-up between two algorithms

圖7 算法總計算時間對比Fig. 7 Comparison of computation time between two algorithms

在串行環(huán)境下，仿真時間大于實(shí)際暫態(tài)過程。在并行環(huán)境下，當(dāng)使用12個核時，算法2加速比達(dá)到了7.01倍。仿真時間是實(shí)際暫態(tài)過程的97.6%，小于實(shí)際暫態(tài)過程，說明了了算法2的有效性和一定的實(shí)用性。

4 結(jié)論

本文在聯(lián)立算法的基礎(chǔ)上，依據(jù)VDHN原理和雙層BBDF算法將整體暫態(tài)穩(wěn)定計算分為3部分：1）動態(tài)元件相關(guān)計算。2）子分區(qū)系統(tǒng)計算。3）邊界系統(tǒng)計算。在高性能GPU計算平臺上，首先將 1）和2）交由多核CPU負(fù)責(zé)的計算。然后，重點(diǎn)研究了邊界系統(tǒng)計算，針對其方程組左邊項為不對稱實(shí)矩陣的特點(diǎn)，引入了BiCGSTAB算法，并進(jìn)一步通過稀疏近似逆預(yù)處理技術(shù)加快收斂速度。最后，通過算例測試表明該算法性能優(yōu)于稀疏近似逆預(yù)處理

CG并行算法，有效的降低了邊界系統(tǒng)的計算時間，獲得了較高的并行加速比，對于12823節(jié)點(diǎn)算例達(dá)到了實(shí)時仿真。

[1] 石恒初. 基于PC機(jī)群的電力系統(tǒng)暫態(tài)穩(wěn)定評估[J]. 電力系統(tǒng)保護(hù)與控制, 2009, 37(10): 5-14.SHI Heng-chu. Personal computer cluster based power system transient stability assessment[J]. Power System Protection and Control, 2009, 37(10): 5-14.

[2] 汪芳宗, 張磊. 基于辛 Radau方法的暫態(tài)穩(wěn)定性并行計算方法[J]. 電力系統(tǒng)保護(hù)與控制, 2011, 39(23):35-38.WANG Fang-zong, ZHANG Lei. Parallel algorithm based on symplectic Radau method for transient stability simulation[J]. Power System Protection and Control,2011, 39(23): 35-38.

[3] 薛巍, 舒繼武, 王心豐, 等. 電力系統(tǒng)暫態(tài)穩(wěn)定仿真并行算法的研究進(jìn)展[J]. 系統(tǒng)仿真學(xué)報, 2002, 14(2):177-182.XUE Wei, SHU Ji-wu, WANG Xin-feng, et al. Advance of parallel algorithm for power system transient stability simulation[J]. Journal of System Simulation, 2002, 14(2):177-182.

[4] Chai J S, Zhu N, Bose A, et al. Parallel Newton type methods for power system stability analysis using local and shared memory multiprocessors[J]. IEEE Trans on Power Systems, 1991, 6(4): 1539-1545.

[5] Chai J S, Bose A. Bottlenecks in parallel algorithms for power system stability analysis[J]. IEEE Trans on Power Systems, 1993, 8(1): 9-15.

[6] 李亞樓, 周孝信, 吳中習(xí). 一種可用于大型電力系統(tǒng)數(shù)字仿真的復(fù)雜故障并行計算方法[J]. 中國電機(jī)工程學(xué)報, 2003, 23(12): 1-5.LI Ya-lou, ZHOU Xiao-xin, WU Zhong-xi. A parallel complex fault computation algorithm for large-scale power system digital simulation[J]. Proceedings of the CSEE, 2003, 23(12): 1-5.

[7] SHU Ji-wu, WEI Xue, ZHENG Wei-min. A parallel transient stability simulation for power systems[J]. IEEE Trans on Power Systems, 2005, 20(4): 1709-1717.

[8] de Pillis L G. A comparison of iterative methods for solving nonsymmetric linear systems[J]. Acta Applicandae Mathematicae, 1998, 51(2):141-159.

[9] Decker I C, Falcao D M, Kaszkurewicz E. Conjugate gradient methods for power system dynamic simulation on parallel computers[J]. IEEE Trans on Power Systems,1996, 11(3):1218-1227.

[10] 韓曉言, 韓禎祥. 預(yù)處理并軛梯度法在電力系統(tǒng)暫態(tài)穩(wěn)定分析并行算法中的應(yīng)用研究[J]. 電力系統(tǒng)及其自動化學(xué)報, 1996, 8(2): 1-6.HAN Xiao-yan, HAN Zhen-xiang. Research on the application of parallel algorithm of preconditioning conjugated gradient of power system transient stability analysis[J]. Proceedings of the CSU-EPSA, 1996, 8(2):1-6.

[11] 林濟(jì)鏗, 仝新宇, 李楊春, 等. 基于預(yù)處理共軛梯度法的電力系統(tǒng)機(jī)電暫態(tài)仿真[J]. 電工技術(shù)學(xué)報, 2008,23(5): 93-99.LIN Ji-keng, TONG Xin-yu, LI Yang-chun, et al.Electrical-machanical transient simulation of power system based on preconditioned conjugate gradient method[J]. Transactions of China Electrotechnical Society, 2008, 23(5): 93-99.

[12] 劉洋, 周家啟, 謝開貴, 等. 基于 Beowulf集群的大規(guī)模電力系統(tǒng)方程并行 PCG求解[J]. 電工技術(shù)學(xué)報,2006, 21(3): 105-111.LIU Yang, ZHOU Jia-qi, XIE Kai-gui, et al. Parallel PCG solution of large scale power system equations based on Beowulf cluster[J]. Transactions of China Electrotechical Society, 2006, 21(3): 105-111.

[13] van der Vorst H A. BI-CGSTAB: a fast and smoothly converging variant of BI-CG for the solution of nonsymmetric linear systems[J]. SIAM Journal on Scientific and Statistical Computing, 1992, 13(2):631-644.

[14] 張浩, 李利軍, 林嵐. GPU的通用計算應(yīng)用研究[J]. 計算機(jī)與數(shù)字工程, 2005(12): 60-62.ZHANG Hao, LI Li-jun, LIN Lan. General-purpose computation on graphics processors[J]. Computer ＆amp;Digital Engineering, 2005(12): 60-62.

[15] Jalili-Marandi V, Dinavahi V. SIMD-based large-scale transient stability simulation on the graphics processing unit[J]. IEEE Trans on Power Systems, 2010, 25(3):1589-1599.

[16] 夏俊峰, 楊帆, 李靜, 等. 基于 GPU 的電力系統(tǒng)并行潮流計算的實(shí)現(xiàn)[J]. 電力系統(tǒng)保護(hù)與控制, 2010,38(18): 100-103.XIA Jun-feng, YANG Fan, LI Jing, et al. Implementation of parallel power flow calculation based on GPU[J].Power System Protection and Control, 2010, 38(18):100-103.

[17] 韓志偉, 劉志剛, 魯曉帆, 等. 基于CUDA的高速并行小波算法及其在電力系統(tǒng)諧波分析中的應(yīng)用[J]. 電力自動化設(shè)備, 2010, 38(1): 98-101.HAN Zhi-wei, LIU Zhi-gang, LU Xiao-fan, et al.High-speed parallel wavelet algorithm based on CUDA and its application in power system harmonic analysis[J].Electric Power Automation Equipment, 2010, 38(1):98-101.

[18] 江涵, 江全元. 一種可變步長的暫態(tài)穩(wěn)定自適應(yīng)修正牛頓組合算法[J]. 中國電機(jī)工程學(xué)報, 2011, 31(34):105-112.JIANG Han, JIANG Quan-yuan. Power system transient stability simulation based on Shamanskii-VDHN algorithm with variable step size[J]. Proceedings of the CSEE, 2011, 31(34): 105-112.

[19] Decker I C, Falcao D M, Kaszkurewicz E. Parallel implementation of a power system dynamic simulation methodology using the conjugate gradient method[J].IEEE Trans on Power Systems, 1992, 7(1): 458-465.

[20] 昆德. 電力系統(tǒng)穩(wěn)定與控制[M]. 北京: 中國電力出版社, 2001.Kundur. Power system stability and control[M]. Beijing:China Electric Power Press, 2001.

[21] Topsakal E, Kindt R, Sertel K, et al. Evaluation of the BICGSTAB(l) algorithm for the finite-element/boundaryintegral method[J]. Antennas and Propagation Magazine,IEEE, 2001, 43(6):124-131.

[22] 谷同祥, 遲學(xué)斌, 劉興平. 稀疏近似逆與多層塊ILU預(yù)條件技術(shù)[J]. 應(yīng)用數(shù)學(xué)和力學(xué), 2004, 25(9): 927-934.GU Tong-xiang, CHI Xue-bin, LIU Xing-ping. AINV and BILUM preconditioning techniques[J]. Applied Mathematics and Mechanics, 2004, 25(9): 927-934.

[23] 汪芳宗, 何一帆, 葉婧. 基于稀疏近似逆預(yù)處理的牛頓-廣義極小殘余潮流計算方法[J]. 電網(wǎng)技術(shù), 2008,32(14): 50-53.WANG Fang-zong, HE Yi-fan, YE Jing. Load flow calculation of Newton-GMRES method with sparse approximate inverse preconditioners[J]. Power System Technology, 2008, 32(14): 50-53.

[24] Benzi Michele, Tuma Miroslav. A sparse approximate inverse preconditioner for nonsymmetric linear systems[J]. SIAM Journal on Scientific Computing, 1998,19(3): 968-994.

[25] Owens J D, Houston M, Luebke D, et al. GPU Computing[J]. Proceedings of the IEEE, 2008, 96(5):879-899.