段長宏, 李 升, 馬思想

(1.南京工程學院 電力工程學院,南京 211167; 2.中廣核新能源南通有限公司,江蘇 南通 210019)

近年來,電力系統(tǒng)的隨機性明顯增加,在“雙碳”的目標下,新型電力系統(tǒng)將會呈現(xiàn)“雙高”與“雙隨機”等特點,電力系統(tǒng)中的隨機性必不可忽視[1-2]。隨著科學的發(fā)展,太陽能、風能等新能源發(fā)電與電動汽車并網(wǎng)給電力系統(tǒng)帶來不斷的隨機擾動,現(xiàn)代互聯(lián)電網(wǎng)規(guī)模增加,由負荷、故障帶來的隨機擾動也更加普遍,傳統(tǒng)的電力系統(tǒng)確定性分析法已不再適用,因此就必須從隨機穩(wěn)定性方面去分析隨機激勵對電力系統(tǒng)的影響[3]。

電力系統(tǒng)的穩(wěn)定包括功角、電壓和頻率的穩(wěn)定。功角穩(wěn)定又分為小擾動與大擾動穩(wěn)定。其中,當小擾動足夠小的時候,可以直接采用Lyapunov第一法進行分析。對于大擾動穩(wěn)定分析主要用數(shù)值方法與能量函數(shù)法。文獻[4]用擴展面積法將多機系統(tǒng)等效為兩群系統(tǒng),建立起隨機微分方程,基于擬哈密頓隨機平均法(quasi-Hamiltonian system and stochastic averaging,QHSSA)得出多機系統(tǒng)電力系統(tǒng)能量函數(shù),分析隨機電力系統(tǒng)暫態(tài)穩(wěn)定性,并與蒙特卡洛模擬法相比有較好的一致性。關(guān)于頻率的穩(wěn)定,文獻[5]提出以域內(nèi)概率作為隨機激勵對系統(tǒng)頻率動態(tài)安全性影響的量化評估指標,并通過后向Kolmogorov方程,提出了快速求解域內(nèi)概率的半解析法。隨機擾動又在隨機理論中稱為“隨機激勵”,初值與參數(shù)在電力系統(tǒng)運行中是固定不變的,外部激勵則是時變的。文獻[6]基于單機系統(tǒng),構(gòu)建了帶有隨機激勵功率項的隨機微分方程,定量分析不同的隨機激勵強度,激勵步長與仿真計算步長對系統(tǒng)功角的影響。文獻[7]基于隨機微積分理論,提出了一種新的隨機暫態(tài)穩(wěn)定性評估框架。文獻[8]在確定型非線性系統(tǒng)的基礎(chǔ)上引入隨機激勵,利用EM數(shù)值算法進行迭代,并通過均值均方穩(wěn)定性證明了隨機小擾動下的系統(tǒng)功角與轉(zhuǎn)速的穩(wěn)定性。文獻[9]提出隨機動力學模型的線性化隨機激勵強度最優(yōu)閾值模型,得出了線性化最優(yōu)閾值的條件。文獻[10]提出改進EEAC法,將多機系統(tǒng)的隨機微分方程模型通過同調(diào)分群等值為單機無窮大模型,推導出含伊藤積分的加減速面積,求出不同隨機激勵強度的極限切除時間。文獻[11]分析考慮在隨機激勵影響下單機無窮大的系數(shù)矩陣特征值的3種不同情況,并證明了這3種情況下均值均差的界與隨機激勵強度與系統(tǒng)相應(yīng)參數(shù)之間的聯(lián)系。

前述文獻都是將外部激勵近似為單一的高斯型激勵,但電力系統(tǒng)的隨機性是復(fù)雜的。本文將連續(xù)且平穩(wěn)的隨機激勵視為Gauss過程,將負荷沖擊,機組投切等離散型偶發(fā)隨機事件近似為復(fù)合泊松過程,利用確定型非線性系統(tǒng)的基本思路和方法,在此基礎(chǔ)上增加高斯與泊松混合激勵而確立為帶跳隨機微分方程,利用Heun數(shù)值方法進行求解,并通過隨機穩(wěn)定性分析方法證明高斯與泊松混合激勵下的電力系統(tǒng)的均值、均方穩(wěn)定性,給出了均值、均方的界與跳躍激勵強度和隨機激勵強度及系統(tǒng)參數(shù)之間的關(guān)系,討論相同隨機激勵強度下不同的跳躍激勵強度與相同跳躍激勵強度下不同隨機激勵強度對多機電力系統(tǒng)的功角的影響。最后,通過多機系統(tǒng)(4機11節(jié)點系統(tǒng))的功角和轉(zhuǎn)速響應(yīng)曲線與系統(tǒng)均值、均方的界進行分析比較,驗證了本文提出的理論的合理性。

1 帶跳隨機微分方程模型的建立

Gauss白噪聲一般用來近似平穩(wěn)、連續(xù)的外部隨機激勵,但電力系統(tǒng)中發(fā)生的功率突變、機組投切等沖擊突發(fā)事件時用Poisson白噪聲來近似更加的合適?？紤]Gauss白噪聲與Poisson白噪聲擾動的n維向量動態(tài)系統(tǒng),其帶跳隨機微分運動方程[12]如下

dX(t)=F[X(t),t]dt+G[X(t),t]dB(t)+
H[X(t),t]dC(t),X(t0)=X0

(1)

式(1)的積分形式為

(2)

式中:X(t0),B(t)與C(t)在正文中,是相互獨立的,X(t)=[X1(t),X2(t),…,Xn(t)]T為n維矢量隨機狀態(tài)變量,在{B(t),t∈T}隨機過程中,在時間段t1

(3)

復(fù)合Poisson過程形式導數(shù)具有Poisson白噪聲特征即dC(t)=WH(t)dt。{N(t),t∈T}為服從參數(shù)λ的泊松過程,表示0～T時間內(nèi)發(fā)生事件的次數(shù),任意s,t≥0,N(t+s)-N(t)～P(λt)。C(t)為獨立增量過程,{Yk,k=1,2,…,n}為一獨立同分布的隨機變量,且與{N(t),t≥0}為單位階躍函數(shù)。Yk與tk分別為隨機事件發(fā)生時的脈沖強度與發(fā)生事件的時刻。

2 Heun數(shù)值算法

隨機微分方程的解析解很難得出,一般都是在特殊情形下才會有解析解,如今的隨機微分方程的研究大部分使用數(shù)值解法。如今的隨機微分方程比較成熟的數(shù)值解法有EM(Euler-Maruyama)法、Milstein法、Heun法、和Runge-Kutta法[14],EM法是其中最簡單的,也是使用最多的,但其收斂階較低,數(shù)值穩(wěn)定性較低,本文采用收斂階較高的Heun數(shù)值算法是基于梯形公式對EM算法進行改進。對式(1)進行時間離散化可近似為

(4)

Heun數(shù)值迭代格式

(5)

圖1 標準Gauss白噪聲與Poisson白噪聲路徑Fig.1 Standard Gauss white noise and Poisson white noise paths

圖2 布朗運動路徑與包含復(fù)合泊松的布朗運動路徑Fig.2 Brownian motion path and Brownian motion path with compound Poisson

3 混合激勵下多機系統(tǒng)的隨機穩(wěn)定性分析

在確定性的情況下,第i臺發(fā)電機的轉(zhuǎn)子運動方程為

(6)

(7)

(8)

(9)

δij=δi-δj

(10)

式中:Mi為第i臺發(fā)電機的慣性;常數(shù)δij為發(fā)電機功角;ωi為發(fā)電機轉(zhuǎn)速;ω0為初始角速度;t為時間;Di為阻尼系數(shù);Pei為電磁功率;n為發(fā)電機得總臺數(shù);δij0為第i臺發(fā)電機與第j臺發(fā)電機之間的初值功角差;Gii為第i臺發(fā)電機的自電導;Yij為第i臺與第j臺發(fā)電機之間的互導納;αij第i臺與第j臺發(fā)電機之間阻抗角的余角;Pmi為機械功率數(shù)值等于Pei的穩(wěn)態(tài)值[15]。

在式(6)的基礎(chǔ)上增加隨機激勵(由Gauss過程驅(qū)動)和泊松激勵(由復(fù)合Poisson過程驅(qū)動),則運動方程變?yōu)?/p>

(11)

(12)

式中:σi為第i臺發(fā)電機受到高斯激勵WG的強度;γi為第i臺發(fā)電機受到泊松激勵的強度。式(12)則為受高斯和泊松混合激勵的發(fā)電機轉(zhuǎn)子運動方程。

對式(9)減去式(8)可得

(13)

(14)

(15)

其中

式(15)可以寫成

(16)

式(16)的解為

(17)

證明式(16)的均值穩(wěn)定性

(E‖X‖)2≤E[‖X‖2]=E[XT(t)X(t)]

(18)

將式(17)代入式(18)得

根據(jù)范數(shù)的定義與隨機過程的數(shù)字特征推導(具體推導過程見附錄A)可得

(19)

由上述關(guān)于系統(tǒng)式(16)的均值、均方穩(wěn)定性的證明可以判斷,當滿足

(20)

(21)

4 仿真分析

選用4機11節(jié)點系統(tǒng)來驗證電力系統(tǒng)在高斯與泊松混合激勵下的穩(wěn)定性。系統(tǒng)參數(shù)為M1=13 s,M2=13 s,M3=12.5 s,M4=12.5 s,D1=0.2,D2=0.25,D3=0.3,D4=0.35。其他具體參數(shù)參考文獻[17]。通過計算可得λA=-0.008 79,C0=2.608 25。首先取系統(tǒng)跳躍事件發(fā)生的概率λ=0.02,即每秒發(fā)生跳躍事件的概率為0.02,則150 s內(nèi)發(fā)生3次跳躍事件, 0～150 s內(nèi)發(fā)生3次跳躍的過程,如圖3所示。為便于之后的研究,統(tǒng)一取隨機激勵強度σ=0.01與跳躍激勵強度γ=0.01,將以上參數(shù)代入式(21)時,可得均方的界E[XT(t)X(t)]<1.461 64×10-4。在此激勵強度下的功角曲線與轉(zhuǎn)速曲線,如圖4、圖5所示。通過功角與轉(zhuǎn)速的仿真曲線可以看出,雖然系統(tǒng)受到高斯與泊松混合激勵的影響,但由于強度較小,同時轉(zhuǎn)子具有慣性,功角波動不超過0.5°,轉(zhuǎn)速波動不超0.005,驗證系統(tǒng)具有均值均方穩(wěn)定性。

圖3 跳躍事件發(fā)生時刻與次數(shù)Fig.3 Occurrence time and number of jumping events

圖4 四機系統(tǒng)在σi=0.01,γi=0.01的功角曲線Fig.4 Power-angle curve of the four-machine system at σi=0.01,γi=0.01

圖5 四機系統(tǒng)在σi=0.01,γi=0.01的轉(zhuǎn)速曲線Fig.5 Rotation speed curve of the four-machine system at σi=0.01,γi=0.01

下面定量分析隨機激勵相等的情況下不同跳躍激勵強度對系統(tǒng)功角的影響,為了避免隨機激勵對系統(tǒng)的影響,統(tǒng)一取隨機激勵σi=0.01。分別取8組不同的跳躍激勵強度來計算均值、均方的界限,具體結(jié)果如表1所示。一號機在相同隨機激勵強度不同跳躍激勵強度的功角響應(yīng)曲線圖,如圖6所示。由圖6可知,在跳躍激勵強度取表1第1組數(shù)據(jù)時,即γ=[0.01,0.02,0.03,0.04],跳躍激勵強度較小,功角狀態(tài)曲線波動較小,基本上可以運行在擾動前的平衡狀態(tài)。圖6中點線“-·”是在跳躍激勵強度取第4組數(shù)據(jù)時即γ=[0.01,0.15,0.30,0.45]時的功角狀態(tài)曲線圖,可以看出由于跳躍激勵強度的增加,在75 s,90 s與120 s發(fā)生跳躍事件時,系統(tǒng)功角的峰值出現(xiàn)了比較明顯的增大,隨著跳躍激勵強度的增加,在發(fā)生跳躍事件時功角曲線波動也逐漸增加。而在跳躍激勵強度γ=[0.01,0.7,1.4,2.1]時,仿真結(jié)果如圖6實線“-”所示,在跳躍激勵強度較大時,發(fā)生跳躍事件后功角響應(yīng)曲線波動幅度出現(xiàn)了非常明顯的增大且出現(xiàn)了嚴重的波動,即在跳躍激勵強度增大時,系統(tǒng)均值、均方的上界增大,系統(tǒng)無法穩(wěn)定運行在擾動前的平衡狀態(tài),最后可得系統(tǒng)失穩(wěn)。

表1 相同隨機激勵強度下不同跳躍激勵強度情況時系統(tǒng)式(16)的均值均方的界Tab.1 The bound of mean square of system (16) under the same random excitation intensity and different jumping excitation intensity

圖6 1號機分別取表1中3組數(shù)據(jù)時的功角曲線對比Fig.6 Comparison of power angle curves of machine 1 when three groups of data in Tab.1 are respectively taken

以下分析對單個發(fā)電機跳躍激勵強度發(fā)生變化時對系統(tǒng)狀態(tài)量的影響,取γ=[0.01,0.01,0.01,0.01]與γ=[0.02,0.01,0.01,0.01],僅改變1號機的跳躍激勵強度,將上述跳躍激勵強度代入式(21)可得在單個發(fā)電機激勵強度增加時,式(16)的均值、均方臨界值也會增加。1號機本身跳躍激勵強度增加時的功角曲線對比圖與3號機在僅1號機跳躍激勵強度發(fā)生改變時的功角曲線對比圖,如圖7、圖8所示。由仿真可以驗證,隨機激勵強度不改變時兩種跳躍激勵強度的功角曲線在發(fā)生跳躍事件之前保持一致,在任意一臺發(fā)電機跳躍激勵強度發(fā)生改變時,本身及其他發(fā)電機功角響應(yīng)曲線波動幅度都會增大。

圖8 3號機在1號機跳躍激勵強度γ=0.01與γ=0.02下的功角曲線Fig.8 Power-angle curve curves of machine No.3 under the jumping excitation intensity γ=0.01 and γ=0.02 of machine No.1

下面分析相同的跳躍激勵強度不同的隨機激勵強度下的攻角響應(yīng)曲線,4臺發(fā)電機跳躍激勵強度統(tǒng)一取0.01,同理取8組不同的隨機激勵強度來計算系統(tǒng)式(16)的均值、均方穩(wěn)定的界限,結(jié)果如表2所示。

表2 相同跳躍激勵強度下不同隨機激勵強度情況時系統(tǒng)式(16)的均值均方的界Tab.2 Bounds of mean square of system (16) for different random excitation intensities under the same jump excitation intensities

圖9為相同跳躍激勵強度下不同隨機激勵強度的1號機功角響應(yīng)曲線圖,可以看出在取第1組隨機激勵強度時,仿真曲線如圖9中虛線“--”所示,功角響應(yīng)曲線波動幅度較小,基本上穩(wěn)定在平衡狀態(tài)。隨著隨機激勵強度的增加,仿真如圖9中點線“-·”與實線“-”所示,功角響應(yīng)曲線逐步向下遠離擾動前的平衡點,并逐漸失穩(wěn)。

圖9 1號機分別取表1中3組數(shù)據(jù)時的功角曲線對比Fig.9 Comparison of power angle curves of machine 1 when three groups of data in table 2 are respectively taken

在跳躍激勵強度保持不變的同時,僅改變1號機的隨機激勵強度,即取σ=[0.01,0.01,0.01,0.01]與σ=[0.02,0.01,0.01,0.01]兩組數(shù)據(jù)進行仿真,結(jié)果如圖10～圖11所示。在75 s第一次跳躍事件發(fā)生前,兩種激勵強度下的1號機與3號機功角曲線相差較小,但隨著跳躍事件的發(fā)生,即擾動的疊加,使兩種激勵強度的功角曲線開始逐漸偏離,相差變大。同時,由仿真也可得出,在僅一臺發(fā)電機隨機激勵強度發(fā)生改變時,其他發(fā)電機的功角響應(yīng)曲線變化與此臺發(fā)電機功角響應(yīng)曲線變化相似。

圖10 1號機在跳躍激勵強度σ=0.01與σ=0.02下的功角曲線Fig.10 Work angle curves of machine 1 under jump excitation intensity σ=0.01 and σ=0.02

圖11 3號機在1號機跳躍激勵強度σ=0.01與σ=0.02下的功角曲線Fig.11 Power angle curves of machine 3 under the jumping excitation intensity σ=0.01 and σ=0.02 of machine 1

5 結(jié) 論

(1)在確定型非線性電力系統(tǒng)模型的基礎(chǔ)上,加入高斯與泊松激勵,建立更加符合電力系統(tǒng)的帶跳隨機微分方程系統(tǒng)模型,然后從系統(tǒng)的均值穩(wěn)定與均方穩(wěn)定性方面來驗證系統(tǒng)的穩(wěn)定性,并給出了跳躍激勵強度γi與隨機激勵強度σi與系統(tǒng)參數(shù)之間的關(guān)系,即在任意一個激勵強度增大時,系統(tǒng)均值、均方的界也會相應(yīng)增加。

(2)對相同跳躍激勵強度下不同隨機激勵強度與相同隨機激勵強度下不同跳躍激勵強度兩種情況對系統(tǒng)進行仿真分析。仿真結(jié)果表明:隨機小激勵強度與跳躍小激勵強度會對系統(tǒng)狀態(tài)變量產(chǎn)生一定影響,但不會造成系統(tǒng)失穩(wěn)。在跳躍激勵與隨機激勵中任一激勵強度增大到一定值時,系統(tǒng)狀態(tài)曲線無法恢復(fù)到平衡狀態(tài),即系統(tǒng)失穩(wěn)。

附錄A

由于Brown運動期望為0,且初值X(t),B(t),C(t)相互獨立。則

由隨機積分的性質(zhì)可得

復(fù)合Poisson性質(zhì)可得

非隨機變量的期望等于本身

E[(eAtX0)TeAtX0]=‖eAtX0‖2,

h[X(t)]=

其中

因此‖h[X(t)]‖2≤T1,T1為大于零的常數(shù)。

由矩陣2范數(shù)可得‖eAt‖2≤C0e2λAt,λA為矩陣A最大特征值實部,C0為大于0的常數(shù)。

式中:‖X0‖2=T0,‖K‖2=T2,‖Q‖2=T3;T0,T2,T3均為大于0的常數(shù)。

在系統(tǒng)穩(wěn)定時,系數(shù)矩陣A所有特征值實部均小于0,因此λA<0,所以