李泳泉

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

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

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

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

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

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

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

1 Nataf轉(zhuǎn)換

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

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

F(x)=Φ(z)

(1)

則可表示為

x=F-1[Φ(z)]

(2)

將兩個(gè)相關(guān)的隨機(jī)變量x1、x2分別由z1、z2表示，為確保經(jīng)由式(2)轉(zhuǎn)換得到的隨機(jī)變量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)。

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

(4)

(5)

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

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

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

(6)

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

Z=LU?U=L-1Z

(7)

RZ=LLT

(8)

L為下三角矩陣。

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

2 擬蒙特卡洛

P-OPF的實(shí)質(zhì)可定性為一個(gè)隨機(jī)非線性規(guī)劃問(wèn)題：

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

(9)

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

假定電力系統(tǒng)中各隨機(jī)因素組成的向量表示成X=(x1，…，xi，…xm)T。不考慮P-OPF的具體運(yùn)算過(guò)程，輸出變量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)過(guò)式(1)的邊際轉(zhuǎn)換后，可以得到一個(gè)標(biāo)準(zhǔn)均勻分布空間：

ui=Φ-1(si)

(12)

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

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

(13)

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

(14)

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

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

(15)

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

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

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

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

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

x=F-1[W(T)]

(16)

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

(17)

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

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

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

(18)

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

(19)

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

(20)

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

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

(21)

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

(22)

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

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

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

(23)

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

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

(24)

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

4 算法步驟

本文求解P-OPF問(wèn)題的具體步驟如下。

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

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

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

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

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

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

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

5 算例分析

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

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

表1 QMCS與MCS的誤差比較

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

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

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

輸出變量的概率分布通過(guò)多項(xiàng)式模型可以精確描述。

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

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

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

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

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

6 結(jié)語(yǔ)

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

參考文獻(xiàn)：

[1] 楊歡，鄒斌．含相關(guān)性隨機(jī)變量的概率最優(yōu)潮流問(wèn)題的蒙特卡羅模擬方法[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)潮流的點(diǎn)估計(jì)算法[J]．中國(guó)電機(jī)工程學(xué)報(bào)，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ī)變量的概率潮流三點(diǎn)估計(jì)法[J]．電力系統(tǒng)自動(dò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．