用于迭代法潮流計算的改進Jacobi預(yù)處理方法

唐 燦, 董樹鋒, 任雪桂, 尹 璐, 鞠 力

(1. 浙江大學(xué)電氣工程學(xué)院, 浙江省杭州市 310027; 2. 北京電力經(jīng)濟技術(shù)研究院, 北京市 100055)

0 引言

牛頓—拉夫遜法是常見的電力系統(tǒng)交流潮流計算方法之一。這一方法需要多次迭代運算,每一次迭代過程中都需要求解線性方程組。當(dāng)方程組的規(guī)模較大時,求解方程組非常耗時,潮流計算的效率受到較大影響。

直接法[1]和迭代法[2-5]是求解線性方程組的兩類方法。其中直接法主要利用矩陣分解技術(shù)求解,例如LU分解法等,雖然直接法能通過有限步驟算出精確解,但計算復(fù)雜度較高,且計算過程不利于并行處理,不適合求解大規(guī)模的線性方程組。相較于直接法,迭代法易于并行計算,在求解大規(guī)模線性方程組時具有明顯的優(yōu)勢[6]。隨著圖形處理器(GPU)的飛速發(fā)展,中央處理器(CPU)與GPU異構(gòu)協(xié)同的計算體系使得串行計算與并行計算協(xié)調(diào)運作,顯著提高了計算能力[7-10],因此可以考慮利用這一計算體系實現(xiàn)迭代法的潮流計算方法。

但是,迭代法也有不足,相比于直接法具有不穩(wěn)定性。迭代法的收斂速度與系數(shù)矩陣的條件數(shù)和譜分布緊密相關(guān),當(dāng)譜分布較為分散時,迭代法的收斂速度明顯降低甚至?xí)l(fā)生不收斂的情況。為保證迭代法線性方程組求解時的穩(wěn)定性,提高求解速度,需要將原方程組轉(zhuǎn)換為等價的、易于求解的線性方程組。針對這一問題,通常采用預(yù)處理技術(shù),即通過對系數(shù)矩陣進行預(yù)處理,改善系數(shù)矩陣的譜分布,使其分布地更為集中,從而提高迭代法的穩(wěn)定性及收斂速度。

在進行預(yù)處理時,往往存在預(yù)處理效果與預(yù)處理時間的矛盾,要想取得較好的預(yù)處理效果往往需要耗費更多時間。本文提出了一種適合于超大規(guī)模電力系統(tǒng)潮流計算中Jacobi矩陣的預(yù)處理方法,這一方法基于適于并行的Jacobi預(yù)處理法,并充分利用Jacobi矩陣自身特點對傳統(tǒng)預(yù)處理方法進行改進。并行求解其預(yù)處理子,不僅擁有高效的求解速度,同時擁有超過其他輕量級預(yù)處理子的效果。經(jīng)算例測試,所提的預(yù)處理方法能有效集中譜分布,改善條件數(shù),提高方程組求解速度。

1 Krylov子空間迭代法的預(yù)處理

Krylov子空間迭代法是一種高效求解大型稀疏線性方程組的有效方法[11-13]。利用Krylov子空間迭代法求解線性方程組時,其收斂速度取決于其譜分布[14-15]。譜分布越集中,其收斂速度越快。

為了提高收斂速度,可以對系數(shù)矩陣進行預(yù)處理。一般來說,預(yù)處理包括以下三種形式。

M-1Ax=M-1b

(1)

(2)

(3)

式中:A為線性方程組的系數(shù)矩陣;b為常數(shù)項;x為待求向量;y為中間向量;M為預(yù)處理子,M1和M2由預(yù)處理子分解得到。

式(1)中分別對Ax和b左乘M-1,然后求解新的線性方程組,其結(jié)果與原方程的解相等；式(2)中對A右乘M-1y,可看作是AM-1(Mx)=b,求解得到Mx,對結(jié)果左乘M-1即可得到與原方程相同的解；式(3)是式(1)和式(2)的結(jié)合,同樣可以得到與原方程同樣的解。對于本文所選取的預(yù)處理子M,既可以使用式(1)與式(2)進行預(yù)處理,也可以將M分解為M1M2,并采用式(3)進行預(yù)處理。

預(yù)處理子的選擇除了應(yīng)當(dāng)使得處理后的方程組更易求解外,也要注意保證原系數(shù)矩陣的稀疏性,盡可能控制非零元注入的數(shù)量。此外,所選的預(yù)處理子的構(gòu)造過程的計算開銷應(yīng)盡可能低。

常用的預(yù)處理方法包括不完全分解預(yù)處理[16-21]、近似逆預(yù)處理[22-23]和多項式預(yù)處理[24-26]等。不完全分解方法的實現(xiàn)過程難以并行,不適用于系數(shù)矩陣規(guī)模較大的情形;近似逆預(yù)處理容易破壞稀疏矩陣逆矩陣的稀疏性;多項式預(yù)處理方法雖具有自然可并行性,但是這一方法得到的預(yù)處理子的預(yù)處理效果較差。

在實際計算中,Jacobi預(yù)處理法是較為常用的預(yù)處理方法。Jacobi預(yù)處理法直接取矩陣的對角元作為預(yù)處理子,其主要優(yōu)勢在于其預(yù)處理的速度非?？?缺點是僅適用于對角元占主導(dǎo)的矩陣。電力系統(tǒng)的潮流計算的Jacobi矩陣是一種類對角元占主導(dǎo)的矩陣,采用Jacobi預(yù)處理法能夠在一定程度上提高潮流計算的速度。但是,潮流計算的Jacobi矩陣并非真正意義上的對角元占主導(dǎo)的矩陣,為了進一步提高預(yù)處理的效果,本文基于適于并行的Jacobi預(yù)處理,充分利用Jacobi矩陣自身特點,設(shè)計了新的預(yù)處理子。

2 改進Jacobi預(yù)處理方法

在潮流計算中,極坐標(biāo)表示的節(jié)點功率方程為[27]:

(4)

式中:Pi和Qi分別為節(jié)點i注入的有功功率和無功功率；Ui和Uj分別為節(jié)點i和節(jié)點j的電壓幅值；δij為節(jié)點i和j的相角差;Gij和Bij分別為節(jié)點i和節(jié)點j間的互電導(dǎo)與互電納。

由此可得到描述電力系統(tǒng)的非線性方程。根據(jù)牛頓—拉夫遜算法,線性化潮流方程可以得到修正方程組[27]為:

(5)

式中:ΔP和ΔQ分別為節(jié)點實際注入功率與迭代過程中值的差額;H,N,J,L分別為迭代過程中Jacobi矩陣的4個子矩陣;ΔV和Δθ分別為相鄰兩次迭代中電壓幅值與相角之間的差值。其中,H為n-1階方陣；N為(n-1)×(n-1-r)階矩陣;J為(n-1-r)×(n-1)階矩陣；L為n-1-r階方陣。H,N,J,L這4個子矩陣都可以通過補充若干全零元素組成的行和列得到n階方陣。

若系統(tǒng)的節(jié)點數(shù)為n,其中PV節(jié)點數(shù)為r，則修正方程組的系數(shù)矩陣為n-r-1階的結(jié)構(gòu)對稱矩陣。

Jacobi矩陣中各元素的計算表達(dá)式為[27]:

(6)

(7)

(8)

(9)

式中:Hii=?ΔPi/?θi,Hij=?ΔPi/?θj,Nii=?ΔPi/?Vi,Nij=?ΔPi/?Vj,Jii=?ΔQi/?θi,Jij=?ΔQi/?θj,Lii=?ΔQi/?Vi,Lij=?ΔQi/?Vj；i和j分別為系統(tǒng)的節(jié)點編號,而不表示元素在矩陣中的行列號。

由式(6)至式(9)可知,Hii,Nii,Jii,Lii分別在H,N,J,L這4個子矩陣中占優(yōu)。因此,利用這些元素,可以構(gòu)造預(yù)處理子。即

(10)

取H,N,J,L中占主導(dǎo)的元素,分別得到Aij,Bij,Cij,Dij并滿足式(11)至式(14)。即

(11)

(12)

(13)

(14)

預(yù)處理子M具有以下性質(zhì)。

1)M為高度稀疏矩陣,其子矩陣A,B,C,D每行的非零元和每列的非零元數(shù)量均不超過1個。同時,A和D為(n-1)×(n-1)階對角矩陣,BT和C為(n-1-r)×(n-1)階矩陣,且非零元的位置相同。

2)對M進行求逆無非零元注入,證明過程如下。

令

(15)

則

(16)

當(dāng)H,L,H-NL-1J,L-JH-1N都可逆時, 解得:

(17)

以M-1的子矩陣X1為例。由于D為n-1階對角矩陣,而BT和C為(n-1-r)×(n-1)階矩陣,且具有相同的非零元素分布,故BD-1C為對角矩陣。則X1=(A-BD-1C)-1與A同為(n-1)×(n-1)階對角矩陣。同理,X2,X3,X4分別與子矩陣B,C,D具有相同的非零元位置。故M-1不引入任何非零元注入。

在計算方面,由于A,D,A-BD-1C，D-CA-1B均為對角矩陣,B與C的非零元位置互為轉(zhuǎn)置且每行至多只有一個非零元?？梢圆捎锰厥獾姆椒焖偾蠼?①稀疏矩陣求逆時,只需對每個對角元分別求逆;②稀疏矩陣相乘時,將其中一個矩陣轉(zhuǎn)置,轉(zhuǎn)置后矩陣的非零元與另一個矩陣中對應(yīng)的非零元分別相乘即可。這一過程具有自然可并行性,可以利用GPU通用計算加速實現(xiàn)。

3 算例分析

3.1 Jacobi矩陣譜分布的改善效果分析

為驗證本文所提出的預(yù)處理方法能夠有效改善系數(shù)矩陣的譜分布,選取IEEE標(biāo)準(zhǔn)系統(tǒng)和PSS/E軟件自帶的BENCH算例測試,對比不同預(yù)處理方法下在預(yù)處理前后的特征值分布。由表1可知,利用本文方法進行預(yù)處理后,各算例譜分布集中于0～2之間,相比較于優(yōu)化前有明顯改善。Jacobi預(yù)處理也能改善矩陣譜分布,但效果不如本方法。

由圖1可知,IEEE 39,IEEE 118,IEEE 300算例和BENCH算例經(jīng)預(yù)處理后,特征值分布變得集中。相比于Jacobi預(yù)處理,本文所提方法效果更為顯著,尤其是虛軸方向上的特征值分布,處理后的最大特征值的虛部約為0。

表1 預(yù)處理前后Jacobi矩陣最大特征值Table 1 Maximum eigenvalue of Jacobi matrix before and after pre-treatment

圖1 預(yù)處理前后Jacobi矩陣特征值分布Fig.1 Distribution of eigenvalue of Jacobi matrix before and after pre-treatment

3.2 M-1與A-1的近似程度

為了比較采用上述方法求取預(yù)處理子時A-1和M-1的近似程度,選取‖Err‖F(xiàn)/n2,‖Err‖1/n,max(Err)作為參考標(biāo)準(zhǔn)。其中Err為誤差矩陣,等于A-1和M-1差值的絕對值；‖Err‖F(xiàn)為Err的Frobenius范數(shù),‖Err‖1為Err的1-范數(shù)?！珽rr‖F(xiàn)/n2為A-1與M-1每個元素的平均偏差,‖Err‖1/n為A-1與M-1每行平均值的最大偏差,max(Err)為A-1與M-1每個元素之間的最大差值。表2統(tǒng)計了用不同算例時，A-1和M-1的近似程度。

表2 預(yù)處理后A-1和M-1的近似程度Table 2 Degree of approximation of A-1 and M-1 after pre-treatment

從表2的結(jié)果可以看出,采用該方法進行預(yù)處理后,M-1和A-1在數(shù)值上十分接近；從‖Err‖F(xiàn)/n2的值可以看出，預(yù)處理后每個元素的平均偏差小于0.001；從‖Err‖1/n可以看出，預(yù)處理后每行平均值最大偏差小于0.1；從max(Err)可以看出,預(yù)處理后對應(yīng)元素之間差額小于1.96。由此可見,該預(yù)處理方法擁有較好的效果。

3.3 預(yù)處理方法的有效性驗證

為了進一步驗證本文提出的預(yù)處理方法的有效性,選取多個算例測試預(yù)處理方法對迭代法迭代次數(shù)的影響。其中,IEEE 8971算例由30個完全相同的IEEE 300系統(tǒng)拼接平衡節(jié)點得到。迭代法采用了穩(wěn)定化的雙共軛梯度(BiCGSTAB)法、雙共軛梯度(BiCG)法、共軛梯度平方(CGS)法、最小殘差法共準(zhǔn)(QMR)這4種常用的Krylov子空間迭代法。表3 對比了無預(yù)處理、Jacobi預(yù)處理和本文預(yù)處理3種情形下求解線性方程組所需的平均迭代次數(shù)。其中,迭代法的計算精度取0.001。

從表3可以看出,將本文的預(yù)處理方法應(yīng)用于不同的迭代法進行線性方程組求解均可以大幅減少迭代次數(shù),從而減少計算時間。當(dāng)系統(tǒng)規(guī)模較大時,如BENCH迭代法,預(yù)處理效果更加明顯,能夠減少90%以上的迭代次數(shù)。相比較于常見的Jacobi預(yù)處理方法,具有顯著的優(yōu)勢。

采用預(yù)處理時能夠減少潮流計算中求解線性方程組時的迭代次數(shù)與迭代時間,但預(yù)處理過程本身需要占據(jù)時間。為了權(quán)衡預(yù)處理方法本身計算的耗時和提升迭代收斂性帶來的計算效率提升之間的關(guān)系,同時考慮到ILU0預(yù)處理子也是商業(yè)軟件常用的預(yù)處理方法,表4統(tǒng)計了本文方法、Jacobi預(yù)處理方法和 ILU0預(yù)處理方法應(yīng)用于計算潮流時所需要的預(yù)處理時間及總時間。計算采用的CPU型號為Intel i7-4710MQ,主頻為2.5 GHz,內(nèi)存為16 GB,線性方程組求解部分采用CULA計算包中的BiCGSTAB方法求解,采用的顯卡為NVDIA GTX860M。其中IEEE 8971算例和IEEE 1496算例由多個IEEE 300算例系統(tǒng)拼接平衡節(jié)點得到。

表3 預(yù)處理與無預(yù)處理時潮流計算中求解線性方程組時的迭代次數(shù)Table 3 Iterations of solving linear equations when doing power flow calculation with and without pre-treatment

由表4可以看出,采用本文的預(yù)處理方法可以大大減少線性方程組的求解時間,并增加求解的穩(wěn)定性。在以上算例中,Jacobi預(yù)處理子耗時均小于1 ms,可忽略不計；ILU0預(yù)處理子耗時相對較大,占總耗時的6.6%～16%；本文方法耗時介于二者之間,小于總耗時的3%。從總耗時來看本方法遠(yuǎn)好于Jacobi預(yù)處理子,與ILU0預(yù)處理子耗時接近。

3.4 算法總體效率測試

本文的預(yù)處理方法能夠有效地加速潮流計算,本文選取了CULA軟件包進行對比測試。在采用CULA軟件包進行潮流計算線性方程組求解時,很多求解器與預(yù)處理子的組合都不能有效地收斂,其中較為穩(wěn)定的組合是BiCGSTAB求解器與ILU0預(yù)處理子或Jacobi預(yù)處理子,但是由于ILU0預(yù)處理子速度更快,于是CULA的數(shù)據(jù)采用了BiCGSTAB+ILU0進行計算。具體數(shù)據(jù)見附錄A表A1。本文提出的預(yù)處理方法應(yīng)用于潮流計算,總體性能優(yōu)于CULA軟件包,可滿足在線計算要求。

表4 不同預(yù)處理方法下進行潮流計算時求解線性方程組耗時Table 4 Time of solving linear equations when doing power flow calculation with different pre-treatment methods

4 結(jié)語

本文提出了一種用于迭代法潮流計算的改進Jacobi預(yù)處理方法。經(jīng)算例測試表明,所提的預(yù)處理方法構(gòu)建預(yù)處理子速度快,同時能夠有效地集中特征值分布、改善條件數(shù)、顯著減少迭代次數(shù),滿足大規(guī)模電網(wǎng)在線潮流計算的需求,具有工程應(yīng)用的潛力和價值。

由于本方法基于牛頓—拉夫遜法潮流計算中的Jacobi矩陣的特性來構(gòu)建預(yù)處理子,尚無法直接應(yīng)用于其他的電力系統(tǒng)計算,在今后的研究中,應(yīng)當(dāng)考慮在其他領(lǐng)域的應(yīng)用。另一方面,改進Jacobi預(yù)處理法預(yù)處理后的系數(shù)矩陣仍有優(yōu)化空間,可考慮多個預(yù)處理子結(jié)合的多步預(yù)處理,進一步加速Krylov子空間法的收斂速度。

附錄見本刊網(wǎng)絡(luò)版(http://www.aeps-info.com/aeps/ch/index.aspx)。

參考文獻(xiàn)

[1] 徐曉飛,曹祥玉,姚旭,等.一種基于Doolittle LU分解的線性方程組并行求解方法[J].電子與信息學(xué)報,2010,32(8):2019-2022.

XU Xiaofei, CAO Xiangyu, YAO Xu, et al. Parallel solving method of linear equations based on Doolittle LU decomposition[J]. Journal of Electronics & Information Technology, 2010, 32(8): 2019-2022.

[2] HOU G, WANG L. A generalized iterative method and comparison results using projection techniques for solving linear systems[J]. Applied Mathematics and Computation, 2009, 215(2): 806-817.

[3] AHAMED A C, MAGOULES F. Iterative methods for sparse linear systems on graphics processing unit[C]// 2012 IEEE 14th International Conference on High Performance Computing and Communications & 2012 IEEE 9th International Conference on Embedded Software and Systems, June 25-27, 2012, Liverpool, United Kingdom: 836-842.

[4] LEON F D, SERNLYEN A. Iterative solvers in the Newton power flow problem: preconditioners, inexact solutions and partial Jacobi updates[J]. IEEE Proceedings: Generation, Transmission and Distribution, 2002, 149(4): 479-484.

[5] KULKARNI A Y, PAI M A, SAUER P W. Iterative solver techniques in fast dynamic calculations of power systems[J]. International Journal of Electrical Power & Energy Systems, 2001, 23(3): 237-244.

[6] 蔡大用,陳玉榮.用不完全LU分解預(yù)處理的不精確潮流計算方法[J].電力系統(tǒng)自動化,2002,26(8):11-14.

CAI Dayong, CHEN Yurong. Solving power flow equations with inexact newton methods preconditioned by incomplete LU factorization with partially fill-in[J]. Automation of Electric Power Systems, 2002, 26(8): 11-14.

[7] 王海峰,陳慶奎.圖形處理器通用計算關(guān)鍵技術(shù)研究綜述[J].計算機學(xué)報,2013,36(4):757-772.

WANG Haifeng, CHEN Qingkui. General purpose computing of graphics processing unit: a survey[J]. Chinese Journal of Computers, 2013, 36(4): 757-772.

[8] 郭春輝.基于GPU的電力系統(tǒng)并行計算的研究[D].濟南:山東大學(xué),2013.

[9] 張逸飛,嚴(yán)正,趙文愷,等.基于GPU的分塊約化算法在小干擾穩(wěn)定分析中的應(yīng)用[J].電力系統(tǒng)自動化,2015,39(22):90-97.DOI:10.7500/AEPS20150126012.

ZHANG Yifei, YAN Zheng, ZHAO Wenkai, et al. Application of block reduction algorithm in power system small-signal stability analysis[J]. Automation of Electric Power Systems, 2015, 39(22): 90-97. DOI: 10.7500/AEPS20150126012.

[10] 周挺輝,趙文愷,嚴(yán)正,等.基于圖形處理器的電力系統(tǒng)稀疏線性方程組求解方法[J].電力系統(tǒng)自動化,2015,39(2):74-80.DOI:10.7500/AEPS20131119005.

ZHOU Tinghui, ZHAO Wenkai, YAN Zheng, et al. A method for solving sparse linear equations power systems based on GPU[J]. Automation of Electric Power Systems, 2015, 39(2): 74-80. DOI: 10.7500/AEPS20131119005.

[11] ZHANG J. Preconditioned Krylov subspace methods for solving nonsymmetric matrices from CFD applications[J]. Computer Methods in Applied Mechanics and Engineering, 1998, 189(3): 825-840.

[12] SAAD Y, VORST H. Iterative solution of linear systems in the 20th century[J]. Journal of Computational and Applied Mathematics, 2000, 123(1): 1-33.

[13] 陳穎,沈沉,梅生偉,等.基于改進Jacobi-Free Newton-GMRES(m)的電力系統(tǒng)分布式潮流計算[J].電力系統(tǒng)自動化,2006,30(9):5-8.

CHEN Ying, SHEN Cheng, 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-8.

[14] AXELSSON O, LINDSKOG G. On the eigenvalue distribution of a class of preconditioning methods[J]. Numerische Mathematik,1986, 48(5): 479-498.

[15] JIA Z. The convergence of Krylov subspace methods for large unsymmetric linear systems[J]. Acta Mathematica Sinica, New Series, 1998, 14(4): 507-518.

[16] ZHANG J. A grid-based multilevel incomplete LU factorization preconditioning technique for general sparse matrices[J]. Applied Mathematics and Computation, 2001, 43(4): 483-500.

[17] SAAD Y, ZHANG J. Enhanced multi-level block ILU preconditioning strategies for general sparse linear systems[J]. Journal of Computational and Applied Mathematics, 2001, 130(1): 99-118.

[18] WILLE S O, STAFF O, LOULA A. Block and full matrix ILU preconditioners for parallel finite element solvers[J]. Computer Methods in Applied Mechanics and Engineering, 2002, 191(13): 1381-1394.

[19] SHEN C, ZHANG J. Parallel two level block ILU preconditioning techniques for solving large sparse linear systems[J]. Parallel Computing, 2002, 28(10):1451-1475.

[20] LEE J, ZHANG J, LU C. Incomplete LU preconditioning for large scale dense complex linear systems from electromagnetic wave scattering problems[J]. Journal of Computational Physics, 2003, 185(1): 158-175.

[21] 柳建新,蔣鵬飛,童孝忠,等.不完全LU分解預(yù)處理的BICGSTAB算法在大地電磁二維正演模擬中的應(yīng)用[J].中南大學(xué)學(xué)報(自然科學(xué)版),2009,40(2):484-491.

LIU Jianxin, JIANG Pengfei, TONG Xiaozhong, et al. Application of BICGSTAB algorithm with incomplete LU decomposition preconditioning to two-dimensional magnetotelluric forward modeling[J]. Journal of Central South University (Science and Technology), 2009, 40(2): 484-491.

[22] KAEBI A, KERAYECHIAN A, TOUTOUNIAN F. Approximate inverse preconditioner by computing approximate solution of Sylvester equation[J]. Applied Mathematics and Computation, 2005, 170(2): 1067-1076.

[23] ZHANG J. A sparse approximate inverse preconditioner for parallel preconditioning of general sparse matrices[J]. Applied Mathematics and Computation, 2002, 130(1): 63-85.

[24] NOTAY Y. Polynomial acceleration of iterative schemes associated with subproper splittings[J]. Journal of Computational and Applied Mathematics, 1988, 24(1): 153-167.

[25] CERDAN J, MAR N, MART N A. Polynomial preconditioners based on factorized sparse approximate inverses[J]. Applied Mathematics and Computation, 2002, 133(1): 171-186.

[26] ZHANG J, ZHANG L. Efficient CUDA polynomial preconditioned conjugate gradient solver for finite element computation of elasticity problems[J]. Mathematical Problems in Engineering, 2013(6): 1-12.

[27] 陳德?lián)P,李亞樓,江涵,等.基于道路樹分層的大電網(wǎng)潮流并行算法及其GPU優(yōu)化實現(xiàn)[J].電力系統(tǒng)自動化,2014,38(22):63-69.DOI:10.7500/AEPS20131014009.

CHEN Deyang, LI Yalou, JIANG Han, et al. A parallel power flow algorithm for large-scale grid based on stratified path trees and its implementation on GPU[J]. Automation of Electric Power Systems, 2014, 38(22): 63-69. DOI: 10.7500/AEPS20131014009.