基于兩點估計法的有源配電網(wǎng)概率潮流計算方法研究*

柳永妍，樂健，高鵬，劉開培

（1.武漢大學(xué) 電氣工程學(xué)院，武漢430072；2.國家電網(wǎng)北京市電力公司，北京100031）

0 引 言

為滿足日益增長的電力需求，高滲透率分布式電源接入配電網(wǎng)運行將是今后配電網(wǎng)發(fā)展的必然趨勢。分布式電源的接入使得原有的輻射型配電網(wǎng)變?yōu)楹卸嚯娫吹木W(wǎng)絡(luò)結(jié)構(gòu)，潮流也從單一流向轉(zhuǎn)變?yōu)榭赡艹霈F(xiàn)回流的情況，為評估配電網(wǎng)中各種不確定因素對節(jié)點電壓、線路網(wǎng)損和諧波等問題的影響，研究相應(yīng)的概率潮流（Probabilistic Load Flow，PLF）計算方法，考慮配電網(wǎng)中的不確定因素，得出配電網(wǎng)狀態(tài)量的統(tǒng)計性結(jié)果將具有重要的實用價值。

現(xiàn)有配電網(wǎng)絡(luò)主要為輻射狀，前推后代法［1-3］因其不存在矩陣計算且收斂性好而得到廣泛的應(yīng)用。其它潮流計算方法包括有回路阻抗法［4］、Z bus法［2，5］、改進(jìn)快速解耦法［6］等。

文獻(xiàn)［6］提出一種基于逆流拓?fù)浞謱铀阉鞯慕怦畛绷魉惴?，對含有多種網(wǎng)絡(luò)結(jié)構(gòu)的復(fù)雜電網(wǎng)進(jìn)行計算。利用逆流拓?fù)浞謱铀阉鞯姆椒▽?fù)雜配電網(wǎng)分解成輻射網(wǎng)和主環(huán)網(wǎng)兩個部分，再分別針對兩部分求解，這樣的處理辦法在很大程度上提高了潮流計算速度、收斂性以及適應(yīng)性，特別是在處理環(huán)網(wǎng)率較低的網(wǎng)絡(luò)結(jié)構(gòu)時，相較單一使用PQ分解法，分層解耦潮流計算方法所需時間要少。

傳統(tǒng)潮流計算方法的優(yōu)化始終都是針對配電網(wǎng)的一組數(shù)據(jù)進(jìn)行潮流分析，分布式電源的接入，需要利用統(tǒng)計性的潮流計算結(jié)果來針對低電壓治理等一系列的問題進(jìn)行分析。概率潮流（隨機(jī)潮流）的計算方法，其注入變量不再表示成單一量的形式，而是針對這一變量取多組數(shù)據(jù)，通過將概率密度函數(shù)或者變量分布所具有的數(shù)字特征引入潮流計算。將概率論與數(shù)理統(tǒng)計的知識與潮流計算方法相結(jié)合，綜合考慮配電網(wǎng)中狀態(tài)量的波動，得到潮流計算結(jié)果的統(tǒng)計特性。

現(xiàn)有的概率潮流計算方法主要分為模擬法和解析法。模擬法中最為主要的是蒙特卡羅（Monte Carlo）法［7］，蒙特卡羅模擬法是分析概率問題最常用的數(shù)學(xué)方法，與不同的采樣方法相結(jié)合可以提高其計算結(jié)果的精確度。但理論與實踐證明，使用蒙特卡羅法需要大量的采樣數(shù)據(jù)才能獲得較高的精度，并且對每一次采樣的數(shù)據(jù)進(jìn)行一次確定性潮流計算，所花費的時間長，實用性不佳。

但當(dāng)該方法抽樣規(guī)模足夠大時，其結(jié)果的準(zhǔn)確性和精確度即可得到很好的保證，因此常常作為其它概率潮流計算方法的驗證標(biāo)準(zhǔn)。文獻(xiàn)［8-9］提出一種能處理變量相關(guān)性的基于拉丁超立方（Latin Hypercube Sampling，LHS）采樣的蒙特卡羅概率潮流計算方法，在盡可能保證精度的同時減少計算所需的時間。同樣的，也可以將蒙特卡羅概率潮流計算方法與其它采樣方法相結(jié)合，提高計算效率。

解析法則分為很多種，其分類依據(jù)是計算采用節(jié)點注入量的數(shù)字特征不同，主要包括有快速傅里葉變換法 （FFTM）、半不變量法［10］、點估計法（PEM）［11］等。

提出了基于兩點估計法的配電網(wǎng)概率潮流計算方法，根據(jù)負(fù)荷與分布式電源的概率分布模型，在選取較少采樣計算點的情況下，獲得與蒙特卡羅潮流計算方法同等精度的潮流計算結(jié)果。提高了配電網(wǎng)含波動性負(fù)荷潮流計算結(jié)果的準(zhǔn)確性和可靠性。

1 負(fù)荷概率模型

本概率潮流計算中需要用統(tǒng)計特性來表征配電系統(tǒng)的相關(guān)狀態(tài)量，因此需要對配電網(wǎng)中具有不確定性對象之一的負(fù)荷建立概率模型。

1.1 一般負(fù)荷概率模型

目前認(rèn)為負(fù)荷預(yù)測所得結(jié)果為正態(tài)分布形式，分別用 μP和表示有功功率的期望值和方差，用μQ和表示無功功率的期望值和方差，則負(fù)荷有功和無功功率的概率密度函數(shù)［12］為：

1.2 光伏發(fā)電概率模型

一般認(rèn)為太陽能輻照度變化曲線［13］近似于正態(tài)分布，即：

式中 r為太陽輻照度（W／m2），服從 N（μ，σ2）分布。考慮到干擾因素，可認(rèn)為一定時間段內(nèi)的太陽輻照度近似服從Beta分布，當(dāng)Beta分布取值區(qū)間為［0，1］時，即為標(biāo)準(zhǔn)Beta分布。其概率密度函數(shù)為：

式中 rmax（W／m2）為該時段太陽最大輻照度，α，β為Beta分布的形狀參數(shù)。

由于光伏發(fā)電系統(tǒng)一般僅向電網(wǎng)提供有功功率，對于無功功率不予考慮。對于一個包含有M個電池組件的太陽能陣列，其總的輸出功率為：

式中 A為光伏陣列的總面積，η0為陣列的光電轉(zhuǎn)換效率，CEE為連接效率，IE為逆變器轉(zhuǎn)換效率。設(shè)定η＝η0·CEE·IE，可將上式簡化為：

根據(jù)連續(xù)型隨機(jī)變量概率分布的定理，對于連續(xù)型隨機(jī)變量 X，設(shè) fx（x）為 X的概率密度，若 y＝g（x）是嚴(yán)格單調(diào)的連續(xù)函數(shù)，且反函數(shù) x＝h（y）有連續(xù)導(dǎo)數(shù)，則Y＝g（X）為連續(xù)型隨機(jī)變量，且概率密度為：

式中α＝min（g（-∞），g（∞））；β＝max（g（-∞），g（∞））。設(shè)定 X＝r／rmax，Y＝PSOLAR，代入式（6）可得：

式中RSOLAR＝rmaxAη為太陽能陣列的最大輸出功率。

1.3 風(fēng)力發(fā)電概率模型

利用威布爾分布雙參數(shù)曲線［14］來對風(fēng)速進(jìn)行統(tǒng)計描述。假設(shè)某風(fēng)電場的風(fēng)速序列（v1，v2，…，vN）服從兩參數(shù)威布爾分布，其概率密度函數(shù)和累計概率函數(shù)可分別表示為：

式中v為風(fēng)速；K和C分別為兩參數(shù)威布爾分布的形狀參數(shù)和尺度參數(shù)。通過改變形狀參數(shù)K，可以使分布曲線變?yōu)椴煌男问?。?dāng)0＜K＜1時，分布的眾數(shù)為0，分布密度為x的函數(shù)；當(dāng)K＝1時，分布曲線呈指數(shù)型；當(dāng)K＝2時，分布即為瑞利分布；當(dāng)K＝3.5時，威布爾分布已經(jīng)逼近于正態(tài)分布了。

兩參數(shù)威布爾分布的均值和標(biāo)準(zhǔn)差［15］可以寫為：

式中μ為平均風(fēng)速；σ為采樣風(fēng)速的標(biāo)準(zhǔn)差；vi為各采樣時刻的實際風(fēng)速；fi為vi對應(yīng)風(fēng)速出現(xiàn)的概率。兩參數(shù)威布爾分布的尺度參數(shù)C和形狀參數(shù)K，可以采用平均風(fēng)速和標(biāo)準(zhǔn)差估計法獲得：

式中μ為平均風(fēng)速；σ為采樣風(fēng)速的方差。

風(fēng)電場的出力是風(fēng)速的函數(shù)［16］，可將其表示為：

式中Pr、vr為風(fēng)輪機(jī)的額定功率和額定風(fēng)速；vci、vco為風(fēng)輪機(jī)的切入、切除風(fēng)速。a、b可表示為：

如果計及尾流效應(yīng)，即考慮風(fēng)向、地形地貌、機(jī)組間距等相關(guān)因素的影響，造成的典型損失值為10%，可以將該典型系數(shù)歸入Pr之中。經(jīng)統(tǒng)計，風(fēng)速一般處于vci＜v＜vr這個區(qū)間范圍內(nèi)，單一風(fēng)機(jī)組的輸出功率是風(fēng)速的函數(shù)，結(jié)合概率論與數(shù)理統(tǒng)計中的定理公式得：

式中Pwind為單臺風(fēng)電機(jī)組輸出功率；K為形狀參數(shù)；C為尺度參數(shù)。

風(fēng)機(jī)通常不會單獨的連入電網(wǎng)，而是以風(fēng)電場的形式輸出功率，假設(shè)風(fēng)電場包含有N臺風(fēng)機(jī)，Pfarm＝N·Pwind，則：

風(fēng)電場的建立，要求具備一定的無功補(bǔ)償容量，因此可近似認(rèn)為風(fēng)電場的功率因數(shù)cosφ保持不變，故有：

式中φ為風(fēng)電場功率因數(shù)角；Pfarm和Qfarm分別為風(fēng)電場有功和無功功率。

2 基于前推回代法的概率潮流計算方法

2.1 蒙特卡洛算法

蒙特卡羅模擬法（MCS）的基本思想是為了求解數(shù)學(xué)、物理工程技術(shù)以及生產(chǎn)管理方面的問題，對待測模型進(jìn)行抽樣實驗來計算待求參數(shù)的統(tǒng)計特征，最后求得待求量的近似值，所得解的精確度可以用估計值的標(biāo)準(zhǔn)誤差來表示。

蒙特卡羅模擬法可以精確地獲得狀態(tài)電壓和潮流的概率描述并且計算模式簡單易行，然而為了獲得精確的結(jié)果需要成千上萬次的模擬，需要較長的計算時間，故通常作為評價各種算法優(yōu)劣的標(biāo)準(zhǔn)。

文中采用的蒙特卡羅模擬法的計算步驟為：

（1）首先根據(jù)配電網(wǎng)的一組負(fù)荷參數(shù)，假設(shè)其服從正態(tài)分布，隨機(jī)產(chǎn)生若干組數(shù)據(jù)；

（2）對于含有分布式電源的節(jié)點，根據(jù)分布式電源的概率分布特性，建立服從該分布的若干組隨機(jī)數(shù)；

（3）利用各組隨機(jī)數(shù)據(jù)分別進(jìn)行確定性潮流計算，得到多組各節(jié)點電壓和幅值、支路有功功率和無功功率的計算結(jié)果；

（4）對所得結(jié)果進(jìn)行統(tǒng)計分析，計算輸出變量的數(shù)學(xué)期望和概率分布情況。

2.2 2n點估計算法

兩點估計法（2PEM）通過在每個不確定變量的均值兩側(cè)確定兩個值，將對節(jié)點注入量的求解和支路功率的求解問題分成若干個子問題，在每一個不確定量處用均值兩側(cè)的值來代替，同時其他不確定量取均值。由此可見，若系統(tǒng)含有k個不確定量，那么需要對該系統(tǒng)進(jìn)行2k次確定性潮流計算。

將節(jié)點注入向量表示為（m個隨機(jī)變量）：

支路潮流表示為節(jié)點注入量的函數(shù)：

兩點估計法通過使用均值兩側(cè)的兩個變量xi1和xi2來匹配隨機(jī)量xi的前三階矩，以取代注入量的概率密度函數(shù)。將xi1和xi2定義為：

式中aver_xi和sd_xi分別表示xi的均值和標(biāo)準(zhǔn)差，spik為位置度量，其表達(dá)式如下：

式中 λi，3為 xi的偏度系數(shù)，E［（xi-aver_xi）3］為 xi的三階中心矩。某一節(jié)點注入隨機(jī)量xi用xi1、xi2代替，其它隨機(jī)量均設(shè)為均值，用2×m組數(shù)據(jù)來進(jìn)行確定性潮流計算：

得到的結(jié)果包括有節(jié)點電壓的統(tǒng)計值和支路潮流的兩個估計 CLr（i，1）和 CLr（i，2）（r＝1，2，…b），其中b為支路數(shù)。用wi，k表示xik的概率集中度，即：

式中xik處位置集中的權(quán)重：

兩點估計法由多組輸入變量的值計算得出其數(shù)字特征，再利用其數(shù)字特征確定估計點，對各組估計數(shù)據(jù)進(jìn)行確定性潮流計算，近似的得到輸出變量的統(tǒng)計分布特性，計算速度較快，但輸出變量的高階矩可能誤差會相對較大。

在沒有給定的概率模型時，將不確定量的個數(shù)選為負(fù)荷節(jié)點的個數(shù)，又由于所討論的為輻射網(wǎng)模型，所以不確定量的個數(shù)為節(jié)點數(shù)N-1＝支路數(shù)b。

文中采用兩點估計法的算法步驟為：

（1）以IEEE標(biāo)準(zhǔn)輻射型配電網(wǎng)為模型，將其中的節(jié)點注入量的值作為節(jié)點的注入量的期望值，默認(rèn)節(jié)點注入量服從正態(tài)分布，自己設(shè)定方差值；

（2）由模型的數(shù)字特征產(chǎn)生服從該分布的隨機(jī)數(shù)；

（3）由這些隨機(jī)數(shù)求得兩點估計法需要的其他數(shù)字特征，例如偏度系數(shù)等；

（4）確定需要進(jìn)行確定性潮流計算的2b組節(jié)點注入量數(shù)據(jù)，然后利用輻射網(wǎng)前推回代潮流計算方法得出2b組計算結(jié)果；

（5）對計算結(jié)果進(jìn)行統(tǒng)計分析，得到支路功率和節(jié)點電壓的相關(guān)統(tǒng)計數(shù)據(jù)。

2.3 2n+1點估計法

2n+1點估計法在2n點估計法的基礎(chǔ)上對每一個隨機(jī)變量xi增加一個估計點，增加的估計點取為隨機(jī)變量的均值aver_xi，其對應(yīng)的位置度量spi3＝0，每個估計點的位置度量 spik，概率集中度 wi，k如下［17］：

式中 λi，4為 xi的峰度系數(shù)，E［（xi-aver_xi）3］為 xi的四階中心矩。這里構(gòu)造了3n個估計點，其中在2n點估計法基礎(chǔ)上增加的n個點是對應(yīng)同一估計點向量（avers_x1，aver_x2，…，aver_xi，…，aver_xm），因此認(rèn)為是2n+1方案。由于2n+1點估計法考慮到了隨機(jī)變量的峰度系數(shù)，因此認(rèn)為2n+1方案比2n方案的點估計精度高。

3 算例分析

以標(biāo)準(zhǔn)IEEE-33節(jié)點系統(tǒng)為對象，負(fù)荷節(jié)點注入量給定值作為期望值，將給定值的10%作為方差構(gòu)造一組服從正態(tài)分布的隨機(jī)數(shù)作為概率模型，假設(shè)在11號節(jié)點存在光伏發(fā)電設(shè)備，其輸入系統(tǒng)的有功功率服從均值為60 kW，方差為10%的正態(tài)分布；24號節(jié)點設(shè)有風(fēng)機(jī)，風(fēng)機(jī)型號為VESTAS RRB 225，其輸入系統(tǒng)的有功功率均值為225 kW，切入風(fēng)速為3.5 m／s、額定風(fēng)速為13.5 m／s、切出風(fēng)速為 25 m／s，因風(fēng)電場通常含有足夠的無功補(bǔ)償裝置，可將功率因數(shù)控制在一定值，本文設(shè)功率因數(shù)為1。分別采用蒙特卡羅模擬法、2n點估計法和2n+1點估計法進(jìn)行概率潮流計算。線路拓?fù)淙鐖D1所示，取系統(tǒng)基準(zhǔn)電壓為12.66 kV，基準(zhǔn)容量為10 kVA，線路參數(shù)及節(jié)點注入功率參數(shù)參考文獻(xiàn)［18］。

圖1 IEEE 33節(jié)點電路拓?fù)銯ig.1 IEEE33 node circuit topology

本文采用的風(fēng)速模型參考文獻(xiàn)［19］，求得尺度參數(shù) C＝7.451 9，形狀參數(shù) K＝3.422 5。繪制服從上述參數(shù)的威布爾分布風(fēng)速見圖2～圖3。

圖2 風(fēng)速概率密度分布Fig.2 Probability density distribution ofwind speed

圖3 風(fēng)速累積概率密度分布Fig.3 Cumulative probability density distribution of wind speed

通過MATLAB仿真得到各支路功率的均值和方差以及各節(jié)點電壓的期望值。

表1 三種算法所得部分支路有功功率均值Tab.1 Part of branch circuit active power average of three algorithms

表1～表4中的數(shù)據(jù)顯示，以蒙特卡羅所得結(jié)果為標(biāo)準(zhǔn)，2n點估計法和2n+1點估計法所得結(jié)果支路有功功率和無功功率的均值誤差在0.01%范圍內(nèi)，方差誤差在0.1×10-5，兩種點估計法所得結(jié)果差距不大。

表2 三種算法所得部分支路無功功率均值Tab.2 Part of branch circuit reactive power average of three algorithms

表3 三種算法所得部分支路有功功率方差Tab.3 Part of branch circuit active power variance of three algorithms

表4 三種算法所得部分支路無功功率方差Tab.4 Part of branch circuit reactive power variance of three algorithms

表5 三種算法所得部分節(jié)點電壓期望值Tab.5 Part of node voltage expectation of three algorithms

表5中的數(shù)據(jù)顯示，2n點估計法和2n+1點估計法與蒙特卡羅模擬法計算所得各節(jié)點電壓期望值相比，其差小于0.05%。

在該算例中，針對給定的負(fù)荷模型和分布式電源模型，蒙特卡羅法產(chǎn)生了10 000組隨機(jī)數(shù)，并對這10 000組隨機(jī)數(shù)分別進(jìn)行了潮流計算，統(tǒng)計其支路功率和節(jié)點電壓的結(jié)果；而兩點估計法在32個含有不確定輸入量節(jié)點上分別取兩個估計點，針對每一個節(jié)點所取的兩個估計點進(jìn)行潮流計算，共進(jìn)行32×2＝64次潮流計算，其結(jié)果與蒙特卡羅模擬法所得結(jié)果只存在很小的誤差。因此，針對含有分布式電源的配電網(wǎng)采用兩點估計概率潮流方法所得的計算結(jié)果有較高的準(zhǔn)確性和精確度。部分節(jié)點電壓分布函數(shù)如圖4所示。

圖4 部分節(jié)點電壓分布函數(shù)Fig.4 Part of node voltage distribution function

采用兩點估計法和蒙特卡羅法計算所得的各節(jié)點電壓分布函數(shù)基本相同，24號節(jié)點因接有服從Weibull分布的風(fēng)電模型而使得電壓分布函數(shù)存在較大的區(qū)別，但是總的來說其走向基本相同，兩點估計法的計算結(jié)果精度高。三種算法計算耗時見表6。

表6 3種算法耗時比較Tab.6 Elapsed time comparison between three algorithms

由表6中的數(shù)據(jù)可以看出2n點估計法的計算時間約為蒙特卡羅法的1.8%，而2n+1點估計法的計算時間約為蒙特卡羅的2.3%，在計算效率上2n點估計法明顯優(yōu)于2n+1點估計法。通過對配電網(wǎng)上接有的負(fù)荷和分布式電源模型數(shù)據(jù)進(jìn)行處理，確定其估計點，可以大大減少采樣計算點的數(shù)量，提高概率潮流的計算效率。

4 結(jié)束語

文章對含有風(fēng)電、光伏和一般負(fù)荷的配電網(wǎng)進(jìn)行了基于前推回代法的概率潮流計算。依據(jù)風(fēng)速和太陽輻照度的統(tǒng)計特性，分別建立服從兩參數(shù)Weibull分布和正態(tài)分布的模型，用以模擬風(fēng)速和太陽輻照度的概率特性。利用兩種點估計法，在選取較少計算點的情況下對含分布式電源的配電網(wǎng)進(jìn)行概率潮流計算。對比兩種方法所得各節(jié)點電壓、各支路有功和無功功率的期望值和方差。仿真結(jié)果表明兩種點估計概率潮流計算方法所得結(jié)果具有較好的準(zhǔn)確性和收斂性，兩種點估計方法所得結(jié)果基本相同，但是2n法的計算速度優(yōu)于2n+1點估計法。文章采用的基于前推回代的2n點估計概率潮流計算方法計算速度遠(yuǎn)高于蒙特羅模擬法，更適用于具有隨機(jī)性、波動性和多樣性的含分布式電源的配電網(wǎng)。