基于GPU計(jì)算平臺(tái)的大規(guī)模電力系統(tǒng)暫態(tài)穩(wěn)定計(jì)算

江 涵，江全元

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

0 引言

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

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

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

基本的并行計(jì)算硬件平臺(tái)主要基于PC集群和共享內(nèi)存的多核計(jì)算機(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)計(jì)算，在包含GPU的多核計(jì)算平臺(tái)上（本文簡(jiǎn)稱為GPU計(jì)算平臺(tái)）使用直接法計(jì)算電力系統(tǒng)暫態(tài)穩(wěn)定性。我國學(xué)者同時(shí)也展開了相應(yīng)研究：文獻(xiàn)[16-17]分析了GPU的特性，并使用其計(jì)算電力系統(tǒng)潮流和進(jìn)行諧波分析，提升了計(jì)算速度，提高了分析的實(shí)用性。

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

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

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

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

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

定義

2）聯(lián)立計(jì)算

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

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

當(dāng)計(jì)及發(fā)電機(jī)凸極效應(yīng)時(shí)，如圖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)的計(jì)算可簡(jiǎn)化為

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

1.2 雙層BBDF并行算法

由電力系統(tǒng)特點(diǎn)，可對(duì)上述聯(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臺(tái)發(fā)電機(jī)；Ag ,Bg,Cg，，分別為A,B,C，，F(xiàn)tk+1的子矩陣。定義 B?g為 B g中非零向量，C?g為 C g中非零向量。

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

第2層：由于整個(gè)電網(wǎng)為大型稀疏網(wǎng)絡(luò)，可將電網(wǎng)分為q個(gè)子分區(qū)系統(tǒng)和一個(gè)邊界系統(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)計(jì)算

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

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

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

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

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

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

步驟1)給定系數(shù)矩陣Y,向量I，及預(yù)處理矩陣M（Y?1的近似矩陣）。設(shè)初始值 v0，最大迭代次數(shù)imax，最大容許誤差εmax，計(jì)算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；

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

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

l l l

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

2 GPU計(jì)算平臺(tái)簡(jiǎn)介及算法實(shí)現(xiàn)

2.1 計(jì)算平臺(tái)架構(gòu)

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

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

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

圖3 CUDA計(jì)算模式Fig. 3 Hierarchy of computing structure of CUDA

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

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

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

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

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

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

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

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

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

3 算例分析

3.1 測(cè)試方案

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

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

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

表 1 算法列表Table 1 Lists of algorithms

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

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

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

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

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

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

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

表 2 邊界系統(tǒng)計(jì)算結(jié)果對(duì)比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í)際計(jì)算時(shí)間的角度，當(dāng)邊界系統(tǒng)增大到131節(jié)點(diǎn)時(shí)，PBiCGSTAB的時(shí)間僅為PCG法的0.62。

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

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

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

圖7 算法總計(jì)算時(shí)間對(duì)比Fig. 7 Comparison of computation time between two algorithms

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

4 結(jié)論

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

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

[1] 石恒初. 基于PC機(jī)群的電力系統(tǒng)暫態(tài)穩(wěn)定評(píng)估[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ì)算方法[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é)報(bào), 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]. 中國電機(jī)工程學(xué)報(bào), 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)及其自動(dòng)化學(xué)報(bào), 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é)報(bào), 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é)報(bào),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的通用計(jì)算應(yīng)用研究[J]. 計(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)并行潮流計(jì)算的實(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]. 電力自動(dòng)化設(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é)報(bào), 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ì)算方法[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.