考慮高斯激勵的多機電力系統(tǒng)的p階隨機穩(wěn)定性分析

秦子璇， 盧占會， 魏軍強

(華北電力大學 數(shù)理學院，北京 102206)

0 引 言

當前全球能源系統(tǒng)面臨的主要問題有能源短缺、環(huán)境惡化以及能源利用效率較低等。以可再生能源逐步替代化石能源，實現(xiàn)能源轉(zhuǎn)型，建設(shè)清潔低碳、安全高效的新一代能源系統(tǒng)是全世界新一輪能源革命的核心目標[1]。與廣泛使用的傳統(tǒng)能源相比，新能源是指基于新技術(shù)開發(fā)和利用的非常規(guī)能源，包括風能、太陽能、生物質(zhì)能、地熱能和海洋能等。近年來，新能源的使用量一直在快速增長，由于新能源發(fā)電的間歇性和隨機波動性，以及新能源接入系統(tǒng)可能出現(xiàn)的低頻和高頻問題等，在將新能源大規(guī)模集成到電網(wǎng)中時，都可能會對電網(wǎng)的穩(wěn)定性產(chǎn)生重大影響[3-5]。例如，風力與光伏發(fā)電，都與風和光的大小和強弱有關(guān)，當大規(guī)模風電站集成到電網(wǎng)中時，由于風能的隨機性，風電將隨著風速的變化而隨時變化，這將影響電能的實時平衡，進而影響電網(wǎng)的安全穩(wěn)定運行。

國內(nèi)外許多學者對電力系統(tǒng)的隨機穩(wěn)定性進行了一些研究，研究成果見諸于不同文獻[6-16]。其中，文獻[6]基于隨機微積分理論，提出了一種新的隨機暫態(tài)穩(wěn)定性評估框架。文獻[7]對擬哈密頓系統(tǒng)利用隨機平均法，對阻尼系數(shù)、激勵強度等對隨機系統(tǒng)的暫態(tài)穩(wěn)定性的影響進行分析。文獻[8]針對單機無窮大系統(tǒng)下的擬哈密頓系統(tǒng)的隨機平均法的一次推廣，同時將單機系統(tǒng)的理論推廣到了多機系統(tǒng)，得到了隨機激勵下的多機系統(tǒng)能量函數(shù)的一維擴散過程公式。文獻[9]通過隨機Lyapunov穩(wěn)定性分析，將電力系統(tǒng)瞬態(tài)能量函數(shù)法用于隨機電力系統(tǒng)的穩(wěn)定性。文獻[10]提出了電力系統(tǒng)的均值和均方穩(wěn)定性的概念，1階穩(wěn)定性為均值穩(wěn)定，2階穩(wěn)定性為均方穩(wěn)定，均值穩(wěn)定意味著，一個隨機輸入所激發(fā)的系統(tǒng)響應(yīng)的均值是有界的；均方穩(wěn)定意味著，一個隨機輸入所激發(fā)的系統(tǒng)響應(yīng)的均方差是有界的，并采用Euler-Maruyama(EM)數(shù)值方法得到電力系統(tǒng)功角曲線的軌跡，同時分析了不同激勵強度下功角曲線的不同。文獻[11] 利用Heun數(shù)值方法對單機無窮大系統(tǒng)進行了仿真模擬，分別分析了計算步長、激勵強度以及激勵步長對功角曲線的影響。文獻[12]采用It公式和It等距法推導出系統(tǒng)狀態(tài)方程的解析解，同時提出利用矩穩(wěn)定性的定義來描述隨機穩(wěn)定性，給出了隨機激勵下系統(tǒng)小擾動的穩(wěn)定性判據(jù)，并通過雙機無窮大系統(tǒng)的仿真驗證了該方法的有效性。文獻[13-14]研究了區(qū)間穩(wěn)定性，利用Lyapunov函數(shù)和矩陣不等式等數(shù)學工具，分析了異步風力機的小干擾區(qū)間模型，給出了系統(tǒng)的區(qū)間穩(wěn)定性。文獻[15]利用Lyapunov方法和其它相關(guān)理論，證明了一類隨機線性系統(tǒng)的p階穩(wěn)定性定理，其中3階穩(wěn)定性和4階穩(wěn)定性的意義與偏度和峰度有關(guān)，峰度值越大，概率密度曲線的尾端越厚，這意味著極值出現(xiàn)的概率更大；相反，如果峰度值很小，則概率分布集中，這意味著極值出現(xiàn)的概率很小，因此隨機系統(tǒng)的穩(wěn)定性更好。但基于隨機微分方程理論研究多機電力系統(tǒng)的p階穩(wěn)定性問題尚未見報道，本文將對多機電力系統(tǒng)在隨機激勵下的p階穩(wěn)定性問題進行研究擴展，擴展思路是先將多機系統(tǒng)進行等效轉(zhuǎn)換，將其等效為單機系統(tǒng)，利用已有理論知識對其穩(wěn)定性進行檢驗。

本文利用隨機微分方程理論，對原多機系統(tǒng)以及由擴展等面積法得到的多機電力系統(tǒng)的等效模型進行線性化處理。選取四機兩區(qū)系統(tǒng)作為實際算例，利用EM數(shù)值方法研究電力系統(tǒng)功角、轉(zhuǎn)速等隨時間變化的趨勢，通過蒙特卡羅法做出等效轉(zhuǎn)換前后電力系統(tǒng)的功角的概率密度函數(shù)曲線，利用其分析等效轉(zhuǎn)換前后系統(tǒng)穩(wěn)定性是否發(fā)生變化。將單機無窮大系統(tǒng)的p階穩(wěn)定性擴展到多機系統(tǒng)的p階穩(wěn)定性。

1 多機電力系統(tǒng)的隨機建模

1.1 多機電力系統(tǒng)的隨機模型

考慮具有高斯型白噪聲的矢量隨機微分方程：

dX(t)=AX(t)dt+QdB(t),t∈[0,+∞)

(1)

式中：X(t)=(X1(t),X2(t),…,Xn(t))T表示n維隨機變量，B(t)=(B1(t),B2(t),…,Bn(t))T表示n維維納過程，dB(t)/dt=W(t)，W(t)為高斯過程，A和Q為常數(shù)矩陣。

(2)

1.2 多機系統(tǒng)的雙機等效隨機模型

(3)

根據(jù)擴展等面積法，將多機系統(tǒng)的失穩(wěn)看作是一臺或是多臺發(fā)電機形成的一個機群相對系統(tǒng)的其余部分失去同步，因此將受擾部分機群稱作S群，其余部分稱為A群，推導出了多機系統(tǒng)的雙機等效模型[7]：

(4)

根據(jù)隨機微分方程理論，對(4)式進行線性化處理，可得：

(5)

2 多機系統(tǒng)的p階穩(wěn)定性

2.1 p階穩(wěn)定性基礎(chǔ)知識

關(guān)于隨機激勵下電力系統(tǒng)的p階穩(wěn)定性，有如下結(jié)論。

定義1[15].如果隨機微分方程(1)的解過程X(t)滿足以下條件：

(6)

則系統(tǒng)(1)具有p階穩(wěn)定性。對于所有的p> 0，其中c> 0，c是常數(shù)。

在上述定義下，可得出系統(tǒng)(1)的穩(wěn)定性條件。

引理1[15]. 對于給定的高斯型隨機激勵下的線性常系數(shù)隨機微分系統(tǒng)(1)，初始值X(0)=X0是有界的并且與維納過程B(t)獨立。 當矩陣A的所有特征值的實部小于零時，隨機微分系統(tǒng)(1)具有p階穩(wěn)定性。

2.2 Euler-Maruyama數(shù)值方法

對于大多數(shù)隨機微分方程來說，其解析解是無法求得的，可以求出解析解的隨機微分方程大多為線性隨機微分方程。一些隨機微分方程雖可求得解析解，但解析解形式復雜，處理起來非常不方便。因此，在實際應(yīng)用中，往往通過數(shù)值積分的方法獲得解過程的軌跡，利用解的軌跡無限逼近精確解即可。

EM 算法是目前最簡單的求解隨機微分方程的數(shù)值解法。設(shè)隨機微分方程的解過程是X(t)，對于某個正整數(shù)N，有Δt=(T-t0)/N，τj=jΔt，j=0,1,2,…,N，Xj=X(τj)。EM數(shù)值方法對應(yīng)的差分迭代格式為

Xj+1=Xj+f(Xj,τj)Δt+G(Xj,τj)ΔWj

(7)

3 案例研究

3.1 四機兩區(qū)系統(tǒng)

四機兩區(qū)電力系統(tǒng)的示意圖如圖1所示[19]。

通過對四機電力系統(tǒng)進行簡化，可得四機兩區(qū)系統(tǒng)的伊藤隨機微分方程為[8]

圖1 四機兩區(qū)電力系統(tǒng)示意圖Fig.1 Schematic diagram of four-machine two-area power system

式中：第一臺發(fā)電機的方程為

(8)

對(8)式進行線性化處理，可得：

(9)

則(9)式的系數(shù)矩陣為

E1E3B13cos(δ10-δ30)-E1E4B14cos(δ10-δ40))。

同理，可以寫出四機兩區(qū)系統(tǒng)中其余三臺發(fā)電機對應(yīng)的伊藤隨機微分方程，同時可寫出其系數(shù)矩陣A2、A3、A4，則四機兩區(qū)系統(tǒng)的系數(shù)矩陣為

3.2 數(shù)值分析

通過對電力系統(tǒng)四機兩區(qū)模型進行EM數(shù)值分析，觀察發(fā)電機功角、發(fā)電機轉(zhuǎn)速在給定初始值的情況下，隨時間的變化的趨勢，對電力系統(tǒng)的p階穩(wěn)定性進行一系列的分析驗證。四機兩區(qū)系統(tǒng)的各個參數(shù)值分別為[8]：σ1=0.05，σ2=0.06，σ3=0.07，σ4=0.08，M1=13 s，M2=13 s，M3=12.35 s，M4=12.35 s，E1=1.03，E2=1.01，E3=1.03，E4=1，P1=0.998 119 479 653 3，P2=0.319 842 206 470 3，P3=0.241 493 636 334 7，P4=1.076 468 049 788 9，B12=B21=1.49，B13=B31=0.62，B14=B41=0.7，B23=B32=0.7，B24=B42=0.76，B43=B34=1.49，D1=0.5，D2=0.55，D3=0.6，D4=0.65，ωN=2·π·60，系統(tǒng)的穩(wěn)定平衡點為δ1 s-δ4 s=37.2°，δ2s-δ4 s=27.5°，δ3s-δ4 s=10.2°。

利用已知參數(shù)，對等效轉(zhuǎn)換后的多機電力系統(tǒng)進行數(shù)值模擬，得到圖2-圖4的三幅圖，分別是系統(tǒng)功角曲線、轉(zhuǎn)速以及500條系統(tǒng)功角模擬曲線的均值與時間t的關(guān)系圖。

通過應(yīng)用EM數(shù)值計算方法對系統(tǒng)(5)進行分析，在xt0=(0,0)T的初始條件下，30 s內(nèi)的運算結(jié)果如圖2和圖3所示。

圖3 ωSA隨時間t的波動圖Fig.3 Movement of ωSA along the time t

從圖4可以看出，在隨機激勵的影響下，功角曲線在初值δSA=0 rad附近波動。在統(tǒng)計學中常用樣本均值來估計總體均值，選擇500條樣本路徑模擬的均值來估計總體的變動情況，其中路徑1和路徑2是500條樣本模擬中的兩條。單獨的每條樣本模擬曲線隨時間的波動都是比較大的，但是500條路徑的樣本均值則波動較小，基本接近穩(wěn)態(tài)初始值δSA=0 rad。

圖4 500條樣本路徑模擬的均值Fig.4 Average angle over 500 paths

利用已有參數(shù)，可以分別計算出等效轉(zhuǎn)換前后兩個模型的系數(shù)矩陣A的特征值。其中，原四機兩區(qū)系統(tǒng)的系數(shù)矩陣A的特征值為

λ1=-0.000 192 307 7+0.882 978 3617j

λ4=-0.000 192 307 7-0.882 978 3617j

λ3=-0.000 211 584 6+0.917 693 818 5j

λ4=-0.000 211 584 6-0.917 693 8185j

λ5=-0.000 242 949 8+0.923 161 262 2j

λ6=-0.000 242 949 8-0.923 161 262 2j

λ7=-0.052 631 587 9

λ8=-0.037 699 111 2

根據(jù)引理1的描述，當隨機微分方程系統(tǒng)的系數(shù)矩陣A的特征值實部均小于零時，可以推斷出系統(tǒng)是p階穩(wěn)定的。因此原四機兩區(qū)系統(tǒng)是p階穩(wěn)定的。

同理可以求出等效轉(zhuǎn)換后的系數(shù)矩陣A的特征值:

λ1=-0.000 226 824 5+0.820 252 071 9j

λ2=-0.000 226 824 5-0.820 252 071 9j

由計算結(jié)果可以看出，系數(shù)矩陣A的特征值實部均小于零。則等效轉(zhuǎn)換前后系統(tǒng)的p階穩(wěn)定性沒有發(fā)生變化。對于多機系統(tǒng)的p階穩(wěn)定性，均可以采用這種方式進行驗證。首先將需要檢驗是否具有p階穩(wěn)定性的多機系統(tǒng)，利用擴展等面積法等[20]，通過進行等效轉(zhuǎn)換的方式，將多機系統(tǒng)轉(zhuǎn)換成雙機或單機系統(tǒng)，其次利用已有定理對等效轉(zhuǎn)換后的系統(tǒng)進行判定，最后根據(jù)等效轉(zhuǎn)換不改變系統(tǒng)穩(wěn)定性這一結(jié)論，可以推斷出原多機系統(tǒng)是否為p階穩(wěn)定。

圖5是等效轉(zhuǎn)換后t=1 s時多機電力系統(tǒng)的功角δ的概率密度函數(shù)圖像，可以看出δ的概率密度函數(shù)曲線是基本符合正態(tài)分布的趨勢。對于3階、4階穩(wěn)定性的判定，通常也可以使用偏度和峰度這一定義[15]。根據(jù)表1中仿真模擬的結(jié)果可以看出，盡管功角δ在各時刻的偏度和峰度不完全接近0和3，但總體的偏度和峰度的模擬效果還是很接近理想狀態(tài)的。這也驗證了前文對于多機系統(tǒng)穩(wěn)定性的判定。

圖5 δ的概率密度函數(shù)Fig.5 Probability density function of δ

表 1 t時刻下δ的偏度和峰度值

3.3 等效轉(zhuǎn)換前后系統(tǒng)功角的概率密度

通過利用蒙特卡羅方法，對等效轉(zhuǎn)換前后的兩個系統(tǒng)在不同時刻下的概率密度函數(shù)的大小進行對比。得到如圖6-圖9所示的，4個不同時刻下的基于等效轉(zhuǎn)換后的系統(tǒng)和原多機系統(tǒng)的四幅圖。對于原始多機系統(tǒng)而言，功角δ是第一臺發(fā)電機功角，對于等效轉(zhuǎn)換后的系統(tǒng)而言為ΔδSA。

圖6 時間t為1 s時δ的概率密度函數(shù)Fig.6 Probability density function of δ when time t is 1 s

由圖6-圖9可以看出，不同時刻下基于原多機系統(tǒng)和轉(zhuǎn)換后系統(tǒng)的功角δ的概率密度函數(shù)曲線的一致性比較好，由此可以判定出等效轉(zhuǎn)化前后系統(tǒng)的穩(wěn)定性基本不會發(fā)生變化，那么多機電力系統(tǒng)的p階穩(wěn)定性問題，可以通過等效轉(zhuǎn)換的方式轉(zhuǎn)換為單機或雙機系統(tǒng)進行驗證。

圖7 時間t為2 s時δ的概率密度函數(shù)Fig.7 Probability density function of δ when time t is 2 s

圖8 時間t為5 s時δ的概率密度函數(shù)Fig.8 Probability density function of δ when time t is 5 s

圖9 時間t為10 s時δ的概率密度函數(shù)Fig.9 Probability density function of δ when time t is 10 s

4 結(jié) 論

本文基于隨機微分方程理論、EM數(shù)值方法以及蒙特卡羅法對多機電力系統(tǒng)的p階穩(wěn)定性進行了分析。在多機電力系統(tǒng)的系數(shù)矩陣形式復雜或參數(shù)不易獲取等情況下，可以將其等效成單機或雙機模型，利用p階穩(wěn)定性定理對等效之后的模型的系數(shù)矩陣進行分析。通過對實際四機兩區(qū)系統(tǒng)的EM數(shù)值模擬的結(jié)果來看，電力系統(tǒng)是p階穩(wěn)定的。根據(jù)蒙特卡羅法，對電力系統(tǒng)的功角的概率密度函數(shù)曲線分析，發(fā)現(xiàn)等效轉(zhuǎn)換前后電力系統(tǒng)的功角穩(wěn)定性不發(fā)生變化。未來的研究工作或?qū)⒓杏诳紤]乘性噪聲下的電力系統(tǒng)的p階穩(wěn)定性。