計及隨機變量相關性的多點線性化概率潮流計算

葉 希，王彥灃，黃 楊，廖幫昆，歐陽雪彤，文云峰

（1.國網四川省電力公司，成都 610041；2.輸變電新技術教育部工程研究中心（湖南大學），長沙 410082）

在“碳達峰/碳中和”發展目標和新型電力系統構建的背景下，風電、光伏等新能源發電裝機容量快速增長，源?網?荷各環節電力電子化程度顯著加深。以同步電源為主的電力系統正逐漸演進為“高比例新能源”并網和“高比例電力電子裝備”饋入的“雙高”電力系統[1]。在“雙高”電力系統中，源荷雙側的功率不確定性顯著增強，電力電量平衡以及容量充裕度的概念與方法將由目前確定性的思路向概率性思路轉變[2]。

概率潮流計算能夠計及不確定性因素并分析其對系統運行特性的影響，可為“雙高”電力系統安全穩定分析提供更為可靠的理論支撐[3?4]。自1974年Borkowska[5]首次提出概率潮流概念后，經過國內外學者近50年的大量研究，形成了多種計算方法，大致可以歸納為3 類：模擬法[6?7]、解析法[8?10]和近似法[11?13]。

模擬法以蒙特卡洛模擬法MCSM（Monte Carlo simulation method）為代表，通過簡單隨機抽樣SRS（simple random sampling）法對輸入隨機變量進行抽樣得到樣本，并對樣本中每組數據進行獨立的確定性潮流計算，最后統計輸出結果得到各節點狀態與支路潮流的概率分布。該方法的優點是簡單易行，適用于任意節點規模的系統，且在樣本規模足夠大的情況下具有較高的精度；缺點是當樣本規模以及節點規模較大時，大量的潮流計算相當耗時，效率低下，難以應用于快速分析或實時監測評估領域，常用于評判其他計算方法的優劣。

為了提升概率潮流的計算效率，國內外學者進行了大量研究。例如：解析法中的半不變量法[10]，即結合了數值特征分析手段的概率潮流分析方法；文獻[14]采用線性化的潮流方程并與半不變量法相結合，降低了計算時間；文獻[15]提出了一種利用最大熵原理改進的C 型Gram?Charlier 即CGC（C?type Gram?Charlier）級數展開法，提升了傳統CGC 方法的計算效率，但該方法中尚未考慮輸入隨機變量之間的相關性，且采用Weibull 分布描述風電場出力概率模型，對于分布特征較為復雜的概率模型可能存在擬合效果較差的情況。點估計法是近似法中的代表，其基本原理與MCSM 相同，但樣本規模相對于MCSM 更少，結合數值統計方法可使效率更高。文獻[16]提出了一種可用于計算含服從非正態分布負荷的概率潮流問題的改進型兩點估計法，其精度高于傳統的兩點估計法，但點估計法的缺點是當系統節點規模十分龐大時，潮流計算次數也隨之增加，可能出現計算效率較低的情況。

由于MCSM的計算精度較高且原理較為簡單，近年來國內外學者針對其抽樣技術以及潮流求解方法進行優化改進，以解決計算耗時長的問題。文獻[6]利用最大似然估計建立混合高斯模型，在新能源機組出力擾動及負荷需求靈活多變的情況下對其概率分布進行精確描述，建立源荷雙側統一概率模型；文獻[17]提出了擴展拉丁超立方采樣技術，改進了MCSM的抽樣手段，在減少抽樣規模的同時保證了采樣覆蓋效果；文獻[18]基于Copula理論構建多維隨機變量聯合概率分布以獲得具有相關性的數據樣本，但其仍基于傳統蒙特卡羅模擬法進行概率潮流計算，未能進一步有效地提高潮流計算效率。

針對蒙特卡洛模擬法計算概率潮流效率低下以及輸入隨機變量概率分布與相關性建模不足的問題，本文提出一種計及隨機變量相關性的多點線性化概率潮流計算方法。對于有實際監測數據的新能源場站或負荷，采用核密度估計法建立其概率分布模型，而無監測數據時則采用參數估計法獲得其概率分布。通過構建Copula 函數描述輸入隨機變量間的相關性，獲取具有相關性的隨機變量樣本。在MCSM 基礎上將非線性潮流方程轉換為線性化方程，并采用多點線性化技術減少線性化帶來的截斷誤差，以適應“雙高”電力系統中考慮強不確定性因素的計算需求。在IEEE?30 節點系統上進行算例分析，驗證所提方法的準確性和有效性。

1 源荷雙側概率分布模型的建立

由于新能源出力不確定性以及負荷隨機擾動性，“雙高”電力系統中節點注入功率具有較強的時變特征，難以應用單一的標準分布模型準確擬合其概率密度函數PDF（probability density function）以及累積分布函數CDF（cumulative distribution func?tion）[10]，而源荷雙側概率分布模型的擬合精確程度直接影響到概率潮流計算結果的精確性，因此需要運用更加準確的擬合方法對源荷概率模型進行刻畫。核密度估計KDE（kernel density estimation）[19]是一種非參數估計方法，完全根據數據本身的特點和性質來擬合分布，能夠比參數估計法得到更好的模型。本文將它應用于源荷雙側有功概率分布模型的建立。

設P1，P2，…，Pn是取自某一新能源機組輸出功率P的n個樣本，則該機組所在節點的輸出功率的概率密度函數f（P）的核密度估計定義為

核函數通常選用以0 為中心的單峰對稱概率密度函數，結合大部分常規風電、光伏機組輸出功率分布特征，可選擇高斯核函數，即

基于高斯核函數的核密度估計函數為

當核函數K定義在[?1，1]上對稱且滿足核函數約束條件[19]，則當h→0、nh→+∞時，有

特別地，當使用高斯核函數進行核密度估計時，則最優窗寬為

圖1 風電場實際輸出功率的核密度估計Fig.1 Kernel density estimation of real power output from wind farm

2 輸入隨機變量間的相關性處理

大部分概率潮流計算方法一般將源荷雙側隨機變量視為相互獨立[15]。然而實際系統運行中，由于氣候、地理等因素影響，新能源機組之間或負荷之間存在著一定相關性。若忽略這些輸入隨機變量之間的相關性，將會給概率潮流計算結果帶來誤差。Copula理論[18]能夠把隨機變量之間的相關關系與變量的邊際分布分開進行研究，不僅能完整描述輸入隨機變量間的相互關系，而且不受輸入隨機變量的邊緣分布類型及是否線性相關的限制，導出的相關性度量在線性變換和嚴格單增變換下都不會改變。為此，本文采用Copula理論描述源荷雙側輸入隨機變量之間的相關性。

本文采用橢球Copula 函數中的正態Copula 函數描述多維隨機變量間的相關性[18]。以構建二維高斯Copula函數為例，相關系數為ρvi,vj的隨機變量vi與vj的Copula函數C的表達式為

式中：Φ?1和Φρvi,vj分別為標準正態分布的累積分布逆函數以及二維正態聯合累積分布函數；ui與uj為隨機變量vi與vj的累積量。多維隨機變量間的Copula 函數可通過參數辨識求解或結合各隨機變量邊緣累積分布函數進行直接求解，詳細求解過程可參考文獻[18]。

步驟3根據等概率轉換原則，生成滿足給定相關性條件的m個具有相關性的隨機變量數據樣本Vm×N，其等概率轉換公式為

3 計及隨機變量相關性的多點線性化概率潮流計算

基于蒙特卡洛模擬的傳統概率潮流計算方法，將規模龐大的數據樣本分別進行多次獨立潮流計算，存在計算耗時、效率低下的難題，無法有效應用于實際大規模電力系統[7]。為此，本文提出了一種計及隨機變量相關性的多點線性化概率潮流計算方法，通過引入潮流方程線性化處理技術，顯著降低了蒙特卡洛模擬過程中的非線性潮流計算量，避免了單點線性化潮流計算無法適應節點注入功率過于偏離基準運行點的情況。對于“雙高”電力系統新能源出力高度不確定性和負荷側功率日內大幅度變化的情況，該方法能夠在兼顧計算精度的同時大幅提升計算效率。

3.1 潮流方程的線性化

牛頓?拉夫遜法潮流方程的節點功率方程以及支路潮流方程簡化為

式中：W為由各節點注入功率構成的矩陣；X為節點狀態矩陣；Z為支路潮流功率矩陣；f(X)和g(X)分別為節點電壓方程與支路潮流方程。

將節點注入功率W視為系統的輸入隨機變量，則潮流計算結果X與Z可視為輸出隨機變量，三者均可寫為基準值+隨機擾動值的形式，即

式中：W0、X0和Z0分別為節點注入功率、節點狀態和支路潮流的基準值；ΔW、ΔX和ΔZ為節點注入功率、節點狀態和支路潮流的隨機擾動值。

將潮流方程在系統運行的基準點進行泰勒展開，并忽略高次項，則可得潮流方程的線性化表達式[14]為

3.2 潮流方程的多點線性化

上述潮流方程的線性化處理屬于單點線性化，在系統基準運行點進行泰勒展開，在各節點的注入功率擾動范圍不大時，具有較好的精度。但高比例新能源并網后源側出力不確定性增強，使得并網節點的注入功率處于大幅度變化狀態，新能源輸出功率偏離給定期望值較大時，會導致單點泰勒展開的潮流方程與實際方程產生較大偏離[14]。此時單點線性化的潮流方程已無法滿足計算精度的需求，需要對其進行優化改進。

此時，若節點i風電場輸出功率的某一樣本數據的擾動處于第k1個區域（k1=1，2，…，n）之內，且其輸出功率全體樣本數據的均值為Pimean，處于第k2個區域（k2=1，2，…，n）之內，則可得到由該節點注入功率擾動引起的節點狀態與支路潮流擾動量的表達式為

應指出，若系統中同時存在多個需要進行多點線性化處理的節點，在對其中單個節點進行線性化處理時，其余節點注入功率按各自的基準值處理。設系統中節點i與j需要進行多點線性化處理，則先對節點i注入功率進行分段劃分，結合式（13）得到節點i引起的ΔXi與ΔZi，此時節點j處注入功率取其基準值；而后對節點j進行處理時，節點i處注入功率取其基準值，得到節點i引起的ΔXj與ΔZj。將上述節點產生的擾動量疊加即可得到多點線性化處理下節點i、j注入功率對系統各節點狀態和各支路潮流的擾動量。綜上所述，多個節點的多點線性化處理原則為：各節點的線性化處理過程相互獨立且依次進行，在針對單個節點的處理過程中不同時考慮其他節點的擾動性。

當系統中存在較多的功率處于大范圍擾動的節點，則在多點線性化處理過程中需要進行多次潮流計算，將會影響計算效率。針對該情況，可適當擴大各節點的注入功率劃分的區間寬度，通過減少分段數量來減少潮流計算次數，且對于系統中功率擾動范圍相對較小的節點，可進一步減少其節點注入功率的分段數量，或對其采用單點線性化處理。顯而易見，上述方法將提升計算效率，但計算誤差略有增大。

3.3 算法流程

本文所提出的概率潮流計算方法，主旨是對蒙特卡洛模擬法中的大量非線性潮流計算進行改進并轉化為線性化潮流計算，進一步針對單點線性化計算所存在的誤差問題進行優化，利用多點線性化處理降低線性化誤差。該算法主要流程如圖2所示。

步驟1源荷雙側概率分布模型的建立。對于有監測數據的新能源場站或負荷，通過核密度估計得到各節點注入功率的累積分布函數F及逆函數F?1；對于無監測數據的新能源場站或負荷，則采用Weibull分布或正態分布等概率分布模型。

步驟3統計各節點注入功率數據，確定系統節點注入功率基準值W0，基于基準值進行潮流計算，得到節點狀態基準X0、支路潮流基準Z0以及雅可比矩陣J0與靈敏度矩陣S0。對功率擾動范圍大的節點進行區域劃分，求解各區域對應的靈敏度矩陣Sk。式（11）中ΔW由上一步獲得的樣本Vn×N減去W0得到。

步驟4對于功率擾動范圍大的節點，先確定其節點注入功率數據所在區域，結合式（13）進行多點線性化計算；其余節點則結合式（12）進行單點線性化計算。兩者疊加得到輸出隨機變量擾動值ΔXn×N與ΔZn×N。

步驟5將X0、Z0與ΔXn×N、ΔZn×N按式（11）進行疊加，得到輸出隨機變量數據Xn×N與Zn×N，最終獲得節點電壓以及支路潮流的概率分布。

4 算例分析

4.1 測試系統

基于Matlab 平臺搭建IEEE?30 節點系統，對所提出的概率潮流計算方法可行性與有效性進行驗證。測試系統的網絡結構以及線路參數等，詳見文獻[20]。在節點6、22、27、28 處各接入30 臺額定出力為2.5 MW 的風電機組，監測數據來自中國東南沿海某4 個風電場2020 年1 月1 日—12 月31 日的8 784組有功出力數據。風電機組控制手段為恒功率因數控制。假設測試系統中各負荷服從正態分布，其期望值為原始系統數據的負荷基準值，標準差為期望值的30%。對測試系統接入的風電場輸出功率概率模型進行核密度估計擬合。以節點6與節點22 風電場輸出功率PW1與PW2概率分布為例，其PDF擬合結果如圖3所示。

圖3 核密度估計與混合高斯模型PDF 曲線比較Fig.3 Comparison of PDF curve between kernel density estimation and Gaussian mixture model

圖3 中，將監測數據實際概率密度直方圖與核密度估計所得的概率密度函數以及由最大似然估計法確定的混合高斯模型三者進行對比，可以看出，核密度估計法更好地刻畫了數據真實的分布情況，相比于混合高斯模型更為精確；且通過改變窗寬可任意調節擬合曲線的平滑度，相對于混合高斯模型有更好的泛用性。

經計算得到PW1與PW2的監測數據相關系數為0.817 4。對2 個風電場的監測數據進行等距抽樣獲得800 組離散數據，繪制出數據直方圖與散點圖，如圖4所示。圖中柱狀圖表示風電場實際輸出功率的概率密度直方圖。分析圖4 中散點分布情況可看出，2個風電場監測數據之間具有相關性，且屬于高度相關。

圖4 風電場實際輸出功率散點圖與直方圖Fig.4 Scatter plot and histogram of real power output from wind farm

結合抽樣技術獲得節點6與節點22的1 000組滿足相關性條件的風電場輸出功率樣本數據PWS1與PWS2，具體數據分布情況如圖5所示。

圖5 風電場模擬輸出功率散點圖與直方圖Fig.5 Scatter plot and histogram of simulated power output from wind farm

圖5 中柱狀圖表示風電場實際輸出功率的概率密度直方圖。分析圖5中數據分布并與圖4進行對比可看出，基于Copula 理論所獲得的節點6與節點22 風電場樣本數據間也具有與原始數據基本一致的相關性。進一步對樣本數據進行統計，得到其相關系數為0.838 8，與實際監測數據的相關系數接近，驗證了基于Copula理論生成具有相關性隨機變量樣本數據方法的有效性。

4.2 仿真計算與誤差分析

采用本文所提出的計算方法對測試系統進行仿真計算。首先對各節點風電場功率進行分段劃分。設定每段功率區間寬度為5 MW，共劃分15個區域，基準容量取100 MV·A，功率區間上、下限取標幺值，分別為[0，0.05），[0.05，0.10），…，[0.70，0.75]。以60 000次MCSM計算結果作為參考值，進行結果比較。

以節點15 電壓幅值、支路24?25 有功功率、支路6?28無功功率累積分布函數為例，2種方法所得曲線對比結果如圖6所示，圖中縱坐標RCDF表示概率的累積量。分析圖6中曲線可知，本文所提出的方法與MCSM具有十分接近的計算精度，兩者所得到的曲線結果十分接近。

圖6 節點電壓幅值或支路潮流CDF 曲線比較Fig.6 Comparison of CDF curve of node voltage magnitude or branch flow power

圖7所示為將MCSM、單點線性化（方法1）和多點線性化（方法2）3 種方法所求得的節點3 電壓PDF 與CDF 曲線同時進行對比的曲線。由圖7 可知，單點線性化方法由于僅在系統運行基準點進行泰勒展開線性化處理，當輸入隨機變量相對于基準點偏離程度增大時，所求得曲線偏離程度也隨之增大。而得益于對非線性潮流方程的多點泰勒展開，本文所提方法求得的曲線與MCSM結果基本吻合，有效地降低了由于新能源機組出力擾動范圍過大而引發的截斷誤差，提升了計算精度。

圖7 不同方法下節點電壓幅值CDF 與PDF 曲線比較Fig.7 Comparison of CDF and PDF curves of node voltage magnitude for different methods

進一步，采用方差和的根均值ARMS（average root mean square）[21]來衡量不同方法下概率潮流計算的誤差。ARMS定義為

式中：CMCi與Ci分別為MCSM與待分析計算方法所求累積分布函數上第i個數據點的值；N為數據規模。ARMS的值越小說明計算結果誤差越小。表1中列舉了部分節點電壓以及支路潮流的ARMS，并將方法1 與方法2 的精確度進行了比較。分析表1數據可見，所提出的多點線性化計算方法在高比例新能源并網節點功率大幅度變化的場景下，仍能保持較高的準確度，如：使用方法2得到的節點3電壓幅值U3的ARMS 為0.16，相較方法1 的1.82 提升了91%；節點14、19 電壓幅值與支路12?25、6?28 有功功率等計算精度提升幅度均在80%以上；多個節點電壓以及支路功率的ARMS 均能控制在0.5%以內。可見，本文所提方法相比于單點線性化方法，精度有大幅度提升。

表1 部分節點電壓幅值或支路功率的ARMSTab.1 ARMS of part node voltage magnitude or branch power

將60 000 次MCSM 所需時間與相同抽樣規模的方法1 和方法2 的進行對比，結果如表2 所示。分析表2 中數據可知，相對于MCSM 的計算時間684.83 s，本文所提方法的計算所需時間縮短至3.20 s，僅為MCSM 的0.4%，有效地提升了計算效率。相對于單點線性化的計算方法，本文所提方法在高比例新能源出力強不確定性場景下，能更有效地降低線性化所帶來的誤差，而代價則是在多點線性化處理過程中需要額外進行少量非線性潮流計算，以確定各區域所對應的靈敏度矩陣，計算時間略有增加，由方法1的2.63 s增加至3.20 s。

表2 不同方法的計算時間對比Tab.2 Comparison of computational time for different methods

5 結 論

針對蒙特卡洛模擬法概率潮流計算效率低下以及輸入隨機變量相關性考慮不足的問題，本文提出了一種考慮隨機變量相關性的多點線性化概率潮流計算方法。該方法主要有以下特點：

（1）基于核密度估計建立輸入隨機變量的概率分布模型，對其分布特征刻畫更為準確；

（2）構建Copula 函數描述多維隨機變量間的相關性，為潮流計算生成更加符合系統運行實際情況的樣本數據；

（3）對傳統蒙特卡洛模擬法中的非線性潮流計算方法進行改進，結合多點線性化技術對潮流方程進行處理，大幅度提高計算效率的同時有效減小了線性化誤差。

IEEE?30節點系統算例測試驗證了該方法相對于MCSM的計算效率有大幅提高，同時與單點線性化方法相比，計算精度也有較大提升，可適用于“雙高”電力系統考慮源荷雙側強不確定性因素下的安全運行與規劃調度等分析計算工作，為規劃調度人員提供了更全面和更準確的輔助決策信息，具有一定的工程實際應用價值。