胡建軍，郭 飛*，王 波，吳佳靜

(1.國網(wǎng)寧夏電力有限公司，銀川 750001；2.國網(wǎng)寧夏電力有限公司吳忠供電公司，吳忠 751100；3.寧夏信通網(wǎng)絡科技有限公司，銀川 750001)

隨著全球能源資源配置的逐步調(diào)整，化石能源和可再生能源供應成本總體呈現(xiàn)出“一升一降”的趨勢[1-2]，風電技術(shù)的發(fā)展和突破使其在可再生能源中所占比重逐漸加大，未來電力系統(tǒng)新場景已全面開啟[3]。但風能具有強波動性和間歇性，并網(wǎng)后使得系統(tǒng)的源端呈現(xiàn)強隨機性，給系統(tǒng)運行穩(wěn)定性帶來嚴峻挑戰(zhàn)。

概率潮流(probabilistic power flow，PPF)可考慮系統(tǒng)隨機信息，研究成果可分為模擬法(monte carlo simulation，MCS)、解析法和點估計法三類[4]。模擬法不受系統(tǒng)規(guī)模和非線性的約束，計算精度高，但比較耗時。拉丁超立方(Latin hypercube sampling，LHS)作為一種改進方法，通過分層抽樣改善了隨機變量在空間中的覆蓋程度，采樣效率提高顯著[5]。解析法常用半不變量法，突出特點是計算速度快，但其需要對潮流方程線性化處理，因而當輸入隨機變量波動范圍較大時，計算結(jié)果誤差較大。點估計法根據(jù)輸入隨機變量的前幾階矩信息，得到輸出隨機變量的均值和方差，但無法描述整體概率特性。

描述風電功率隨機波動特性的概率模型已有廣泛研究，文獻[6-7]采用Weibull分布和常規(guī)正態(tài)分布分別擬合風電出力和負荷概率模型，其理論上具有普適性，但無法正確捕捉變量多峰性和強隨機性。文獻[8]采用高斯混合模型(Gaussian mixture model，GMM)刻畫負荷隨機性，并通過歷史數(shù)據(jù)驗證了該模型的準確性。文獻[9]探究了采用GMM不同維數(shù)情形下擬合風電場和集群風電場的可行性。事實上，GMM可以更加準確地描述變量強隨機波動，其擬合曲線可呈現(xiàn)不對稱、多峰分布的概率特性。

在構(gòu)建變量的概率模型基礎上，如何有效解決計及新能源出力強隨機性的 PPF 是現(xiàn)有潮流計算研究的難點?，F(xiàn)有PPF方法難以計及狀態(tài)變量的先驗信息和當前觀測數(shù)據(jù)的內(nèi)在聯(lián)系，貝葉斯理論將狀態(tài)變量視為變動的參數(shù)，隨著觀測數(shù)據(jù)的不斷累積而改變[10]。近似貝葉斯計算(approximate Bayesian computation，ABC)方法通過對比模擬數(shù)據(jù)集和觀測數(shù)據(jù)集來估計變量后驗概率[11]，可用作解決電力系統(tǒng)似然函數(shù)未知的問題。此外，已有研究中處理隨機變量相關(guān)性較廣的方法包括Cholesky分解、隨機排序等[12-13]。但在實際應用中，描述隨機變量兩兩之間關(guān)系的相關(guān)系數(shù)通常是非理想條件的，采用奇異值分解的方法得到考慮相關(guān)性的狀態(tài)變量樣本，可靈活處理相關(guān)系數(shù)矩陣非正定的情形。

基于上述分析，現(xiàn)提出一種基于序貫重要性采樣(sequential importance sampling，SIS)改進ABC的方法解決PPF問題。首先采用GMM建立風電出力不確定模型，利用奇異值分解結(jié)合Nataf變換靈活處理相關(guān)系數(shù)矩陣非正定的情形。其次，將具有先驗信息的狀態(tài)變量轉(zhuǎn)化為模擬數(shù)據(jù)集，通過ABC對比觀測數(shù)據(jù)集以獲取后驗估計。最后，采用樣本加權(quán)方式和擾動核來避免出現(xiàn)低接受率情況，通過算例仿真證明所提方法的高精度和高效性。

1 隨機潮流不確定模型

1.1 風電場概率模型建立

風電場并網(wǎng)概率模型的建立方法主要分為參數(shù)方法和非參數(shù)方法兩種。文獻[14]采用GMM建立風電出力不確定模型，該模型概念上簡單易懂，融合了參數(shù)模型和非參數(shù)模型各自的優(yōu)勢，可以靈活高效地考慮多維隨機變量之間的相關(guān)性變化。采用混合高斯分布表示風電場出力的概率函數(shù)fGMM(x)為

(1)

(2)

文獻[9]驗證了二維GMM可有效表示風電場出力。因此用二分量GMM表示風電場有功出力的概率分布為

(3)

式(3)中：α為權(quán)重。

采用最大期望(EM)算法求解GMM中未知參數(shù)值，其分作兩步迭代更新隱藏變量和模型分布參數(shù)，最終收斂后，保證得到全局最大值。從而降低計算最大似然估計的難度。

給定風電場出力數(shù)據(jù)X={x1,x2,…，xN}，GMM的似然函數(shù)為

(4)

式(4)中：Θ=(θ1,θ2,…,θN)。對式(4)取對數(shù)得

(5)

使用EM算法求解風電場出力概率函數(shù)中未知參數(shù)結(jié)果為

(6)

式(6)中：m為求解最大化對數(shù)似然函數(shù)次數(shù)；n為參數(shù)變量個數(shù)；p為參數(shù)的迭代次數(shù)；Gm為后驗密度函數(shù)。

1.2 線性化交流潮流模型

交流潮流方程一般可表示為

(7)

式(7)中：W為節(jié)點注入功率，即輸入隨機變量，包括負荷、風電節(jié)點注入功率；Z為支路潮流，即輸出隨機變量；x為系統(tǒng)的狀態(tài)變量，包括電壓幅值和相角。

概率潮流的目的是根據(jù)節(jié)點注入功率W的概率分布，通過W與狀態(tài)變量x的線性化關(guān)系g，得到x的數(shù)字特征，進而求得支路有功、無功功率的數(shù)學特征。

2 考慮相關(guān)隨機變量概率潮流計算

2.1 輸入變量相關(guān)性控制

評估新能源并網(wǎng)對系統(tǒng)運行的影響，應同時考慮出力的強隨機性和相關(guān)性。具有相關(guān)性的隨機功率注入對系統(tǒng)產(chǎn)生潮流擾動，常用隨機排序、Cholesky分解和遺傳算法的處理方法來獲得具有相關(guān)性的樣本序列。為了可靈活處理相關(guān)系數(shù)矩陣非正定的情況，在此采用奇異值分解結(jié)合Nataf變換的方法處理多維輸入隨機變量。假定存在n維輸入變量X=[X1,X2,…,Xn]T，其相關(guān)系數(shù)矩陣ρX=(ρXij)n×n，累積分布函數(shù)Fi(Xi)。通過式(8)等概率變換為標準正態(tài)隨機變量Z=[Z1,Z2,…,Zn]T。

(8)

式(8)中：Φ(·)和Φ-1(·)分別為Zi的累積分布函數(shù)及其反函數(shù)；Fi(Xi)和Φ(Zi)服從U[0,1]。

變換后ρZ=(ρZij)n×nρX和ρX的非對角元素由式(9)得出，即

ρZij=T(ρXij)ρXij

(9)

式(9)中：T(ρXij)為ρXij和ρZij之間的轉(zhuǎn)換系數(shù)，該表達式與Xi、Xj的概率分布有關(guān)。文獻[15]給出了已知ρXij時，如何快速求得ρZij的方法。

上述過程完成后，需再次變換到獨立標準的正態(tài)空間。將ρZ進行奇異值分解，可處理非正定或非滿秩矩陣等情況。對于任意矩陣Am×n，都存在式(10)分解：

(10)

式(10)中：Σ為由Am×n奇異值構(gòu)成的對角陣，對角元素降序排列；Um×n為列正交矩陣；Vn×n為正交陣。當m=n時，U=V。

奇異值分解具有以下性質(zhì)，假設H為n維獨立標準正態(tài)分布向量，對ρZ進行奇異值分解生成矩陣L為

(11)

式(11)中：UρZ為酉矩陣；ΣρZ為由矩陣奇異值組成的對角陣，對角元素按大小排列；則式(12)定義相關(guān)系數(shù)矩陣Z*為

(12)

2.2 ABC序貫蒙特卡洛算法

2.2.1 ABC簡介

ABC的基本理念是在已知輸入隨機變量先驗概率分布的基礎上，利用貝葉斯定理匹配觀測數(shù)據(jù)，從而快速得到輸出變量的概率特性。傳統(tǒng)貝葉斯公式為

(13)

式(13)中：P(θ|Y)為參數(shù)θ基于樣本Y的后驗分布；f(Y|θ)為條件概率，即似然函數(shù)；P(θ)為基于已知經(jīng)驗提供的θ的先驗概率分布；f(Y)為總條件概率，常歸一化為常值。

在無法獲取似然函數(shù)的條件下，ABC的思想是先確立θ的先驗概率分布P(θ)，用式(14)來估計后驗分布。

P[d(Y,Y′)≤ε]P(θ)

(14)

式(14)中：條件分布Pε(Y|Y′)用作衡量樣本差距；P(Y′|θ)為由概率模型生成的模擬量；誤差函數(shù)d(Y,Y′)和閾值ε用于判斷采樣點是否符合觀測值。

通過這一過程，θ先驗分布的采樣點通過觀測值Y的作用下轉(zhuǎn)化為后驗分布。具體是，從先驗分布P(θ)中隨機采樣得到θ*，從而構(gòu)成模擬數(shù)據(jù)集Y′，考慮到計算效率，通常不應采取完整數(shù)據(jù)集進行比較，而是選取適當?shù)母乓y(tǒng)計量來替代完整數(shù)據(jù)集，通過概要統(tǒng)計量S(Y)、S(Y′)比較模擬數(shù)據(jù)集與觀測數(shù)據(jù)集，計算d[S(Y),S(Y′)]≤ε的結(jié)果，則可以判斷是否將θ*作為參數(shù)的后驗概率估計值。

然而，ABC算法具有接受率低的不足，主要有以下兩類原因，一是當先驗、后驗分布概率特征差異較大時，無信息性的先驗樣本會降低接受率；二是閾值ε的選取會產(chǎn)生不同的后驗估計結(jié)果。這將ABC局限在原始數(shù)據(jù)集是低維、離散的情況，對于高維、連續(xù)的原始數(shù)據(jù)集，往往很難應用。

2.2.2 基于SIS改進ABC

(15)

以上方法在此記為ASMC，具體步驟如下：

(2)前一代粒子集θi以權(quán)重wi抽取一個樣本θ*后，通過擾動核得到θ**=q(θ**|θ*)。

(3)將θ**代入P(θ)，P(θ**)=0成立則代入實際模型計算模擬數(shù)據(jù)集Y′，否則回到步驟(2)。

(4)判斷d[S(Y),S(Y′)]≤ε,成立則保留該樣本θ**，否則回到步驟(2)。

(6)回到步驟(1)，令t=t+1，直到迭代次數(shù)達到T，得到最終樣本序列θT。

2.3 基于ASMC的概率潮流計算流程

將本文ASMC方法結(jié)合奇異值分解的Nataf變換用于風電并網(wǎng)系統(tǒng)概率潮流計算當中，流程如圖1所示，具體如下。

圖1 PPF計算流程圖Fig.1 Flow chart of PPF calculation

(2)觀測數(shù)據(jù)集生成。對于系統(tǒng)中存在的n個狀態(tài)變量(電壓和相角)，設定先驗分布P(θn)，基于1.2節(jié)建立實際模型g(·|θn)，根據(jù)2.1節(jié)生成具有相關(guān)性的隨機變量觀測的樣本數(shù)據(jù)集Y。

(3)ASMC模型生成。建立誤差函數(shù)d[S(Y),S(Y′)]，擾動核q(·|θ*)，概要統(tǒng)計量S(·)，閾值ε={ε1,ε2,…,εT}。根據(jù)2.2.2節(jié)完成模擬數(shù)據(jù)集Y′和觀測數(shù)據(jù)集Y的差異程度對比。

(4)概率潮流計算。將得到的各組狀態(tài)變量樣本集代入潮流計算式中，利用概率統(tǒng)計方法得到輸出變量的數(shù)學特征。

3 算例分析

3.1 算例介紹

在改造的IEEE-30節(jié)點系統(tǒng)中進行仿真測試，系統(tǒng)結(jié)構(gòu)如圖2所示。

圖2 改造的IEEE-30系統(tǒng)結(jié)構(gòu)圖Fig.2 System structure of modified IEEE-30

在節(jié)點17、節(jié)點19、節(jié)點21、節(jié)點24接入 30 MW 的風電場，新能源滲透率水平約為20%，恒定節(jié)點功率因數(shù)0.9。此外，假定系統(tǒng)中的負荷均服從正態(tài)分布，以標準系統(tǒng)中給定的負荷值作為期望值，取標準差為期望值的10%。假設風功率相關(guān)系數(shù)矩陣為非正定陣，表示為

(16)

以下采用本文方法進行潮流計算，由于ASMC是通過對粒子集賦予相應權(quán)重，再經(jīng)過閾值迭代后，使樣本逐步符合真實后驗概率分布，先驗分布的選取不再是影響樣本接受率的重要指標。在此根據(jù)歷史經(jīng)驗對狀態(tài)變量先驗分布進行假設，假定系統(tǒng)中的電壓服從N(1,0.12)，相角服從N(μθ,0.12)，采用傳統(tǒng)直流潮流算法解出μθ。選擇正態(tài)分布η～N(0,σ2)作為擾動核q(·|θ*)，疊加η作為樣本點的偏移，將上一代形成的樣本集方差作為σ。為避免選擇閾值ε的主觀性對后驗概率估計精度造成影響，取ε={ε1,ε2,…,εT}，根據(jù)迭代次數(shù)其取值逐漸減小。距離函數(shù)采用經(jīng)典的均方根誤差作為判定公式，即

(17)

式(17)中：XY,m和XY′,m分別為第m個觀測數(shù)據(jù)集Y和模擬數(shù)據(jù)集Y′計算得到的概要統(tǒng)計量。

3.2 風電出力概率擬合

從德國某電力公司獲取的風電出力數(shù)據(jù)，將其作為觀測數(shù)據(jù)代入EM算法計算公式中，求解出GMM中的未知參數(shù)。表1給出了采用二維 GMM 擬合實際風功率得到的各子成分參數(shù)。圖3比較了采用Weibull函數(shù)、GMM擬合的風電功率波動特性，可以看出，相比于Weibull分布，采用二維加權(quán)高斯混合概率模型得到的概率曲線更加貼合原始樣本概率曲線，具有更好地擬合精度，可以用作描述風電出力(pu為標幺值)的波動特性。

表1 二維GMM子成分參數(shù)Table.1 Two dimensional GMM subcomponent parameters

圖3 實際風電出力GMM擬合Fig.3 GMM approximation of wind power output

3.3 潮流計算

在3.1節(jié)所述系統(tǒng)進行潮流計算，在相同的條件下，分別采用常規(guī)ABC、LHS和本文方法進行PPF計算，以MCS方法對輸入隨機變量采樣104次結(jié)果作為真實后驗分布。得到節(jié)點17、節(jié)點19電壓如圖4所示?？梢钥闯?，ASMC方法較常規(guī)ABC而言，其計算獲得的概率密度曲線更接近與基于MCS方法得到的真實后驗分布，即常規(guī)ABC在先驗分布與真實后驗分布相差較大時，具有容易出現(xiàn)接受率低和準確性不足的情況，圖4中ABC計算得到的概率密度曲線明顯方差過小，不符合實際運行情況。表2 給出了不同方法下PPF計算的耗時比較，其中基準MCS、LHS、ABC和ASMC仿真時間分別為24.125、12.743、32.378、6.489 s。結(jié)合圖4和表2可知，常規(guī)ABC方法耗時最長，原因在于其具有接受率低下的缺點，選擇較低的閾值ε就會帶來計算開支時間長的問題。ASMC與LHS方法二者的計算精度相仿，但ASMC方法具有更高的計算效率，在迭代過程中通過擾動核和賦予樣本權(quán)重來產(chǎn)生樣本數(shù)據(jù)集，可以進一步提高效率。

圖4 負荷節(jié)點電壓的概率分布Fig.4 Probability distribution for voltage of bus

表2 各算法耗時比較Table 2 Time consuming comparison of algorithms

圖5為節(jié)點10電壓的累積概率分布函數(shù)，可以看出，該節(jié)點電壓存在越上限(>1.05 pu)的情況，采用基準MCS、LHS、ABC和ASMC所得的概率分別為5.18%、5.35%、9.46%和5.33%。ABC計算得到的結(jié)果偏大，高估了系統(tǒng)潛在運行風險。圖6為支路16、17的有功潮流概率分布，ASMC得到的有功功率概率密度分布函數(shù)(probability density function,PDF)、累積分布函數(shù)(cumulative distribution function,CDF)都和基準MCS非常接近。可以說明本文方法可以獲取更準確的概率特性，計算結(jié)果可為系統(tǒng)運行提供一定的參考。

圖5 節(jié)點10電壓累積分布函數(shù)Fig.5 Cumulative distribution function for voltage of bus 10

圖6 支路16、17有功潮流概率分布Fig.6 Probability distribution for active power flow of line 16，17

圖7 不同εT條件下電壓概率密度曲線Fig.7 Voltage probability density curve under different εT conditions

對于采用ASMC方法計算PPF而言，閾值ε={ε1,ε2,…,εT}的最小取值εT可能會產(chǎn)生不同的后驗估計結(jié)果，采用不同的εT取值分別進行概率潮流計算，得到如圖7所示的不同結(jié)果?？梢钥闯?，不同閾值εT取值情況下，得到的狀態(tài)變量概率密度差異較大，在總體趨勢上，隨著εT的增大，概率分布越來越分散，即方差在增大。這說明了閾值越小，抽取的樣本將全部來自后驗分布，對于后驗的估計精度更高，產(chǎn)生的近似后驗分布差異很小。因此，需綜合考慮實際處理問題的復雜性，選擇合適的閾值下限。

4 結(jié)論

(1)針對新能源風電出力，建立的二維GMM具有更好的擬合精度，可以用作描述風電出力的波動特性，模型的精度和建模效率得到提升。

(2)奇異值分解結(jié)合Nataf變換可準確計及輸入變量相關(guān)系數(shù)矩陣非正定的情形，更符合實際應用場景，可適用于新能源風力發(fā)電并網(wǎng)及電力負荷具有波動性的未來電網(wǎng)的概率潮流分析。

(3)提出的ASMC方法用于PPF計算效率高，其既能夠利用ABC計及狀態(tài)變量先驗信息的思想，又采用樣本加權(quán)方式和擾動核來避免低接受率情況的發(fā)生，最終將狀態(tài)變量概率分布求解問題轉(zhuǎn)變?yōu)橐环N參數(shù)反演問題，為PPF注入了一種新思路。