王玥嬌，郭俊山，王 輝，田珊也，蒿天衢，張興友

（1.國網(wǎng)山東省電力公司電力科學(xué)研究院，山東 濟南 250003；2.山東大學(xué)，山東 濟南 250016）

0 引言

“碳達(dá)峰、碳中和”是我國政府對世界的莊嚴(yán)承諾。預(yù)計到2050 年，我國70%的電量將來自風(fēng)光發(fā)電。2020 年我國新能源新增裝機1.2 億kW，預(yù)計到“十四五”期間仍將以每年新增1.3 億kW 的速度發(fā)展［1］。分布式發(fā)電（Distributed Generations，DGs）的容量一般在幾十千瓦到數(shù)兆瓦之間，以接入35 kV及更低電壓等級配電網(wǎng)為主［2-4］?！笆奈濉逼陂g，山東規(guī)劃了500 萬kW 城鎮(zhèn)光伏發(fā)電、1 000 萬kW 鄉(xiāng)村分布式光伏發(fā)電和50 萬kW 分散式風(fēng)力發(fā)電項目。分布式發(fā)電因其一次能源的不確定性具有顯著的波動性和間歇性。多能互補的區(qū)域能源系統(tǒng)成為進(jìn)一步擴大新能源發(fā)電規(guī)模同時穩(wěn)定電網(wǎng)潮流的有效形式。在此背景下，不同形式的DGs 對配電網(wǎng)影響越來越顯著，為了有效規(guī)劃電源分布，研究和開發(fā)含不同類型DGs的新的配網(wǎng)潮流算法尤為重要［5］。

“隨機潮流”的概念由BORKOWSKA 在1974 年提出［6］，傳統(tǒng)的確定性潮流算法中，輸入變量和輸出變量皆為確定值，而隨機潮流（Probabilistic Load Flow，PLF）計算中，負(fù)荷、發(fā)電機的有功功率、無功功率等變量不再是確定值，而是服從一定概率分布的變化量，并由此計算出狀態(tài)變量和支路潮流的概率分布及所需的概率統(tǒng)計指標(biāo)。PLF 計算在實踐中依據(jù)不同的建模、求解方法，演化出多種研究方法。目前，PLF 算法大體可歸類為解析法和模擬法兩大類，根據(jù)采用的計算方法又可分為卷積法、蒙特卡洛模擬法、點估計法和半不變量法等［7-17］。

文獻(xiàn)［7］基于分布式發(fā)電的隨機性特征，通過Gram-Charlier 級數(shù)展開，分析輸出的隨機變量的概率分布。文獻(xiàn)［12］對微電網(wǎng)潮流進(jìn)行了預(yù)測，考慮了光伏發(fā)電和風(fēng)力發(fā)電的相關(guān)性，應(yīng)用聯(lián)合概率密度分布及馬爾科夫鏈和拉丁超立方抽樣方法進(jìn)行PLF 預(yù)測研究。文獻(xiàn)［13］和［14］將馬爾科夫鏈、蒙特卡洛法同時結(jié)合應(yīng)用，僅考慮到風(fēng)力發(fā)電的隨機性，對不同運行模式下微電網(wǎng)的潮流進(jìn)行分析預(yù)測。目前，國內(nèi)外學(xué)者對蒙特卡洛模擬法、點估計法［15-16］、一次二階矩法［17］和半不變量法等PLF 算法進(jìn)行了多方面的研究。

受地理微環(huán)境因素影響，同一配電網(wǎng)區(qū)域內(nèi)不同光伏電站的光照強度、不同風(fēng)場的風(fēng)速總體趨勢存在一致性，具體量測數(shù)值存在差異性，導(dǎo)致相關(guān)的光伏發(fā)電和風(fēng)力發(fā)電等分布式發(fā)電出力也具有相似性。針對同一區(qū)域內(nèi)分布式發(fā)電的相關(guān)性，提出一種基于三階多項式正態(tài)變換的計及多種分布式電源相關(guān)性的PLF 算法，改變了傳統(tǒng)PLF 只能利用獨立的新能源發(fā)電系統(tǒng)有功出力概率空間的限制，本方法將相關(guān)隨機性有功出力進(jìn)行轉(zhuǎn)換，使PLF方法應(yīng)用范圍得到推廣。將研究中具有相關(guān)性的考慮多種影響因素的多種新能源發(fā)電出力轉(zhuǎn)換為傳統(tǒng)PLF 可以利用的獨立非正態(tài)分布的有功出力概率，利用蒙特卡洛或三點估計等方法求出估計值，由確定性潮流進(jìn)行計算，進(jìn)一步通過Gram-Charlier 級數(shù)展開擬合出曲線。以高比例分布式發(fā)電配電網(wǎng)為例，仿真驗證考慮風(fēng)光發(fā)電相關(guān)性的潮流分布情況。

1 三階多項式正態(tài)變換

多項式正態(tài)變換方法（Polynomial Normal Transformation，PNT）于1978 年 由FLEISHMAN 提出［17］，可以用于求取非正態(tài)分布的新能源發(fā)電系統(tǒng)隨機有功出力。三階多項式正態(tài)變換（Thirdorder Polynomial Normal Transformation，TPNT）的原理及其在多維隨機變量解析中的應(yīng)用如下。

1.1 單個隨機變量TPNT

設(shè)非正態(tài)分布變量P，TPNT 是利用服從正態(tài)分布的隨機變量Q對P表示出來，表達(dá)式為

式中：P為隨機變量；Q為服從正態(tài)分布的隨機變量；a,b,c,d均為三階多項式系數(shù)。

式（1）成立需要滿足式（2）和式（3）的條件。

各系數(shù)a,b,c,d的值可由矩法（Product Moment，PM）、Fisher-Conish（FC）級數(shù)法或L-Moment Estimation（LM）法求得［18］。

1.2 多維隨機變量三階PNT

對于具有非正態(tài)分布形式的多維隨機變量矩陣X，X的相關(guān)系數(shù)矩陣RX=［ρXij］m×m在被求得后，進(jìn)一步處理得到符合正態(tài)分布的隨機變量矩陣Z的相關(guān)系數(shù)矩陣RZ=［ρZij］m×m。其中，ρXij和ρZij為系數(shù)矩陣對應(yīng)元素，存在函數(shù)關(guān)系，且滿足條件為：

需要將滿足正態(tài)分布的隨機變量矩陣Z的相關(guān)系數(shù)矩陣RZ去相關(guān)性。設(shè)服從標(biāo)準(zhǔn)正態(tài)分布且獨立的多維隨機變量矩陣S，對矩陣RZ進(jìn)行Cholesky分解，得出下三角矩陣L，其關(guān)系式滿足：

由式（8）將多維相關(guān)性隨機變量矩陣X進(jìn)一步變換為服從正態(tài)分布隨機變量矩陣S，其中A、B、C、D為對應(yīng)系數(shù)矩陣。通過TPNT 去除了隨機變量的相關(guān)性，使輸入變量變?yōu)楠毩?、服從?biāo)準(zhǔn)正態(tài)分布的隨機變量。再經(jīng)過反變換，變換為去相關(guān)性的非正態(tài)分布隨機變量，為PLF 計算提供符合條件的隨機輸入變量。

2 Gram-Charlier級數(shù)展開

在PLF 計算結(jié)果的處理中，通過Gram-Charlier級數(shù)展開分析支路輸出功率、相應(yīng)節(jié)點電壓，描述其概率分布并繪制圖形。Gram-Charlier 級數(shù)展開可在函數(shù)f(x)表達(dá)式未知時，利用特定的數(shù)學(xué)計算流程處理，相應(yīng)可求出f(x)的前k階矩陣Mk［7］，即

式中：n(x)為連續(xù)正函數(shù)；Ck為第k項系數(shù)；m0(x)，m1(x)，···，mk(x)為n(x)的正交多項式，各多項式需要滿足的正交條件如式（10）所示。

式中：a=0，1，2，…，k；b=0，1，2，…，k。

f(x)每一項的系數(shù)Ck為

綜上，任一隨機變量x的概率密度函數(shù)f（x）都可以通過對應(yīng)隨機變量前k階矩陣的求解來獲取。實際工程應(yīng)用中，隨機變量的概率密度函數(shù)和累積分布函數(shù)通常是通過符合標(biāo)準(zhǔn)正態(tài)分布的概率密度函數(shù)及累積分布函數(shù)求得。因此，推廣至隨機變量矩陣X，求得其前k階矩陣，即可得到其概率密度的函數(shù)f（X）表達(dá)式。

設(shè)μX為隨機變量矩陣X的均值，標(biāo)準(zhǔn)差為σX，標(biāo)準(zhǔn)化后X′=(X-μX)/σX，概率密度函數(shù)f（X）和概率分布函數(shù)F（X）的表達(dá)式為：

式中：Φ(X′)為標(biāo)準(zhǔn)正態(tài)分布函數(shù)；φ(X′)為標(biāo)準(zhǔn)正態(tài)分布概率密度函數(shù)。

在PLF 計算中，常用到前六項系數(shù)，如式（14）所示，σ為標(biāo)準(zhǔn)差、β為中心矩。

3 計及分布式發(fā)電相關(guān)性的配電網(wǎng)潮流計算

地理位置相近的分散式風(fēng)力發(fā)電或分布式光伏發(fā)電，由于風(fēng)速、光照強度之間存在一定的關(guān)聯(lián)，因此其出力也存在相關(guān)性，其相關(guān)系數(shù)矩陣可通過歷史數(shù)據(jù)進(jìn)行擬合［19］?；诙囗検秸龖B(tài)變換處理含相關(guān)性非正態(tài)分布隨機變量的概率問題可參考文獻(xiàn)［20-21］。通過TPNT 變換，同一配電網(wǎng)區(qū)域內(nèi)各風(fēng)力發(fā)電、光伏發(fā)電之間的相關(guān)性可以去除。基本步驟為，將具有相關(guān)性的光伏發(fā)電和風(fēng)力發(fā)電隨機出力進(jìn)行相關(guān)性處理，由PNT 變換到獨立正態(tài)空間。然后對獨立正態(tài)空間采樣，依據(jù)相對應(yīng)的系數(shù)矩陣映射關(guān)系進(jìn)行逆變換，返回原變量空間，最終得到去相關(guān)性的考慮多種影響因素的非正態(tài)隨機發(fā)電出力以應(yīng)用于PLF計算中。

基于以上計算結(jié)果，可利用蒙特卡洛、三點估計等方法求取各新能源節(jié)點功率的估計值，然后應(yīng)用確定性潮流算法，計算得到支路功率、節(jié)點電壓的點估計值，并進(jìn)一步計算其中心矩，利用Gram-Charlier級數(shù)展開繪制出分布曲線，從而實現(xiàn)計及具有隨機性的分布式發(fā)電相關(guān)性的配電網(wǎng)潮流計算。具體步驟如圖1所示。

圖1 考慮分布式發(fā)電隨機相關(guān)性的配電網(wǎng)PLF計算方法

常見DGs 有風(fēng)力發(fā)電、光伏發(fā)電、燃料電池和儲能系統(tǒng)。在分析時，須注意其不同的出力概率特性及因其接入電網(wǎng)的方式不同表現(xiàn)出的節(jié)點特性。

光伏發(fā)電系統(tǒng)的有功出力，如果不考慮棄光，直接由光伏發(fā)電系統(tǒng)所在地區(qū)的光照強度決定，因此其概率密度為光照強度的概率密度乘以代表其光電轉(zhuǎn)換效率的系數(shù)。分析數(shù)個地區(qū)歷史光照的數(shù)據(jù)發(fā)現(xiàn)，光照強度的概率密度曲線大致呈現(xiàn)貝塔分布［20］。光伏發(fā)電系統(tǒng)所發(fā)電力需要經(jīng)電壓源型逆變器轉(zhuǎn)化為交流電能，在進(jìn)行電力系統(tǒng)分析時通常把其并網(wǎng)點作為PV節(jié)點考慮。

風(fēng)力發(fā)電系統(tǒng)輸出的有功出力與風(fēng)速直接相關(guān)。但與光伏發(fā)電系統(tǒng)出力不同，光伏發(fā)電系統(tǒng)幾乎可以利用所有光照強度下的一次能源產(chǎn)生電能，而風(fēng)力發(fā)電機組無法利用所有風(fēng)速下的風(fēng)能，小于切入風(fēng)速的風(fēng)能因摩擦損耗、電機啟動等原因無法被利用；大于切出風(fēng)速的風(fēng)能受風(fēng)力發(fā)電機組機械應(yīng)力與機組電氣容量的限制，無法全部利用，當(dāng)風(fēng)速大于額定風(fēng)速、小于切出風(fēng)速時，風(fēng)力發(fā)電機組出力不變。在以上因素的影響下，風(fēng)力發(fā)電機組的出力概率曲線可由威布爾分布曲線進(jìn)行描述［21］。在進(jìn)行電力系統(tǒng)分析時，風(fēng)力發(fā)電機組的并網(wǎng)點通常作為PQ節(jié)點考慮。

4 仿真測試與結(jié)果

研究采用MATLAB 進(jìn)行編程并進(jìn)行仿真測試。首先建立仿真模型，模型基于一個含有多種DGs 的配電網(wǎng)系統(tǒng)，各支路中分別包含多種DGs 和負(fù)載的組合形式，各DGs 所發(fā)電能經(jīng)電力電子變換器以交流形式并網(wǎng)，電壓等級為380 V，各線路長度以及線路連接示意如圖2 所示。在提取到獨立非正態(tài)分布新能源電源出力估計值后，確定性潮流算法應(yīng)用牛頓-拉夫遜法。

圖2 含分布式發(fā)電的配電網(wǎng)

仿真中，由于3 個風(fēng)力發(fā)電單元DG1、DG2 和DG3 的空間距離較近，可知其風(fēng)速具有相關(guān)性，DG6、DG7 同為光伏電站，與風(fēng)力發(fā)電單元出力之間的關(guān)系類似，同理可知其光照強度也具有相關(guān)性。相關(guān)性系數(shù)處理上均采用LM 法。根據(jù)相關(guān)研究結(jié)果，光強概率密度的曲線形狀與貝塔分布曲線最為接近，光照強度概率分布曲線的形狀與威布爾分布最為接近，通過調(diào)節(jié)貝塔分布和威布爾分布中決定其最終形狀的參數(shù)，可建立隨機的光照強度和風(fēng)速模型［20-21］?；谝陨涎芯拷Y(jié)論，本研究開展的仿真驗證中，光伏出力和風(fēng)力發(fā)電機組出力分別服從貝塔分布和威布爾分布。

3 個風(fēng)場風(fēng)速相關(guān)系數(shù)以式（15）中矩陣R11表示，轉(zhuǎn)換為服從正態(tài)分布空間的相關(guān)系數(shù)矩陣可得R21，如式（16）所示；同理，兩個光伏電站的光照強度相關(guān)系數(shù)矩陣R12，如式（17）所示，可轉(zhuǎn)換為正態(tài)分布空間的相關(guān)系數(shù)矩陣R22，如式（18）所示。

光照強度和風(fēng)速的三階多項式系數(shù)如表1所示。表2 給出了含相關(guān)性的風(fēng)光發(fā)電出力處理后的數(shù)學(xué)期望。

表1 光照強度和風(fēng)速三階多項式系數(shù)

表2 含相關(guān)性輸入變量的各節(jié)點有功功率期望值 單位：pu

以電網(wǎng)末端的節(jié)點6和節(jié)點7為例，分析配電網(wǎng)的電壓質(zhì)量。節(jié)點6 具有光伏發(fā)電和風(fēng)力發(fā)電的相關(guān)性，圖3 顯示了其電壓概率密度曲線，電壓波動范圍為0.93～1.02 pu，波動值為0.09 pu。節(jié)點7 的電壓概率密度曲線如圖4所示，電壓范圍為0.97～1.02 pu，波動值為0.05 pu。分析可知因節(jié)點7 中接入了燃料電池發(fā)電機組，其有功出力可控性高，對比節(jié)點6 的電壓情況，其電壓質(zhì)量也較高。

圖3 節(jié)點6電壓概率密度曲線

圖4 節(jié)點7電壓概率密度曲線

圖5、圖6 給出了經(jīng)TPNT 變換后，利用三點估計法和蒙特卡洛模擬法得到的節(jié)點6、節(jié)點7 電壓概率分布曲線對比圖。由曲線可以看出，考慮輸入隨機變量的相關(guān)性后，TPNT 變換得到的獨立非正態(tài)隨機變量可以很好地應(yīng)用到潮流算法中，在三點估計法和蒙特卡洛模擬法中都可有效應(yīng)用，針對具備相關(guān)性輸入隨機變量的配電網(wǎng)PLF 是可行有效的，可以推測出考慮分布式發(fā)電區(qū)域相關(guān)性更貼近電網(wǎng)實際運行情況。

圖5 節(jié)點6電壓概率分布曲線

圖6 節(jié)點7電壓概率分布曲線

圖7、圖8給出了利用Gram-Charlier 級數(shù)展開后繪制的支路7—1、7—8 潮流概率密度曲線。從圖7可以看出，由于DGs 出力隨機性、波動性的影響，支路7—1的功率流向大多數(shù)情況下由節(jié)點7到節(jié)點1，但當(dāng)DGs 出力不足時，有功潮流可能也會反向變?yōu)橛晒?jié)點1 到節(jié)點7。圖8 說明支路7—8 有功功率主要在小于0.35 pu 的范圍之內(nèi)波動，功率流向為從節(jié)點7流向節(jié)點8。

圖7 支路7—1潮流概率密度曲線

圖8 支路7—8潮流概率密度曲線

5 結(jié)語

針對同一區(qū)域配電網(wǎng)中不同分布式發(fā)電單元具有的隨機性、波動性以及相關(guān)性特點，提出一種考慮分布式發(fā)電單元出力相關(guān)性的變換算法，并將其應(yīng)用到隨機潮流計算。該算法將具有相關(guān)性的不同新能源出力利用TPNT 變換為考慮多因素影響的獨立非正態(tài)隨機出力，進(jìn)一步通過隨機潮流的計算，獲得相應(yīng)的概率統(tǒng)計特征。仿真結(jié)果表明，所提出的潮流算法靈活、簡潔、應(yīng)用范圍廣，且具有較好的準(zhǔn)確度，可應(yīng)用于隨機潮流計算的數(shù)據(jù)處理中。經(jīng)分析驗證，所提方法能夠合理反映高比例分布式發(fā)電配電網(wǎng)的穩(wěn)態(tài)運行特性，對高比例分布式發(fā)電配電網(wǎng)運行和規(guī)劃有一定的指導(dǎo)意義。