黃強，張劉冬，陳兵

(國網(wǎng)江蘇省電力有限公司電力科學(xué)研究院，南京 211103)

0 引 言

對于含大規(guī)模風(fēng)電場并網(wǎng)的電力系統(tǒng)，風(fēng)電場出力和負(fù)荷的變化不僅存在不確定性，而且分別具有相關(guān)性。忽視這兩類相關(guān)性，可能會導(dǎo)致系統(tǒng)潮流概率分布出現(xiàn)較大誤差，從而造成系統(tǒng)規(guī)劃和運行風(fēng)險評估偏離實際。因此，概率潮流計算需要具備同時處理負(fù)荷相關(guān)性以及風(fēng)速相關(guān)性的能力[1-9]。

針對同時處理風(fēng)速相關(guān)性和負(fù)荷相關(guān)性的概率潮流算法，現(xiàn)有研究大致可歸納為三類方法：

(1)蒙特卡羅法[2-4]；(2)點估計法[5-6]；(3)累積量法(半不變量法)[7-9]。其中，文獻[2-4]在概率潮流計算中，采用蒙特卡羅法抽樣生成大量具有相關(guān)性的風(fēng)速及負(fù)荷樣本。但該方法缺點在于為提高計算精度，往往需要提高系統(tǒng)抽樣規(guī)模，從而導(dǎo)致計算時間過長，一般用來作為基準(zhǔn)方法進行比較，衡量其他方法的準(zhǔn)確性。點估計法是一種近似方法，一般利用輸入隨機變量的統(tǒng)計信息來逼近系統(tǒng)狀態(tài)變量的數(shù)字特征，即可根據(jù)已知隨機變量的概率分布，求得待求隨機變量的各階矩。該方法避開了大規(guī)模的重復(fù)抽樣，因而求解速度較快。文獻[5-6]在點估計法中，使用矩陣變換的方法處理風(fēng)速相關(guān)性和負(fù)荷相關(guān)性，計算速度很快，但獲得支路潮流的概率密度函數(shù)不夠精確，一般只用于求取其均值和方差。累積量法又稱半不變量法，其思想是只通過解析計算就能求出系統(tǒng)狀態(tài)變量的概率分布。該方法要求系統(tǒng)狀態(tài)變量是輸入變量的線性表達式，利用累積量的特性，只通過簡單的代數(shù)運算就可以求得狀態(tài)變量的各階累積量，并結(jié)合Gram-Charlier級數(shù)獲得支路潮流的概率密度函數(shù)，大大地節(jié)省了計算時間。文獻[7-9]在線性化的交流潮流模型中，利用采樣生成風(fēng)速和負(fù)荷相關(guān)性樣本的方法計算累積量，并使用累積量法計算支路潮流的概率分布。目前該算法的研究重點是非正態(tài)分布的變量相關(guān)性處理，改善各階矩和累積量的計算精度與計算效率，獲得更為準(zhǔn)確的概率密度函數(shù)。

上述同時處理負(fù)荷相關(guān)性和風(fēng)速相關(guān)性的概率潮流計算方法，實質(zhì)上忽視了兩種隨機變量概率特性的差異，采用相同的方法處理具有不同概率分布特性的隨機變量的相關(guān)性，沒有充分利用一般呈正態(tài)分布的負(fù)荷在計算相關(guān)性時具有的便利。文獻[10]利用負(fù)荷的正態(tài)分布特性，首次采用解析法快速求取具有相關(guān)性的多個負(fù)荷對應(yīng)總負(fù)荷的正態(tài)分布函數(shù)，并應(yīng)用了累積量法計算支路潮流的概率密度函數(shù)。文獻[11]首次采用卷積法求解概率潮流問題，并考慮了負(fù)荷相關(guān)性為線性的情況。文獻[12]將蒙特卡羅法與卷積法相結(jié)合，采用蒙特卡羅法模擬出風(fēng)電場出力的概率密度函數(shù)，然后代入解析的直流概率潮流計算中，結(jié)果具有足夠精度的情況下盡量減少計算時間。

為此，提出基于分類處理風(fēng)速相關(guān)性和負(fù)荷相關(guān)性的直流概率潮流算法。它延用蒙特卡羅法計算相關(guān)風(fēng)電場總出力的概率密度曲線；利用負(fù)荷一般呈正態(tài)分布的特性，采用解析法快速求得總負(fù)荷的正態(tài)分布函數(shù)；最后采用卷積法將兩類概率計算結(jié)果進行綜合。由于在處理具有相關(guān)性負(fù)荷時避免了負(fù)荷樣本的生成與采樣，僅需要進行一次卷積計算，可以大幅提高概率潮流的計算效率。文章對含多風(fēng)電場的IEEE RTS-96系統(tǒng)和IEEE 108系統(tǒng)進行了仿真分析，并與蒙特卡羅法、點估計法及累積量法進行比較，驗證了本文方法的有效性和優(yōu)越性。

1 直流潮流模型

1.1 基于直流潮流的線路有功功率線性表達式

采用直流潮流計算模型或者線性化的交流潮流計算模型。由直流潮流模型[12-13]可得：

(1)

式中B為節(jié)點導(dǎo)納矩陣的虛部矩陣；B′為B去掉平衡節(jié)點后的矩陣，所以B′存在逆矩陣；Pij為線路i-j的有功功率；δij為線路i-j首末節(jié)點的相角差，δ為相角差矩陣；xij為線路i-j的電抗；PN為節(jié)點注入功率矩陣。

綜合式(1)中的3個等式，可以獲得線路i-j的有功功率線性表達式：

(2)

式中αn為節(jié)點n對線路i-j的轉(zhuǎn)移分布因子；Pn為節(jié)點n的注入功率。節(jié)點注入功率可以分為發(fā)電機組出力和負(fù)荷需求2部分，其中發(fā)電機組又包括傳統(tǒng)機組和風(fēng)電機組。因此，式(2)可以分為3部分，即:

(3)

式中M為傳統(tǒng)機組數(shù)目；Pgm為傳統(tǒng)機組m的出力；K為風(fēng)電場數(shù)目；Pwk為風(fēng)電場k的出力；L為負(fù)荷數(shù)目；Pdl為負(fù)荷l的大小。

1.2 風(fēng)電場出力模型

(4)

式中系統(tǒng)接入的K個風(fēng)電場中包含n個具有相關(guān)性的風(fēng)電場組，每組風(fēng)電場分別有K1,K2, …,Kn個風(fēng)電場，風(fēng)電場組之間相互獨立，組內(nèi)風(fēng)電場之間具有相關(guān)性；其余(K-K1-K2-…-Kn)個風(fēng)電場為獨立風(fēng)電場。

1.2.1 風(fēng)電場出力的概率分布

風(fēng)速的概率分布一般采用威布爾分布模型來模擬[12-13]。

(5)

式中k為形狀參數(shù);c為尺度參數(shù)。

風(fēng)電機組的有功功率輸出特性[12-13]如下：

(6)

式中vci為切入風(fēng)速;vrate為額定風(fēng)速;vco為切出風(fēng)速；Prate為風(fēng)電機組額定出力。

由式(5)和(6)可以獲得風(fēng)電場出力的概率分布：

(7)

1.2.2 風(fēng)電場出力的相關(guān)性描述

風(fēng)電場出力相關(guān)性一般采用線性相關(guān)系數(shù)來描述，多個風(fēng)電場的出力相關(guān)性常用線性相關(guān)系數(shù)矩陣描述[14-16]。在式(4)中，第m個相關(guān)風(fēng)電場組包含Km個具有相關(guān)性的風(fēng)電場，其相關(guān)系數(shù)矩陣Cm如式(8)所示：

(8)

式中ρij為風(fēng)電場i和風(fēng)電場j的線性相關(guān)系數(shù)，i,j=1,2,…,Km；Vi和Vj為對應(yīng)風(fēng)電場的風(fēng)速隨機變量；σi和σj為對應(yīng)風(fēng)速的標(biāo)準(zhǔn)差。

1.3 負(fù)荷模型

(9)

式中系統(tǒng)L個負(fù)荷中包含m個具有相關(guān)性的負(fù)荷組，每組負(fù)荷分別有L1,L2, …,Lm個負(fù)荷，負(fù)荷組之間相互獨立，組內(nèi)負(fù)荷之間具有相關(guān)性；其余(L-L1-L2-…-Ln)個負(fù)荷為獨立負(fù)荷。

1.3.1 負(fù)荷的概率分布

一般認(rèn)為負(fù)荷是服從正態(tài)分布的[2-12]，其概率密度函數(shù)如下：

(10)

式中μ是均值;σ2是方差。

1.3.2 負(fù)荷的相關(guān)性描述

負(fù)荷是服從正態(tài)分布的，正態(tài)分布變量的相關(guān)性可以用線性相關(guān)系數(shù)描述，因此多個負(fù)荷的相關(guān)性用線性相關(guān)系數(shù)矩陣描述。在式(9)中，第k個相關(guān)負(fù)荷組包含Lk個具有相關(guān)性的負(fù)荷，其相關(guān)系數(shù)矩陣Ck如式(11)所示。

(11)

2 基于卷積的分類處理風(fēng)速相關(guān)性和負(fù)荷相關(guān)性的直流概率潮流算法

2.1 直流概率潮流算法的流程

鑒于風(fēng)速和負(fù)荷的不同概率分布特性[2-12]，本文提出了基于卷積的分類考慮負(fù)荷相關(guān)性和風(fēng)速相關(guān)性的概率潮流計算方法。算法流程圖如圖1所示。

具體步驟如下：

步驟1：參照第1節(jié)方法，根據(jù)直流潮流模型或者線性化的交流潮流模型，獲得線路的有功功率線性表達式，并將其分為傳統(tǒng)機組出力、風(fēng)電場出力及負(fù)荷3個部分;

步驟2：獲得風(fēng)電總出力的概率密度函數(shù);

步驟2.1：對具有風(fēng)速相關(guān)性的風(fēng)電場組及獨立的風(fēng)電場，采用拉丁超立方采樣 (Latin Hypercube Sampling, LHS) 方法獲得風(fēng)速樣本;

圖1 直流概率潮流算法流程圖

步驟2.2：根據(jù)式(6)，由風(fēng)速樣本轉(zhuǎn)化為風(fēng)電功率樣本，根據(jù)式(4)，將風(fēng)電功率樣本相加獲得風(fēng)電總出力樣本；

步驟2.3：通過離散點擬合風(fēng)電總出力樣本，獲得風(fēng)電場總出力的概率密度曲線。

步驟3：總負(fù)荷的概率密度函數(shù);

步驟3.1：對具有負(fù)荷相關(guān)性的負(fù)荷組，采用解析法獲得其線性組合的概率密度函數(shù)表達式；

步驟3.2：對獨立的負(fù)荷，采用解析法獲得其累積和的概率密度函數(shù)表達式。

步驟4：采樣卷積方法，將風(fēng)電和負(fù)荷的概率密度函數(shù)以及常規(guī)機組出力結(jié)合起來，獲得線路有功功率的概率密度函數(shù)。

2.2 基于蒙特卡羅法的風(fēng)電總出力的概率密度函數(shù)求取

2.2.1 風(fēng)速樣本的LHS采樣

分別對式(4)中的n組具有相關(guān)性的風(fēng)電場組及其它獨立風(fēng)電場進行風(fēng)速的LHS采樣。

(1)對風(fēng)電場組進行相關(guān)風(fēng)速模擬。

針對式(4)中n組具有相關(guān)性的風(fēng)電場組，采用LHS方法進行相關(guān)風(fēng)速模擬。如其中第m組具有相關(guān)性的風(fēng)電場組，包含Km個具有相關(guān)性的風(fēng)電場，首先進行分層抽樣生成初始風(fēng)速樣本數(shù)據(jù)，詳細(xì)步驟如下：

然后根據(jù)風(fēng)電場的相關(guān)系數(shù)矩陣對風(fēng)速樣本進行重新排序，詳細(xì)步驟如下：

第m組風(fēng)電場組的相關(guān)系數(shù)矩陣為Cm，將Cm進行Cholesky分解Cm=BBT，獲得變換矩陣B。對Km個相互獨立的標(biāo)準(zhǔn)正態(tài)分布隨機變量，進行采樣規(guī)模為N的隨機抽樣，獲得Km×N的樣本矩陣W，由Z=BW獲得樣本矩陣Z，根據(jù)Z和W的順序變化獲得順序矩陣Ls(Ls為Km×N的矩陣，每一行為整數(shù)1到N的一個排列，對應(yīng)于Z中相應(yīng)行元素的大小順序)。對分層抽樣得到的初始風(fēng)速樣本數(shù)據(jù)按照順序矩陣Ls進行排序，即可獲得最終的風(fēng)速樣本。

(2)對獨立風(fēng)電場進行風(fēng)速模擬。

對于式(4)中獨立的風(fēng)電場，采用LHS方法進行獨立風(fēng)速模擬。具體流程與相關(guān)風(fēng)速模擬類似，首先進行分層抽樣生成初始風(fēng)速樣本數(shù)據(jù)；然后根據(jù)風(fēng)電場的相關(guān)系數(shù)矩陣對風(fēng)速樣本進行重新排序，順序矩陣直接由對標(biāo)準(zhǔn)正態(tài)分布隨機變量進行采樣規(guī)模為N的隨機抽樣獲得；最后對分層抽樣得到的樣本按照順序矩陣進行排序，即可獲得最終的風(fēng)速樣本。

若將獨立風(fēng)速的相關(guān)系數(shù)視為0，則可將所有風(fēng)電場視為一個相關(guān)風(fēng)電場組，相關(guān)風(fēng)速和獨立風(fēng)速采樣可以同時進行。

2.2.2 風(fēng)電總出力樣本轉(zhuǎn)化

根據(jù)式(6)，由風(fēng)速樣本轉(zhuǎn)化為風(fēng)電功率樣本。給定風(fēng)機參數(shù)，將每個風(fēng)速樣本帶入式(6)中，即可得到相應(yīng)的風(fēng)電功率。再根據(jù)式(4)，將所有風(fēng)電場的風(fēng)電功率樣本帶入式中，可以得到風(fēng)電總出力樣本，樣本規(guī)模為N。

2.2.3 風(fēng)電總出力的概率密度曲線擬合

通過離散點擬合風(fēng)電總出力樣本可以獲得風(fēng)電總出力的概率密度函數(shù)。即根據(jù)風(fēng)電總出力的取值范圍，得到其取值落在相應(yīng)范圍內(nèi)的頻數(shù)，進而以頻率代替概率，從而獲得風(fēng)電總出力的概率密度函數(shù)。

2.3 基于解析法的總負(fù)荷的概率密度函數(shù)求取

2.3.1 具有相關(guān)性負(fù)荷組的正態(tài)分布函數(shù)求取

已知正態(tài)分布函數(shù)的均值和方差就可以描述其概率密度函數(shù)。對于具有相關(guān)性的正態(tài)分布負(fù)荷組，其線性組合也服從正態(tài)分布，因此可以采用解析法快速計算出總負(fù)荷的均值和方差，得出總負(fù)荷的正態(tài)分布函數(shù)的表達式。

在式(9)中，對于相關(guān)性負(fù)荷組中的第k組負(fù)荷組，可以采用解析法快速得出總負(fù)荷的正態(tài)分布函數(shù)的均值和方差。計算公式如下：

(12)

式中k為負(fù)荷組號，第k個負(fù)荷組包含Lk個相關(guān)性負(fù)荷，k=1, 2, …,m；Ak為系數(shù)矩陣Ak=[α1,α2, …,αLk]；Pdk為負(fù)荷變量矩陣Pdk= [Pd1,Pd2, …,PdLk]；μ為對應(yīng)的Pdk的均值矩陣，C為Pdk的協(xié)方差矩陣。

2.3.2 總負(fù)荷的正態(tài)分布函數(shù)求取

在式(9)中，具有相關(guān)性的負(fù)荷組通過式(12)轉(zhuǎn)換為獨立的負(fù)荷，其總負(fù)荷的正態(tài)分布函數(shù)可以通過式(13)求出。

(13)

對于式(9)，也可以將獨立負(fù)荷視為相關(guān)性負(fù)荷，其與其它負(fù)荷間相關(guān)性系數(shù)為0，即協(xié)方差矩陣對應(yīng)項為0，則解析法可以同時處理變量相關(guān)和獨立的情況。

2.4 基于卷積的線路功率的概率密度函數(shù)求取

在式(3)中，根據(jù)2.2節(jié)和2.3節(jié)分別求出的風(fēng)電總出力以及總負(fù)荷的概率密度函數(shù)，采用卷積計算可以獲得線路功率的概率密度函數(shù)。

(14)

3 算例仿真與分析

本節(jié)分別對含多風(fēng)電場的IEEE RTS-96系統(tǒng)和IEEE 108系統(tǒng)[4,10]進行仿真分析，并與蒙特卡羅法、點估計法[5-6,17]及累積量法[9,18]進行比較，驗證文中方法的有效性和優(yōu)越性。

3.1 IEEE RTS-96系統(tǒng)算例

3.1.1 系統(tǒng)參數(shù)

本節(jié)對IEEE RTS-96系統(tǒng)進行了仿真計算。系統(tǒng)相關(guān)參數(shù)設(shè)置如下[10]：節(jié)點1、2、16、21為風(fēng)電場接入節(jié)點，其相關(guān)系數(shù)矩陣為：

節(jié)點16、18、19、20的負(fù)荷具有相關(guān)性，其相關(guān)系數(shù)矩陣為：

系統(tǒng)中風(fēng)速模型、風(fēng)機出力模型參數(shù)如下：

風(fēng)速參數(shù)：k=2.0，c=7.5；風(fēng)機參數(shù)：vci=4 m/s，vrate=10 m/s，vco=22m/s

3.1.2 各類方法比較分析

文中方法與累積量法以及點估計法進行了比較，以蒙特卡羅法作為基準(zhǔn)，蒙特卡羅法采樣次數(shù)為10 000。

(1)概率密度曲線比較。

對文中方法與累積量法以及點估計法分別獲得的線路有功功率的概率密度函數(shù)曲線進行比較，選取其中3條支路，分別如圖2～圖4所示。對3種方法獲得的系統(tǒng)支路的概率密度函數(shù)誤差進行比較，如圖5所示。

圖2 支路1-2有功功率的概率密度函數(shù)曲線的比較

圖3 支路1-5有功功率的概率密度函數(shù)曲線的比較

圖4 支路3-9有功功率的概率密度函數(shù)曲線的比較

圖5 不同方法的概率密度誤差比較

由圖2～圖4可以看出，文中方法獲得的概率密度曲線精度最高，其他兩種方法尤其是點估計法在高峰和低谷時差別較大。具體總結(jié)如下：

(a)對于極少數(shù)線路，如圖2所示，點估計法的誤差特別大，這是由于點估計法使用Gram-charlier展開求取概率密度函數(shù)的特性造成的，因此一般不使用點估計法求概率密度函數(shù)；

(b)對于多數(shù)線路，其有功功率的概率密度曲線為明顯的雙峰，如圖3所示，本文方法精度最高，累積量法和卷積法得到概率密度函數(shù)在峰值時誤差較大;

(c)對于少數(shù)線路，其有功功率的概率密度曲線為接近單峰，如圖4所示，仍然是本文方法精度最高，累積量法和卷積法得到概率密度函數(shù)誤差略小;

(d)圖5計算的概率密度誤差取的是絕對值，由圖5可以看出，本文方法的準(zhǔn)確性遠(yuǎn)高于累積量法和點估計法。并由圖2～圖4可以看出，本文方法的誤差來源其實是蒙特卡羅法的波動，其與蒙特卡羅法誤差很小，證明了本文方法的準(zhǔn)確性。

(2) 計算時間比較。

本文方法與蒙特卡羅法對IEEE RTS-96系統(tǒng)進行仿真計算。系統(tǒng)中風(fēng)電場個數(shù)為4，節(jié)點數(shù)為24，計算時間如表1所示。此外，為了比較系統(tǒng)中風(fēng)電場節(jié)點數(shù)量對計算時間的影響，本文調(diào)整算例中的風(fēng)電場個數(shù)，計算時間統(tǒng)計如表1所示。

表1 風(fēng)電場數(shù)量對計算時間的影響

從表1可以看出，隨著風(fēng)電場個數(shù)增加，計算時間增加，這是因為風(fēng)電場增加導(dǎo)致采樣規(guī)模增加，采樣時間增加。風(fēng)電場數(shù)目本中文方法和蒙特卡羅法的影響是相同的，因此文中方法可以在保證準(zhǔn)確性的前提下，大幅減少計算時間，證明了文中方法的優(yōu)越性。

3.2 IEEE 108系統(tǒng)算例

3.2.1 系統(tǒng)參數(shù)

對IEEE 108系統(tǒng)進行了仿真計算如表2所示。

表2 風(fēng)電場數(shù)據(jù)

IEEE 108系統(tǒng)規(guī)模比IEEE RTS-96系統(tǒng)更大，且系統(tǒng)中包含3組風(fēng)電場，各組間相互獨立，組內(nèi)風(fēng)電場間具有相關(guān)性。系統(tǒng)中風(fēng)電場節(jié)點的參數(shù)及相關(guān)性設(shè)置如下[4]：

節(jié)點2、3、5、13、14、16為風(fēng)電場接入節(jié)點，其相關(guān)系數(shù)矩陣為：

節(jié)點44、50、52、53為風(fēng)電場接入節(jié)點，其相關(guān)系數(shù)矩陣為：

節(jié)點82、83、84、86為風(fēng)電場接入節(jié)點，其相關(guān)系數(shù)矩陣為：

節(jié)點106、108、109、110的負(fù)荷具有相關(guān)性，其相關(guān)系數(shù)矩陣為：

系統(tǒng)中風(fēng)速模型、風(fēng)機出力模型參數(shù)如下：風(fēng)速參數(shù)：k=2.0，c=7.5；風(fēng)機參數(shù)：vci=4 m/s，vrate=10 m/s，vco=22 m/s

3.2.2 各類方法比較分析

(1)概率密度曲線比較

對文中方法與累積量法以及點估計法分別獲得線路有功功率的概率密度函數(shù)曲線進行比較，選取其中2條支路，如圖6、圖7所示。對3種方法獲得的支路的概率密度函數(shù)誤差進行比較，如圖8所示。由圖6～圖8可以看出，文中方法獲得的概率密度曲線和蒙特卡羅法擬合精度最高，而累積量法和點估計法在概率密度曲線峰值時誤差較大。

圖6 支路1-2有功功率的概率密度函數(shù)曲線的比較

圖7 支路3-5有功功率的概率密度函數(shù)曲線的比較

圖8 不同方法的概率密度誤差比較

(2)計算時間比較

計算時間比較如表3所示。

表3 計算時間的比較

4 結(jié)束語

提出了基于卷積的分類處理風(fēng)速相關(guān)性和負(fù)荷相關(guān)性的直流概率潮流算法。該方法采用基于拉丁超立方采樣的蒙特卡羅法處理風(fēng)速相關(guān)性，解析法處理負(fù)荷相關(guān)性，并使用卷積方法將兩者結(jié)合起來。解析法的使用避免了大規(guī)模采樣，縮減了計算時間。文中對IEEE RTS-96和IEEE 108系統(tǒng)進行了仿真計算，結(jié)果表明文中使用的方法可以在保證準(zhǔn)確性的前提下，大幅減少計算時間，可以應(yīng)用于大規(guī)模電力系統(tǒng)的概率潮流計算。