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

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

0 引言

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

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

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

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

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

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

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

M-1Ax=M-1b

(1)

(2)

(3)

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

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

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

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

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

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

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

(4)

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

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

(5)

式中:ΔP和ΔQ分別為節(jié)點(diǎn)實(shí)際注入功率與迭代過程中值的差額;H,N,J,L分別為迭代過程中Jacobi矩陣的4個(gè)子矩陣;Δ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個(gè)子矩陣都可以通過補(bǔ)充若干全零元素組成的行和列得到n階方陣。

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

Jacobi矩陣中各元素的計(jì)算表達(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é)點(diǎn)編號(hào),而不表示元素在矩陣中的行列號(hào)。

由式(6)至式(9)可知,Hii,Nii,Jii,Lii分別在H,N,J,L這4個(gè)子矩陣中占優(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個(gè)。同時(shí),A和D為(n-1)×(n-1)階對(duì)角矩陣,BT和C為(n-1-r)×(n-1)階矩陣,且非零元的位置相同。

2)對(duì)M進(jìn)行求逆無(wú)非零元注入,證明過程如下。

令

(15)

則

(16)

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

(17)

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

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

3 算例分析

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

為驗(yàn)證本文所提出的預(yù)處理方法能夠有效改善系數(shù)矩陣的譜分布,選取IEEE標(biāo)準(zhǔn)系統(tǒng)和PSS/E軟件自帶的BENCH算例測(cè)試,對(duì)比不同預(yù)處理方法下在預(yù)處理前后的特征值分布。由表1可知,利用本文方法進(jìn)行預(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ù)處理子時(shí)A-1和M-1的近似程度,選取‖Err‖F(xiàn)/n2,‖Err‖1/n,max(Err)作為參考標(biāo)準(zhǔn)。其中Err為誤差矩陣,等于A-1和M-1差值的絕對(duì)值；‖Err‖F(xiàn)為Err的Frobenius范數(shù),‖Err‖1為Err的1-范數(shù)?！珽rr‖F(xiàn)/n2為A-1與M-1每個(gè)元素的平均偏差,‖Err‖1/n為A-1與M-1每行平均值的最大偏差,max(Err)為A-1與M-1每個(gè)元素之間的最大差值。表2統(tǒng)計(jì)了用不同算例時(shí)，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é)果可以看出,采用該方法進(jìn)行預(yù)處理后,M-1和A-1在數(shù)值上十分接近；從‖Err‖F(xiàn)/n2的值可以看出，預(yù)處理后每個(gè)元素的平均偏差小于0.001；從‖Err‖1/n可以看出，預(yù)處理后每行平均值最大偏差小于0.1；從max(Err)可以看出,預(yù)處理后對(duì)應(yīng)元素之間差額小于1.96。由此可見,該預(yù)處理方法擁有較好的效果。

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

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

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

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

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

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

3.4 算法總體效率測(cè)試

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

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

4 結(jié)語(yǔ)

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

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

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

參考文獻(xiàn)

[1] 徐曉飛,曹祥玉,姚旭,等.一種基于Doolittle LU分解的線性方程組并行求解方法[J].電子與信息學(xué)報(bào),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ì)算方法[J].電力系統(tǒng)自動(dò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] 王海峰,陳慶奎.圖形處理器通用計(jì)算關(guān)鍵技術(shù)研究綜述[J].計(jì)算機(jī)學(xué)報(bào),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)并行計(jì)算的研究[D].濟(jì)南:山東大學(xué),2013.

[9] 張逸飛,嚴(yán)正,趙文愷,等.基于GPU的分塊約化算法在小干擾穩(wěn)定分析中的應(yīng)用[J].電力系統(tǒng)自動(dò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)自動(dò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] 陳穎,沈沉,梅生偉,等.基于改進(jìn)Jacobi-Free Newton-GMRES(m)的電力系統(tǒng)分布式潮流計(jì)算[J].電力系統(tǒng)自動(dò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é)報(bào)(自然科學(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)化實(shí)現(xiàn)[J].電力系統(tǒng)自動(dò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.