基于自適應(yīng)冪級數(shù)初始點的電力系統(tǒng)全純嵌入潮流并行計算

姜 濤，張道遠(yuǎn)，李 雪，張 勇，陳厚合

（東北電力大學(xué) 現(xiàn)代電力系統(tǒng)仿真控制與綠色電能新技術(shù)教育部重點實驗室，吉林 吉林 132012）

0 引言

潮流計算是電力系統(tǒng)規(guī)劃與運行的基礎(chǔ)。目前，牛頓-拉夫遜（Newton-Raphson，NR）法已在電力系統(tǒng)潮流計算中廣泛應(yīng)用［1-2］，但采用NR 法求解潮流存在以下不足［3-4］：NR 法對初值敏感，不合理的初值將導(dǎo)致潮流計算收斂速度慢甚至不收斂；高維潮流雅可比矩陣易奇異，從而導(dǎo)致潮流計算結(jié)果不收斂。隨著高比例可再生能源的并網(wǎng)以及電力電子裝備的大規(guī)模應(yīng)用［5］，電力系統(tǒng)運行環(huán)境更加復(fù)雜多變，這使得采用NR 法求解電力系統(tǒng)潮流時的初值選取愈加困難。此外，在電力供應(yīng)日益緊張的背景下，當(dāng)電網(wǎng)運行處于臨界狀態(tài)［6］時，采用NR 法計算潮流會使雅可比矩陣奇異，從而導(dǎo)致潮流求解失敗的風(fēng)險增大［7］。為適應(yīng)未來電力系統(tǒng)運行場景復(fù)雜多變的趨勢［8］以及提高電網(wǎng)安全穩(wěn)定分析效率，亟需探索適用性更好的潮流計算新方法［9］。

全純嵌入法（holomorphic embedding method，HEM）是一種高維非線性方程組求解方法，不同于NR 法，該方法基于復(fù)分析理論，采用遞歸思想求解非線性方程組，無須設(shè)定初值，不形成雅可比矩陣，具有全局收斂性，在求解電力系統(tǒng)潮流時具有諸多優(yōu)勢［10］。文獻(xiàn)［11］首次將HEM 引入電力系統(tǒng)潮流計算中，提出基于全純嵌入的電力系統(tǒng)潮流計算構(gòu)想，構(gòu)建PQ 節(jié)點的HEM 潮流模型；文獻(xiàn)［12］探討電力系統(tǒng)潮流計算中的PV 節(jié)點HEM 建模方法；文獻(xiàn)［13］提出電力系統(tǒng)的全純嵌入潮流計算方法（holomorphic embedding load flow method，HELM），并從算法收斂性的角度論證該方法在系統(tǒng)潮流有解時可準(zhǔn)確收斂至系統(tǒng)真實解，在潮流無解時可通過潮流解振蕩發(fā)出無解信號［14］；文獻(xiàn)［15］計及換流站的控制特點，將HELM應(yīng)用到交直流電力系統(tǒng)的潮流計算中。

然而，HELM 仍存在以下問題［16-17］：①HELM 在求解電力系統(tǒng)潮流時參考平啟動方式確定電壓冪級數(shù)初始點，當(dāng)系統(tǒng)規(guī)模較大，潮流解的實際電壓值偏離該初始點時，雖然HELM 能求解出系統(tǒng)潮流解，但其收斂所需的冪級數(shù)階數(shù)較高，計算量較大，從而導(dǎo)致計算耗時增加；②在求解各節(jié)點電壓的過程中，基于Padé 近似的節(jié)點電壓解析延拓法的計算量較大，計算耗時較長。為解決問題①，文獻(xiàn)［18］結(jié)合迭代法與HEM，將高斯-賽德爾法或快速解耦法迭代3次的計算結(jié)果作為HELM 的電壓冪級數(shù)初始點，在一定程度上改善了HELM 的收斂性，但計算過程復(fù)雜，收斂速度也受高斯-賽德爾法或快速解耦法計算結(jié)果的影響。為解決問題②，文獻(xiàn)［19］計及高階Padé近似的Froissart doublets 效應(yīng)，采用潮流重啟計算策略來縮短HELM 求解電力系統(tǒng)潮流過程中Padé近似的計算耗時，在一定程度上提高了Padé 近似的計算效率，但其根本仍是基于中央處理器（central processing unit，CPU）的串行架構(gòu)進(jìn)行計算。隨著系統(tǒng)節(jié)點規(guī)模的增大，HELM 中Padé 近似的總計算量將急劇增加。由于系統(tǒng)各節(jié)點電壓的Padé 近似計算間相互獨立，存在并行潛力，因此，文獻(xiàn)［20］采用矩陣拼接的方式將各節(jié)點電壓Padé 近似計算中的線性方程組進(jìn)行合并，并利用圖形處理器（graphics processing unit，GPU）對其進(jìn)行并行求解，實現(xiàn)了各節(jié)點電壓Padé 近似的并行計算，但合并后的矩陣維度隨著系統(tǒng)節(jié)點規(guī)模的增大而急劇增大，且矩陣元素也隨著冪級數(shù)階數(shù)的增加而變化，這導(dǎo)致CPU 與GPU 間的數(shù)據(jù)交互頻繁且交互量增大，從而使數(shù)據(jù)交互耗時顯著增長，降低了GPU并行計算效率。

本文提出一種基于自適應(yīng)冪級數(shù)初始點的全純嵌入潮流并行計算方法。首先，根據(jù)節(jié)點導(dǎo)納矩陣和平衡節(jié)點電壓自適應(yīng)調(diào)整系統(tǒng)全純嵌入潮流方程中節(jié)點電壓冪級數(shù)系數(shù)的初始點，進(jìn)而構(gòu)建基于自適應(yīng)冪級數(shù)初始點的全純嵌入潮流方程；其次，基于冪級數(shù)系數(shù)遞歸關(guān)系式建立冪級數(shù)系數(shù)求解遞歸方程組，并求解冪級數(shù)系數(shù)，實現(xiàn)全純隱函數(shù)顯式化；然后，根據(jù)Padé 近似的計算特點，基于Numba 的多核CPU 并行計算架構(gòu)實現(xiàn)各節(jié)點電壓Padé 近似的并行化，以提高潮流計算效率；最后，通過節(jié)點數(shù)在4～13 802 范圍內(nèi)的不同規(guī)模測試系統(tǒng)算例對所提全純嵌入潮流并行計算方法進(jìn)行驗證和分析。

1 電力系統(tǒng)潮流的全純嵌入計算方法

不同于NR 法采用數(shù)值迭代思想將非線性潮流方程近似轉(zhuǎn)化為線性方程進(jìn)行求解，HELM 采用遞歸思想求解電力系統(tǒng)潮流，基于復(fù)分析理論構(gòu)造節(jié)點電壓和PV節(jié)點無功功率的隱式全純函數(shù)，利用全純函數(shù)冪級數(shù)的展開特性，通過求解全純函數(shù)的冪級數(shù)將隱式全純函數(shù)顯式化，實現(xiàn)對潮流方程的求解。HELM 將非線性潮流方程求解問題轉(zhuǎn)化為隱式電壓函數(shù)的顯式化求解問題，因此，無須預(yù)設(shè)初值，也不形成雅可比矩陣，理論上該方法有良好的收斂性。

采用HELM 求解潮流包括全純嵌入潮流模型的構(gòu)建和全純嵌入潮流模型的求解2個部分。

1.1 電力系統(tǒng)全純嵌入潮流模型的構(gòu)建

在直角坐標(biāo)系下，節(jié)點數(shù)為N的電力系統(tǒng)潮流方程［21］可表示為：

式中：Yij為連接節(jié)點i、j的支路導(dǎo)納；Vi為節(jié)點i的電壓相量；S～i為節(jié)點i的注入功率；“*”為共軛運算符。

根據(jù)電力系統(tǒng)不同類型節(jié)點的特點［13］，在式（1）中嵌入復(fù)變量α，構(gòu)造關(guān)于α的各節(jié)點電壓的全純函數(shù)和PV 節(jié)點無功功率全純函數(shù)，參考文獻(xiàn)［15］，傳統(tǒng)電力系統(tǒng)潮流的全純嵌入方程為：

式中：l、p、w分別為PQ 節(jié)點、PV 節(jié)點、平衡節(jié)點編號；Ylk，tr、Ypk，tr、Yl，sh、Yp，sh分別為串聯(lián)支路導(dǎo)納矩陣中第l行第k列元素、串聯(lián)支路導(dǎo)納矩陣中第p行第k列元素、PQ 節(jié)點l對地并聯(lián)支路導(dǎo)納、PV 節(jié)點p對地并聯(lián)支路導(dǎo)納；Vk(α)、Vl(α)、Vp(α)、Vw(α)分別為節(jié)點k、l、p、w電壓的全純函數(shù)；Pp為節(jié)點p的注入有功功率；Qp(α)為節(jié)點p無功功率的全純函數(shù)；Vpm為節(jié)點p的電壓幅值；Vwm為節(jié)點w的電壓幅值。

全純嵌入相關(guān)理論如附錄A所示。

1.2 電力系統(tǒng)全純嵌入潮流模型的求解

不同于傳統(tǒng)NR 法采用數(shù)值迭代思想求解非線性潮流方程，HELM 采用遞歸思想求解電力系統(tǒng)潮流。HELM 的求解過程可總體分為冪級數(shù)初始點的獲取、冪級數(shù)系數(shù)遞歸求解方程的構(gòu)建、全純嵌入潮流函數(shù)的顯式化、Padé 近似4 個部分，限于篇幅，詳細(xì)求解過程如附錄B 所示。由式（B3）可知，HELM通過求解初始狀態(tài)得到電壓冪級數(shù)初始點（電壓冪級數(shù)常數(shù)項）。在任意系統(tǒng)中，由Ylk，tr、Ypk，tr構(gòu)成的導(dǎo)納矩陣對稱且每行元素之和為0，此外，平衡節(jié)點電壓冪級數(shù)初始點Vw［0］和PV節(jié)點電壓冪級數(shù)初始點Vp［0］均為1+j0 p.u.，因此，當(dāng)且僅當(dāng)所有PQ 節(jié)點電壓冪級數(shù)初始點Vl［0］同時為1+j0 p.u.且所有PV 節(jié)點的Qp［0］同時為0 時，才可求得HELM 電壓冪級數(shù)初始點為1+j0 p.u.?？芍?，傳統(tǒng)HELM 電壓冪級數(shù)初始點只能固定為1+j0 p.u.。然而，當(dāng)系統(tǒng)規(guī)模較大時，網(wǎng)絡(luò)拓?fù)涓鼮閺?fù)雜，實際電壓解偏離1+j0 p.u.較多，這將導(dǎo)致傳統(tǒng)HELM 需較高的冪級數(shù)階數(shù)才能實現(xiàn)潮流收斂。由Padé 近似部分可知，基于斯塔爾理論，對角／近對角Padé近似函數(shù)可生成冪級數(shù)最大解析延拓［22］，以保證電壓全純函數(shù)在α=1 處可靠收斂，本文采用Padé 近似實現(xiàn)各節(jié)點電壓解Vx（1）（x=l，p，w）的準(zhǔn)確求取。

2 基于自適應(yīng)冪級數(shù)初始點的全純嵌入潮流計算方法的潮流并行求解

針對傳統(tǒng)HELM 在求解電力系統(tǒng)潮流時參考平啟動模式確定節(jié)點電壓冪級數(shù)的初始點，可能導(dǎo)致潮流收斂所需冪級數(shù)階數(shù)增大，從而導(dǎo)致計算耗時增加的問題，本文提出一種基于自適應(yīng)冪級數(shù)初始點的全純嵌入潮流計算方法（adaptive power series initial point based holomorphic embedding load flow method，A-HELM）。A-HELM 根據(jù)節(jié)點導(dǎo)納矩陣和平衡節(jié)點電壓自適應(yīng)調(diào)整電壓冪級數(shù)初始點，以提高潮流收斂性和計算效率。同時，考慮到求解較大規(guī)模系統(tǒng)全純嵌入潮流時，在采用基于Padé 近似的解析延拓法求解各節(jié)點電壓的過程中，各節(jié)點電壓的Padé 近似計算彼此獨立，具有良好的并行潛質(zhì)特點，借助多核CPU的多線程計算優(yōu)勢，將A-HELM 中逼近各節(jié)點電壓真實解的Padé 近似計算過程并行化，以提升全純嵌入潮流的求解效率。

2.1 電力系統(tǒng)的A-HELM潮流計算模型

首先，構(gòu)建電力系統(tǒng)的A-HELM 潮流計算模型，然后，利用待定系數(shù)法推導(dǎo)A-HELM 潮流計算模型的冪級數(shù)系數(shù)遞歸求解方程，實現(xiàn)A-HELM 潮流計算模型的求解。

2.1.1 A-HELM潮流計算模型的構(gòu)建

將式（1）所示潮流方程采用矩陣形式表示為：

式中：Y為節(jié)點導(dǎo)納矩陣；V為節(jié)點電壓向量；S為節(jié)點注入復(fù)功率向量。

由于平衡節(jié)點的電壓幅值和相角已知，參考NR法計算潮流的方法，不計及平衡節(jié)點功率方程，則式（3）變?yōu)榻惦A潮流方程，如式（4）所示。

式中：Yreduced、Vreduced、Sreduced分別為不計及平衡節(jié)點的節(jié)點導(dǎo)納矩陣、節(jié)點電壓向量、注入功率向量；Islack為平衡節(jié)點的等效注入電流向量，如式（5）所示。

式中：Yslack為由平衡節(jié)點與其余節(jié)點間互導(dǎo)納元素組成的列向量；Vslack為系統(tǒng)平衡節(jié)點電壓。

將Yreduced拆分為串聯(lián)支路導(dǎo)納矩陣Yred，tr和對地并聯(lián)支路導(dǎo)納矩陣Yred，sh兩部分，如式（6）所示。

為便于理解，以包括1個平衡節(jié)點、2個PQ 節(jié)點和2個PV節(jié)點的5節(jié)點系統(tǒng)為例，簡要描述式（3）—（6）中矩陣的變化過程，如附錄C圖C1所示。

將式（6）代入式（4）得：

在式（7）中嵌入復(fù)變量α，可得：

根據(jù)不同類型節(jié)點的特點，由式（8）得到PQ 節(jié)點和PV節(jié)點的A-HELM潮流計算方程為：

式中：Ylw、Ypw分別為PQ節(jié)點l與平衡節(jié)點w間的互導(dǎo)納、PV節(jié)點p與平衡節(jié)點w間的互導(dǎo)納。

為使A-HELM 嵌入方程可解，應(yīng)保證線性方程組數(shù)量與待求變量數(shù)量相等?？紤]到PV 節(jié)點的電壓幅值Vpm已知，電壓相角未知，應(yīng)補充A-HELM 的PV節(jié)點電壓幅值約束方程，如式（10）所示。

當(dāng)式（9）和式（10）中復(fù)變量α=0 時，將式（3）代入式（9）和式（10），得到A-HELM的初始狀態(tài)為：

為求解式（11），將式（3）代入式（7），并令α=0，得到：

將式（5）代入式（12），則式（11）中A-HELM 潮流方程的電壓冪級數(shù)初始點可由式（13）求得。

式中：Vreduced[]0 為所有電壓冪級數(shù)初始點組成的列向量；Vw為平衡節(jié)點w的電壓幅值組成的列向量。

由式（11）—（13）可知，A-HELM 的節(jié)點電壓冪級數(shù)初始狀態(tài)是由系統(tǒng)參數(shù)決定的。與式（4）所示的HELM 電壓冪級數(shù)初始狀態(tài)對比可知，所提A-HELM 不是將各節(jié)點電壓冪級數(shù)初始點強(qiáng)制設(shè)置為1+j0 p.u.，而是根據(jù)系統(tǒng)參數(shù)獲取冪級數(shù)初始點。隨著系統(tǒng)運行狀態(tài)的改變，A-HELM 可根據(jù)系統(tǒng)節(jié)點導(dǎo)納矩陣生成的Yred，tr、Yslack和平衡節(jié)點電壓Vslack實現(xiàn)對節(jié)點電壓冪級數(shù)初始點的自適應(yīng)調(diào)整。

相較于傳統(tǒng)HELM，所提A-HELM 具有如下特點：在潮流計算方程中無須對平衡節(jié)點進(jìn)行建模；嵌入方程中的節(jié)點電壓冪級數(shù)初始狀態(tài)根據(jù)系統(tǒng)參數(shù)特征進(jìn)行選取，以加快電壓冪級數(shù)的收斂速度。

為確保所構(gòu)建的A-HELM 潮流計算方程具有可行性，應(yīng)檢驗所構(gòu)建的A-HELM 全純嵌入潮流方程是否滿足以下4 個條件［2］：①當(dāng)嵌入的復(fù)變量α=0時，系統(tǒng)各節(jié)點電壓的冪級數(shù)初始點真實存在且易于求取；②當(dāng)α=1 時，對應(yīng)實際潮流解，此時全純嵌入潮流方程應(yīng)與原潮流方程完全等價；③所構(gòu)造的待求隱式函數(shù)具有全純性；④根據(jù)斯塔爾定理，在到達(dá)鞍節(jié)分岔點前的α路徑上，全純嵌入潮流方程具有唯一解。

對于條件①，由式（13）可得到A-HELM 冪級數(shù)初始點，因此，所提方法滿足條件①。對于條件②，當(dāng)α=1 時，式（13）完全等價于式（1）的原潮流方程，因此，所提方法滿足條件②。對于條件③和條件④，采用HELM 求解電力系統(tǒng)潮流的合理性和實用性已被驗證［23］，A-HELM 僅是對HELM 的嵌入方程進(jìn)行恒等變換，而未破壞嵌入方程的全純性和解的唯一性，因此，所提方法滿足條件③和條件④。

綜上可知，所構(gòu)建的潮流計算模型合理且可行。

2.1.2 A-HELM潮流計算模型的遞歸求解

借鑒HELM求解電力系統(tǒng)潮流的基本思路，采用A-HELM求解式（9）、（10）的冪級數(shù)系數(shù)過程如下。

為便于描述，令Xx(α)=1/V*x(α*)(x=l，p，w)，將式（3）代入式（9）、（10）所示的A-HELM 潮流計算方程，基于同次冪級數(shù)系數(shù)的待定系數(shù)法，分別構(gòu)建PQ節(jié)點的電壓冪級數(shù)系數(shù)遞歸求解方程以及PV 節(jié)點的節(jié)點電壓和無功功率的冪級數(shù)系數(shù)遞歸求解方程，分別如式（14）、（15）所示。

為便于計算，將式（14）、（15）中復(fù)數(shù)方程的實部和虛部拆分為純實數(shù)型線性方程。由于復(fù)數(shù)型遞歸關(guān)系式兩側(cè)實部、虛部分別相等，可得對應(yīng)的實數(shù)型遞歸方程為：

式中：Glk、Blk分別為串聯(lián)支路導(dǎo)納矩陣中第l行第k列導(dǎo)納元素的實部和虛部；Gpk、Bpk分別為串聯(lián)支路導(dǎo)納矩陣中第p行第k列導(dǎo)納元素的實部和虛部；“（Re）”“（Im）”分別表示相應(yīng)復(fù)數(shù)的實部、虛部，后同。

根據(jù)式（16）所示的冪級數(shù)系數(shù)遞歸求解方程組，可由電壓冪級數(shù)初始點Vx[0] (x=l，p，w)逐階獲取各階冪級數(shù)系數(shù)Vx[k]和Qx[k]，實現(xiàn)A-HELM 中各節(jié)點電壓全純函數(shù)Vx(α)和PV 節(jié)點無功功率全純函數(shù)Qx(α)的顯式化。

2.2 Padé近似并行化

目前，在電力系統(tǒng)潮流并行計算領(lǐng)域，關(guān)于交流潮流并行化的研究主要聚焦于加速修正方程組的并行求解。然而，對于全純嵌入潮流的計算，Padé近似部分的計算耗時在其潮流求解總時間中占比較高，其原因在于，雖然采用近對角Padé 近似法可提高單個節(jié)點電壓的解析延拓計算效率，但仍需對系統(tǒng)中每個節(jié)點的電壓逐一進(jìn)行Padé 近似計算，即串行計算系統(tǒng)各待求節(jié)點的電壓，具體過程如附錄C 圖C2所示，隨著系統(tǒng)節(jié)點規(guī)模的增加，A-HELM 中Padé近似的總計算量也將顯著增大。由于A-HELM 需對各節(jié)點電壓依次進(jìn)行Padé 近似且近似過程相互獨立，因此，所提A-HELM 的Padé 近似過程具有較高的潛在并行能力，可根據(jù)節(jié)點電壓Padé 近似的計算特點設(shè)計相應(yīng)的并行計算架構(gòu)，并采用合適的并行計算技術(shù)將各節(jié)點電壓解析延拓的Padé 近似計算過程并行化，以提升Padé近似的總體計算效率。

針對數(shù)值并行計算，Anaconda 公司推出了Python 環(huán)境下的Numba 高性能計算庫，以解決工程科學(xué)計算領(lǐng)域中的數(shù)值計算密集型問題。Numba采用即時編譯技術(shù)，基于Numba 并行架構(gòu)可對Python原代碼進(jìn)行編譯優(yōu)化，加快運行效率，且可釋放Python 中的全局解釋器鎖，以多線程并行計算的方式實現(xiàn)自動并行化，為此，本文借助Numba 高性能計算庫實現(xiàn)所提A-HELM潮流求解的并行化。

所提A-HELM 潮流方程節(jié)點電壓的Padé 近似過程中涉及大量的數(shù)組運算，結(jié)合Numba 并行架構(gòu)的優(yōu)點，可借助Numba 對A-HELM 中的Padé 近似函數(shù)代碼進(jìn)行即時編譯，并在此基礎(chǔ)上進(jìn)行并行化處理，以提高Padé近似的效率，進(jìn)而提升所提A-HELM的潮流求解效率。Numba的CPU自動并行化技術(shù)可自動識別具有并行計算潛力的數(shù)值計算部分并對其進(jìn)行優(yōu)化處理，以生成最優(yōu)并行方案，再將計算任務(wù)分割為若干個獨立的子任務(wù)，以并行執(zhí)行計算任務(wù)。在所提A-HELM 的節(jié)點電壓Padé 近似過程中，由于各節(jié)點電壓的Padé 近似計算互不依賴，不受計算順序的影響，Numba可識別該特點，進(jìn)而將各節(jié)點電壓Padé 近似的串行計算過程拆分為多個獨立的子任務(wù)，并映射至各CPU 核中進(jìn)行計算，在不影響計算準(zhǔn)確性的基礎(chǔ)上實現(xiàn)各節(jié)點電壓Padé 近似的并行化?；贜umba 并行架構(gòu)，所提A-HELM 的節(jié)點電壓Padé近似并行計算實現(xiàn)流程如圖1所示。

圖1 Padé近似并行計算示意圖Fig.1 Schematic diagram of parallel calculation of Padé approximation

圖1 中，在對Padé 近似函數(shù)代碼進(jìn)行即時編譯的基礎(chǔ)上，Numba開辟并行區(qū)域，以劃分計算塊的方式將各節(jié)點電壓Padé 近似的計算任務(wù)分解為若干個獨立的子任務(wù)，并為并行區(qū)域的各線程分配1 個Padé 計算塊?？紤]到每個節(jié)點電壓Padé 近似的計算量基本相同，為使各并行線程的計算負(fù)載均衡，可采用類似于OpenMP 中的線程靜態(tài)調(diào)度方式來劃分Padé計算塊［4］，即各Padé計算塊大小（子任務(wù)包含的節(jié)點數(shù)）相等，從而均勻分配計算任務(wù)。以包含N個待求節(jié)點的系統(tǒng)為例，Numba 自動根據(jù)CPU 的核數(shù)安排t個線程進(jìn)入并行區(qū)域的迭代塊，每個Padé 計算塊的大小為N/t，即每個線程負(fù)責(zé)執(zhí)行N/t個節(jié)點電壓的Padé 近似計算，因此，在不同的CPU 核中將同時執(zhí)行并行區(qū)域的t個線程，進(jìn)而將A-HELM 中逐次求解各節(jié)點電壓的Padé 近似計算任務(wù)均衡分配至各處理器核中，實現(xiàn)多核CPU 平臺下節(jié)點電壓Padé近似的并行化。

2.3 計算流程

綜上，所提基于A-HELM 的并行計算方法流程如附錄C圖C3所示，具體步驟如下。

1）根據(jù)式（13），由系統(tǒng)節(jié)點導(dǎo)納矩陣與平衡節(jié)點電壓計算系統(tǒng)潮流的冪級數(shù)初始點。

2）利用式（16）所示的遞歸關(guān)系，由冪級數(shù)初始點開始逐階計算各階電壓冪級數(shù)系數(shù)。

3）對所求得的電壓冪級數(shù)系數(shù)進(jìn)行近似計算，采用式（B14）所示的Padé 近似方法計算各節(jié)點電壓近似值，并將此計算過程交由多核CPU進(jìn)行并行處理。

4）由式（17）和式（18）計算得到各節(jié)點功率組成的列向量Scalc及其最大功率不平衡量max_mis。

式中：VPade為所有節(jié)點電壓近似值組成的列向量；“·”表示向量點乘運算；Scalc_PQ為由VPade計算得到的PQ 節(jié)點的復(fù)功率列向量；SPQ為系統(tǒng)已知的PQ 節(jié)點注入的復(fù)功率列向量；Pcalc_PV為由VPade計算得到的PV 節(jié)點的有功功率列向量；PPV為系統(tǒng)給定PV 節(jié)點注入的有功功率列向量；mPQ為PQ 節(jié)點復(fù)功率不平衡量列向量；mPV為PV 節(jié)點有功功率不平衡量列向量；mPQ(a)、mPV(b)分別為mPQ、mPV的第a、b個元素。

5）將max_mis與給定潮流收斂精度ε進(jìn)行比較：若max_mis＜ε，則輸出系統(tǒng)潮流解；否則返回執(zhí)行步驟2）— 4），直至冪級數(shù)階數(shù)達(dá)到階數(shù)上限，若此時max_mis仍不小于ε，則程序終止，輸出“算法不收斂”。

3 算例分析

為驗證所提基于A-HELM 的并行計算方法的準(zhǔn)確性和可行性，通過不同規(guī)模的測試系統(tǒng)對該方法進(jìn)行分析和驗證，測試系統(tǒng)相關(guān)信息如附錄C 表C1 所示。計算平臺硬件配置為CPU AMD Ryzen5 4 600 H，6核12線程，主頻為3.00 GHz，內(nèi)存為16 GB；軟件平臺采用Anaconda 的Spyder 平臺（Python 版本為3.8.8）?；谠撥浖脚_開發(fā)所提A-HELM 并行計算方法，并將其潮流計算結(jié)果與基于NR 法的開源潮流計算軟件PYPOWER 的計算結(jié)果進(jìn)行對比［24］。算例中潮流收斂精度均為10-3p.u.。

3.1 方法準(zhǔn)確性與通用性驗證

3.1.1 方法準(zhǔn)確性驗證

本節(jié)以Case 6ww 測試系統(tǒng)為例，詳細(xì)說明采用所提方法求解電力系統(tǒng)潮流的基本過程。

首先，根據(jù)本文基于自適應(yīng)冪級數(shù)初始點的全純嵌入潮流計算流程，生成Case 6ww 測試系統(tǒng)的節(jié)點導(dǎo)納矩陣Y，并從Y中提取Yslack和Yreduced，進(jìn)一步根據(jù)式（6）求取Yred，tr，矩陣變化過程如圖C3所示。

然后，將Case 6ww 測試系統(tǒng)的Yred，tr、Yslack、Vslack代入式（13），求取節(jié)點電壓冪級數(shù)初始點。根據(jù)式（16）的冪級數(shù)系數(shù)求解遞歸方程，利用所提方法得到的冪級數(shù)初始點V1［0］—V5［0］求取各待求節(jié)點電壓的1階冪級數(shù)系數(shù)V1［1］—V5［1］，此時，式（B14）所示電壓冪級數(shù)表達(dá)式為Vx（α）=Vx［0］+Vx［1］α，對應(yīng)的潮流最大功率不平衡量為0.104 3 p.u.，大于潮流收斂精度10-3p.u.，不滿足潮流收斂條件，因此，需返回式（16），由冪級數(shù)系數(shù)遞推關(guān)系式繼續(xù)求取各節(jié)點電壓的2 階冪級數(shù)系數(shù)V1［2］—V5［2］，再通過Padé近似計算獲得新的電壓解，并判定是否滿足潮流收斂條件，如此循環(huán)，直至滿足潮流收斂判定條件。電壓冪級數(shù)初始點以及1 階和2 階冪級數(shù)系數(shù)如附錄C表C2所示。圖2展示了潮流求解過程中電壓冪級數(shù)階數(shù)與最大功率不平衡量的關(guān)系（圖中ΔPmax為系統(tǒng)最大功率不平衡量，為標(biāo)幺值，后同）。

圖2 電壓冪級數(shù)階數(shù)與最大功率不衡量的關(guān)系Fig.2 Relation between voltage power series order and maximum power unbalance

由圖2 可知，隨著電壓冪級數(shù)階數(shù)的增加，系統(tǒng)最大功率不平衡量減小，即潮流收斂精度與冪級數(shù)階數(shù)密切相關(guān)，冪級數(shù)階數(shù)越高，潮流計算精度越高，但計算耗時也隨之增加。Case 6ww 測試系統(tǒng)僅需3 階電壓冪級數(shù)即可實現(xiàn)潮流快速收斂，利用AHELM 得到的各節(jié)點電壓3 階冪級數(shù)系數(shù)和各節(jié)點實際電壓解如表C2所示。

表1為Case 6ww測試系統(tǒng)中A-HELM 串行計算方法與NR 法潮流計算結(jié)果的對比（表中2種方法的電壓幅值為標(biāo)幺值，相對誤差為A-HELM 串行計算方法結(jié)果相對NR 法結(jié)果的誤差，后同）。由表可知，A-HELM 串行計算方法與NR 法的潮流計算結(jié)果幾乎一致，最大相對誤差小于0.2%，由此驗證了所提A-HELM 串行計算方法可實現(xiàn)電力系統(tǒng)潮流的準(zhǔn)確求解。

表1 Case 6ww測試系統(tǒng)中A-HELM串行計算方法與NR法潮流計算結(jié)果的對比Table 1 Load flow calculation results comparison between A-HELM serial calculation method and NR method in Case 6ww test system

Case 6ww 測試系統(tǒng)中A-HELM 并行計算方法與NR 法潮流計算結(jié)果的對比如附錄C 表C3 所示。由表可知，相較于NR 法，A-HELM 并行計算方法的最大相對誤差也小于0.2%，驗證了所提A-HELM 并行計算方法可實現(xiàn)潮流的準(zhǔn)確求解。

A-HELM 串行計算方法、A-HELM 并行計算方法、傳統(tǒng)HELM、NR 法求解Case 6ww 系統(tǒng)潮流所需時間的對比如附錄C 表C4 所示。由表可知：相較于傳統(tǒng)HELM 和NR 法，所提A-HELM 計算方法的計算效率有所提高；相較于所提A-HELM 串行計算方法，所提A-HELM并行計算方法的計算效率進(jìn)一步提高，但由于Case 6ww 測試系統(tǒng)的節(jié)點數(shù)較少，潮流求解的總計算量相對較小，因此，對Padé 近似的并行化并未使計算效率顯著提升，但隨著系統(tǒng)規(guī)模的擴(kuò)大，所提A-HELM并行計算方法的優(yōu)勢將進(jìn)一步凸顯。

3.1.2 方法通用性驗證

本節(jié)進(jìn)一步通過節(jié)點數(shù)為4～13 802的不同規(guī)模測試系統(tǒng)對所提A-HELM 串行和并行計算方法進(jìn)行測試，計算過程與Case 6ww 測試系統(tǒng)的潮流求解過程大體一致。以NR 法的計算結(jié)果為基準(zhǔn)，不同規(guī)模測試系統(tǒng)中A-HELM 串行計算方法與A-HELM 并行計算方法潮流計算結(jié)果的最大誤差如附錄C 表C5 所示。由表可知，所提A-HELM 串行和并行計算方法的節(jié)點電壓幅值及相角最大相對誤差均小于0.1 %，有效驗證了所提方法可實現(xiàn)不同規(guī)模電力系統(tǒng)潮流的準(zhǔn)確計算，具有良好的通用性。

3.2 方法收斂性驗證

為驗證所提方法的收斂性，本節(jié)對比所提AHELM串行計算方法、傳統(tǒng)HELM、NR法在不同規(guī)模測試系統(tǒng)下的計算耗時，結(jié)果如表2所示。

表2 A-HELM串行計算方法、傳統(tǒng)HELM、NR法的計算耗時對比Table 2 Calculation time comparison among A-HELM serial calculation method，traditional HELM and NR method

由表2 可知，隨著測試系統(tǒng)規(guī)模的擴(kuò)大，系統(tǒng)網(wǎng)絡(luò)拓?fù)涓訌?fù)雜，NR 法不收斂的問題愈加嚴(yán)重，這主要是由于NR 法采用平啟動初值，無法滿足潮流收斂條件，而對于NR 法不收斂的測試系統(tǒng)，所提AHELM串行計算方法均可收斂至系統(tǒng)潮流解。

A-HELM 串行計算方法與傳統(tǒng)HELM 冪級數(shù)初始點潮流計算耗時的對比如附錄C 表C6 所示。由表可知：雖然傳統(tǒng)HELM 能收斂，但由于冪級數(shù)初始點與真實潮流解偏差相對較大，計算速度較慢；相較于傳統(tǒng)HELM，本文所提A-HELM 串行計算方法的自適應(yīng)冪級數(shù)初始點更接近真實潮流解，因此，潮流收斂所需冪級數(shù)階數(shù)更少，驗證了所提A-HELM 串行計算方法具有較高的計算效率。

A-HELM 與傳統(tǒng)HELM 初始最大功率不平衡量的對比如附錄C 表C7 所示。由表可知：相較于傳統(tǒng)HELM，采用所提A-HELM可有效減少系統(tǒng)功率不平衡量，這是由于傳統(tǒng)HELM節(jié)點電壓冪級數(shù)的初始值通常選取平啟動方式下的電壓幅值，而所提A-HELM根據(jù)系統(tǒng)節(jié)點導(dǎo)納矩陣和平衡節(jié)點電壓自適應(yīng)調(diào)整節(jié)點電壓的冪級數(shù)初始點，這使得節(jié)點電壓冪級數(shù)初始點的電壓值更加接近實際潮流解；相較于傳統(tǒng)HELM，當(dāng)系統(tǒng)節(jié)點數(shù)較少時，所提A-HELM 的優(yōu)勢不明顯，其最大功率不平衡量未顯著減小，而當(dāng)系統(tǒng)規(guī)模變大，網(wǎng)絡(luò)拓?fù)涓鼜?fù)雜時，所提A-HELM 的優(yōu)勢凸顯，其最大功率不平衡量明顯減小。

為進(jìn)一步分析冪級數(shù)初始點對潮流收斂的影響，圖3 以Case 4gs、Case 300、Case 2 746wp、Case 3 012wp、Case 3 375wp、Case 10 790 測試系統(tǒng)為例，給出了A-HELM 與傳統(tǒng)HELM 的系統(tǒng)最大功率不平衡量隨電壓冪級數(shù)階數(shù)增加的變化趨勢。由圖可知：當(dāng)系統(tǒng)規(guī)模較小時，A-HELM 與傳統(tǒng)HELM 的收斂趨勢較接近；隨著系統(tǒng)規(guī)模的增大，A-HELM 的系統(tǒng)最大功率不平衡量均小于傳統(tǒng)HELM。在潮流收斂過程中，傳統(tǒng)HELM 的節(jié)點電壓冪級數(shù)系數(shù)的初始點均固定為1+j0 p.u.，這使得潮流計算起始的系統(tǒng)功率不平衡量較大，進(jìn)而影響電壓高階次冪級數(shù)系數(shù)的求取，導(dǎo)致求解過程中出現(xiàn)較大的收斂振蕩。與傳統(tǒng)HELM 相比，自適應(yīng)初始點可減小A-HELM在計算初始時的潮流功率不平衡量，因此，減小了冪級數(shù)初始點對電壓高階次冪級數(shù)系數(shù)求解的影響。

圖3 不同測試系統(tǒng)中A-HELM與傳統(tǒng)HELM最大功率不平衡量變化趨勢對比Fig.3 Variation trend comparison of maximum power unbalance between A-HELM and traditional HELM in different test systems

由上述分析可知，與傳統(tǒng)HELM 相比，A-HELM基于節(jié)點導(dǎo)納矩陣和平衡節(jié)點電壓自動生成更適用的冪級數(shù)初始點，使得A-HELM 的潮流收斂趨勢更加平滑，從而減少了潮流收斂所需的電壓冪級數(shù)階數(shù)，因此，自適應(yīng)冪級數(shù)初始點的選取方法使得AHELM比傳統(tǒng)HELM具有更好的潮流收斂性。

3.3 計算效率對比

本節(jié)進(jìn)一步驗證所提A-HELM 串行計算方法和A-HELM 并行計算方法的計算效率。在內(nèi)存需求方面，在求解全純嵌入潮流模型的過程中，A-HELM 冪級數(shù)系數(shù)遞歸求解方程組的系數(shù)矩陣元素均為系統(tǒng)網(wǎng)絡(luò)節(jié)點導(dǎo)納值，在較大規(guī)模的實際算例中該系數(shù)矩陣占用內(nèi)存的情況如附錄C 表C8 所示。由表可知：該系數(shù)矩陣為常數(shù)矩陣且具有很高的稀疏度，因此，在編程時可采用稀疏格式存儲矩陣，并以該格式進(jìn)行科學(xué)計算；實際求解時，A-HELM 系數(shù)矩陣占用的內(nèi)存略高于傳統(tǒng)NR 法的雅可比矩陣占用的內(nèi)存，但由于A-HELM 系數(shù)矩陣是系統(tǒng)節(jié)點導(dǎo)納矩陣的一種變形，因此，在系統(tǒng)網(wǎng)絡(luò)拓?fù)浣Y(jié)構(gòu)不發(fā)生變化時，該系數(shù)矩陣數(shù)值不變，相較于傳統(tǒng)NR 法的雅可比矩陣，無須對該系數(shù)矩陣進(jìn)行反復(fù)更新和存儲。

A-HELM 串行計算方法、A-HELM 并行計算方法以及傳統(tǒng)HELM 在求解不同規(guī)模測試系統(tǒng)時的潮流計算耗時對比如附錄C 表C9 所示。由表可知：A-HELM 串行計算方法的計算耗時短于傳統(tǒng)HELM的計算耗時，A-HELM 并行計算方法的計算耗時顯著短于A-HELM 串行計算方法和傳統(tǒng)HELM 的計算耗時；與傳統(tǒng)HELM 相比，A-HELM 串行計算方法的計算效率平均提升28.86 %，其中，對于Case 2 746wp測試系統(tǒng)，A-HELM 串行計算方法的計算效率提升了40 % 以上。由于A-HELM 的冪級數(shù)初始點比傳統(tǒng)HELM 更有利于潮流收斂，A-HELM 達(dá)到收斂所需電壓冪級數(shù)階數(shù)更低，計算量更小，縮短了整體用時，此外，A-HELM 簡化了傳統(tǒng)HELM 的算法邏輯，這使得求取每階冪級數(shù)系數(shù)的耗時更短，從而在整體上提升了算法的計算速度。

為更直觀地體現(xiàn)A-HELM 并行計算方法在求解電力系統(tǒng)潮流時的高效性，本文將A-HELM 串行計算方法的計算耗時與A-HELM 并行計算方法的計算耗時之比定義為加速比，加速比越大說明A-HELM并行計算方法的計算效率越高。表3 給出了不同測試系統(tǒng)潮流計算中的加速比。

表3 不同測試系統(tǒng)潮流計算中的加速比Table 3 Acceleration ratios in load flow calculation of different test systems

由表3 可知：對于不同規(guī)模的測試系統(tǒng)，加速比均穩(wěn)定在4～5 之間，這說明A-HELM 并行計算方法具有較好的加速效果；測試系統(tǒng)規(guī)模越大，A-HELM并行計算方法的加速效果越明顯，這是由于隨著系統(tǒng)規(guī)模的增大，采用Padé 近似求解節(jié)點電壓的總計算量也急劇增加，而對各節(jié)點電壓Padé 近似的并行化可提升A-HELM 算法的計算效率，加快電力系統(tǒng)潮流的求解速度。

4 結(jié)論

本文提出一種基于自適應(yīng)冪級數(shù)初始點的電力系統(tǒng)潮流全純嵌入并行計算方法，并通過節(jié)點數(shù)在4～13 802間的不同規(guī)模測試系統(tǒng)對所提方法進(jìn)行分析和驗證，得到結(jié)論如下：

1）所提A-HELM 的潮流計算模型合理且可行，可實現(xiàn)電力系統(tǒng)潮流的準(zhǔn)確求解；

2）與傳統(tǒng)HELM 確定節(jié)點電壓冪級數(shù)初始點的方法不同，所提A-HELM 可根據(jù)節(jié)點導(dǎo)納矩陣與平衡節(jié)點電壓自適應(yīng)調(diào)整系統(tǒng)潮流冪級數(shù)初始點，從而加快電壓冪級數(shù)的收斂速度，降低潮流收斂所需冪級數(shù)階數(shù)，因此，所提A-HELM 具有較好的潮流收斂性能和計算效率；

3）基于多核CPU 的Numba 并行架構(gòu)，所提方法實現(xiàn)了A-HELM 中Padé 近似計算的并行化，在保證計算結(jié)果準(zhǔn)確的前提下提升了潮流求解效率。

附錄見本刊網(wǎng)絡(luò)版（http：//www.epae.cn）。