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

孟代江 張淑榮

同濟大學 電子與信息工程學院 上海 201804

1 研究背景

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

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

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

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

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

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

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

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

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

2 隨機潮流數學模型

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

U=f(W)

(1)

U=(U1,U2,…,UD)

W=(W1,W2,…,Wd)

式中:U為節點電壓、相角等系統狀態隨機變量;D為輸出隨機變量數量;W為風光電場出力及負荷等輸入隨機變量;d為輸入隨機變量數量;f表示潮流方程隱含的關于系統狀態與節點注入功率的隱函數關系。

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

U=f(v,r,PL)

(2)

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

3 自適應稀疏偽譜法計算時間復雜度

對于隨機潮流模型,自適應稀疏偽譜法構造的單輸出多項式逼近函數為:

(3)

(4)

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

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

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

(5)

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

表1 總積分點數與輸入維數關系

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

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

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

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

低階多項式逼近法構建的單輸出逼近函數表達式為:

(6)

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

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

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

(7)

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

t1s=Drmaxdtsolve

(8)

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

tErr=DrmaxNd(pmax+1)tbasis

(9)

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

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

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

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

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

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

(3) 采用自適應稀疏偽譜法求解隨機潮流。若全局誤差η滿足小于TOL,且抽樣次數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是基準問題參數,tci,0為對應組參數的tci,這些參數由已知系統先行計算得出,以便在計算t1、t2時更加快捷。

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

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

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

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

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

6 算例

通過調整設置不同規模大小的系統和隨機輸入變量數量,檢驗綜合算法的逼近精度和計算時間,以及驗證方法的有效性和優越性。所采用系統參數全部來自Matpower 6.0。

6.1 算例1 IEEE-30節點系統

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

這一隨機潮流問題的輸入隨機變量維數為9,可以認為是低維問題。采用基于綜合算法求解這一問題時,使用了自適應稀疏偽譜法。其中,總積分點數為229,計算時間為0.175 s,總指標數為60,誤差估計為0.004 189。由于誤差估計小于給定誤差收斂閾值0.05,因此采用自適應稀疏偽譜法的計算結果作為最終輸出結果。計算結果輸出期望對比如圖2所示,輸出標準差對比如圖3所示?；鶞世〕⒎椒ㄓ嬎銜r間為3.596 s。

圖2 算例1輸出期望對比

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

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

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

圖4 算例1采樣次數對比

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

6.2 算例2 IEEE-85節點系統

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

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

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

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

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

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

7 結束語

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