交叉熵蒙特卡洛法在含風(fēng)電場電力系統(tǒng)小擾動穩(wěn)定性評估中的應(yīng)用

蔣 慧

(淮南聯(lián)合大學(xué) 智能制造學(xué)院, 安徽 淮南 232038)

0 引 言

風(fēng)電作為最具有開發(fā)潛力的可再生能源在世界范圍內(nèi)迅速發(fā)展。根據(jù)全球風(fēng)能理事會《2022年全球風(fēng)能報告》發(fā)布的數(shù)據(jù),2021年,全球新增風(fēng)電裝機(jī)容量93.6 GW,累計風(fēng)電裝機(jī)容量達(dá)837 GW,同比增長12%。2021年,中國風(fēng)電新增并網(wǎng)裝機(jī)容量47.57 GW,占全球當(dāng)年新增裝機(jī)容量的50.91%。截至2022年底,中國風(fēng)電累計裝機(jī)容量395.6 GW,占全球風(fēng)電總裝機(jī)容量三分之一以上。隨著風(fēng)電場裝機(jī)規(guī)模的不斷增加,產(chǎn)生的風(fēng)功率波動致使風(fēng)電輸出功率不穩(wěn)定,甚至造成風(fēng)力發(fā)電機(jī)組脫網(wǎng),特別是隨著滲透率的不斷提升,風(fēng)電功率隨機(jī)性和波動性被視為隨機(jī)擾動注入電網(wǎng),從而對電網(wǎng)的穩(wěn)定性產(chǎn)生影響[1-2]。

電力系統(tǒng)小擾動穩(wěn)定性分析方法主要包括確定分析法和概率法兩類,目前多采用概率分析方法[2-4],隨機(jī)理論中的隨機(jī)響應(yīng)面法、點估計法、蒙特卡洛法相繼被應(yīng)用于分析隨機(jī)電力系統(tǒng)的隨機(jī)小擾動穩(wěn)定性問題[5-9]。文獻(xiàn)[5]考慮風(fēng)速相關(guān)性的風(fēng)電場群模型,運用隨機(jī)響應(yīng)面法分析了系統(tǒng)的小擾動穩(wěn)定性。文獻(xiàn)[6]將風(fēng)速作為威布爾分布的隨機(jī)變量,基于2m+1點估計法和Cornish Fisher展開法,采用正交變換技術(shù)處理風(fēng)電場之間的相關(guān)性,提出一種對具有相關(guān)風(fēng)源的電力系統(tǒng)進(jìn)行小信號穩(wěn)定性分析的概率方法。其中,由于蒙特卡洛仿真(Monte Carlo Simulation, MCS)對系統(tǒng)輸入概率函數(shù)模型無特殊要求,且收斂性與系統(tǒng)維數(shù)無關(guān),所以被廣泛應(yīng)用。文獻(xiàn)[7]將系統(tǒng)的發(fā)電方式、負(fù)荷水平和網(wǎng)絡(luò)拓?fù)涞牟淮_定性作為隨機(jī)擾動,基于蒙特卡洛狀態(tài)抽樣方法,通過分析系統(tǒng)特征參數(shù)的概率特性和失穩(wěn)概率研究了電力系統(tǒng)小擾動概率穩(wěn)定問題。文獻(xiàn)[8]針對小擾動穩(wěn)定概率評估問題,利用特征值分析法結(jié)合蒙特卡洛模擬法對系統(tǒng)進(jìn)行評估。通過建立負(fù)荷水平和網(wǎng)絡(luò)拓?fù)潆S機(jī)變量的概率模型和小干擾穩(wěn)定概率指標(biāo),對新英格蘭系統(tǒng)進(jìn)行了評估。但是傳統(tǒng)的蒙特卡洛算法在計算精度和計算時間兩個方面存在矛盾,為了克服這一矛盾,常采用分層抽樣法、重要抽樣法、等分散抽樣法、拉丁超立方抽樣法等減方差技術(shù)提高收斂速度和仿真效率[9]。其中,重要抽樣法最為常用,基于交叉熵的蒙特卡洛法(Cross Entropy Monte Carlo Simulation, CEMCS)作為一種高效的重要抽樣法,可以克服傳統(tǒng)蒙特卡洛法對稀有時間的不敏感問題,引入到電力系統(tǒng)可靠性分析中,有效地克服了傳統(tǒng)蒙特卡洛仿真收斂滿的缺點[10-11]。文獻(xiàn)[10]提出一種基于交叉熵的蒙特卡洛法,并將其應(yīng)用于發(fā)電系統(tǒng)充裕度評估中,使用重要抽樣密度函數(shù),通過求解最優(yōu)問題獲得該函數(shù)的最優(yōu)參數(shù),從而提高傳統(tǒng)蒙特卡洛法的抽樣效率。文獻(xiàn)[11]提出一種計及多狀態(tài)離散型和連續(xù)型隨機(jī)變量的蒙特卡洛擴(kuò)展交叉熵法,融合預(yù)抽樣階段和最優(yōu)抽樣階段的樣本信息得到綜合可靠性指標(biāo),進(jìn)一步提升了電網(wǎng)可靠性評估擴(kuò)展交叉熵法的計算效率。

文中將風(fēng)電功率波動和負(fù)荷變化作為隨機(jī)小擾動源,構(gòu)建基于特征值的小擾動穩(wěn)定性概率指標(biāo),通過交叉熵非序貫蒙特卡洛法對系統(tǒng)小擾動穩(wěn)定性進(jìn)行仿真。

1 小擾動穩(wěn)定性概率指標(biāo)

小擾動穩(wěn)定分析即在系統(tǒng)穩(wěn)定運行點時受到隨機(jī)小干擾后的系統(tǒng)穩(wěn)定性。

綜合文獻(xiàn)[3-8,12]對基于特征值的小擾動穩(wěn)定性進(jìn)行論述。

設(shè)電力系統(tǒng)描述一階微分-代數(shù)方程組為

(1)

式中:x----系統(tǒng)狀態(tài)向量;

u----系統(tǒng)輸入向量;

y----系統(tǒng)輸出向量。

在穩(wěn)定運行點(x0,y0)處對上式進(jìn)行線性化可得

(2)

當(dāng)零輸入狀態(tài)時,式(2)變換為

(3)

其特征方程為

|λI-A|=0。

(4)

特征方程的每個解為共軛復(fù)數(shù)

λi=σi±jωi,

(5)

式中:σi----特征值的實部,反映了阻尼振蕩模式;

ωi----振蕩頻率,ωi=2πfi。

由于阻尼比與特征值存在以下關(guān)系,

(6)

可知阻尼比和特征值的實部互為異號,根據(jù)李雅普諾夫穩(wěn)定性第一定理可知,狀態(tài)矩陣A的特征值有負(fù)實部,則系統(tǒng)小擾動穩(wěn)定性,同時,在電力系統(tǒng)中,阻尼比反應(yīng)的是系統(tǒng)振蕩模式,綜合兩個參數(shù)的關(guān)系以及系統(tǒng)穩(wěn)定狀態(tài)見表1。

表1 小擾動穩(wěn)定性對照表

由表1可知,阻尼比小于零為系統(tǒng)小擾動失穩(wěn)的條件,即小擾動失穩(wěn)指數(shù)(Index of Small Signal Instability, ISSI)為

(7)

式中:f(ξ)----阻尼比ξ的概率密度函數(shù)。

當(dāng)ISSI>0時,系統(tǒng)不穩(wěn)定。

實際情況下,只有當(dāng)ξ>0.03時,才能保證系統(tǒng)具有足夠的低頻振蕩穩(wěn)定性,所以,當(dāng)ξ∈(0,0.03)時,屬于弱阻尼振蕩模式,存在不穩(wěn)定風(fēng)險,確定不穩(wěn)定風(fēng)險指標(biāo),

(8)

當(dāng)IR>0時,存在不穩(wěn)定風(fēng)險。

2 交叉熵非序貫蒙特卡洛算法

2.1 非序貫蒙特卡洛法

非序貫蒙特卡洛法是一種大量隨機(jī)試驗法,在電力系統(tǒng)穩(wěn)定性分析應(yīng)用中,非序貫蒙特卡洛法不考慮故障時間和維修時間,主要采用按照某種概率密度函數(shù)隨機(jī)抽樣系統(tǒng)各元件的狀態(tài),繼而組合成系統(tǒng)的狀態(tài),通過系統(tǒng)分析和穩(wěn)定性判定[7-8]。

設(shè)變量x的期望值μ=E(x),x按照某種概率密度分布函數(shù)隨機(jī)抽樣獲得互相獨立的樣本空間為x=(x1,x2,…,xN),計算平均值為

(9)

根據(jù)大數(shù)定律可知

(10)

(11)

由中心極限定律可知,隨著樣本數(shù)N的增加,樣本均值的概率分布趨于正態(tài)分布N(μ,δ2),即樣本均值在期望值附近的概率最大。則樣本均值方差為

(12)

定義方差系數(shù)β作為收斂精度和條件為

(13)

2.2 交叉熵算法

在大多數(shù)情況下,小擾動穩(wěn)定指標(biāo)ISSI=P{ξ<0}屬于稀有事件的概率,采用傳統(tǒng)非序貫蒙特卡洛算法進(jìn)行無偏估計需要大量樣本才能獲得足夠的計算精度。

交叉熵法作為重要抽樣法的一種,其核心是構(gòu)造一個重要抽樣概率密度函數(shù)g(x)代替原概率密度函數(shù)f(x),計算均值與期望值的距離,以交叉熵距離最小為收斂依據(jù),表示g(x)與最優(yōu)gopt(x)越相似,g(x)具體構(gòu)造過程如下[9-11]：

理論上,系統(tǒng)失穩(wěn)指標(biāo)的期望值為

(14)

式中:f(x;α)----隨機(jī)抽樣概率密度函數(shù);

L(x)----狀態(tài)變量x對應(yīng)的系統(tǒng)失穩(wěn)測度函數(shù),滿足以下關(guān)系:

(15)

構(gòu)造的重要抽樣密度分布函數(shù)g(x;ν)需滿足,

Eg[L(x)W(x;α,ν)],

(16)

則式(16)的無偏估計為

(17)

式中:xi----通過概率密度函數(shù)g(x;ν)隨機(jī)抽取的獨立同分布樣本。

(18)

然而,實際值l是未知量,無法直接得到最優(yōu)重要抽樣密度函數(shù),因此在尋優(yōu)階段,通過迭代減小待更新的重要抽樣概率密度函數(shù)g(x;ν)和理論上最優(yōu)的重要抽樣概率密度函數(shù)gopt(x;ν)的交叉熵來尋找到近似最優(yōu)的重要抽樣概率密度函數(shù)。具體為

(19)

(20)

將式(18)代入式(20)可得

(21)

在重要抽樣概率密度的交叉熵尋優(yōu)階段,按照式(20)進(jìn)行隨機(jī)抽樣,式(17)進(jìn)行小擾動失穩(wěn)指標(biāo)的無偏估計。

3 基于交叉熵蒙特卡洛的含風(fēng)電場電力系統(tǒng)小擾動失穩(wěn)判定

通過上述分析可知,基于交叉熵蒙特卡洛算法的核心是構(gòu)造以交叉熵距離最小為依據(jù)構(gòu)造重要抽樣概率密度函數(shù)g(x)替代原初始抽樣概率密度函數(shù)f(x),則文中的隨機(jī)小擾動考慮風(fēng)電功率波動和負(fù)荷變化,兩者的初始隨機(jī)抽樣概率密度函數(shù)構(gòu)造有以下幾點。

3.1 風(fēng)電功率隨機(jī)抽樣

考慮因風(fēng)速隨機(jī)波動引起風(fēng)電功率輸出隨機(jī)波動作為隨機(jī)小擾動。由于雙參數(shù)Weibull分布的曲線形狀與現(xiàn)實狀況很匹配,被廣泛用來描述風(fēng)速的分布。風(fēng)速波動可以采用雙參數(shù)Weibull分布模擬風(fēng)速條件概率密度為[13]

(22)

式中:v----風(fēng)速;

k----形狀系數(shù);

c----尺度系數(shù)。

風(fēng)力發(fā)電機(jī)組一般工作在停車、最大風(fēng)能跟蹤和恒功率輸出三個階段。為了簡化計算,風(fēng)力發(fā)電機(jī)在不同風(fēng)速下的有功輸出采用近似功率曲線直接轉(zhuǎn)換風(fēng)速值的方法進(jìn)行計算,近似計算模型為[14]

(23)

式中:ν、νin、νout、νN、PN----分別為風(fēng)力發(fā)電機(jī)的實時、切入、切出、額定風(fēng)速以及額定功率。

綜合式(1)和式(2)可得每臺發(fā)電機(jī)的風(fēng)電有功輸出的抽樣值?？紤]同一風(fēng)電場的風(fēng)速基本相同,將每臺發(fā)電機(jī)的有功輸出值求和,可得該風(fēng)電場的有功輸出隨機(jī)抽樣值。

3.2 負(fù)荷變化隨機(jī)抽樣

實際系統(tǒng)中節(jié)點負(fù)荷是一個時刻波動的不確定變量,負(fù)荷的變化對系統(tǒng)潮流影響比較大,而電力系統(tǒng)小擾動穩(wěn)定與網(wǎng)絡(luò)潮流的初始解有關(guān),因此在小干擾穩(wěn)定評估中負(fù)荷隨機(jī)性越大,系統(tǒng)的穩(wěn)定性也就越差。

節(jié)點負(fù)荷變化的隨機(jī)性也是影響電力系統(tǒng)小干擾穩(wěn)定的主要因素之一,所以文中將負(fù)荷變化作為隨機(jī)小擾動。負(fù)荷功率概率模型一般采用服從正態(tài)分布的概率密度函數(shù)進(jìn)行模擬,

(24)

式中:u----負(fù)荷功率的期望值,取統(tǒng)計數(shù)據(jù)的均值;

σ2----方差,取均值的2%。

3.3 小擾動失穩(wěn)判定計算流程

設(shè)定初始抽樣的樣本數(shù)目Nk,收斂精度γ,按照以收斂精度結(jié)束仿真的條件,得出基于交叉熵蒙特卡洛的含風(fēng)電場電力系統(tǒng)小擾動概率穩(wěn)定性判定的具體流程,如圖1所示。

圖1 小擾動失穩(wěn)判定計算流程

4 算例仿真與分析

為了驗證本算法風(fēng)電功率隨機(jī)波動對系統(tǒng)穩(wěn)定性影響的有效性,在3機(jī)9節(jié)點系統(tǒng)的基礎(chǔ)上并上一個風(fēng)電場[15]作為穩(wěn)定性測試系統(tǒng),其總裝容量為100 MW,其余為300 MW的常規(guī)發(fā)電機(jī),采用標(biāo)幺值計算且基準(zhǔn)值均取100 MVA,如圖2所示。

圖2 含風(fēng)電場的3機(jī)9節(jié)點電力系統(tǒng)

具體Matlab仿真結(jié)果與分析見表2。

表2 各常規(guī)發(fā)電機(jī)有功輸出概率模型

兩臺常規(guī)發(fā)電機(jī)有功輸出均采用表2中降額概率模型,輸電線路采用兩狀態(tài)概率模型。設(shè)定切入風(fēng)速、額定風(fēng)速以及切除風(fēng)速分別為3.5 m/s,13 m/s和30 m/s,按照式(22)和式(24)分別對風(fēng)速和負(fù)荷波動進(jìn)行隨機(jī)抽樣。特征值采用QR算法,潮流計算采用牛頓拉夫算法。

通過風(fēng)速隨機(jī)抽樣概率密度函數(shù)計算得出風(fēng)電功率抽樣的直方圖及核密度估計如圖3所示。

圖3 風(fēng)電輸出功率概率直方圖

從圖3可以看出,風(fēng)電場零功率輸出和滿額輸出的概率最大,且風(fēng)電場的功率輸出多集中在0.2～0.4 pu之間,能起到隨機(jī)小擾動的作用。

CEMCS與MCS在小擾動失穩(wěn)指數(shù)(ISSI)和不同收斂精度下的結(jié)果見表3。

表3 小擾動失穩(wěn)指數(shù)計算及仿真性能表

由表3可知,隨著設(shè)定方差系數(shù)的減小,兩個算法計算的小擾動失穩(wěn)指數(shù)均在減小,說明ISSI的計算精度都有所提高,且CEMCS算法的計算精度略高于MCS算法。隨著收斂方差系數(shù)的減小,傳統(tǒng)MCS算法需要大量的抽樣次數(shù)才可以獲得相應(yīng)的ISSI計算精度,而CEMCS算法的抽樣次數(shù)和計算時間變化不大。橫向比較可知,同一方差系數(shù)下,CEMCS算法的抽樣次數(shù)和計算時間均優(yōu)于MCS,達(dá)到減方差的效果。

設(shè)定抽樣次數(shù)為1 000次,CEMCS和MCS在此抽樣次數(shù)下的小擾動失穩(wěn)指標(biāo)及其方差系數(shù)收斂過程分別如圖4和圖5所示。

圖4 小擾動失穩(wěn)指標(biāo)收斂過程

圖5 小擾動失穩(wěn)指標(biāo)方差系數(shù)收斂過程

從圖4可以看出,MCS得到的小擾動失穩(wěn)指標(biāo)在抽樣400次以內(nèi)的波動非常大,而CEMCS得到小擾動失穩(wěn)指標(biāo)快速趨于平穩(wěn)。圖5中CEMCS的方差系數(shù)最小值是0.017,而MCS的方差系數(shù)最小值為0.081。由此可知,在相同的抽樣次數(shù)下,CEMCS的方差系數(shù)小于MCS的方差系數(shù),即在收斂精度方面CEMCS較優(yōu)。在收斂速度方面,CEMCS在抽樣400次左右已趨于最優(yōu)值,而傳統(tǒng)算法近1 000次迭代才達(dá)到最優(yōu)值。

5 結(jié) 語

將風(fēng)速引起的風(fēng)電功率波動和負(fù)荷波動作為隨機(jī)擾動源,利用特征值時域分析法確定小擾動失穩(wěn)指標(biāo),采用非序貫蒙特卡洛分析方法進(jìn)行概率模擬評估。引入交叉熵降低失效事件的稀疏性,提高蒙特卡洛算法的收斂精度和計算效率。通過算例仿真結(jié)果顯示,交叉熵非序貫蒙特卡洛算法在小擾動穩(wěn)定性評估中具有一定的有效性和優(yōu)越性。