基于綜合算法的隨機潮流計算方法

孟代江 張淑榮

同濟大學(xué) 電子與信息工程學(xué)院 上海 201804

1 研究背景

隨著風(fēng)能和太陽能等清潔能源在電力系統(tǒng)中發(fā)電占比快速增長,影響電力系統(tǒng)運行狀態(tài)的不確定性大大增加,電力系統(tǒng)穩(wěn)定運行面臨著更大的挑戰(zhàn)。隨機潮流計算能夠在系統(tǒng)各節(jié)點功率注入,系統(tǒng)參數(shù)具有不確定性條件下,獲得電壓和相角等系統(tǒng)狀態(tài)量信息的概率分布,是解決電力系統(tǒng)不確定性這一難題的有效方法。隨機潮流計算方法的研究近年來受到廣泛關(guān)注。

目前,隨機潮流的計算方法主要有模擬法、解析法、代理模型類方法。

模擬法[1-4]只要保證足夠的采樣次數(shù),即可達到很高的精度,通常作為其它方法的基準方法,但收斂速度較慢。

解析法包括半不變量法[5-6]、點估計法[7-8]等,不需要抽樣,通過對原函數(shù)的局部線性化或利用某些特征量的性質(zhì),直接求取輸出變量的期望和方差,雖然計算效率高,但是精度較低,適用范圍有限。

隨機潮流問題的原函數(shù)通常是隱含在復(fù)雜非線性潮流方程中的隱函數(shù),求解速度相對較慢。對此,代理模型類方法用某種特定結(jié)構(gòu)的正交多項式逼近函數(shù),構(gòu)造隨機潮流模型的代理模型,進而基于逼近函數(shù)進行大量抽樣計算,輸出變量的概率特征,也可以通過基函數(shù)系數(shù)直接計算輸出變量的期望、方差等概率特征值。代理模型類方法又可以分為隨機伽遼金方法和隨機配置點法。

隨機伽遼金方法用逼近函數(shù)代替各輸出變量,代入潮流方程,通過正交多項式性質(zhì)求解方程,獲得各基函數(shù)系數(shù)。這種方法通過選擇合理的基函數(shù)及展開階次,可以獲得很高的精度,但是基函數(shù)項由于隨機輸入維數(shù)的增加而迅速增加,方程也更加復(fù)雜,不適用高維系統(tǒng)的求解[9-10]。

隨機配置點法可以分為偽譜類方法和插值類方法。偽譜類方法基于數(shù)值積分及二次方逼近理論,能夠精確計算任意項基函數(shù)系數(shù),但是由于基函數(shù)項數(shù)隨維度指數(shù)的上升而會導(dǎo)致產(chǎn)生維數(shù)災(zāi)難。文獻[11-12]提出自適應(yīng)稀疏偽譜法,在基函數(shù)項擴展時,優(yōu)先擴展對于提高精度作用大的基函數(shù)項,大幅減少基函數(shù)項的數(shù)量。文獻[13]提出更精確的逼近誤差估計,提高計算精度,并提出多輸出同時擴展計算的策略。以上文獻雖然在一定程度上改善了維數(shù)災(zāi)難問題,但是并未根本解決,在高維場景中仍無法高效率應(yīng)用。

文獻[14-16]提出一種插值型方法——低階多項式逼近法,通過構(gòu)造由若干階次的低次多維正交多項式基函數(shù)的線性組合來逼近原函數(shù),通過交替最小二乘法逐維求解基函數(shù)系數(shù)。這一方法的基函數(shù)項數(shù)與隨機輸入變量維數(shù)呈現(xiàn)線性關(guān)系,即所需抽樣點與輸入維數(shù)呈線性增長關(guān)系,避免維數(shù)災(zāi)難問題,但在低維時的計算精度與自適應(yīng)稀疏偽譜法相差甚遠。

筆者針對代理模型類方法展開研究,分析自適應(yīng)稀疏偽譜法和低階多項式逼近法的計算時間復(fù)雜度,基于兩種方法的優(yōu)勢提出一種基于綜合算法的隨機潮流計算方法。這一方法在低、中、高維度場景中都有較高的計算精度和計算效率,通過算例驗證了這一方法的有效性和可行性。

2 隨機潮流數(shù)學(xué)模型

隨機潮流模型可以簡寫為:

U=f(W)

(1)

U=(U1,U2,…,UD)

W=(W1,W2,…,Wd)

式中:U為節(jié)點電壓、相角等系統(tǒng)狀態(tài)隨機變量;D為輸出隨機變量數(shù)量;W為風(fēng)光電場出力及負荷等輸入隨機變量;d為輸入隨機變量數(shù)量;f表示潮流方程隱含的關(guān)于系統(tǒng)狀態(tài)與節(jié)點注入功率的隱函數(shù)關(guān)系。

考慮風(fēng)光出力的不確定性通常直接由風(fēng)速、光照強度等因素決定,隨機潮流模型又可以寫為:

U=f(v,r,PL)

(2)

式中:v為各風(fēng)電場的風(fēng)速向量,一般認為服從韋布爾分布;r為各光伏場的光照強度向量,一般認為服從貝塔分布;PL為各節(jié)點負載有功功率,一般認為服從正態(tài)分布。

3 自適應(yīng)稀疏偽譜法計算時間復(fù)雜度

對于隨機潮流模型,自適應(yīng)稀疏偽譜法構(gòu)造的單輸出多項式逼近函數(shù)為:

(3)

(4)

式中:Qk(〈φn,f〉)為待計算的高斯積分;φn為正交基函數(shù);ck為網(wǎng)格系數(shù)。

自適應(yīng)稀疏偽譜法的計算時間主要花費在對原函數(shù)的抽樣上,對自適應(yīng)稀疏偽譜法計算時間復(fù)雜度進行分析,主要取決于自適應(yīng)稀疏偽譜法積分點數(shù)量的分析。

K={k:d≤k1+k2+...+kd≤l+d}

(5)

這一形式的指標集包含了所有總階次不超過l的正交基函數(shù)。文獻[17]簡單推導(dǎo)了當階次l為0～7時指標集所需總積分點數(shù)n與輸入隨機變量維數(shù)d的函數(shù)關(guān)系,見表1。

表1 總積分點數(shù)與輸入維數(shù)關(guān)系

根據(jù)文獻[18],對于隨機潮流問題,當階次l為3時,使用經(jīng)典Smolyak指標集的稀疏偽譜法即達到最高逼近精度。繼續(xù)增大l,稀疏偽譜法的逼近精度幾乎不變。

由此,可以估算基于自適應(yīng)稀疏偽譜法的隨機潮流計算時間復(fù)雜度為O(Dd3),輸出隨機變量數(shù)量一般是輸入隨機變量數(shù)量的一到兩倍,即D為cd,c為1～2。綜合以上分析,自適應(yīng)稀疏偽譜法的隨機潮流計算時間復(fù)雜度為O(cd4)。

另外,對于隨機潮流問題,通常認為10 000點拉丁超立方法可以獲得隨機潮流輸出變量各階矩的精確值。根據(jù)自適應(yīng)稀疏偽譜法所需采樣點數(shù)量與基函數(shù)階次、輸入維度的關(guān)系,以及基準拉丁超立方法所需采樣點數(shù),采取如下策略獲得高維問題和低維問題的分界值。由于自適應(yīng)稀疏偽譜法的指標集形式靈活,指標集可以位于階次為2和階次為3的經(jīng)典指標集之間?？紤]到反映潮流方程的非線性特性,基函數(shù)總次數(shù)一般需大于2。由此,根據(jù)表1,并令自適應(yīng)稀疏偽譜法的積分點數(shù)下限2d2+2d+1的值為10 000,可求得輸入維數(shù)為70.2。據(jù)此,可認為輸入維數(shù)不小于70時屬于高維輸入情況,此時自適應(yīng)稀疏偽譜法計算效率低于拉丁超立方法,應(yīng)采用無維數(shù)災(zāi)難的求解方法求解隨機潮流。反之,輸入維數(shù)小于70時屬于低維輸入情況,此時自適應(yīng)稀疏偽譜法的計算效率可能更高,且由于自適應(yīng)稀疏偽譜法計算精度上限較低階多項式逼近法更高,可嘗試先使用自適應(yīng)稀疏偽譜法進行隨機潮流求解,若計算精度不滿足要求,則再轉(zhuǎn)向低階多項式逼近法或拉丁超立方法進行求解。

4 低階多項式逼近法計算時間復(fù)雜度

低階多項式逼近法構(gòu)建的單輸出逼近函數(shù)表達式為:

(6)

低階多項式逼近法的計算分為兩部分,一部分是對原函數(shù)的抽樣,另一部分是采用交替最小二乘算法求解基函數(shù)系數(shù)。假設(shè)抽樣一次的時間為teval,則抽樣時間tsample為Nteval,N為抽樣次數(shù)。低階多項式逼近法采用交替最小二乘方法求解基函數(shù)系數(shù)的過程還可以細分為三部分。

(1) 計算凍結(jié)系數(shù),這部分的總計算時間tci為:

tci=DrmaxdN(d-1)(pmax+1)tbasis

(7)

(2) 計算最小二乘,這部分的總計算時間t1s為:

t1s=Drmaxdtsolve

(8)

(3) 計算收斂判據(jù),這部分的總計算時間tErr為:

tErr=DrmaxNd(pmax+1)tbasis

(9)

式中:tbasis為計算一個單維正交基函數(shù)在一個點處的值所需要的計算時間;tsolve為求解一次最小二乘問題所需要的時間。

主要計算時間還是在最小二乘計算上,rmax一般取3～10,輸出隨機變量數(shù)量一般是輸入隨機變量數(shù)量的一到兩倍,因此低階多項式逼近法的隨機潮流計算復(fù)雜度約為O(d2)。

5 基于綜合算法的隨機潮流計算過程

基于綜合算法的隨機潮流計算步驟如下。

(1) 給定隨機潮流系統(tǒng)函數(shù)X=f(W),若d<70,進行步驟(2);反之,進行步驟(4)。

(2) 設(shè)定最大抽樣次數(shù)Nmax,誤差收斂閾值為TOL。

(3) 采用自適應(yīng)稀疏偽譜法求解隨機潮流。若全局誤差η滿足小于TOL,且抽樣次數(shù)N滿足小于Nmax,則進行步驟(8);反之,進行步驟(4)。

(4) 比較低階多項式逼近法估計計算時間t1和t2,若t1≤t2,則進行步驟(5);反之,進行步驟(7)。t1、t2分別為:

+Nsteval

(10)

t2=Nmaxteval

(11)

D0、rmax,0、d0、Ns,0、d0-1、pmax,0是基準問題參數(shù),tci,0為對應(yīng)組參數(shù)的tci,這些參數(shù)由已知系統(tǒng)先行計算得出,以便在計算t1、t2時更加快捷。

(5) 設(shè)定抽樣次數(shù)Ns為drmax(pmax+1)2,使用低階多項式逼近法求解隨機潮流,計算低階多項式逼近法逼近誤差ε。

(6) 若ε小于3TOL,則進行步驟(8);反之,進行步驟(7)。

(7) 使用Nmax點拉丁超立方法求解隨機潮流。

(8) 輸出電壓幅值、相角的期望和標準差。

基于綜合算法的隨機潮流計算流程如圖1所示。

6 算例

通過調(diào)整設(shè)置不同規(guī)模大小的系統(tǒng)和隨機輸入變量數(shù)量,檢驗綜合算法的逼近精度和計算時間,以及驗證方法的有效性和優(yōu)越性。所采用系統(tǒng)參數(shù)全部來自Matpower 6.0。

6.1 算例1 IEEE-30節(jié)點系統(tǒng)

在IEEE-30節(jié)點系統(tǒng)基礎(chǔ)上,添加五個服從韋布爾分布的風(fēng)電功率注入和四個服從貝塔分布的光伏注入。其中,五個風(fēng)電場分別配置在1、5、6、9、10節(jié)點處,參數(shù)λ為8,k為3,切入風(fēng)速為2 m/s,額定風(fēng)速為12 m/s,切出風(fēng)速為25 m/s,功率因數(shù)為0.2,額定功率為100 MW。四個光伏電場分別配置在12、13、14、15節(jié)點處,參數(shù)α、β為2.222 2,光照強度峰值為1 000 W/m2,光伏電場輻射點輻射為150 W/m2,標準太陽輻射為1 000 W/m2,額定功率為50 MW。自適應(yīng)稀疏偽譜法設(shè)置最大抽樣次數(shù)為10 000次,誤差收斂閾值為0.005。以10 000點拉丁超立方法計算結(jié)果為基準。

這一隨機潮流問題的輸入隨機變量維數(shù)為9,可以認為是低維問題。采用基于綜合算法求解這一問題時,使用了自適應(yīng)稀疏偽譜法。其中,總積分點數(shù)為229,計算時間為0.175 s,總指標數(shù)為60,誤差估計為0.004 189。由于誤差估計小于給定誤差收斂閾值0.05,因此采用自適應(yīng)稀疏偽譜法的計算結(jié)果作為最終輸出結(jié)果。計算結(jié)果輸出期望對比如圖2所示,輸出標準差對比如圖3所示。基準拉丁超立方法計算時間為3.596 s。

圖2 算例1輸出期望對比

圖3 算例1輸出標準差對比

由圖2和圖3可以看出,自適應(yīng)稀疏偽譜法對系統(tǒng)各輸出變量的計算結(jié)果幾乎與拉丁超立方法計算結(jié)果重合,說明自適應(yīng)稀疏偽譜法對于低維場景具有很高的精度。

由于低階多項式逼近法和拉丁超立方法的采樣點數(shù)量與自適應(yīng)稀疏偽譜法每次計算時按需所得不同,而是先驗給定的,因此為了體現(xiàn)筆者所提方法在求解問題時的有效性,做出采樣次數(shù)為200、230、260時的低階多項式逼近法和拉丁超立方法實際逼近誤差曲線,實際逼近誤差用各輸出變量的標準差相對于拉丁超立方法計算出的基準值的平均相對誤差來衡量,比較結(jié)果如圖4。

圖4 算例1采樣次數(shù)對比

拉丁超立方法的運行時間分別為0.121 s、0.133 s、0.135 s,低階多項式逼近法的運行時間分別為0.319 s、0.322 s、0.354 s。無論是低階多項式逼近法、自適應(yīng)稀疏偽譜法,還是拉丁超立方法,計算時間都極短,采用計算精度最高的方法對隨機潮流進行求解,無疑是最優(yōu)的選擇,因此從算例1可以看出,綜合算法充分利用了自適應(yīng)稀疏偽譜法在低維時的高計算精度,展現(xiàn)了方法的有效性。

6.2 算例2 IEEE-85節(jié)點系統(tǒng)

在IEEE-85節(jié)點系統(tǒng)的基礎(chǔ)上,添加三個服從韋布爾分布的風(fēng)電注入、八個服從貝塔分布的光伏注入、59個服從高斯分布的負荷噪聲,風(fēng)電場位于1、2、3節(jié)點處,光伏電場位于4、5、6、7、10、11、14、15節(jié)點處。有84個負荷噪聲在其它節(jié)點處,均值為原負荷量,標準差為5%均值。風(fēng)光電場參數(shù)及自適應(yīng)稀疏偽譜法參數(shù)同算例1。

此時,輸入維度為70,屬于中高維問題。綜合算法在求解這一問題時,先采用自適應(yīng)稀疏偽譜法進行求解,然后采用10 000點拉丁超立方法進行求解。自適應(yīng)稀疏偽譜法的求解結(jié)果為采樣次數(shù)10 025、指標數(shù)1 906、全局誤差指示器0.125 01,計算時間19.675 s。拉丁超立方法與自適應(yīng)稀疏偽譜法計算結(jié)果輸出標準差對比如圖5所示。采用拉丁超立方法進行這一系統(tǒng)隨機潮流的求解,精度遠遠高于自適應(yīng)稀疏偽譜法。因此,筆者所提方法可以在很大程度上保證對各種規(guī)模系統(tǒng)隨機潮流求解的精度和計算時間。

圖5 算例2拉丁超立方法與自適應(yīng)稀疏偽譜法輸出標準差對比

綜合算法在求解這一問題時,首先采用自適應(yīng)稀疏偽譜法進行求解,在求解精度不能滿足要求時,進而判斷低階多項式逼近法與拉丁超立方法的計算時間,選擇時間較短的進行計算。根據(jù)自適應(yīng)稀疏偽譜法的抽樣時間和抽樣點數(shù),可以計算出單次采樣時間為1.097 9×10-3s。再根據(jù)已知系統(tǒng)計算tbasis為4.916 5×10-8s。設(shè)置低階多項式逼近法采樣點數(shù)為1 000。計算出低階多項式逼近法預(yù)計所需時間為75 s,而拉丁超立方法預(yù)計所需時間為10.9 s。可見,綜合算法最終采用拉丁超立方法進行計算。

綜合算法與低階多項式逼近法計算結(jié)果輸出標準差對比如圖6所示,可以看出此時低階多項式逼近法與綜合算法精度近似,采取綜合算法計算時間更短,充分體現(xiàn)綜合算法的合理性和有效性。

圖6 算例2綜合算法與低階多項式逼近法輸出標準差對比

7 結(jié)束語

基于綜合算法的隨機潮流計算方法在隨機潮流問題求解的過程中,充分結(jié)合了三種各具優(yōu)勢的隨機潮流計算方法,并且通過計算精度、計算時間等條件,將三種方法銜接,發(fā)揮三種方法的長處,極大避免自適應(yīng)稀疏偽譜法在高維計算量時的維數(shù)災(zāi)難問題。通過算例,驗證了所提方法的有效性。