大規(guī)模電力系統(tǒng)潮流計算的分布式GESP算法

謝開貴 張懷勛 胡 博 曹 侃 吳 韜

（重慶大學輸配電裝備及系統(tǒng)安全與新技術國家重點實驗室 重慶 400044）

1 引言

區(qū)域大電網(wǎng)互聯(lián)的運行行為分析以及電力市場環(huán)境下實時安全控制和潮流跟蹤計算，都迫切需要快速地進行大規(guī)模甚至超大規(guī)模電力潮流計算[1]。在潮流計算中，需對形如 Ax=b的高維稀疏線性修正方程進行反復求解，以得到收斂的潮流解。文獻[2]指出：在用牛頓法計算電力系統(tǒng)潮流時，超過80%的時間用于求解高維修正方程組，而 LU分解更占據(jù)了大規(guī)模電力系統(tǒng)潮流單步迭代的絕大部分時間。在求解超大規(guī)模電力系統(tǒng)潮流問題時，已有的潮流算法呈現(xiàn)出計算機內(nèi)存不足、收斂速度慢等維數(shù)災問題。

高性價比可擴展集群并行系統(tǒng)的逐步成熟和應用，使大規(guī)模電力系統(tǒng)潮流并行計算和分布式仿真成為可能。采用并行技術處理 LU分解和三角方程的求解，已成為大規(guī)模電力系統(tǒng)潮流計算的新思路。并行計算的本質(zhì)是將多任務映射到多處理機中進行處理。電力系統(tǒng)潮流并行計算的主要研究內(nèi)容正是設計具有較高并行性、較低數(shù)據(jù)相關性的高效并行算法。

現(xiàn)有的潮流并行計算方法主要有3類：分解協(xié)調(diào)計算法、稀疏矢量法、Krylov子空間迭代法。分解協(xié)調(diào)計算法[3-4]通過對大規(guī)模電力網(wǎng)絡進行物理意義上的區(qū)域分塊，并將各子塊通過協(xié)調(diào)層聯(lián)系起來。這樣便將潮流計算轉(zhuǎn)換為各個子塊潮流的并行計算和協(xié)調(diào)層的計算。但是，對大電網(wǎng)進行合理區(qū)域劃分的策略和依據(jù)，以及對協(xié)調(diào)層進行平衡計算已成為該算法的難點。稀疏矢量法[5]需提前判斷前代和回代過程中所必需的運算，然后根據(jù)運算過程中的因子表路徑找出其中可以并行執(zhí)行的部分。文獻[6]將消去樹理論與稀疏向量法相結(jié)合使得LU分解過程也具有了可并行性。Krylov子空間方法[7-8]在潮流計算過程中將高維非線性方程組劃分為內(nèi)、外雙層迭代，在求解內(nèi)層線性化方程組時采用Krylov子空間方法并行求解。但是，文獻[9]指出Krylov子空間方法的收斂速度和穩(wěn)定性嚴重依賴于系數(shù)矩陣的條件數(shù)和特征值分布，所以必須對其進行恰當?shù)念A條件處理才能取得良好的收斂特性。

GESP（Gaussian Elimination with Static Pivoting）算法[10]已在 20世紀末提出，其適用于分布式并行處理大規(guī)模稀疏方程組。該算法通過采用靜態(tài)選主元技術來代替?zhèn)鹘y(tǒng)稀疏高斯消去過程中的部分選主元技術。文獻[11]指出，GESP算法已在電路仿真、石油工程、流體力學等領域取得了初步的應用，并且在分布式存儲的并行計算平臺上顯示出良好的穩(wěn)定性。但是到目前為止，GESP算法在電力系統(tǒng)潮流計算中的設計、實現(xiàn)及應用還很鮮見。

已有的電力系統(tǒng)潮流并行算法大都基于共享存儲的并行計算平臺，使得潮流并行計算嚴重依賴于昂貴的共享存儲并行計算機，以至于在實際工程應用中潮流并行計算很難實現(xiàn)。本文在分布式集群環(huán)境下，結(jié)合電力系統(tǒng)潮流方程固有的稀疏性特點，利用GESP算法對潮流計算中修正的大型稀疏方程組進行并行LU分解以及 LU分解后的三角方程進行并行求解。

2 GESP算法及其在潮流計算中的應用

傳統(tǒng)牛頓潮流算法將大量非線性方程組轉(zhuǎn)換成線性修正方程組進行反復迭代、串行求解。但是，隨著網(wǎng)絡規(guī)模不斷擴大，串行計算面臨著越來越嚴重的維數(shù)災問題。利用潮流問題修正方程系數(shù)矩陣固有的超稀疏性，對其進行快速并行求解成為解決問題的新思路。

2.1 GESP算法

稀疏線性方程組的直接求解常需進行系數(shù)矩陣的LU分解。進行矩陣分解時，可能會引入填充元，從而增加了計算復雜性和儲存需求。另外，分解過程中的主元為0或者主元的模太小也會導致分解無法進行或者數(shù)值不穩(wěn)定。

傳統(tǒng)的部分選主元 LU分解算法在計算過程中通過對主元進行修正，在保證矩陣的非零元結(jié)構(gòu)情況下取得較高的數(shù)值穩(wěn)定性。但是，部分選主元算法必須動態(tài)改變數(shù)據(jù)的存儲結(jié)構(gòu)，由此所采用的細粒度通信模式使其難以在分布式存儲的集群系統(tǒng)上并行執(zhí)行。

靜態(tài)選主元相對于部分選主元的優(yōu)點在于可以提前決定高斯消去過程中的數(shù)據(jù)結(jié)構(gòu)和通信模式?；陟o態(tài)選主元技術的GESP算法在LU數(shù)值分解之前對矩陣A進行對稱置換和選主元操作，減少了分解過程中非零元的填入，提高了分解過程的并行性。這使得GESP算法可以提前確定填入模式、分布式的數(shù)據(jù)結(jié)構(gòu)和通信模式，便于對A矩陣進行塊狀劃分的粗粒度的并行計算。

設 A=[aij]為線性方程組 Ax=b的系數(shù)矩陣，ε為收斂容差（可根據(jù)實際問題需要選取合適的值），eb、elb分別表示當前回代誤差和最小回代誤差。GESP算法[11]求解大規(guī)模稀疏方程組的計算步驟如下：

（1）通過雙向加權匹配算法對A執(zhí)行行置換，以實現(xiàn)其對角占優(yōu)。

（2）對 A進行符號分解并執(zhí)行列置換，以保證其在分解過程中的稀疏性，同時對A進行靜態(tài)選主元。

（3）因式分解·=AL U并控制對角元絕對值的大?。?/p>

（4）用矩陣L和U進行三角求解。如有必要，進行迭代修正：

從GESP算法流程可以看出，第2步的靜態(tài)選主元操作可能使所得解中包含較大的誤差，但通過（4）的迭代修正保證了最終解的精度。GESP算法作出這種改進的目的正是為了適應大規(guī)模線性方程組分布式存儲、求解的需要。

2.2 基于GESP算法的牛頓潮流算法

牛頓潮流算法的突出優(yōu)點是算法具有二次方收斂特性，且算法的迭代次數(shù)與網(wǎng)絡規(guī)?；緹o關，對某些病態(tài)問題也具有較可靠的收斂特性[12]。但牛頓法每次迭代計算時間過長以及所需的存儲量較大，所以很難直接應用于超大規(guī)模電力系統(tǒng)潮流求解。由于牛頓潮流修正方程的系數(shù)矩陣是一個超稀疏矩陣，且系數(shù)矩陣的非零元主要分布于對角帶上。如果能夠充分利用系數(shù)矩陣的稀疏性，進行各塊都較小的某種分塊劃分時，只有少量塊中存在非零元素，而其他塊中的元素都是 0。利用 GESP算法對牛頓法潮流的修正方程組進行求解，可以充分利用修正方程系數(shù)矩陣的稀疏性，降低計算過程中的存儲需求。此外，GESP算法所具有的可并行性，使得對修正方程進行直接并行求解成為可能，從而大大降低了程序所需的計算時間。根據(jù)牛頓法潮流計算的步驟，可形成基于GESP算法的牛頓潮流分布式算法的流程，如圖1所示。

3 基于MPI的分布式潮流計算

下面將根據(jù)潮流修正方程系數(shù)矩陣的結(jié)構(gòu)和分布的特點，設計GESP算法的數(shù)據(jù)結(jié)構(gòu)以及稀疏LU分解的分布式處理算法，算法流程如圖1所示。

本文利用MPI（Message Passing Interface）進行數(shù)據(jù)通信和進程同步。MPI是目前國內(nèi)外高性能集群中最廣泛使用的并行通信接口，采用消息傳遞的方式實現(xiàn)并行程序間通信[13]。在分布式計算過程中數(shù)據(jù)復制的速度和延遲是影響消息傳遞效率的關鍵因素，因而如何減少進程之間的數(shù)據(jù)交換是提高并行效率的關鍵。

3.1 分布式數(shù)據(jù)結(jié)構(gòu)的設計

如 2.2節(jié)所述：大規(guī)模電力系統(tǒng)牛頓潮流修正方程的系數(shù)矩陣為一超稀疏矩陣且非零元主要集中于對角帶。由此可見，對系數(shù)矩陣的分塊存儲可以在很大程度上減少存儲量與計算量。

本文將系數(shù)矩陣基于超節(jié)點（超節(jié)點即是基于L或U陣的非零元結(jié)構(gòu)劃分的方陣）[12]的邊界劃分為若干個二維分塊矩陣。由于系數(shù)矩陣的非零元主要集中于對角帶上，所以每個超節(jié)點中只存在少量的非零元。如果一個N×N矩陣里面有M（M?N）個超節(jié)點，那么該矩陣即被劃分為M2個子塊。

如圖2所示：圖2a為一個由6個處理器節(jié)點（即0，1，2，3，4，5）排列成 2×3的 2維處理器網(wǎng)格；圖 2b表示系數(shù)矩陣分塊劃分后，分塊存儲在由圖2a中的處理器網(wǎng)格循環(huán)排列成的處理器網(wǎng)格上，即方格中的數(shù)字表示該子塊對應的處理器編號。圖2b中位于系數(shù)矩陣對角線上的陰影子塊皆為超節(jié)點，而其他子塊的維數(shù)則取決于其所處行塊與列塊中超節(jié)點的維數(shù)?；谥苯訉Ψ謮K進行存儲的思想，系數(shù)矩陣的每個子塊被看作一個元素存儲在2維處理器網(wǎng)格的某個節(jié)點上。設Pr、Pc分別表示2維處理器網(wǎng)格中行、列所含處理節(jié)點的個數(shù)，則在這種存儲方式中，子塊A(i, j)（0≤i, j＜M-1）將被存儲在坐標為（i mod Pr, j mod Pc）的處理器網(wǎng)格節(jié)點上。如圖2中，設圖2b中A(3，5)表示A陣中第4行、第6列的子塊，則其被存儲在圖2a中坐標為（1，2）的5號處理器上。

圖2 對應于處理器網(wǎng)格的矩陣分塊存儲Fig.2 Distributed storage mode for the matrix based on the processor grid

對于稀疏矩陣，盡管1維劃分更常用且更易于執(zhí)行，但是采用上述2維分塊劃分存儲方式后，系數(shù)矩陣中下三角子塊 L(i, j)在分解過程中僅被存儲了行塊i的那些處理器用到。同樣，上三角子塊U(i,j)也僅被存儲了列塊 j的那些處理器用到。由此可見，2維分塊存儲可以帶來更好的負載平衡和更低的通信量，其也更利于算法的升級。

3.2 基于流水線技術的并行LU分解

在矩陣分塊存儲中，系數(shù)矩陣的每個列塊、行塊都存儲在多個處理器節(jié)點上。記 PR(i)、PC(j)分別為存儲了行塊i、列塊j的處理器集合。如圖2中，PR(2)={3,4,5}，PC(1)={0，3}。

并行稀疏 LU分解算法在各個超節(jié)點之間進行迭代分解。其第K步迭代主要由3步組成：

（1）處理節(jié)點列 PC(K)參與列塊 L(K∶M, K)（即列塊K中位于超節(jié)點下方的所有子塊）的因子分解。

（2）處理節(jié)點行 PR(K)參與分解行塊 U(K,K+1∶K)。

（3）對所有節(jié)點用 L(K+1∶M, K)和 U(K,K+1∶M)進行更新。

上述過程中第3步為并行LU分解的主要工作，同時相對于其他兩步也具有更大的可并行性。設dia為當前迭代過程所處理的超節(jié)點，并行稀疏 LU分解算法的基本流程如下：

并行稀疏LU分解單步迭代時的第3步工作耗時最長，可以與下一個單步迭代的第 1、第 2步工作并行執(zhí)行。本文采用流水線技術來實現(xiàn)其并行計算，即在處理器列PC(K)完成第K步迭代列塊K+1的分解后，處理器列 PC(K+1)可以同時并行執(zhí)行第K+1步迭代的前兩步計算，而不必等到處理器列PC(K)完成第 K步迭代的所有列塊分解后才開始第K+1步迭代的計算。其過程如圖 3所示，圖中 Tij表示第 i步迭代的第 j步計算工作。由圖可見，基于流水線技術的并行算法相對于串行算法，可以使相鄰的兩步迭代在計算時間上得到重疊，從而達到了減少程序整體計算時間的目的。

圖3 基于流水線技術的并行處理過程Fig.3 Parallel processing based on the pipeline technique

在并行計算過程中，本文使用MPI標準非阻塞通信接口mpi_irecv和mpi_isend完成節(jié)點間的數(shù)據(jù)接收和發(fā)送操作。非阻塞通信可以解決節(jié)點間數(shù)據(jù)同步時等待時間過長的問題，實現(xiàn)計算與通信的時間重疊。因為在此通信模式下，計算節(jié)點可以在等待接收所需數(shù)據(jù)的同時向其他節(jié)點發(fā)送數(shù)據(jù)。而在傳統(tǒng)阻塞通信模式下，節(jié)點在一個通信時間段內(nèi)只能執(zhí)行發(fā)送操作或者接收操作，從而導致節(jié)點間耗費大量的相互等待時間，極端情況下甚至導致進程間發(fā)生死鎖，使程序的計算效率大大降低?；贛PI非阻塞通信模式的流水線技術大大減少了進程間的等待、數(shù)據(jù)同步時間，提高了程序的并行效率。

4 算例測試及分析

本文利用曙光 TC1700服務器和 Cisco千兆以太網(wǎng)交換機搭建了具有若干處理節(jié)點的分布式計算平臺。每個節(jié)點處理器為Xeon 2.0GHz/512KB，具有1GB內(nèi)存存儲空間，操作系統(tǒng)為Linux v9.0，節(jié)點間的并行通信接口為：MPICH-1.2.5。

由于對系數(shù)矩陣采用了基于超節(jié)點的分塊存儲，如果不同的超節(jié)點之間規(guī)模差異過大將不利于節(jié)點間負載平衡，同時單個超節(jié)點的規(guī)模過大會使其失去稀疏性從而增加計算量、降低了單個處理器節(jié)點的計算效率，所以本文選定單個超節(jié)點的最高維數(shù)不超過20。潮流計算中，收斂容差為10-4（基準功率100MVA）。

測試系統(tǒng)包括由 IEEE標準測試系統(tǒng)（30和118節(jié)點系統(tǒng)）合成的多個大規(guī)模系統(tǒng)。合成系統(tǒng)[15]將IEEE標準系統(tǒng)按N×N網(wǎng)格結(jié)構(gòu)通過聯(lián)絡線互聯(lián)起來，從而得到更大規(guī)模的測試系統(tǒng)。合成系統(tǒng)的節(jié)點導納矩陣與同等規(guī)模的實際輸電網(wǎng)絡的節(jié)點導納矩陣具有接近的條件數(shù)。表1給出了8個算例系統(tǒng)的網(wǎng)絡規(guī)模和Jacobi矩陣條件數(shù)。

表1 測試系統(tǒng)的規(guī)模及Jacobi矩陣條件數(shù)Tab.1 The scale of test systems and condition number of Jacobi matrix

表2為本文分布式GESP法與分布式牛頓法、串行牛頓法潮流計算結(jié)果的對比，其中分布式GESP法與分布式牛頓法為4節(jié)點4進程的計算結(jié)果，串行牛頓法為單節(jié)點的串行求解結(jié)果。收斂時間指10次潮流計算的平均時間。從表2可以看出：分布式GESP潮流算法與分布式牛頓法的收斂特性相同，且收斂速度明顯快于分布式牛頓法。這表明GESP算法已對直接牛頓法做出了重要改進，使其利于分布式并行計算。

表2 三種潮流計算方法的速度比較Tab.2 Comparison of calculation speeds for three power flow calculation methods

由表 2可見，隨著系統(tǒng)規(guī)模的增大，分布式GESP求解潮流的收斂時間呈超線性增長。此時，增加計算節(jié)點成為解決問題的思路之一。表3給出3000和12000節(jié)點系統(tǒng)在采用不同數(shù)目計算節(jié)點時的加速比與并行效率。由表3可看出：隨著計算節(jié)點的增加，并行加速比的增大會出現(xiàn)瓶頸并呈逐步下降趨勢，并行效率亦將趨于降低。其原因為：一方面，隨著并行度的增大，并行計算部分的計算負荷將減小而串行計算負荷將增大，串行計算所引起的開銷在整個計算過程中所占的比例越來越高；另一方面，GESP算法涉及全局通信操作，其耗用的時間將隨著參與通信操作的并行處理單元的數(shù)目以及通信量的增加而增加。因此，對大規(guī)模電力系統(tǒng)方程進行分布式計算時，隨著計算節(jié)點的增多，并行加速比在達到某一值后，便會隨著并行度的繼續(xù)增大而下降，并行效率亦趨于下降。這符合Amdahl并行加速比定律[14]。

表3 逐步增加計算節(jié)點時的加速比與并行效率Tab.3 Comparison of parallel speedup and efficiency with different number of processors

表4給出了節(jié)點數(shù)為3000的系統(tǒng)在采用不同并行算法時加速比和并行效率的變化情況。從表中可以看出：在相同計算節(jié)點的情況下，GESP算法比分布式牛頓法具有更高的加速比和并行效率。這說明本文采用的稀疏矩陣分塊存儲技術可以有效降低并行計算過程中進程間的通信量，而基于流水線技術的 LU并行分解則進一步減少了進程間的等待和同步時間，以上兩者使得本文方法的加速比和并行效率相對傳統(tǒng)并行算法得到了明顯的提升。

表4 采用不同并行算法計算SYN3000網(wǎng)絡的加速比與并行效率Tab.4 Comparison of parallel speedup and efficiency using different parallel methods for the SYN3000 system

5 結(jié)論

基于新型數(shù)學模型的潮流快速并行計算越來越得到研究人員的重視，為提高潮流并行計算的加速比和并行效率，本文利用GESP算法對牛頓潮流修正線性方程組的分布式并行求解進行了研究。

根據(jù)牛頓潮流方程修正線性方程組系數(shù)矩陣的稀疏性特點，本文對其進行基于超節(jié)點劃分的分塊存儲，并在此基礎上將流水線技術應用于并行 LU分解過程中。在分布式計算平臺上，采用 Linux+MPI的軟件架構(gòu)，實現(xiàn)了GESP牛頓潮流算法的并行計算。通過對不同規(guī)模的測試系統(tǒng)計算表明：本文方法對較大規(guī)模電力網(wǎng)絡進行潮流計算時，在滿足相同求解精度的條件下，相對于串行潮流算法、傳統(tǒng)并行算法具有較好的速度優(yōu)勢。通過對同一系統(tǒng)在不同計算節(jié)點下測試表明：本文方法相對傳統(tǒng)潮流分布式算法具有可擴展性，且有較高的加速比與并行效率。

[1]盧強，周孝信，薛禹勝，等. 面向21世紀電力科學技術講座[M]. 北京：中國電力出版社，2000.

[2]劉洋，周家啟，謝開貴，等. 基于 Beowulf集群的大規(guī)模電力系統(tǒng)方程并行 PCG求解[J]. 電工技術學報，2006，21 (3)： 105-111.Liu Yang, Zhou Jiaqi, Xie Kaigui, et al. Parallel PCG solution of large scale power system equations based on Beowulf cluster[J]. Transactions of China Electrotechnical Society, 2006, 21 (3)： 105-111.

[3]張海波，張伯明，孫宏斌. 基于異步迭代的多區(qū)域互聯(lián)系統(tǒng)動態(tài)潮流分解協(xié)調(diào)計算[J]. 電力系統(tǒng)自動化，2003，27(24)： 1-9.Zhang Haibo, Zhang Boming, Sun Hongbin. A decomposition and coordination dynamic power flow calculation for multi-area interconnected system based on asynchronous iteration[J]. Automation of Electric Power Systems，2003，27(24)： 1-9.

[4]陳穎，沈沉，梅生偉，等. 基于改進 Jacobian-Free Newton-GMRES(m)的電力系統(tǒng)分布式潮流計算[J].電力系統(tǒng)自動化，2006，30(9)： 5-9.Chen Ying, Shen Chen, Mei Shengwei, et al.Distributed power flow calculation based on an improved Jacobian-Free Newton-GMRES (m)method[J].Automation of Electric Power Systems, 2006, 30(9)： 5-9.

[5]Wu Junqiang, Bose. A parallel solution of large sparse matrix equations and parallel power flow[J].IEEE Transactions on Power System, 1995, 10(3)：1343-1349.

[6]徐得超，李亞樓，郭劍，等. 消去樹理論及其在潮流計算中的應用[J]. 電網(wǎng)技術，2007，31(22)： 12-16.Xu Dechao, Li Yalou, Guo Jian, et al. Elimination tree theory and its application in power flow calculation[J].Power System Technology, 2007, 31(22)： 12-16.

[7]Chen Tsung Hao, Charlie Chung Ping Chen. Efficient large-scale power grid analysis based on preconditioned Krylov-subspace iterative methods[C]. Proceedings of Design Automation Conference, USA, Las Vegas,Nevada, 2001： 559-562.

[8]Alexander Flueck J, Chiang Hsiao Dong. Solving the nonlinear power flow equations with an inexact Newton method using GMRES[J]. IEEE Transactions on Power System, 1998, 13(2)： 267-273.

[9]劉洋，周家啟，謝開貴，等. 預條件處理 CG法大規(guī)模電力系統(tǒng)潮流計算[J]. 中國電機工程學報，2006，26(7)： 89-94.Liu Yang, Zhou Jiaqi, Xie Kaigui, et al. The preconditioned CG method for large scale power flow solution[J]. Proceedings of the CSEE, 2006, 26(7)： 89-94.

[10]Li X S, Demmel J W. Making sparse Gaussian elimination scalable by static pivoting[C]. In Proceedings of SC98： High Performance Networking and Computing Conference (Orlando, FL), 1998： 236-240.

[11]Patrick R Amestoy, Iain S Duff, Xiaoye S Li.Analysis and comparison of two general sparse solvers for distributed memory computers[J]. ACM Transactions on Mathematical Software (TOMS),2001, 27(4)： 388-421.

[12]諸駿偉. 電力系統(tǒng)分析[M]. 北京： 中國電力出版社，1995.

[13]都志輝. 高性能并行計算之并行編程技術——MPI并行程序設計[M]. 北京： 清華大學出版社，2002.

[14]James W Demmel, Stanley C Eisenstat, John R Gilbert, et al. A supernode approach to sparse partial pivoting[J]. SIAM Journal on Matrix Analysis and Applications, 1999, 20(3)： 1-15.

[15]劉洋. 大規(guī)模電力系統(tǒng)并行處理技術及可靠性評估Web計算系統(tǒng)研究[D]. 重慶： 重慶大學, 2006.

[16]陳國良. 并行計算—結(jié)構(gòu)、算法、編程[M]. 北京：高等教育出版社，1999.