計及多個風場預測誤差的電力系統風險快速計算方法

張 沛 田佳鑫 謝 樺

（北京交通大學電氣工程學院 北京 100044）

0 引言

受能源轉型的需求及國家政策影響，新能源并網規模越來越大，但大規模風電的波動性會給電網的調度運行帶來風險，所以需要對電網運行開展量化風險評估。因此，風險的計算方法的研究顯得尤為重要。

電力系統風險計算方法主要有四種：點估計法、狀態枚舉法、蒙特卡洛法和半不變量法。在點估計法和狀態枚舉法的研究方面，文獻[1-3]均采用點估計法來計算系統的運行風險，構建系統風險指標與可用電源輸出隨機變量的函數關系，從隨機變量中選取少量估計點進行函數運算，進而求得風險指標。點估計法雖然在計算風險指標時有較快的速度，但是當系統規模較大時，輸出隨機變量的高階距和概率分布誤差增加，其實用性與準確性會降低，還需進一步考證。文獻[4]中所列舉的狀態枚舉法，通過對系統出現的所有狀態逐一枚舉并進行潮流計算，進而得出風險，使用該方法進行風險計算時，結果固然準確，但當系統中接入多個風場時，基于狀態枚舉的方法計算運行風險時會出現系統狀態組合爆炸的問題。例如，某電力系統中風場數量為n，每一風場的出力區間為Ni，則系統的狀態數為，采用枚舉法計算運行風險時不僅需要進行次潮流計算，還需要在每一次潮流計算結果的基礎上計算運行風險，耗費大量時間，無法滿足調度運行的時效性。

在蒙特卡洛法的研究方面，文獻[5]采用改進的分散抽樣蒙特卡洛法，其中等分散抽樣僅僅在處理多狀態元件上具有較高的抽樣效率，可以加快計算的收斂速度；文獻[6]采用傳統的非序貫蒙特卡洛法來計算風險，其中風場預測誤差運用TLS（T-location scale）擬合，具有好的擬合效果，文獻中提出切負荷、電壓越限、線路有功越限及電壓崩潰等風險指標，全方面地評估系統運行風險；文獻[7]采用蒙特卡洛模擬法評估光伏并網系統的時變風險，通過失負荷概率和期望缺供電量風險指標來評估風險水平；文獻[8-10]均采用非序貫蒙特卡洛法建立電力系統的運行風險評估模型，較序貫蒙特卡洛法而言，該算法具有需求內存少、模型簡單及收斂快的優點。文獻[5-10]都采用蒙特卡洛法或者蒙特卡洛的改進算法，要進行大規模抽樣才能保證其準確度，但會導致計算效率下降。這些方法雖然提高了狀態抽樣的效率，可以加快收斂速度，但是無法避免大量重復的潮流計算，計算周期長，內存消耗大，不利于電力系統的實時調度運行。

半不變量法是根據潮流計算結果構建輸出隨機變量和輸入隨機變量的線性化關系，在此基礎之上，利用半不變量的性質計算出輸出隨機變量的半不變量，結合級數展開理論得到輸出隨機變量的概率分布函數。文獻[11]采用半不變量法計算線路隨機故障下的運行風險，同時改進了分布式電源的隨機模型，但僅應用于小型配電系統，針對大型主電網的風險計算并沒體現；文獻[12]在快速解耦的潮流模型基礎上采用半不變量的改進Von Mises 方法計算運行風險，加快了計算速度，但該方法是建立在快速解耦法之上，因而前提為電抗遠大于電阻，這一條件不是所有電力系統都能滿足，因此存在一定的局限性；文獻[13]先采用基于序列運算理論實現隨機潮流計算，之后將隨機潮流引入安全評估，通過概率與風險兩個角度來實現風險的量化。文獻[11-13]均采用半不變量法或者其改進算法，雖然避免了大量重復的潮流計算，提高了計算效率，但輸入隨機變量波動范圍較大時，計算結果會產生較大線性化誤差。文獻[14]對這種誤差進行了研究和說明。同時由于這種線性化誤差的存在會導致輸出隨機變量的概率分布發生局部畸變，使得風險評估結果準確度不能夠達到實時調度運行的要求。

針對上述不足，本文提出快速計算運行風險的方法，可以快速計算含有多個風場預測誤差的電力系統的運行風險。本文所提方法在半不變量法基礎之上加入了聚類技術。在該方法中，考慮風電預測誤差的節點注入功率（簡稱節點注入功率）作為輸入隨機變量，節點電壓與支路潮流作為輸出隨機變量。基于K-means 聚類技術對輸入隨機變量樣本進行聚類，以降低樣本的波動范圍。針對每一簇類進行潮流計算，根據潮流計算結果構建輸出隨機變量與輸入隨機變量的線性化關系，并結合半不變量的性質計算輸出隨機變量的K階半不變量，進而采用Gram-Charlier 級數展開理論計算輸出隨機變量累積分布函數（Cumulative Distribution Function, CDF）。根據CDF 曲線，對越限部分進行積分運算求得運行風險，并計算其運行風險指標。經過聚類處理后，輸入隨機變量樣本的波動范圍大大減小，在每一類簇中采用半不變量法相比整個樣本采用半不變量法得到的線性化關系更加準確，提高了運行風險評估結果的準確性。

1 考慮風場預測誤差波動的潮流線性化關系構建

本文考慮多個風場于不同節點同時并網的運行工況，計及預測誤差的風場出力作為系統的輸入隨機變量。在該運行工況，由于風電預測誤差的波動范圍較大，所以每個風機并網點的風場出力也具有較大的波動性。由于這種波動性的存在，以風場出力期望值代入系統得到的潮流線性化關系進行風險計算會產生較大的線性化誤差，因此本文提出采用K-means 對輸入隨機變量樣本進行聚類，建立每一簇類樣本對應的潮流線性化關系，可以降低每一簇類內風場出力的波動性，減小線性化誤差。

1.1 輸入隨機變量的概率分布

輸入隨機變量是風場出力，其中風場出力是風場預測功率與預測誤差之和，故輸入隨機變量的概率分布與風場預測誤差的概率分布相同。其中，由于風場之間存在一定的差異，所以不同風場預測誤差的概率分布也不同。本文風場預測誤差的概率分布是根據某地區風場出力實際功率與預測功率的差值進行統計計算而來，因此更符合實際運行工況。

由于輸入隨機變量之間相互獨立是使用半不變量法的前提條件，因此輸入隨機變量之間需要進行相關性處理。根據文獻[15-16]中的Cholesky 分解去除多個輸入隨機變量之間的相關性。根據輸入隨機變量的概率分布，并進行相關性處理后，生成符合該分布的樣本，可表示為

式中，Xwind為輸入隨機變量樣本；xij為第j個隨機變量中第i個樣本的數值；m為輸入隨機變量的數量；N為每個輸入隨機變量中樣本數量。

1.2 基于K-means 技術的樣本聚類

為了減小由風場出力波動性較大而造成的線性化誤差，采用K-means 技術[17-18]將輸入隨機變量聚為多簇類。其聚類過程如下：

（1）選取初始聚類中心。為了保證能夠得到較好的聚類效果，針對輸入隨機變量樣本Xwind而言，隨機選取10%的樣本進行聚類，得到a個聚類中心，即

式中，Mwind為聚類中心矩陣；Mij為第j個輸入隨機變量第i簇類的聚類中心。通過該方式得到的聚類中心可以反映整個樣本的分布情況，且可以提高聚類效率。因此，將這10%樣本的聚類中心作為初始聚類中心。

（2）計算所有點到聚類中心的歐式距離。歐式距離D為

（3）計算簇中樣本均值并更新聚類中心。根據歐式距離大小，更新每一簇類中的點，并計算簇類中樣本均值，更新聚類中心Mwind。

（4）重復步驟（2）、步驟（3）直到簇類中的點和聚類中心不發生改變，將最后一次得到的聚類中心作為輸入隨機變量的聚類中心。

1.3 潮流線性化關系的建立

針對多個風場并網的系統，將輸入隨機變量經過K-means 聚類后形成的聚類中心Mwind代入系統進行a次潮流計算。根據每次潮流計算結果，得到每一簇類樣本下的潮流線性化關系。以某一簇類為例，其計算方法如下。

電力系統的潮流方程[15,19]為

式中，f(·)為潮流方程；S為節點注入功率向量；V為節點電壓向量。在該簇類的中心處，應用泰勒級數對式（4）展開，可得該運行點處的潮流線性化關系為

式中，Z為支路功率向量；V0和Z0分別為當前運行點處的節點電壓和支路功率向量；ΔV和ΔZ分別為節點電壓和支路功率的變化向量；和T0分別為節點電壓和支路功率對節點注入功率的靈敏度矩陣；ΔSwind為節點注入功率的變化向量。

由于輸入隨機變量樣本具有較大的波動性，采用整個樣本的期望值計算潮流線性化關系誤差較大，不能完全滿足整個樣本的波動特性。因此基于Kmeans 技術將輸入隨機變量樣本聚為a類，在每一類中，樣本的波動性會大大降低，其構建的潮流線性化關系更加符合該類樣本的波動特性，使得線性化誤差大大減小，解決了文獻[11-13]中的線性化誤差問題。

2 輸出隨機變量概率分布的計算方法

基于K-means 技術對輸入隨機變量樣本聚類的結果，針對每一簇類而言，通過式（5）可以得到輸入隨機變量與輸出隨機變量之間的線性化關系，根據線性化關系可以計算出輸出隨機變量的概率分布。在這一過程中，對聚類后的每一類簇中應用半不變量法，將復雜的卷積運算簡化為加法運算，依次計算輸入隨機變量和輸出隨機變量的半不變量，并采用級數展開的方式計算輸出隨機變量的概率分布，不僅降低了計算的復雜程度，而且K-means 聚類方法的使用解決了使用傳統半不變量法時由于輸入隨機變量波動大而導致線性化誤差的問題，使輸出隨機變量概率分布更加準確，提高了風險計算的精度。

2.1 輸入隨機變量半不變量的計算方法

根據聚類結果，采用式（6）[20]可以計算每一簇類樣本對應輸入隨機變量的k階半不變量

式中，m0和β分別為經過聚類后簇中樣本的均值和原點矩。β的計算公式為

式中，xj為簇類中的樣本；p(xj)為xj的概率。

2.2 輸出隨機變量半不變量的計算方法

以每一簇類樣本作為研究對象，式（5）中的ΔSwind需要通過各風場預測誤差的概率分布卷積計算得到，運算復雜。利用半不變量的性質可以將卷積運算簡化為加法運算，故每一簇類樣本對應的輸出隨機變量半不變量計算公式為

式中，ΔV(k)為當前簇類樣本對應的節點電壓k階半不變量；ΔZ(k)為當前簇類樣本對應的支路潮流k階半不變量。

將每一簇類樣本對應的輸入隨機變量半不變量代入式（8）中，可以計算該簇類樣本對應的輸出隨機變量半不變量。

針對全樣本而言，計算輸出隨機變量的半不變量方法如下：

（1）計算全概率值。計算各簇類樣本數量占總樣本數量的概率。

（2）計算各簇類樣本的原點矩。根據式（6）計算各簇類樣本對應的各階原點矩。

（3）計算全樣本的原點矩。將各簇類樣本的原點矩與全概率值對應相乘并累加得到全樣本的原點矩。

（4）計算全樣本下的輸出隨機變量半不變量。將全樣本的原點矩代入式（6）中，計算全樣本下的輸出隨機變量半不變量。

2.3 輸出隨機變量累積分布的計算方法

在輸出隨機變量半不變量已知的狀況下，本文采用級數展開法計算輸出隨機變量的累積分布函數。目前，常用的級數展開方法有Gram-Charlier 級數展開、Edgeworth 級數展開和Cornish-Fisher 級數展開。其中，Cornish-Fisher 級數展開對非正態分布變量更加準確，且展開的多項式中是以隨機變量的分位數作為系數，而本文采用的是半不變量，故不采用Cornish-Fisher 級數展開法；Edgeworth 級數展開在隨機變量的CDF 和PDF 尾部存在較大誤差；Gram-Charlier 級數展開卻具有較高的精度[21]，且節點電壓越上限和支路潮流越限一定出現在輸出隨機變量CDF 曲線的尾部，若該部位出現較大誤差，其運行風險的計算結果也會有較大誤差。由此可見，Gram-Charlier 級數展開法比Edgeworth 級數展開法更加合適。因此本文采用Gram-Charlier 級數展開法[22]計算輸出隨機變量的累積分布函數。

Gram-Charlier 級數展開求取輸出隨機變量的概率密度函數f(x)，即

式中，?(·)為標準正態分布的概率密度函數；μ為輸出隨機變量的期望值；σ為輸出隨機變量標準差；gk為輸出隨機變量標準化后的k階半不變量；Hk(·)為Hermite 多項式[23-24]。

對式（9）進行積分運算可得出輸出隨機變量的累積分布函數F(x)，即

式中，1x和2x分別為積分上、下限。

將節點電壓和支路潮流的各階半不變量代入式（10），可以計算相應的概率分布函數。

3 運行風險指標計算和評估流程

3.1 運行風險評估指標的計算方法

系統運行風險是通過風險發生概率與產生影響的乘積來綜合表征。其中，產生的影響在電力系統中體現為兩方面：節點電壓越限和支路潮流越限。因此，本文采用節點電壓越限指標和支路潮流越限指標來定量評估系統的運行風險。在計算運行風險時，本文提出采用積分運算理論，對CDF 曲線中越限部分進行積分運算求得風險，進而計算運行風險指標。

3.1.1 節點電壓越限風險指標

節點電壓越限風險指標是對系統在一定條件下發生越上、下限的度量。根據電力系統運行導則可知，110kV 及其以下的節點電壓的限值范圍為0.9(pu)～1.1(pu)，220kV 及其以上的節點電壓的限值范圍為0.95(pu)～1.05(pu)。以220kV 以上的電壓等級為例，節點電壓的越限風險指標計算如下。

節點電壓的 CDF 曲線如圖 1 所示。若Vimin≥0.95(pu)且Vimax≤1.05(pu)時，節點電壓越限風險為0；若Vimin或Vimax不在0.95(pu)～1.05(pu)之間，節點電壓越限風險為圖1 中陰影面積之和，計算方法為

式中，Vimin為第i節點電壓可能出現的最小值（即圖1 中累計概率剛好大于零所對應的電壓值）；Vimax為第i節點電壓可能出現的最大值（即圖1 中累積概率剛好等于1 時所對應的電壓值）；Ri為節點i的電壓越限風險。針對具有n個節點的電力系統，系統的節點電壓越限風險指標

圖1 節點電壓的CDF 曲線Fig.1 CDF curve of node voltage

3.1.2 支路潮流越限風險指標

與節點電壓越限指標相似，但支路潮流只有越上限的風險。支路潮流上限值為線路允許的最大傳輸功率Pmax，其大小與線路型號有關。支路潮流的越限風險指標計算如下。

若PLmax≤Pmax，支路潮流的越限風險為0；若PLmax＞Pmax，支路潮流的越限風險為

式中，PLmax為某一支路潮流可能出現的最大值；Rij為支路ij的潮流越限的風險。對具有h條支路的電力系統，其支路潮流越限風險指標

3.2 運行風險評估流程

采用本文所提方法進行運行風險評估的具體計算流程如圖2 所示。

圖2 計算流程Fig.2 Compute process

（1）采用K-means 對計及預測誤差的風電出力樣本聚類。將輸入隨機變量樣本進行聚類，得到a個聚類簇和a個聚類中心。

（2）建立每一簇類的潮流線性化關系。將每一簇類的聚類中心分別代入系統進行潮流計算，確定在每一簇類樣本下的潮流線性化關系。

（3）計算輸出隨機變量的半不變量。首先計算每一簇類中輸入隨機變量的k階半不變量，再根據潮流線性化關系計算在每一簇類樣本下的輸出隨機變量的k階半不變量，最后根據全概率公式計算整個樣本下的輸出隨機變量的k階半不變量。

（4）計算累積分布函數。利用Gram-Charlier 級數展開理論分別計算節點電壓和支路潮流的累積分布函數。

（5）計算運行風險指標。根據累積分布函數計算各節點和支路的越限風險，分別累加各節點和支路的風險值，得到運行風險指標。

4 案例分析

本文選取某實際電網作為算例，基于該算例設定了兩種不同場景進行仿真研究。在這兩個場景中，場景1 風電并網點為兩個，風電裝機總量為400MW；場景2 風電并網點為四個，風電裝機總量4 000MW。針對上述兩個場景而言，本文將場景1 的結果詳細列出并進行分析；由于場景2 所采用方法和仿真過程與場景1 一致，且出于篇幅的原因，詳細結果和過程本文不再贅述，其得到的風險指標及計算耗時將和場景1 進行對比分析。

4.1 算例系統說明

在該算例系統中，共有810 個節點和977 條支路，其中有263 個負荷節點及48 臺發電機，總共負荷為17 400.2MW。系統電壓等級有500kV、220kV、110kV、66kV 和35kV 及其以下，圖3 的單線圖展示了500kV 和220kV 主網架。

圖3 電網單線Fig.3 Single-line diagram of power grid

場景1 中，在220kV 等級的母線55 和母線190 處分別接入裝機容量為200MW 的風電場；場景2 中，分別在母線55 和母線190 接入裝機容量為200MW 的風電場，將母線7 和母線16 處的火電機組用等容量的風電場代替，故總的裝機容量為4 000MW。風場預測誤差的概率分布是根據該地區風場實際數據計算而來，選取該地區風電場兩年的歷史數據（實際功率和預測功率）作為數據支撐，基于上述數據進行統計，計算出適合該地區特性的風場預測誤差概率分布。風場預測誤差的概率分布通過對所選數據進行數值擬合得到。場景1 中，采用風場實際運行數據進行擬合，得到風場預測誤差概率分布為正態分布（擬合結果和比較證明見附錄），這也滿足了Gram-Charlier 級數展開法對輸入數據分布特征的要求。基于此，風場預測誤差的概率分布如圖4 所示。圖4 中分別為節點55 和190 處風電場預測誤差的概率分布，風場出力概率分布與其相同。本文將得到的正態分布參與后續的風險評估計算，這樣更能反映該地區的實時運行工況。

圖4 風場預測誤差的概率分布Fig.4 Probability distribution of wind farms forecasting errors

4.2 輸入隨機變量樣本的聚類結果

根據圖4 的概率分布，采用蒙特卡洛法生成20 000 組符合該分布的樣本，使用K-means 對場景1 中兩個風場的20 000 組樣本進行聚類。對輸入隨機變量進行聚類時，聚類數量越多效果越好，但耗時較長。因此聚類數量的合理選取不僅可以保證聚類效果，還能提高計算效率。本文使用加權平均半徑[21]（Weighted Average Radius, WAR）來確定聚類數目，其表達式為

式中，WAR 為聚類后各簇類的加權平均半徑；Pi為第i簇類中樣本數量占總樣本的比例；ri為第i簇類的半徑。通過上述計算得出場景1 中WAR 與聚類數量的關系如圖5 所示。

圖5 加權平均半徑和聚類數量關系Fig.5 Relationship between WAR and number of clusters

根據圖5 可知，隨著聚類數量的增加，加權平均半徑越小，說明聚類的效果越好。從圖中可以看出，聚類數量為10 時，WAR 有一個明顯的轉折；當聚類數量超過20 時，WAR 趨于穩定，且變化較小，因此最佳聚類數量選取應在10～20 之間。基于此，在該區間選取10、12、15、17、19 對比聚類效果，經過對比發現，聚類數量為10 和12 時，對于邊緣點聚類效果較差；聚類數量為15、17 和19 時，聚類效果相差較小，考慮到聚類數量和計算耗時成反比的因素，同時在圖5 中發現在聚類數量為15 時有一個轉折，綜合上述幾種因素，所以選擇聚類數量為15。對輸入隨機變量樣本聚為15 類，其聚類效果如圖6 所示。

圖6 聚類結果Fig.6 Clustering results

在圖6 的聚類結果中，15 個聚類中心坐標Mij見表1。

表1 聚類中心坐標Tab.1 Coordinates of the cluster center

4.3 輸出隨機變量半不變量的計算

根據聚類結果和2.2 節可以計算15 簇類中各類樣本對應輸入隨機變量半不變量、潮流線性化關系和輸出隨機變量半不變量。由于聚類數目和輸出隨機變量數量較多，以第15 個簇類樣本為研究對象，母線55 處風場注入功率的半不變量見表2。

表2 母線55 注入功率的半不變量Tab.2 Cumulant of injected power at bus 55

根據表1 將第15 簇類的聚類中心代入算例系統進行潮流計算，得到該簇類樣本對應的潮流線性化關系。以節點190 的電壓和支路805 的潮流為例，其線性化關系為

結合表2 和式（14）可以計算該簇類樣本對應的節點190 的電壓和支路805 的潮流半不變量，見表3。

針對其他簇類而言，將聚類中心Mij代入系統中進行潮流計算，其潮流線性化關系和輸出隨機變量半不變量均可計算。計算每個簇類樣本數量占全樣本數量概率見表4。

表3 輸出隨機變量的半不變量Tab.3 Cumulant of output random variable

表4 全概率值Tab.4 Total probability value

根據表4 的全概率值計算整個樣本下的輸出隨機變量的半不變量。還以節點190 和支路805 為例，其8 階半不變量見表5。

表5 輸出隨機變量的半不變量Tab.5 Cumulant of output random variable

4.4 輸出隨機變量累計分布的計算

根據2.3 節中的Gram-Charlier 級數展開理論求得輸出隨機變量的累積分布，以節點475 和支路544為例，CDF 曲線如圖7 所示。

圖7 輸出隨機變量的CDF 曲線Fig.7 CDF curve of output random variable

圖7a 中陰影部分面積為節點475 的越限風險值，圖7b 中陰影部分面積為支路544 的越限風險值。

4.5 系統運行風險指標

在求得輸出隨機變量的累積分布的基礎上，根據式（11）可以計算各節點電壓的越限風險值；同理，各支路潮流的越限風險值也可以由式（12）計算得到。經過計算，場景1 中有121 個節點和9 條支路存在越限風險；經過同樣的仿真計算，場景2 中有154 個節點和8 條支路存在越限風險。兩個場景運行風險指標見表6。

表6 運行風險指標Tab.6 Operational risk indices

從表6 中可見，針對場景1 而言，由于風電波動而引起的節點電壓總的風險為0.002 1，支路潮流總的風險為0.012 7；針對場景2 而言，節點電壓的總風險為0.018 7，支路潮流的總風險為0.032 6。由此可見，場景2 的節點電壓越限風險和支路潮流越限風險均比場景1 大。

4.6 準確性和快速性驗證

為了驗證本文所提計算運行風險方法的準確性和快速性，與蒙特卡洛法、未聚類的半不變量法和狀態枚舉法分別進行對比分析。

本文為了保證計算結果的準確性，采用文獻[7]中的蒙特卡洛模擬法作為基準與之比較。蒙特卡洛模擬法是根據樣本滿足相應分布產生大量隨機數，在此基礎上進行大量重復性潮流計算。由于樣本數量足夠大，故計算結果足夠精準，且不受系統規模的影響，是衡量其他方法準確性的重要參考。因此，本文基于輸入隨機變量的概率分布，采用蒙特卡洛模擬法生成20 000 組隨機數，進行了20 000 次仿真結果作為對比的標準。以節點400 為例，本文方法所得CDF 和蒙特卡洛法的對比見圖8。由圖8 可見，本文方法的計算結果與蒙特卡洛法吻合度較高，但仍存在一定誤差。為了避免偶然性，同時采用了文獻[20]中平均方均根（Average Root Mean Square,ARMS）描述兩種方法計算結果的誤差，經過大量計算（計算結果見附錄），結果表明本文提出的方法與蒙特卡洛法相比，誤差較小，可以作為隨機概率潮流結果。

圖8 CDF 曲線對比Fig.8 The comparison of CDF curves

同時，本文還采用狀態枚舉法對系統進行了風險評估。由圖4b 可知，場景1 中每個風場有15個出力區間（風場出力區間的選擇過程見附錄），因此風場出力共有152種組合狀態；同理場景2 中風場出力共有 154種組合狀態。對兩個場景的所有狀態進行潮流計算，之后分別計算各個場景的風險指標。

為避免人為主觀層面的影響，選取K-means 不同聚類數對計算效率的影響進行驗證。分別選取10、12、15、17 和20 個聚類數并使用本文方法進行仿真計算，其耗時結果見表7。

表7 耗時對比Tab.7 Time comparison

由表7 可見，隨著聚類數量增加，其耗時在逐漸增加，計算效率會降低（主要原因是聚類數量的增加導致潮流計算的次數增加）。結合之前聚類數量選取的分析，聚類數量選取15 合理。

在基于相同的服務器環境和仿真平臺下，對蒙特卡洛法、未聚類的半不變量法、狀態枚舉法及本文提出的方法分別進行了編譯計算，針對兩個場景采用上述四種方法的計算用時和風險指標見表8。

表8 綜合對比Tab.8 Comprehensive comparison

表8 中，以蒙特卡洛法的計算結果為標準，在計算的準確度方面，對于場景1，本文方法的準確度為97.4%，半不變量法（未聚類）的準確度為83.6%，狀態枚舉法的準確度為82.9%；對于場景2，本文方法的準確度為98.6%，半不變量法（未聚類）的準確度為89.8%，狀態枚舉法的準確度為88.5%。通過兩個場景均可以看出，本文提出的方法較半不變量法（未聚類）準確度更高，這是由于采用了聚類技術使得線性化誤差大大減小，從而提高了計算精度。綜上，本文方法的準確度最高，計算結果與蒙特卡洛法最為接近，證明了本文方法的準確性。

在計算耗時方面，針對場景1 而言，本文方法用時12.3s，僅為蒙特卡洛法的2.68%，其耗時較少的原因有三方面：①本文方法對輸入變量樣本進行聚類，其聚類數遠小于樣本數量，因此大大減少了潮流計算次數；②采用本文方法在計算輸出隨機變量的概率分布時，運用半不變量的方法將卷積運算轉換為加法運算，大大降低了計算復雜程度；③根據累積分布計算風險時，僅對越限部分進行統計計算，減少了計算量。半不變量法（未聚類）與本文方法相比，由于只進行一次潮流計算，所以耗時更少，但使用該方法在風電出力波動大的點存在較大的線性化誤差，準確度低。狀態枚舉法耗時11.3s，為蒙特卡洛法的2.47%，與本文方法相差較小。通過對比相同方法下兩個場景的耗時情況，發現本文方法、半不變量法（未聚類）和蒙特卡洛法在兩種場景下耗時相差很小，但狀態枚舉法相差非常大，這是由于風電并網點的增加，使系統組合狀態數量從152暴增到154，潮流計算量呈指數增長，突顯出狀態枚舉法在處理該種問題上存在一定的弊端。綜上所述，通過各方法之間的橫向比較和場景之間的縱向比較，可見本文所提方法在計算運行風險時突顯其快速性。

5 結論

針對含多個風場的電力系統，本文提出考慮風電預測誤差的系統運行風險快速計算方法。其中，K-means 聚類技術的運用有效降低了由風電出力波動性大而導致的線性化誤差，在此基礎之上，在每一類簇中通過將半不變量法和Gram-Charlier 級數結合，可以準確快速地計算節點電壓和支路潮流的概率分布，最后通過積分運算計算運行風險。通過算例分析對本文方法進行了驗證，場景1 的結果表明：本文方法計算運行風險比蒙特卡洛法耗時更少，計算效率提高了97.3%；以蒙特卡洛法為標準，與半不變量法（未聚類）和狀態枚舉法相比，本文方法計算的運行風險準確度最高，達97.4%。

綜上所述，采用本文方法進行風險評估，不僅具有高的計算精度，還可大大降低風險評估的耗時，能夠滿足電力系統實時調度運行在精度和時效上的需求。

附 錄

1. 場景1 中，對風場預測誤差數據進行擬合，其擬合結果如附圖1 所示。

為了驗證擬合的精確性，分別計算了附圖1a 和附圖1b 的ARMS。附圖1a 中的ARMS 為7.2×10-4，附圖1b中的ARMS 為6.3×10-4，其值均非常小，因此風場預測誤差概率分布滿足正態分布。

附圖1 風場預測誤差擬合結果App.Fig.1 Fitting result of wind farms forecasting error

2. 由于本文所選算例規模較大，節點和支路數較多，若全部列出，篇幅較長，因此附表1 中只列出部分數據。

附表1 部分CDF 曲線的ARMS 值App.Tab.1 ARMS value of partial CDF curve

3. 采用狀態枚舉法計算系統的運行風險時，選取合適的風電出力區間數量，可以節省大量計算時間。分別選取6、8、10、12、15、20、30、40、50、60、70、80個區間進行計算，將各種情況下的耗時和風險值記錄于附表2。

附表2 不同區間的對比App.Tab.2 Comparison of different discrete intervals

基于上述計算結果可見，隨著區間數量的增加，其支路潮流越限風險也趨于穩定，但計算耗時在逐漸增加。若選擇區間數量較少，計算速度快，但計算結果精度較低；若選擇區間數量較多，計算精度可以滿足要求，但耗時較長。支路潮流越限風險的均值為0.018 1，觀察附表2 中0.017 8 與之最為接近，因此選擇區間數量為15。