含風(fēng)電場概率最優(yōu)潮流的擬蒙特卡洛算法

李泳泉

(國網(wǎng)浙江云和縣供電有限公司，浙江 云和 323600)

風(fēng)力在清潔能源中占有不可替代的位置，云和縣內(nèi)風(fēng)力資源豐富，境內(nèi)黃源風(fēng)電站現(xiàn)已建成并接入電網(wǎng)。在風(fēng)電場接入電網(wǎng)后， 積極收集相關(guān)運行數(shù)據(jù)，開展相關(guān)分析研究工作。相鄰風(fēng)電場在并網(wǎng)運行后，其出力之間往往存在一定的相關(guān)性；風(fēng)電場并入電網(wǎng)運行后，出力的隨機性和相關(guān)性對電力系統(tǒng)影響不可忽略。在電網(wǎng)規(guī)劃與運行中，概率最優(yōu)潮流(Probabilistic Optimal Power Flow，簡稱P-OPF)能評估電力系統(tǒng)中各種不確定因素，有助于潛在薄弱環(huán)節(jié)的提前發(fā)現(xiàn)和解決，降低電網(wǎng)運行風(fēng)險。

P-OPF將電力系統(tǒng)中的不確定量視為服從某種概率分布的隨機變量，通過輸入變量的統(tǒng)計信息來求取輸出變量的統(tǒng)計信息，比如各階統(tǒng)計矩、概率密度函數(shù)(Probability Density Function，簡稱 PDF)、累積分布函數(shù)(Cumulative Distribution Function，簡稱CDF)等?，F(xiàn)存的P-OPF問題解法主要有蒙特卡洛模擬法[1](Monte Carlo Simulation，簡稱MCS)、累積量法[2-5](Cumulant Method)和點估計法[6-7](Point Estimate Method)。

這三種方法中，MCS的適用性最廣，只要建立相關(guān)的非正態(tài)隨機變量的數(shù)學(xué)模型，再生成大量樣本進(jìn)行最優(yōu)潮流(Optimal Power Flow，簡稱OPF)計算，進(jìn)而輸出變量的統(tǒng)計矩、PDF和CDF能夠準(zhǔn)確得到。但MCS最后結(jié)果的準(zhǔn)確、收斂需要龐大的計算量來支撐。

累積量法將P-OPF計算中的優(yōu)化過程視為一種概率映射，通過線性化潮流運算模型來求取輸出量的統(tǒng)計信息。但該方法有兩個不足，由于P-OPF計算中不等式約束的存在，輸出量與輸入量間函數(shù)關(guān)系的非線性很大，線性化該函數(shù)關(guān)系引入的誤差不可忽略；累積量法在相關(guān)性的處理上需繁雜的數(shù)學(xué)計算保證，并且輸入變量之間必須相互獨立，而風(fēng)電場機組出力之間的相關(guān)性必須考慮。

點估計[8]則是采用另一種方法來逼近輸出、輸入量間的函數(shù)關(guān)系。 在概率潮流問題中，能較準(zhǔn)確地求得輸出量的數(shù)學(xué)期望和標(biāo)準(zhǔn)差[9-10]，易于輸入變量相關(guān)性的處理[11-12]。 但其為減小計算量而引入的誤差也是很大的，不能準(zhǔn)確求取輸出量的高階統(tǒng)計矩[13-14]。 在 P- OPF運算中，由于不等式約束的限制，即使輸入隨機變量是正態(tài)的，輸出變量也不是正態(tài)的， 不少輸出變量的概率分布甚至是截斷的，該方法的精度更是不夠。

這三種方法均存在不足之處，將采取擬蒙特卡洛法(Quasi Monte Carlo Simulation，簡稱QMCS)來求取輸出量的統(tǒng)計矩。為達(dá)到加快收斂速度和保證較高的精度目的。將采用Sobol數(shù)列中的低偏差序列(Low Discrepancy Sequence)替代MCS的偽隨機序列(Pseudo Random Sequence)的方法?？紤]到隨機因素的相關(guān)性對概率最優(yōu)潮流的顯著影響[1]，本文用 Nataf轉(zhuǎn)換對輸入變量間的相關(guān)性進(jìn)行處理，用基于 Gauss- Hermite積分的插值法求解原相關(guān)的隨機變量映射到標(biāo)準(zhǔn)正態(tài)空間的相關(guān)系數(shù)，通過 Cholesky分解將問題轉(zhuǎn)換到獨立的正態(tài)空間。同時，利用輸出變量的概率加權(quán)矩(Probability Weighted Moments，簡稱PWM)，通過多項式轉(zhuǎn)換模型來建立輸出量的PDF和CDF。

1 Nataf轉(zhuǎn)換

相鄰區(qū)域的風(fēng)電場風(fēng)速因素可視為相關(guān)的非正態(tài)隨機向量?？紤]到直接得到多變量聯(lián)合概率分布函數(shù)較為困難。聯(lián)合概率分布模型是常用的方法之一，該模型能滿足各個變量之間的相關(guān)性。在概率分布問題處理中，正態(tài)分布擁有難以比擬的易處理特點，通過多維正態(tài)分布的聯(lián)合概率密度函數(shù)建立多維非正態(tài)隨機變量的模型， 即 Nataf變換[15-16]，在工程應(yīng)用中較多使用。 這里介紹一種用Nataf變換來生成風(fēng)電場風(fēng)速樣本的方法。

設(shè)有隨機變量，其CDF為F(x)；z為服從標(biāo)準(zhǔn)正態(tài)分布的隨機變量，其CDF為Φ(z)。根據(jù)等邊際概率原則有：

F(x)=Φ(z)

(1)

則可表示為

x=F-1[Φ(z)]

(2)

將兩個相關(guān)的隨機變量x1、x2分別由z1、z2表示，為確保經(jīng)由式(2)轉(zhuǎn)換得到的隨機變量x1、x2間的相關(guān)系數(shù)為ρx，需確定z1、z2間的相關(guān)系數(shù)ρz。ρx、ρz有如下函數(shù)關(guān)系：

(3)

式中φ(z1，z2，ρz)——二維標(biāo)準(zhǔn)正態(tài)分布的聯(lián)合概率密度函數(shù)；μi、σi——xi的數(shù)學(xué)期望和標(biāo)準(zhǔn)差(i=1，2)。

對于給定的ρx值，直接通過求解式(3)中的積分方程來求得ρz是很困難的。但注意到ρx是關(guān)于ρz的連續(xù)函數(shù)，且相關(guān)系數(shù)ρz的取值位于區(qū)間[1，-1]內(nèi)，因此可以用一個關(guān)于ρz的多項式來逼近兩者的函數(shù)關(guān)系。

(4)

(5)

通過對各類概率分布的測試，一個階數(shù)不超過9階的多項式便可非常精準(zhǔn)地逼近ρx、ρz的函數(shù)關(guān)系(n≤9)。式(4)中的二重積分可由11點的Gauss-Hermite積分準(zhǔn)確求得。

對于給定的ρx值，求解式(5)中的一元n次方程便可求得ρz，所求的ρz應(yīng)滿足如下條件：

|ρz|≤1，ρzρx≥0

(6)

設(shè)隨機向量X=(x1，…，xi，…xm)T的相關(guān)系數(shù)矩陣為Rx，可將式(2)由標(biāo)準(zhǔn)正態(tài)隨機向量Z=(z1，…，zi，…zm)T表示。對X中的任意兩個元素xi、xj(i≠j)，zi、zj間的相關(guān)系數(shù)可由插值法求得，從而得到Z的相關(guān)系數(shù)矩陣Rz。對Rz進(jìn)行Cholesky分解。便可將Z由相互獨立的標(biāo)準(zhǔn)正態(tài)隨機向量U表示：

Z=LU?U=L-1Z

(7)

RZ=LLT

(8)

L為下三角矩陣。

通過式(7)、式(2)，可以將獨立的標(biāo)準(zhǔn)正態(tài)隨機向量U轉(zhuǎn)換成具有指定邊際概率分布和相關(guān)系數(shù)矩陣Rx的非正態(tài)隨機向量X，利用多維正態(tài)分布的聯(lián)合概率密度函數(shù)描述相關(guān)的非正態(tài)多維隨機變量。

2 擬蒙特卡洛

P-OPF的實質(zhì)可定性為一個隨機非線性規(guī)劃問題：

min{f(x)|g(x)=b0，h(x)≤0|}

(9)

其中：g(x)=b0表示P-OPF計算中的等式約束，h(x)≤0表示不等式約束，P-OPF的數(shù)學(xué)模型可參見文獻(xiàn)[1]和[6]。

假定電力系統(tǒng)中各隨機因素組成的向量表示成X=(x1，…，xi，…xm)T。不考慮P-OPF的具體運算過程，輸出變量y和X之間的關(guān)系表達(dá)為

y=H(x1，…，xi，…Xm)

(10)

經(jīng)由Nataf轉(zhuǎn)換，可將xi由ui表示，將y視為ui的函數(shù)：

y=H*(u1，…，ui，…Um)

(11)

經(jīng)過式(1)的邊際轉(zhuǎn)換后，可以得到一個標(biāo)準(zhǔn)均勻分布空間：

ui=Φ-1(si)

(12)

si為標(biāo)準(zhǔn)均勻分布隨機變量(si∈[0，1])。此時，y可表示為

y=h(S)S=(s1，…，si，…sm)

(13)

QMCS的基本思路通過計算輸出變量的數(shù)學(xué)期望為例來展示：

(14)

式(14)的積分區(qū)域是m維單位空間[0，1]m。

在空間[0，1]m內(nèi)隨機選取點列Sj，應(yīng)用MCS，對式(14)求積分，得到y(tǒng)j的數(shù)學(xué)期望：

(15)

通過算法生成偽隨機序點列Sj，為確保所得的E[y]值準(zhǔn)確收斂，點列數(shù)n通常很大，計算時間長。本文采取QMCS來改善計算的收斂速度和精度，即：選取盡可能均勻地分布在m維空間[0，1]m的點列Sj，以低偏差序列來替代MCS的偽隨機序列，從而以較小的計算量得到較高精度的結(jié)果。

Sobol數(shù)列[18]在高維空間內(nèi)的均勻性較好，本文選用Sobol數(shù)列進(jìn)行計算。計算過程中所選點列是確定的，采用一種選取n個點列的點估計方法，為兼顧計算量和所得結(jié)果的精度，選用前(2-2001)列的Sobol數(shù)列進(jìn)行運算。

3 多項式轉(zhuǎn)換模型

對于概率分布未知且只有觀測樣本的隨機變量，需要找到一個合適的模型來表述其分布情況。本文采用基于韋伯分布(α=1，β=4)的多項式轉(zhuǎn)換模型對樣本進(jìn)行分布擬合。

設(shè)隨機變量t服從韋伯分布W(1，4)，由式(1)邊際轉(zhuǎn)換可得：

x=F-1[W(T)]

(16)

用一個10階多項式逼近函數(shù)F-1[W(T)]：

(17)

解出系數(shù)ak，隨機變量x可用韋伯分布對其進(jìn)行表達(dá)。

隨機變量的各階統(tǒng)計矩可以描述其統(tǒng)計特性，只要使得式(17)中多項式模型的統(tǒng)計矩等于x的統(tǒng)計矩，便可得到較佳的模擬效果。本文采用PWM來確定ak的值。

隨機變量的r階PWM(βr)定義[19]如下：

(18)

式(17)多項式模型的βr為

(19)

W(t)、w(t)分別為韋伯分布W(1，4)的CDF和PDF。Mk，r定義如下：

(20)

M0，0-M10，10可通過數(shù)值積分得到。

對于離散系統(tǒng)，可以將變量因素按x1≤…xi≤…≤xn排序，則βr的無偏估計值[20]如下：

(21)

依式(21)求解得到的前11階PWM，利用式(19)，得到如下的線性方程組：

(22)

求解方程得到系數(shù)ak。

需要說明的是，韋伯分布的逆累積分布函數(shù)(Inverse Cumulative Distribution Function，簡稱ICDF)存在閉式解析式：

t=[-ln(1-s)]1/4

(23)

式中s——標(biāo)準(zhǔn)均勻分布隨機變量。

將式(23)代入式(17)就可以得到ICDF：

(24)

使用ICDF的優(yōu)點在于：對于給定的概率點，可通過式(24)直接求得x的分位數(shù)。

4 算法步驟

本文求解P-OPF問題的具體步驟如下。

(1)選取m維Sobol數(shù)列的前(2-2001)點列，m為輸入變量的個數(shù)。通過式(12)，解得轉(zhuǎn)換到標(biāo)準(zhǔn)正態(tài)空間后的向量Uj。

(2)根據(jù)各輸入變量的ICDF，對式(3)用插值法求得任意兩個變量xi、xj在標(biāo)準(zhǔn)正態(tài)空間的相關(guān)系數(shù)ρz，建立相關(guān)系數(shù)矩陣Rz。

(3)按式(8)把Rz實行Cholesky分解，得到矩陣L(下三角)，向量Uj轉(zhuǎn)換為向量Zj，zj=LUj。

(4)將向量Zj中的各元素zj，i依其ICDF按式(2)轉(zhuǎn)換為xj，i，得到向量Xj。

(5)通過Xj代入到OPF，得到樣本的輸出量。

(6)利用所得的2 000個樣本求得各輸出變量的數(shù)學(xué)期望和標(biāo)準(zhǔn)差；根據(jù)式(21)求出其前11階PWM，代入式(22)解得ak，得到式(24)的ICDF，最后得到PDF、CDF。

MCS法的樣本也可先生成維相互獨立的標(biāo)準(zhǔn)正態(tài)隨機向量，再依步驟2-4生成樣本即可。

5 算例分析

利用文獻(xiàn)[1]中算例，其中均值為靜態(tài)平衡點處的有功負(fù)荷值，而有功負(fù)荷服從一般正態(tài)分布，均值的5%取為標(biāo)準(zhǔn)差，風(fēng)力發(fā)電的出力之間還應(yīng)考慮更強的相關(guān)性：

以105次MCS計算結(jié)果作為比較基準(zhǔn)，2 000次QMCS、MCS得到的數(shù)學(xué)期望、標(biāo)準(zhǔn)差的相對誤差平均值如表1所示。

表1 QMCS與MCS的誤差比較

由表1數(shù)據(jù)可知，次數(shù)相等情況下，QMCS所得數(shù)學(xué)期望和標(biāo)準(zhǔn)差的精度值要明顯高于MCS；由于 QMCS使用確定點列計算，其收斂性要優(yōu)于MCS。以有功潮流數(shù)學(xué)期望相對誤差平均值這個指標(biāo)為例，執(zhí)行10次MCS(2000次)運算，其誤差波動范圍為0.383～2.053。

對各輸出變量，基于所得的樣本，依式(24)求出其ICDF，在5%、25%、50%、75%、95%處求得對應(yīng)的分位數(shù)，最后計算其和MCS(105次)所得結(jié)果的相對誤差平均值，如表2所示。

表2 基于多項式模型所得分位數(shù)的誤差

輸出變量的概率分布通過多項式模型可以精確描述。

以線路32-113的無功潮流為例，其對應(yīng)的多項式模型系數(shù)如表3所示；所得 PDF、CDF如圖1、圖2所示，單峰概率分布很難描述線路32-113的無功潮流，本文提出的模型能夠很精確地求解出相對應(yīng)的 PDF、CDF圖。

表3 多項式模型的系數(shù)(線路32-113無功潮流)

圖1 線路32-113無功潮流的概率密度函數(shù)

在圖1的右端，可以看到有一個誤差呈現(xiàn)明顯尖峰狀。對P-OPF所有的輸出變量進(jìn)行模擬測試，此誤差位于部分不規(guī)則概率分布的PDF兩端，但由多項式模型建立的CDF位于5%～95%這一概率區(qū)間段的模擬效果，依舊很準(zhǔn)確。

圖2 線路32-113無功潮流的累積分布函數(shù)

6 結(jié)語

針對MCS法、累積量法、點估計法等P-OPF計算方法的不足，本文基于QMCS提出了一種求取P-OPF輸出變量PDF、CDF的方法。 該方法具有如下特點：只需一次插值就能得到任意兩種概率分布的相關(guān)系數(shù)計算公式，且能對服從任意連續(xù)概率分布隨機變量間的相關(guān)性進(jìn)行處理；較小計算量即可準(zhǔn)確求得輸出變量的數(shù)學(xué)期望、標(biāo)準(zhǔn)差等統(tǒng)計信息；對服從不規(guī)則分布輸出量的PDF、CDF進(jìn)行描述；在相同的計算次數(shù)下，其收斂性和精度要優(yōu)于MCS。

參考文獻(xiàn)：

[1] 楊歡，鄒斌．含相關(guān)性隨機變量的概率最優(yōu)潮流問題的蒙特卡羅模擬方法[J]．電力系統(tǒng)保護(hù)與控制，2012，40(19)：110-115.

YANG Huan，ZOU Bin．A Monte Carlo simulation method for probabilistic optimal power flow with correlated stochastic variables[J]．Power System Protection and Control，2012，40(19)：110-115.

[2]M M，P K，Q V H．Probabilistic optimal power flow[C]//IEEE Canadian Conference on Electrical and Computer Engineering，Waterloo，Canada，1998.

[3]S A，R W，A J．Cumulant- based probabilistic optimal power flow(P-OPF) with Gaussian and Gamma distributions[J]．IEEE Transactions on Power Systems，2005，20(2)：773-781.

[4]S A，R W，A J．Cumulant based probabilistic optimal power flow (P-OPF) [C]//International Conference on Probabilistic Methods Applied to Power Systems，Ames，American，2004.

[5]S A，R W，A J．Introduction to cumulant-based probabilistic optimal power flow (P-OPF)[J]．IEEE Transactions on Power Systems，2005，20(2)：1184-1186.

[6]V G，C C A．Probabilistic Optimal Power Flow in Electricity Markets Based on a Two-point Estimate Method[J]．IEEE Transactions on Power Systems，2006，21(4)：1883-1893.

[7]潘煒，劉文穎，楊以涵．概率最優(yōu)潮流的點估計算法[J]．中國電機工程學(xué)報，2008，28(16)：28-33.

PAN Wei，LIU Wenying，YANG Yihan．Point Estimate Method for Probabilistically Optimal Power Flow Computation[J]．Proceedings of the CSEE，2008，28(16)：28-33.

[8]H H P．An efficient point estimate method for probabilistic analysis[J]．Reliability Engineering and System Safety，1998，59(3)：261-267.

[9]S C L．Probabilistic Load-Flow Computation Using Point Estimate Method[J]．IEEE Transactions on Power Systems，2005，20(4)：1843-1851.

[10]M J M，P R J．Point Estimate Schemes to Solve the Probabilistic Power Flow[J]．IEEE Transactions on Power Systems，2007，22(4)：1594-1601.

[11]M J M，B L，C A J，et al．Probabilistic power flow with correlated wind sources[J]．IET Generation Transmission and Distribution，2010，4(5)：641-651.

[12]楊歡，鄒斌．含相關(guān)性隨機變量的概率潮流三點估計法[J]．電力系統(tǒng)自動化，2012，36(15)：51-56.

YANG Huan，ZOU Bin．A Three-point Estimate Method for Solving Probabilistic Power Flow Problems with Correlated Random Variables[J]．Automation of Electric Power Systems，2012，36(15)：51-56.

[13]U J．Probabilistic load flow with correlated wind power injections[J]．Electric Power systems research，2010，80(5)：528-536.

[14]U J．Probabilistic load flow in systems with wind generation[J]．IET Generation Transmission and Distribution，2009，3(12)：1031-1041.

[15]L P L，D K A．Multivariate distribution models with prescribed marginals and convariances [J]．Probabilistic Engineering Mechanics，1986，1(2)：105-112.

[16]Q X．Evaluating correlation coefficient for Nataf transformation [J]．Probabilistic Engineering Mechanics，2014(37)：1-6.

[17]N J C，S A F M．Applications of a method for the efficient computation of posterior distributions [J]．Journal of the Royal Statistical Society. Series C (Applied Statistics)，1982，31(3)：214-225.

[18]M W J，C R E．Quasi-Monte Carlo integration[J]．Journal of Computational Physics，1995， 122(2)：218-230.

[19]G J A，L J M，M N C．Probability weighted moments: definition and relation to parameters of several distributions expressable in inverse form[J]．Water Resources Research，1979， 15(5)：1049-1054.

[20]H J R M，W J R，W E F．Estimation of the generalized extreme-value distribution by the method of probability-weighted moments [J]．Technometrics，1985， 27(3)：251-261．