電力系統(tǒng)潮流并行計算中的方程組求解方法?

楊淑丹 董方敏

（三峽大學計算機與信息學院 宜昌 443002）

1 引言

電力系統(tǒng)潮流計算在電力系統(tǒng)規(guī)劃預(yù)測、生產(chǎn)運行、調(diào)度管理、科學計算中有著廣泛的應(yīng)用。隨著電力網(wǎng)絡(luò)大規(guī)模、強耦合的特點日益突出，傳統(tǒng)的串行潮流計算已經(jīng)無法滿足大電網(wǎng)仿真和實時模擬的要求，高效的并行潮流算法正成為研究的熱點。從數(shù)學角度來看，潮流計算是求解一組由有功、無功功率、電壓和功角構(gòu)成的非線性方程組，求解時間占達整個潮流計算時間的70%以上，因此對潮流并行計算中的方程組求解方法的研究是實現(xiàn)實時仿真計算的關(guān)鍵。

本文針對電力系統(tǒng)潮流并行計算中的方程組求解這一核心問題，分析了潮流并行計算大型矩陣常用算法，通過對稀疏矩陣壓縮為稠密矩陣的算法以及矩陣分塊算法進行性能比較，給出了不同算法的適用范圍；同時，對常用的潮流計算線性方程組求解并行算法進行探討，對比分析了共軛梯度法、廣義極小殘差法等不同算法的優(yōu)點和局限性；并對進一步的研究方向進行了展望。

2 電力系統(tǒng)潮流計算模型

電網(wǎng)穩(wěn)定狀態(tài)特性通常采用電力系統(tǒng)潮流計算進行描述。設(shè)網(wǎng)絡(luò)節(jié)點數(shù)為n，PQ節(jié)點數(shù)為m-1，PV節(jié)點數(shù)為n-m，則可以將潮流方程表示成為以下直角坐標形式：

式（1）是第i個電力系統(tǒng)有功功率方程式，共有n-1個，式（2）是第i個無功功率方程式，共有n-m-1個。其中，Pi為PQ和PV節(jié)點注入的有功功率，Qi為PV節(jié)點注入的無功功率，ei和 fi分別為迭代過程中節(jié)點電壓的實部、虛部，ej和 fj分別為平衡點電壓的實部、虛部，Gij為節(jié)點的自導(dǎo)，Bij為節(jié)點間的互導(dǎo)。

求解潮流計算的方法有很多，最常用的算法為牛頓拉夫遜法。其修正方程式如下

式中，ΔPi、ΔQi為節(jié)點的注入有功功率和無功功率，Δei、Δfi為節(jié)點電壓不平衡量的實部和虛部，i=1，…，n。修正方程也可簡寫為

其中，Δf為不平衡量的列向量，J為雅可比矩陣，Δx為節(jié)點電壓列向量。根據(jù)上述公式，可以將電力系統(tǒng)潮流計算轉(zhuǎn)化為對線性方程組的求解問題，即

在大規(guī)模電網(wǎng)中，由于節(jié)點之間的連接線遠少于電力系統(tǒng)，大型系數(shù)矩陣 A一般為高度稀疏陣。因此，通常采用不同的大型稀疏陣壓縮算法來節(jié)省存儲空間，以及通過矩陣分塊方法提高潮流并行計算效率；同時，可對線性方程組求解算法進行優(yōu)化，針對不同情況選用合適的并行算法進行求解，提高運行效率。本文將從大型稀疏矩陣處理算法以及方程組并行求解算法兩個方面，來探討分析電力系統(tǒng)潮流并行計算中方程組求解的問題。

3 潮流并行計算矩陣處理算法

3.1 稀疏矩陣壓縮為稠密矩陣算法

為了減少存儲空間，便于方程組并行求解，通常采用稀疏陣壓縮算法改進方程組的求解，下面將依次介紹COO，ELL，HYB三種壓縮方法。

1）COO（Coordinate form）存儲

COO存儲［1～2］是采用三元組對稀疏矩陣的非零元進行存儲，其中row，column，value分別代表元素所在的行、列和元素的值，三元組分別用三個數(shù)組進行存儲。COO存儲不存在數(shù)據(jù)填充的問題，壓縮比較徹底，但是缺乏行、列的規(guī)整性，因此并行時的性能一般。

2）ELL（ELL-pack form）存儲

ELL存儲［3］是使用一個n×k的矩陣進行存儲，k指非零元素最多行的非零元素數(shù)目，ELL存儲不適用于行之間非零元相差很大和非零元很少的稀疏矩陣。

3）HYB（Hybrid form）存儲

HYB存儲［4］是上述兩種存儲方式的綜合，可實現(xiàn)GPU+CPU同時進行運算，來減少運算時間。其實現(xiàn)算法可描述如下：

輸入：稀疏矩陣A，稀疏矩陣行數(shù)N，列數(shù)M，行非零元素的下限值min；

若每行中非零元數(shù)目＜min

則非零元存儲在ELL格式中，在GPU中運算；

否則

將其存儲在COO格式中，在CPU中運算。

但該算法存在任務(wù)分配不均、GPU和CPU的訪問存儲器和調(diào)度模式不同、而容易造成時間浪費等問題，通常可從數(shù)據(jù)分割、數(shù)據(jù)共享、存儲算法等方面進行改進。

表1對上述三種矩陣壓縮存儲方法進行了對比。

表1 稀疏矩陣壓縮為稠密矩陣的算法對比

3.2 矩陣分塊算法

采用矩陣分塊算法，能夠?qū)⒋笮途仃嚪指畛啥鄠€小型矩陣，有利于進行并行運算，減輕運算壓力。下面分別介紹對角加邊法，單、雙對角加邊法這三種分塊算法：

1）對角加邊法矩陣A進行對角加邊變換（block bordered diag?onal form，BBDF）［5～7］的方式大致可以分為三類：節(jié)點撕裂法，支路切割法，統(tǒng)一網(wǎng)絡(luò)分塊法。下面分別介紹這三種分塊方法：

（1）節(jié)點撕裂法

節(jié)點撕裂法［8］是將兩個系統(tǒng)中交集部分取一合適節(jié)點，進行撕裂，使得兩個系統(tǒng)都多一個節(jié)點，使得研究的系統(tǒng)簡單方便，如圖1所示。

圖1 節(jié)點撕裂法

其中S1和S2為兩個區(qū)域電網(wǎng)，SA為由邊界節(jié)點構(gòu)成的上級電網(wǎng)，S1in和S2in為兩個區(qū)域電網(wǎng)的內(nèi)部節(jié)點，A為邊界節(jié)點的集合，A1，A2，A3為同一邊界節(jié)點在不同區(qū)域的分區(qū)，節(jié)點撕裂后平衡狀態(tài)方程為

式中，uA1，uA2為兩個分裂區(qū)域的電壓，θ1，θ2為分裂區(qū)的相角，PA1，QA1，PA2，QA2分別為兩個分裂區(qū)的有功和無功功率。

（2）支路切割法

支路切割法［9～10］即在子網(wǎng)絡(luò)之間的聯(lián)絡(luò)線處加入電流源來代替各子網(wǎng)絡(luò)間的耦合狀態(tài)，電流源矢量作為該方法的協(xié)調(diào)變量，如圖2～3所示。

根據(jù)推導(dǎo)，最終得到第i個子網(wǎng)絡(luò)的注入電流的不平衡方程：

圖2 互聯(lián)電力系統(tǒng)

圖3 互聯(lián)電力系統(tǒng)的分解

（3）統(tǒng)一網(wǎng)絡(luò)分塊法

統(tǒng)一網(wǎng)絡(luò)分塊法是對大規(guī)模電力系統(tǒng)進行網(wǎng)絡(luò)分塊，利用節(jié)點編號優(yōu)先順序的理論來對網(wǎng)絡(luò)進行劃分。通常有靜態(tài)優(yōu)化法，半動態(tài)優(yōu)化法，動態(tài)優(yōu)化法三種編號方法［11］，半動態(tài)優(yōu)化法使用最為廣泛。下面以五節(jié)點網(wǎng)絡(luò)為例依次介紹這三種編號方法：

圖4 有向A圖

①靜態(tài)優(yōu)化法：按照每個節(jié)點出線度由小到大的順序進行編號，則有向A圖編號順序即為靜態(tài)優(yōu)化編號方法。

②半動態(tài)優(yōu)化法：先編號出線度最小的節(jié)點，進行消去，其次為消去節(jié)點后剩余節(jié)點出線度最少的節(jié)點進行編號，依次類推…，圖中1，2節(jié)點編號不用更改，消去這兩個節(jié)點節(jié)點以后，剩余節(jié)點的出線度相同，所以可任意編號。

③動態(tài)優(yōu)化法：對節(jié)點先不進行編號，對任意節(jié)點進行消去，為消去后出現(xiàn)支路最少的節(jié)點先編號，依次類推…，圖中1，2，3，4節(jié)點，由于消去不產(chǎn)生注入元，可任意順序消去。

對角加邊法中三種方法各有特點，表2對這三種方法特點進行比較分析。

2）單、雙邊塊對角法

采用對角加邊分塊法，能夠?qū)⒋笮途仃嚪指畛啥鄠€小塊矩陣，減輕計算壓力，便于并行求解。但由于電力網(wǎng)絡(luò)的不確定性，該方法無法滿足求解要求，因此在該方法基礎(chǔ)上提出單、雙邊塊對角變換法（singly/doubly bordered block diagonal form，SBBD/DBBD）［12］，下面依次介紹這兩種分塊方法。

表2 BBDF中矩陣分塊方法的對比

（1）對A進行雙邊塊對角DDBD變換

對矩陣A進行變換后轉(zhuǎn)置，則DBBD的表現(xiàn)形

式為

其中塊 A11，A22，…，ANN是nN×nN，邊界矩陣C1是n1×p矩陣，R1是 p×n1矩陣，且 p≤n。邊界塊中的每一列叫做耦合列，每一行叫做耦合行，對A進行DDBD變換時，只有當耦合行和列中的 p很小時，才能滿足協(xié)調(diào)矩陣All之間的協(xié)調(diào)性。其中A矩陣中的元素aij需要滿足的邊界條件是， 門 檻 值u∈(0，1)。采用該方法，容易忽略列邊界條件，導(dǎo)致u無法滿足條件，矩陣變換失敗。因此，只有在不消除邊界條件的前提下，通過移動行、列的位置，DBBD的變換方式才有效。

（2）對A進行SBBD變換

采用DBBD變換時，移動矩陣行列會增加矩陣的維數(shù)，增加方程組求解的迭代次數(shù)。為了簡化求解方法，對矩陣進行更進一步變換，即SBBD變換，其表現(xiàn)形式如下：

其中塊 All是 ml×nl矩陣，ml＞nl，且邊界矩陣為ml×p(p≤nl)，矩陣 A中的元素aij需要滿足的邊界條件是，在進行變換時，應(yīng)盡可能減少耦合列的數(shù)量以保證分解塊的穩(wěn)定性。

4 潮流計算線性方程組求解并行算法

線性方程組并行求解最常用的處理方法有共軛梯度法，不完全Cholesky預(yù)處理共軛梯度法，重開始廣義極小殘差法等。下面分別對這些并行算法進行介紹。

4.1 共軛梯度法

共軛梯度法（Conjugate Gradient，CG）［13］也稱共軛斜量法，是求取正定線性方程組首選方法。該方法最大的特點為當求解的線性方程組階數(shù)很高時，通常只需經(jīng)過比階數(shù)小很多的迭代次數(shù)，就可以取得方程組的解。其算法描述如下：

對于 Ax=b線性方程，輸入 A，b，其中 A為n×n矩陣，b為n階向量，輸出x

1）給定初值 x(0)=0，d(0)=0，r(0)=-b，迭代精度 ε＞0；

2）在第k次迭代時，求取x(k)分五步：

（1）計算殘向量：r(k)=Ax(k-1)-b；

（2）計算修正方向向量：

（3）計算迭代步長：

（4）求出新近似量：x(k)=x(k-1)+ξ(k)d(k)；

（5）如不滿足迭代精度要求，則重復(fù)以上步驟。

該算法適用于A為正定矩陣的情況，若A為非對稱非正定矩陣，LDLT分解矩陣不一定存在，使用該算法可能會導(dǎo)致程序中斷。

為了滿足實際需求，求解更大規(guī)模的線性方程組，可在共軛梯度算法上加入不同種類預(yù)處理。下面介紹一種常見的預(yù)處理方法，即不完全Cholesky預(yù)處理共軛梯度法（Modified Incomplete Cholesky Preconditioner Conjugate Gradient，MIC PCG）［14］，算法和PCG基本一致，這里不再贅述。

MIC PCG較PCG相比，迭代次數(shù)顯著減少，在進行向量內(nèi)積運算時，共享內(nèi)存和全局內(nèi)存達到最優(yōu)，能夠提高2倍以上的加速比。

4.2 廣義極小殘差法

重開始廣義極小殘差法（Generalized Minimal Residual Method，GMRES）［15～16］是 krylov子空間法的一種算法，是求解大規(guī)模線性方程組最受歡迎的一種，具有收斂速度快且比較穩(wěn)定的特征。該算法能很好地在GPU上并行求解方程組，通過對顯存中線程塊、線程數(shù)的合理分配，使得每個線程能夠同時并行求解方程，達到性能最優(yōu)。

其算法描述如下：

對于 Ax=b線性方程，輸入 A，b，其中 A為n×n矩陣，b為n階向量，輸出x

1）給定初值 x(0)=0，r(0)=-b，迭代精度ε＞0；

2）選取合適的m，使用Arnoldi迭代方法求?得Vn和，其中m是標準正交基向量的個數(shù)，Vn是由標準正交基組成的m×n矩陣，H~n是(n +1)×n階的上Hessenberg矩陣；

3）極小化殘余量范數(shù)，即 ‖H~nyn-βe1‖，求取使殘余量范數(shù)最小的 yn。其中 β=‖b- Ax0‖，e1=(1 ，0，…0 )；

4）根據(jù) xn=Vnyn，求得 xn；

5）如不滿足迭代要求，則重復(fù)以上步驟。

根據(jù)以上算法，可以得出每次迭代所需要的計算次數(shù)以及它的復(fù)雜度，設(shè)矩陣A為n×n階稀疏陣，m個標準正交基，每行有k個非零元，則可以得出表3。

從表中可以看出，稀疏矩陣與向量乘積、向量內(nèi)積與范數(shù)在整體運算中是最耗時間的。當k比m 2+1大時，稀疏矩陣與向量乘積運算主要影響整個程序的性能，當k比m 2+1小時，向量乘積與范數(shù)運算主要影響整個程序的性能。考慮到在GPU運算過程中，向量內(nèi)積與范數(shù)運算的瓶頸比稀疏矩陣與向量乘積運算少很多，因此加速效果主要受m大小所影響，m越大加速效果越好。但受浮點精度的影響，m不宜取值過大，取值過小又會導(dǎo)致收斂停滯。通常根據(jù)GPU的性能、系數(shù)矩陣A的特征、浮點精度的大小等因素折中選取合適的m值。表4對上述并行算法進行了比較。

表3 每次迭代中不同計算所需的計算次數(shù)和復(fù)雜度

表4 并行算法性能比較

5 結(jié)語

本文討論對比了電力系統(tǒng)潮流并行計算中方程組的常用求解方法。其中，稀疏矩陣壓縮算法適用于節(jié)點間物理連線較少的電力系統(tǒng)，可通過壓縮矩陣來節(jié)省顯存存儲空間，提高并行速度；矩陣分塊算法能進一步降低問題粒度，便于在GPU上實現(xiàn)多線程塊并行求解線性方程組，提高并行化效率；線性方程組并行算法，通過預(yù)處理降低條件數(shù)來提升收斂速度，并在迭代求解過程中實現(xiàn)矩陣與向量內(nèi)積等的并行計算，實現(xiàn)求解效率的有效提升。

目前已有很多較為有效的電力系統(tǒng)潮流并行計算方法，但由于潮流并行計算方程組的求解是一個比較復(fù)雜的問題，也還存在許多尚待進一步完善和解決的問題。如求解線性方程組時，在迭代次數(shù)減少的同時，如何解決收斂性能變差的問題；如何通過線程塊的合理劃分，保證二者之間的平衡性，達到性能最優(yōu)的問題；如何在降低時間復(fù)雜度的情況下，進行網(wǎng)絡(luò)中大型矩陣壓縮算法的合理混合使用的問題；以及對方程組并行求解添加預(yù)處理算法來提高并行效率等等，這些將是未來需要進一步研究的方向。

［1］Hoang-Vu Dang，Bertil Schmidt.The Sliced COO Format for Sparse Matrix Vector Multiplication on CUDA enabled GPUs［J］.Procedia Computer Science，2012，9：57-61.

［2］李佳佳，張秀霞，譚光明，等.選擇稀疏矩陣乘法最優(yōu)存儲格式的研究［J］.計算機研究與發(fā)展，2014，51（4）：882-888.LI Jiajia，ZHANG Xiuxia，TAN Guangming，et al.Study of Choosing the Optimal Storage Format of Sparse Matrix Vector Multiplication［J］.Journal of Computer Research and Development，2014，51（4）：882-888.

［3］商磊.大規(guī)模線性方程組求解的并行算法及應(yīng)用［D］.西安：西北工業(yè)大學，2007，20-23.SHANG Lei.Parallel Algorithms and Applications for Solv?ing Large scale Linear Equations［D］.Xi'an：Northwest Polytechnic University，2007，20-23.

［4］陽王東，李肯立.基于HYB格式稀疏矩陣與向量乘在CPU+GPU異構(gòu)系統(tǒng)中的實現(xiàn)與優(yōu)化［J］.計算機工程與科學，2016，38（2）：205-206.YANG Wangdong，LI Kenli.Implementation and Optimiza?tion of HYB based SPMV on CPU+GPU Heterogeneous Computing Systems［J］.Computer Engineering&Science，2016，38（2）：205-206.

［5］薛巍，舒繼武，王心豐，等.電力系統(tǒng)潮流并行算法的研究進展［J］.清華大學學報（自然科學版），2002，42（9）：1194-1199.XUE Wei，SHU Jiwu，WANG Xinfeng，et al.Advance of parallel algorithms for power flow simulation［J］.J Tsingh ua Univ（Sci&Tech），2001，42（9）：1193-1199.

［6］楊挺，向文平，王洪濤，等.電力系統(tǒng)對角加邊模型的數(shù)據(jù)中心求解方法［J］.中國電機工程學報，2015，35（3）：512-514.YANG Ting，XIANG Wenping，WANG Hongtao，et al.An Algorithm for Solving the Block Bordered Diagonal Form of Electrical Power System in Data Center［J］.Proceedings of the CSEE，2015，35（3）：512-514.

［7］朱永興，張步涵.電力系統(tǒng)潮流分解協(xié)調(diào)并行計算［J］.電力系統(tǒng)及其自動化學報，2010：22（5）.ZHU Yongxing，ZHANG Buhan.Decomposition and Coor?dination Power Flow Parallel Calculation［J］.Proceedings of the CSUEPSA，2010：22（5）.

［8］張伯明，陳壽孫.高等電力系統(tǒng)網(wǎng)絡(luò)分析［M］.北京：清華大學出版社，1996：144-165.ZHANG Boming，CHEN Shousun.Network Analysis of Ad?vanced Power System［M］.Beijing：Tsinghua University Press，1996：144-165.

［9］黃彥全，肖建，劉蘭，等.基于支路切割方法的電力系統(tǒng)潮流并行協(xié)調(diào)算法［J］.電網(wǎng)技術(shù)，2006，30（4）：21-22.HUANG Yanquan，XIAO Jian，LIU Lan，et al.A Coordina?tional Parallel Algorithm for Power Flow Calculation Based on Branch Cutting［J］.Power System Technology，2006，30（4）：21-22.

［10］龔仁敏，張永浩，翟學峰，等.面向整定計算的大型電力網(wǎng)絡(luò)分塊計算的實用化通用計算方法［J］.電力信息化，2007年增刊：42-44.GONG Renmin，ZHANG Yonghao，ZHAI Xuefeng，et al.A Practical General Method of Large Scale Power Net?work Partition Computing for Setting Calculation［J］.Electric power Information，2007 Supplement：42-44.

［11］陳珩.電力系統(tǒng)穩(wěn)態(tài)分析［M］.3版.北京：中國電力出版社，2007：160-163.CHEN Heng.Steady-State Analysis of Power System［M］.The Third Release.Beijing：China Electric Power Press，2007：160-163.

［12］Iain S.Duff，Jennifer A.Scott.Stabilized block diagonal forms for parallel sparse solvers［R］.Council for the Cen?tral Laboratory of the Research Councils，2004.

［13］吳長江.基于CUDA的大規(guī)模線性稀疏方程組求解器的設(shè)計［D］.成都：中國電子科技大學，2013：27-30.WU Changjiang.Design of the of the solver for large scale linear sparse equations based on CUDA［D］.Chengdou：University of Electronic Science and Technol?ogy of China，2013：27-30.

［14］Jiaquan Gao，Ronghua Liang，Jun Wang.Research on the conjugate gradient algorithm with a modified incom?plete Cholesky preconditioner on GPU［J］.J.Parallel Dis?trib.Comput.74，2014.

［15］Kai He，Sheldon X.-D.Tan，Hengyang Zhao，et al.Par?allel GMRES solver for fast analysis of large linear dy?namic systems on GPU platforms［J］.INTEGRATION，the VLSI journal，2016（52）.

［16］張慧，于春肖，白雪婷，等.Krylov子空間E-變換GMRES（m）算法［J］.遼寧工程技術(shù)大學學報（自然科學版），2014，33（9）：1290-1291.ZHANG Hui，YU Chunxiao，BAI Xueting，et al.E-transform GMRES（m）algorithm based on krylov sub?space［J］.Journal of Liaoning Technical University（Natu?ral Science），2014，33（9）：1290-1291.