林 瓏 李建杰

(1. 重慶水利電力職業(yè)技術(shù)學(xué)院， 重慶 402160; 2. 重慶大學(xué) 輸配電裝備及系統(tǒng)安全與新技術(shù)國家重點(diǎn)實(shí)驗(yàn)室， 重慶 400030)

隨著智能電網(wǎng)的進(jìn)一步發(fā)展，我國電網(wǎng)規(guī)模不斷擴(kuò)大，結(jié)構(gòu)日益復(fù)雜，區(qū)域間耦合日益緊密[1-3]．面對以上變化，各級(jí)調(diào)度應(yīng)具備全網(wǎng)統(tǒng)一的網(wǎng)絡(luò)分析能力，傳統(tǒng)的潮流計(jì)算方法面臨挑戰(zhàn)[4]．而國民經(jīng)濟(jì)的迅速發(fā)展使得電網(wǎng)常常重負(fù)荷運(yùn)行，此時(shí)電力系統(tǒng)潮流方程呈現(xiàn)病態(tài)特征，需要特殊的計(jì)算方法．

病態(tài)潮流計(jì)算相關(guān)研究由來已久．最早出現(xiàn)，也是迄今為止應(yīng)用最廣泛的方法是最優(yōu)乘子法[5-9]，在牛頓-拉夫遜法的基礎(chǔ)上進(jìn)行魯棒性變換，添加乘子并優(yōu)選其數(shù)值，從而搜索病態(tài)潮流解．另一種病態(tài)潮流的求解算法是非線性規(guī)劃法[10]，該方法將潮流計(jì)算中非線性方程組求解的問題轉(zhuǎn)化為一組無約束的非線性規(guī)劃問題，優(yōu)勢在于總是能夠給出潮流有解或無解的解答，但計(jì)算量過大且結(jié)果依賴于優(yōu)化求解方法的選擇．此外，文獻(xiàn)[11]構(gòu)造了基于同倫法的病態(tài)潮流算法，該方法發(fā)揮同倫方法的大范圍收斂性，使得病態(tài)潮流易于收斂．文獻(xiàn)[12]則采用隱式Cholesky分解方法進(jìn)行大規(guī)模病態(tài)潮流計(jì)算．文獻(xiàn)[13]對多種病態(tài)潮流算法進(jìn)行了論述以及比較分析．

以上方法能夠完成特定應(yīng)用場景下的病態(tài)潮流計(jì)算，并已被證明在一定規(guī)模的電力系統(tǒng)中是有效的．但其有效性受系統(tǒng)規(guī)模限制的特性明顯，在系統(tǒng)節(jié)點(diǎn)數(shù)目超過1000后，方法的收斂性受到影響，同時(shí)計(jì)算耗時(shí)呈指數(shù)特性上升，對于大規(guī)模系統(tǒng)效率極低，無法滿足全網(wǎng)一體化分析需求．因此有必要研究適用于大規(guī)模電力系統(tǒng)的新型病態(tài)潮流計(jì)算方法．

本文針對病態(tài)潮流方程的數(shù)學(xué)特性，采用連續(xù)牛頓法將非線性方程組轉(zhuǎn)化為常微分方程組，并進(jìn)行雅克比矩陣的相應(yīng)處理，構(gòu)建潮流計(jì)算的微分求解模型．在此基礎(chǔ)上，采用一定的數(shù)值積分方法，能夠準(zhǔn)確高效地完成大規(guī)模病態(tài)電力系統(tǒng)的潮流求解．結(jié)合實(shí)際系統(tǒng)算例，其正確性得以驗(yàn)證．

1 數(shù)學(xué)原理簡述

連續(xù)牛頓法是一種將非線性方程組求解問題轉(zhuǎn)化為常微分方程組積分求解問題的數(shù)學(xué)方法，最初是通過將牛頓-拉夫遜法和顯示歐拉法求解步驟進(jìn)行類比和替換形成的[14]．現(xiàn)對其原理簡述如下：對于如式(1)的非線性方程組

f(x)=0

(1)

牛頓-拉夫遜迭代法[15]在定點(diǎn)進(jìn)行一階泰勒展開，得到用方程雅克比矩陣表示的基本迭代格式如式(2)所示．式中x為待求解變量，J(xn)為在xn處的雅克比矩陣．

xn+1=xn+Δxn

Δxn=-J(xn)-1*f(xn)

(2)

另一方面，對于如式(3)的常微分方程組，

(3)

顯式歐拉法[16]給出的積分公式為：

xn+1=xn+Δxn

Δxn=Δt*g(xn)

(4)

其中Δt為積分步長，一般有Δt=1．對比式(2)和式(4)可發(fā)現(xiàn)，如果做如式(5)的假設(shè)，則兩式完全相同．

g(xn)=-J(xn)-1*f(xn)

(5)

由式(1)，(3)和(5)，可得到下式：

(6)

式(6)即為CNM的基本數(shù)學(xué)模型．結(jié)合一定的數(shù)值積分方法，該方程可以求解．設(shè)積分求解得到的結(jié)果為x*，易得：

g(x*)=-J(x*)-1*f(x*)=0

(7)

只要J(x*)非奇異，必有f(x*)=0，即x*為原非線性方程組的解．

應(yīng)注意的是，式(6)是否具有穩(wěn)定平衡點(diǎn)也需要驗(yàn)證，否則不能保證數(shù)值積分過程可以求解得到穩(wěn)定結(jié)果．設(shè)g(x)=0的解為x=x*，則x*一定為原方程的一個(gè)平衡點(diǎn)，但該平衡點(diǎn)不一定穩(wěn)定．由相關(guān)數(shù)學(xué)知識(shí)可知，x*為穩(wěn)定平衡點(diǎn)的條件為g(x)在x=x*處的梯度矩陣的所有特征值均在左半平面(即右端函數(shù)雅克比矩陣負(fù)定)．對式(6)進(jìn)行梯度矩陣的計(jì)算，

(8)

明顯是負(fù)定矩陣，式(6)具有穩(wěn)定平衡點(diǎn)．

2 病態(tài)潮流微分求解模型

2.1 針對潮流計(jì)算構(gòu)建的微分方程

電力系統(tǒng)潮流方程表征電力系統(tǒng)中的功率平衡關(guān)系，由電力網(wǎng)絡(luò)參數(shù)與節(jié)點(diǎn)電氣參數(shù)構(gòu)成．對于N節(jié)點(diǎn)電力系統(tǒng)，設(shè)其非平衡節(jié)點(diǎn)數(shù)目為n=N-1，其中PV節(jié)點(diǎn)數(shù)目為r，PQ節(jié)點(diǎn)數(shù)目為n-r，則在極坐標(biāo)下其潮流方程形式如式(9)所示．其中ΔPi、ΔQi為節(jié)點(diǎn)i的有功和無功功率不平衡量，在系統(tǒng)穩(wěn)態(tài)時(shí)其值應(yīng)為0；PGi、PLi為節(jié)點(diǎn)i的發(fā)電機(jī)有功出力和負(fù)荷有功注入，QGi、QLi為節(jié)點(diǎn)i的發(fā)電機(jī)無功出力和負(fù)荷無功注入；Vi、Vj為節(jié)點(diǎn)i和j的電壓幅值，θij=θi-θj為節(jié)點(diǎn)i和j的相角差值；Gij和Bij為網(wǎng)絡(luò)節(jié)點(diǎn)電導(dǎo)矩陣和電納矩陣中第i行j列的元素．

(9)

潮流方程本質(zhì)為一組復(fù)雜的非線性方程組，可用式(1)表達(dá)，其中f為不平衡功率計(jì)算對應(yīng)的映射，x為潮流計(jì)算中的未知狀態(tài)變量，一般為節(jié)點(diǎn)電壓幅值和相角．利用前述CNM數(shù)學(xué)方法，可將其轉(zhuǎn)化為式(6)所示的常微分方程組進(jìn)行求解，能夠得到正確的潮流解．

當(dāng)電力系統(tǒng)重負(fù)荷運(yùn)行時(shí)，潮流方程呈現(xiàn)病態(tài)，其顯著特征之一即為潮流雅克比矩陣接近奇異，對應(yīng)于式(6)即為J(x)-1的計(jì)算效率很低，對于大規(guī)模系統(tǒng)該求逆過程更是極為困難，整體求解遇到瓶頸．為充分利用CNM方法對于基態(tài)和病態(tài)潮流的一致收斂性，避免由于雅克比矩陣求逆造成的收斂性和效率問題，本文對CNM方法進(jìn)行了定雅克比矩陣處理，即在數(shù)值計(jì)算過程中采用固定的雅克比矩陣，其因子表或逆矩陣均只需要形成一次．采用定雅克比矩陣的CNM方法得到的微分方程如式(10)所示，其中J0為固定的雅克比矩陣，一般由初始點(diǎn)x0處計(jì)算得到．

(10)

采用定雅克比矩陣處理后，微分求解模型的正確性需要驗(yàn)證．首先，由于J0非奇異，對式(10)進(jìn)行積分求解的結(jié)果一定為原方程的解．至于式(10)穩(wěn)定平衡點(diǎn)的存在性問題，同樣計(jì)算梯度矩陣，可得：

(11)

無法直接判斷其是否負(fù)定．本文研究發(fā)現(xiàn)選擇平凡點(diǎn)處的J0(如初始求解點(diǎn)處)，采用本文算法計(jì)算實(shí)際電力系統(tǒng)病態(tài)潮流時(shí)矩陣-J0-1*J(x*)的特征根均在左半平面，滿足具有穩(wěn)定平衡點(diǎn)的條件．后文算例測試結(jié)果也印證了其正確性．

2.2 數(shù)值積分方法應(yīng)用

對式(10)所示微分方程進(jìn)行求解時(shí)，需采用一定的數(shù)值積分方法．盡管任何有效的數(shù)值積分方法均可應(yīng)用于通過CNM建立的微分方程的求解[17]，但不同方法的收斂性和收斂速度有很大差異．由第1節(jié)中相關(guān)內(nèi)容可知，若采用一階的顯示歐拉法，則本文算法計(jì)算過程與普通的牛頓-拉夫遜法基本相同；為完成大規(guī)模系統(tǒng)病態(tài)潮流的求解，需選用更高階數(shù)、更精確的數(shù)值積分方法，這也正體現(xiàn)了采用本文算法與牛頓法計(jì)算的區(qū)別．

考慮到計(jì)算需求，本文主要采用龍格庫塔積分方法．該方法分為不同階數(shù)的顯式和隱式求解方法，具有不同特性[18]．為保證算法的最優(yōu)性能，研究中將3種不同階數(shù)的顯式龍格庫塔法RK1、RK2、RK4分別應(yīng)用于實(shí)際系統(tǒng)算例的求解中，通過收斂性和計(jì)算效率的對比進(jìn)行優(yōu)選，最終選擇RK2作為具體實(shí)現(xiàn)時(shí)的積分方法．其積分計(jì)算式如式(12)所示，其中xi為第i步積分時(shí)由電壓幅值和相角組成的狀態(tài)變量，ki(1)、ki(2)為計(jì)算的中間向量，ti為虛擬時(shí)間變量序列，Δt為固定步長．

(12)

2.3 大規(guī)模病態(tài)潮流求解流程

圖1 大規(guī)模病態(tài)潮流求解流程

3 算例分析

本文應(yīng)用了不同規(guī)模的電力系統(tǒng)算例，為方便說明，將n節(jié)點(diǎn)系統(tǒng)算例簡稱為casen．文中小規(guī)模算例case39，case118，case300，均為IEEE標(biāo)準(zhǔn)算例；大規(guī)模系統(tǒng)算例case2383，case3120，case9241以及萬節(jié)點(diǎn)以上算例case13659則為來自波蘭電網(wǎng)與歐洲電網(wǎng)的真實(shí)算例．為全面分析算法對于不同規(guī)模、不同病態(tài)水平的系統(tǒng)潮流求解時(shí)的性能，本文采用發(fā)電機(jī)有功和負(fù)荷有功、無功同比例增長的方法構(gòu)造病態(tài)潮流算例，即如式(13)所示．其中λ為負(fù)荷增長因子，PL、PL0分別為病態(tài)算例和標(biāo)準(zhǔn)算例的所有負(fù)荷有功功率，無功功率同理；PG、PG0分別為病態(tài)算例和標(biāo)準(zhǔn)算例中所有發(fā)電機(jī)的有功出力．

PL=λPL0，QL=λQL0，PG=λPG0

(13)

為提高算法效率，本文采用C/C++語言進(jìn)行算法編程實(shí)現(xiàn)，仿真平臺(tái)為Microsoft Visual Studio 2012；同時(shí)采用Matlab的MATPOWER計(jì)算庫采用牛頓-拉夫遜法進(jìn)行潮流計(jì)算結(jié)果正確性對比，仿真平臺(tái)為Matlab R2015b．計(jì)算機(jī)配置為Intel?Xeon(R) CPU E3-1230 v5，8GB內(nèi)存．計(jì)算中收斂精度均為10-4．

不失一般性，以case2383為例，本文算法及牛頓-拉夫遜法迭代次數(shù)隨負(fù)荷因子的變化見表1．在負(fù)荷水平比較低時(shí)，N-R法收斂較快，本文算法由于固定雅可比矩陣的處理而收斂次數(shù)較多但穩(wěn)定于10次左右，此時(shí)還未進(jìn)入病態(tài)；隨著負(fù)荷水平不斷提高，收斂越來越困難，在λ=2.0時(shí)N-R法已經(jīng)無法收斂，本文算法仍能夠求解，計(jì)算的殘差曲線如圖2所示；λ=2.1時(shí)本文算法也無法收斂．

表1 不同負(fù)荷水平下算法收斂性(迭代次數(shù))對比

圖2 潮流病態(tài)求解殘差曲線

以最優(yōu)乘子法計(jì)算相應(yīng)病態(tài)潮流的結(jié)果作為參考值，評(píng)估本文算法的電壓和相角誤差，進(jìn)一步驗(yàn)證算法準(zhǔn)確性．同樣添加變雅克比矩陣的測試結(jié)果，針對不同規(guī)模的病態(tài)潮流算例進(jìn)行驗(yàn)證，其結(jié)果見表2，表中誤差值為各個(gè)節(jié)點(diǎn)的電壓(標(biāo)幺值)和相角(弧度)與參考值作差的絕對值中的最大值．結(jié)果顯示，電壓和相角誤差均不超過5×10-3，在工程可以允許的范圍內(nèi)．

表2 病態(tài)潮流求解正確性

*注：采用最優(yōu)乘子法對case3120與case13659的相應(yīng)病態(tài)潮流進(jìn)行計(jì)算時(shí)無法收斂，因此以變雅克比法作為參考值．

為進(jìn)一步驗(yàn)證建立的微分求解模型的合理性，結(jié)合2.1節(jié)中相應(yīng)理論分析，針對不同規(guī)模的病態(tài)算例，計(jì)算如式(11)所示梯度矩陣在復(fù)平面上最靠右的特征值的實(shí)部a，結(jié)果見表3．可以看到，所有特征根均位于負(fù)半平面，微分方程組具有穩(wěn)定平衡點(diǎn)．

表3 梯度矩陣特征根分析

以上算例分析結(jié)果能夠驗(yàn)證本文算法應(yīng)用于病態(tài)潮流計(jì)算時(shí)的正確性．對于大規(guī)模系統(tǒng)而言，算法計(jì)算效率同樣需要關(guān)注．?dāng)?shù)值積分方法的步長對于常微分方程組求解的收斂速度有較大影響，對于本文方法同樣如此．在應(yīng)用算法時(shí)為達(dá)到效率的最優(yōu)化，應(yīng)選取使得迭代次數(shù)最少的積分步長，但并未有研究能夠?qū)Ψe分步長的選取進(jìn)行理論指導(dǎo)，因此本文對不同步長時(shí)算法的收斂速度進(jìn)行了實(shí)際測試與比較．求解時(shí)在Δt=0～1.5之間設(shè)置間隔為0.01的150個(gè)點(diǎn)，對Δt值在1.0附近波動(dòng)時(shí)的算法迭代次數(shù)進(jìn)行測試，結(jié)果如圖3所示．迭代次數(shù)上限設(shè)置為100次，達(dá)到該次數(shù)即視為不收斂．可以看出對于不同規(guī)模系統(tǒng)具有相似規(guī)律，在積分步長過小或過大時(shí)無法收斂，最優(yōu)積分步長在0.8～1.0之間．

圖3 不同積分步長時(shí)算法收斂速度

為進(jìn)一步提高整體效率，本文對于因子分解及前代回代計(jì)算等步驟均采用了Intel MKL[19]計(jì)算庫中的相應(yīng)函數(shù)，不同規(guī)模系統(tǒng)病態(tài)潮流計(jì)算的整體耗時(shí)結(jié)果見表4．可見算法對于上千節(jié)點(diǎn)大規(guī)模系統(tǒng)的病態(tài)潮流計(jì)算均有較高效率，對于超過萬節(jié)點(diǎn)的系統(tǒng)也在1s內(nèi)能夠給出病態(tài)潮流計(jì)算結(jié)果．

表4 算法整體效率

4 結(jié) 論

本文完成了基于微分求解模型的大規(guī)模電力系統(tǒng)病態(tài)潮流計(jì)算，主要結(jié)論如下：1)基于連續(xù)牛頓法，非線性方程組的求解可以轉(zhuǎn)化為常微分方程組的數(shù)值積分求解．2)建立電力潮流的微分求解模型，并進(jìn)行相應(yīng)簡化處理，使其適用于大規(guī)模電力系統(tǒng)的病態(tài)潮流求解．3)將數(shù)值積分方法應(yīng)用于求解中，得到整體積分計(jì)算式，并根據(jù)電力網(wǎng)絡(luò)特點(diǎn)構(gòu)建整體求解算法流程．4)結(jié)合具體系統(tǒng)算例進(jìn)行分析，通過對N-R法以及最優(yōu)乘子法潮流計(jì)算的結(jié)果對比，本文算法求解病態(tài)潮流時(shí)的有效性和準(zhǔn)確性得以證明；同時(shí)結(jié)果顯示算法對于大規(guī)模系統(tǒng)有較高效率，是可以應(yīng)用于實(shí)際場景的實(shí)用性算法．