王明軒， 陳 穎， 黃少偉， 魏 巍， 常曉青

(1. 清華大學電機工程與應(yīng)用電子技術(shù)系， 北京市 100084； 2. 國網(wǎng)四川省電力公司電力科學研究院， 四川省成都市 610072)

0 引言

隨著電網(wǎng)規(guī)模的不斷擴大和智能電網(wǎng)的進一步發(fā)展，其遠距離、重負荷、大區(qū)域聯(lián)網(wǎng)的特點日益凸顯[1-2]。對于某些重負荷系統(tǒng)、梳狀放射性系統(tǒng)或具有鄰近多解的系統(tǒng)，潮流方程常常呈現(xiàn)病態(tài)。此時，普通潮流算法無法收斂，需要采用病態(tài)潮流算法。

病態(tài)潮流算法的相關(guān)研究由來已久。文獻[3-5]提出并發(fā)展了最優(yōu)乘子法，即在牛頓—拉夫遜法的基礎(chǔ)上進行了魯棒性變換，添加乘子并優(yōu)選其數(shù)值，從而搜索病態(tài)潮流解。文獻[6]采用隱式Cholesky分解方法進行大規(guī)模病態(tài)潮流計算。文獻[7]則利用狀態(tài)空間搜索的方法研究了病態(tài)系統(tǒng)的求解問題。文獻[8]構(gòu)造了基于同倫法的病態(tài)潮流算法，該方法發(fā)揮同倫方法的大范圍收斂性，使得病態(tài)潮流易于收斂。

由于含有特殊處理邏輯，計算效率較低，前述病態(tài)潮流算法一般并不適合求解普通潮流。文獻[9]則提出了基于連續(xù)牛頓法(continuous Newton’s method，CNM)的潮流求解方法，同時適用于普通和病態(tài)潮流求解。該方法將潮流方程組等效為常微分方程組，進而可采用不同的數(shù)值積分方法求取其穩(wěn)態(tài)，對應(yīng)獲得原始潮流方程的解。上述方法雖然通用性較強，但由于算法復(fù)雜度較高，求解大規(guī)模(病態(tài))潮流時效率較低。因此，有必要研究相關(guān)加速方法，提升大規(guī)模電網(wǎng)基于CNM的病態(tài)潮流算法在線應(yīng)用能力。

針對基于CNM的潮流求解方法，本文采用定雅可比矩陣和2階龍格—庫塔(RK)方法實現(xiàn)算法邏輯。進一步，采用中央處理器+圖形處理器(CPU+GPU)協(xié)同計算架構(gòu)對大規(guī)模潮流計算過程進行加速。設(shè)計了GPU計算內(nèi)核，采用海量并發(fā)線程完成潮流求解所涉節(jié)點不平衡功率計算；根據(jù)計算規(guī)模和計算內(nèi)容的不同特性，優(yōu)選軟硬件組合提升其計算效率。最后，以大規(guī)模病態(tài)潮流算例證明所提方法的正確性和實用性。

1 基于CNM的潮流計算方法及其改進

電網(wǎng)潮流可用如下非線性方程組描述：

f(x)=0

(1)

CNM將上式轉(zhuǎn)換為以下常微分方程組：

(2)

式中：J(x)為f(x)在x處的雅可比矩陣。

為減少LU分解計算量，對式(2)進行改造，采用初始點處的定雅可比矩陣J0替換J(x)，形成式(3)。

(3)

文獻[10-11]對采用J0時算法的收斂性進行了理論論證，同時指出該定常處理可能使得收斂速度變慢。但考慮到LU分解在并行架構(gòu)下的巨大計算量，采用該處理對于提高本文算法整體計算效率是有意義的。給定合適的初值，對式(3)進行積分，所得穩(wěn)態(tài)平衡點即為原潮流方程的解。相應(yīng)的，經(jīng)過比選，本文選擇二階顯式龍格—庫塔(RK2)方法[12]對式(3)進行積分計算，具體流程如附錄A圖A1所示。算法中在判定潮流收斂的同時，考慮了待求量改變幅度和積分時的限制，兼顧了效率和準確性。

2 CPU+GPU協(xié)同架構(gòu)下的求解方法

2.1 J0的LU分解

CPU+GPU架構(gòu)中，稀疏矩陣LU分解可采用Right-looking，Left-looking[13]及Crout分解方法[14]等。本文采用了適用于CPU的NICSLU[15]和PARDISO[16]計算庫，以及適用于GPU的cuSOLVER[17]計算庫進行稀疏矩陣的LU分解。后文給出了基于不同數(shù)學庫對典型電網(wǎng)雅可比矩陣進行LU分解的性能對比。通過比選，本文采用NICSLU進行J0的稀疏矩陣LU分解，其特色在于使用Left-looking方法，并可輸出J0因子表，方便連續(xù)積分計算。

2.2 潮流右端項并行計算

經(jīng)過前述改進，所提CNM潮流計算中，潮流右端項(向量)計算是重要的一部分。為了提升其計算效率，設(shè)計了專門的GPU核函數(shù)實現(xiàn)多線程并行處理。

以單節(jié)點有功不平衡量計算為例(其余同理)，相關(guān)計算可以拆分為以下兩個部分，

(4)

考慮到導(dǎo)納矩陣中大部分元素為0，非零元素為與b條線路對應(yīng)的2b個非對角元素及N個對角元素。因而，可設(shè)計(b+N)個線程，每個線程計算一個非零元素對于與其相關(guān)的電流項I產(chǎn)生的遞加量，并將電流項與遞加量相加。這樣實現(xiàn)了排零細粒度并行處理。另一方面，不同線程可能需要同時寫入I的某一元素，從而帶來了資源訪問競爭，阻塞了GPU中的并行計算。為消除此現(xiàn)象，本文采用原子加法操作(atomicAdd函數(shù))實現(xiàn)I并發(fā)修正。

如圖1所示對上述并發(fā)計算過程進行了描述，相關(guān)GPU核函數(shù)偽代碼見附錄A表A1。

由以上方法得到向量I之后，根據(jù)式(4)中的第1個式子，重新安排n個線程，每個線程進行單個元素的計算，即可得到n維潮流右端項(向量)。

2.3 前代回代求解

前代回代過程實際上是稀疏矩陣為上三角和下三角矩陣的線性方程組的求解過程。同樣，在本文潮流算法中需要進行稀疏的前代回代計算，每步積分需要進行兩次前代回代。

有多個提供稀疏前代回代算法實現(xiàn)的CPU與GPU數(shù)學庫，其中適用于CPU的NICSLU庫和適用于GPU的cuSPARSE庫[15]效率較高，且接口友好。后文給出了基于以上兩種數(shù)學庫完成潮流求解的前代回代部分的性能比較。通過對比，本文采用NICSLU計算庫進行前代回代求解。

3 算例測試

3.1 測試平臺

進行測試的軟、硬件平臺如附錄B表B1所示。其中正確性測試僅在CPU上完成。

3.2 測試算例

本文采用了不同節(jié)點數(shù)目的6個算例。將39節(jié)點算例稱為case39，其余同理。case39和case300均為IEEE標準算例，case2383，case3120，case9241，case13659均取自MATPOWER 6.0，其中case2383和case3120分別來自2000年冬季和2008年夏季的波蘭電網(wǎng)數(shù)據(jù)，case9241和case13659則為來自歐洲電網(wǎng)的真實算例。所采用的最大規(guī)模的算例節(jié)點數(shù)目為13 659，方程組維數(shù)為23 225。所計算的均為交流潮流，結(jié)果精度為10-4。涉及的時間測試結(jié)果均為運行1 000次的平均耗時。

3.3 正確性測試

病態(tài)潮流常由于負荷增長引起。負荷增長可能有多種方向，本文測試時采用固定比例同步增長方法。為不失一般性，以case3120為例分析其正確性測試結(jié)果，在標準算例的基礎(chǔ)上不斷增加所有負荷功率和發(fā)電機有功出力。負荷因子值增大時采用本文算法求解所需的迭代次數(shù)如表1所示。系統(tǒng)未進入病態(tài)時迭代次數(shù)為10次左右。隨著負荷水平提高，系統(tǒng)開始進入病態(tài)，算法收斂性變差，迭代次數(shù)隨著病態(tài)程度加重而增加，同時算法求解速度下降。負載和發(fā)電機出力達到基態(tài)的2.3倍時，傳統(tǒng)的牛頓—拉夫遜法已無法收斂。此時采用本文算法，系統(tǒng)殘差r=‖f(xi)‖隨虛擬時間增長變化如圖2所示，殘差逐漸趨于0，病態(tài)潮流得解。負荷因子進一步增長至2.5時，本文算法也無法收斂。

表1 不同負荷因子時的算法迭代次數(shù)Table 1 Iteration times for different load factors

圖2 潮流病態(tài)求解殘差曲線Fig.2 Residual curve of ill-conditioned case

為進一步證明本文方法的正確性，以預(yù)測—校正法進行計算時相應(yīng)病態(tài)潮流斷面的結(jié)果作為參考值，計算本文算法相應(yīng)的電壓和相角誤差值。為證明采用定雅可比矩陣的合理性，同時添加采用變雅可比矩陣的測試結(jié)果作為對比。針對不同規(guī)模的病態(tài)潮流算例進行驗證，結(jié)果如表2所示，表中誤差值為各個節(jié)點的電壓和相角與參考值作差的絕對值中的最大值。電壓誤差(標幺值)和相角誤差(弧度)都不超過0.005，在工程可允許的范圍內(nèi)。

表2 算法正確性測試Table 2 Method accuracy test

3.4 效率測試

3.4.1LU分解計算部分

采用NICSLU，cuSOLVER，PARDISO三種計算庫進行電力系統(tǒng)雅可比矩陣分解，計算效率對比如附錄B圖B1所示。結(jié)果顯示，NICSLU計算庫耗時更短，具有更高的效率。其高效性在矩陣維數(shù)較大時更為明顯。

3.4.2潮流右端項計算部分

下文比較基于GPU和CPU的潮流右端項計算效率。附錄B圖B2展示了不同算例基于CNM的潮流計算中該部分計算總耗時的差異。結(jié)果表明，在算例規(guī)模較小時采用CPU計算潮流右端項效率更高；而在算例規(guī)模較大時，GPU算法可取得明顯加速效果，其加速比變化趨勢如附錄B圖B3所示。

3.4.3前代回代計算部分

采用NICSLU計算庫與cuSPARSE計算庫分別完成前代回代求解，比較不同算例基于CNM的潮流求解中該部分計算總耗時，結(jié)果如附錄B圖B4所示。NICSLU庫中的前代回代計算在各種規(guī)模下均更加高效。

3.4.4整體計算效率

根據(jù)前文測試，設(shè)計切換規(guī)則，針對病態(tài)潮流求解所涉計算內(nèi)容可選擇不同實現(xiàn)方法，使得算法執(zhí)行總體效率最優(yōu)。具體包括：①對于小于1 000節(jié)點的小規(guī)模算例，潮流右端項計算在CPU上完成，求解更大規(guī)模的算例時，則在GPU上完成潮流右端項計算；②前代回代求解部分均在CPU上完成，其中大于1 000節(jié)點的算例求解中涉及CPU與GPU間的通信，但此類通信耗時占比極小，不影響整體算法效率。

根據(jù)以上規(guī)則，對算法進行整體耗時測試。其中病態(tài)潮流算例仍由固定比例同步增長方法得到。值得注意的是，步長Δt的選取對算法收斂性及收斂結(jié)果有較大影響，步長過大可能會造成無法收斂的結(jié)果。對于本文算法，選取Δt=1.0是對于所有算例可行的一個參數(shù)值。為進一步加快迭代收斂的速度，本文根據(jù)不同系統(tǒng)的具體情況，對Δt進行了一定的修正，修正原則為盡量選取使得結(jié)果正確且迭代次數(shù)最少的Δt值。

整體效率測試結(jié)果及步長修正情況如表3所示。表中Δt為不同算例整定的步長參數(shù)值。可以看出，整體算法具有較高的效率，對于10 000節(jié)點以上規(guī)模的系統(tǒng)，可以在100 ms左右完成普通潮流和病態(tài)潮流求解。

表3 整體計算效率測試Table 3 Efficiency test of all calculations

以case13659的普通潮流計算為例，各部分計算占比如附錄B圖B5所示。經(jīng)優(yōu)化，潮流右端項計算占比已經(jīng)非常??；而LU分解計算雖僅一次，仍占據(jù)了大量時間，這也證明了采用定雅可比矩陣的必要性。減少LU分解次數(shù)和耗時對于提高整體算法效率十分重要。

4 結(jié)語

本文設(shè)計并完成了適用于CPU+GPU協(xié)同架構(gòu)的大規(guī)模病態(tài)潮流計算方法。將CNM應(yīng)用于大規(guī)模病態(tài)潮流求解中并進行相關(guān)改進，形成高效大規(guī)模病態(tài)潮流計算方法，并驗證了所提方法的正確性和高效性。本文研究僅針對病態(tài)負荷水平下潮流的單次計算過程，沒有解決在負荷增加過程中的潮流連續(xù)求解問題。在本文算法基礎(chǔ)上，后續(xù)工作將開展大規(guī)模連續(xù)潮流GPU加速并行計算的相關(guān)研究。

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

參 考 文 獻

[1] 李立浧,張勇軍,陳澤興,等.智能電網(wǎng)與能源網(wǎng)融合的模式及其發(fā)展前景[J].電力系統(tǒng)自動化,2016,40(11):1-9.DOI:10.7500/AEPS20150912002.

LI Licheng, ZHANG Yongjun, CHEN Zexing, et al. Merger between smart grid and energy-net: mode and development prospects[J]. Automation of Electric Power Systems, 2016, 40(11): 1-9. DOI: 10.7500/AEPS20150912002.

[2] 于宏文,鄭春偉,汪洋,等.智能電網(wǎng)調(diào)度控制系統(tǒng)中歷史數(shù)據(jù)服務(wù)優(yōu)化方案[J].電力系統(tǒng)自動化,2016,40(19):113-118.DOI:10.7500/AEPS20150928001.

YU Hongwen, ZHENG Chunwei, WANG Yang, et al. Historical data service optimization scheme for smart grid dispatching and control systems[J]. Automation of Electric Power Systems, 2016, 40(19): 113-118. DOI: 10.7500/AEPS20150928001.

[3] IWAMOTO S, TAMURA Y. A load flow calculation method for ill-conditioned power systems[J]. IEEE Transactions on Power Apparatus and Systems, 1981, 100(4): 1736-1743.

[4] 杜正春,周佃民,董繼民.考慮負荷電壓靜特性的最佳乘子牛頓潮流算法[J].中國電機工程學報,2002,22(1):102-105.

DU Zhengchun, ZHOU Dianmin, DONG Jimin. Optimal multiplier Newton method of load flow with static load characteristics[J]. Proceedings of the CSEE, 2002, 22(1): 102-105.

[5] 胡澤春,王錫凡.基于最優(yōu)乘子潮流確定靜態(tài)電壓穩(wěn)定臨界點[J].電力系統(tǒng)自動化,2006,30(6):6-11.

HU Zechun, WANG Xifan. Determination of static voltage collapse critical point based on load flow method with optimal multiplier[J]. Automation of Electric Power Systems, 2006, 30(6): 6-11.

[6] 張道天,嚴正,徐瀟源,等.采用隱式Cholesky分解的大規(guī)模病態(tài)潮流計算[J].電網(wǎng)技術(shù),2016,40(4):1197-1203.

ZHANG Daotian, YAN Zheng, XU Xiaoyuan, et al. Large-scale ill-conditioned power flow calculation using implicit Cholesky factorization method[J]. Power System Technology, 2016, 40(4): 1197-1203.

[7] SHAHRIARI A, MOKHLIS H, BAKAR A H A, et al. The calculation of low voltage solution based on state space search method in ill-conditioned system[J]. Corrosion Science, 2012, 20(8/9): 1059-1066.

[8] 周佃民,廖培金.電力系統(tǒng)病態(tài)潮流的同倫方法求解[J].電力系統(tǒng)及其自動化學報,1999,11(5/6):67-71.

ZHOU Dianmin, LIAO Peijin. Homotopy method for ill-conditioned power system load flow calculation[J]. Proceedings of the CSU-EPSA, 1999, 11(5/6): 67-71.

[9] MILANO F. Continuous Newton’s method for power flow analysis[J]. IEEE Transactions on Power Systems, 2009, 24(1): 50-57.

[10] RAMM A G, SMIRNOVA A B, FAVINI A. Continuous modified Newton’s-type method for nonlinear operator equations[J]. Annali Di Matematica Pura Ed Applicata, 2003, 182(1): 37-52.

[11] AIRAPETYAN R. Continuous Newton method and its modification[J]. Applicable Analysis, 1999, 73(3/4): 463-484.

[12] BUTCHER J C. The numerical analysis of ordinary differential equations: Runge-Kutta and general linear methods[J]. Mathematics of Computation, 1987, 51(183): 693.

[13] SCHENK O, GRTNER K, FICHTNER W. Efficient sparse LU factorization with left-right looking strategy on shared memory multiprocessors[J]. Bit Numerical Mathematics, 2000, 40(1): 158-176.

[14] FORD W. Chapter 11—Gaussian elimination and the LU decomposition[M]// Numerical Linear Algebra with Applications, 2015: 205-239.

[15] CHEN X, WANG Y, YANG H. NICSLU: an adaptive sparse matrix solver for parallel circuit simulation[M]. Piscataway, USA: IEEE Press, 2013.

[16] SCHENK O, GRTNER K. PARDISO[M]. Boston: Springer, 2011: 1458-1464.

[17] NVIDIA. CUDA toolkit documentation v8.0[EB/OL]. [2017-06-23]. http://docs.nvidia.com/cuda/index.html.