基于分段Copula函數和高斯混合模型的多段線性化概率潮流計算

江雪辰，袁越，吳涵，徐蘊岱，黃阮明，王躍峰

(1. 河海大學能源與電氣學院，南京市 211100；2.國網上海市電力公司，上海市 200120；3.中國電力科學研究院有限公司，北京市 100192)

0 引 言

近年來，為緩解能源危機、治理環境污染，以風電為代表的可再生能源大規模集群和高滲透率分散并網，電動汽車、儲能等多樣化負荷的使用也持續增加，源荷的不確定性給電力系統的安全運行帶來巨大挑戰[1-2]。概率潮流(probabilistic power flow, PPF)[3]作為電力系統運行分析和決策的重要工具，能夠考慮各種不確定因素。因此，尋找一種能夠適應電網新環境的PPF計算方法具有重要意義。

準確建立各類輸入變量的數學概率模型是PPF計算的基礎。針對含有風電場的電力系統，文獻[4]假設風速服從雙參數威布爾分布，可快速模擬生成風速樣本，適用于缺乏實際數據的場景，但直接假定相鄰風電場的風速相關性易偏離客觀事實。已有文獻基于風速或風功率的實際歷史數據，采用單一Copula函數建立風電相關性模型[5]，方法簡單但精度有限；為提高建模精度，文獻[6]采用混合Copula函數建立風功率相關性模型，然而同組風電場在不同季節呈現的相關性程度是不同的，上述文獻采用一整段長期歷史數據構建的模型掩蓋了風電相依結構的季節性變化。基于此，本文提出一種風功率分段Copula函數，既能保證模型的準確性，又能呈現風電的季節特性。對于負荷模型，為簡化計算通常假定負荷服從標準正態分布[4,7]，而實際負荷的多樣性導致正態分布擬合存在較大誤差。文獻[8]通過建立負荷高斯混合模型(Gaussian mixture model，GMM)，有效擬合具有不對稱、多峰概率分布特性的負荷，但人工設置GMM的子成分個數限制了模型結構，且須重復進行擬合度檢驗，影響建模效率。因此，本文提出一種改進的期望最大化(expectation maximization，EM)算法建立GMM，能夠避免多次的EM計算和檢驗，在保證模型精度的同時有效提高建模效率。

在準確構建輸入變量的概率模型基礎上采用合適的PPF計算方法有助于提高計算精度與效率。現有PPF求解算法主要包括模擬法[9]、點估計法[10]及解析法[11-12]，解析法中的半不變量法(cumulant method，CM)在系統基準運行點處對潮流方程線性化以簡化計算，提高求解效率，因而被廣泛運用，但實際風功率和負荷波動范圍較廣，直接運用解析法中線性化模型求解將產生較大的誤差。文獻[11]針對離散分布的發電機注入功率進行多點線性化，但其模型缺乏普遍適用性。文獻[12]將風電出力分段線性化以減小線性化誤差，但假定風速服從相互獨立的威布爾分布且未計及實際負荷的強波動性。因此，現有的PPF算法很難同時對變量自身波動性和變量間的相關性進行有效處理。

基于上述分析，本文提出一種基于風功率的分段Copula函數和負荷GMM的多段線性化PPF計算方法。采用分段Copula函數對相鄰風電場出力進行計及季節性變化的相關性建模；通過改進EM算法確定負荷GMM的子成分個數及參數，快速準確地建立負荷GMM。考慮風電與負荷大范圍波動對PPF計算精度的影響，在所提模型基礎上，采用計及相關性的多段線性化CM進行PPF計算。以蒙特卡羅模擬法(Monte Carlo simulation method，MCSM)計算結果為基準，通過算例仿真證明所提方法計算準確且高效。

1 風電相關性模型

Copula理論[13]能夠有效處理非線性、非對稱等各類相關性問題，廣泛運用于電力系統的新能源出力建模領域。然而，像風電、光伏等受氣候影響嚴重的新能源，其出力具有更加復雜的時空相關性。采用單一Copula函數難以綜合描述多樣化的相關性結構，而混合Copula函數缺乏對相關性隨時間變化的考慮。因此，本文提出采用分段Copula對相鄰風電場出力進行計及季節性變化的相關性建模。

假設{Xi,Yi}(i=1,2,…,n)為二維隨機變量(X,Y)的樣本，分段Copula函數C(u,v)定義如下：

(1)

式中：u、v分別為隨機變量X和Y的邊緣分布函數；Ωk為原始數據按四季劃分后的樣本集合；Ik為示性函數，若{xi,yi|xi∈X,yi∈Y}∈Ωk，Ik=1，否則Ik=0；Ck(u,v,θk,Ωk)為常用的Copula函數，其類型及參數θk根據Ωk計算選取；k為季節參數，k=1、k=2、k=3、k=4分別對應春、夏、秋、冬。

構建季節相關性模型時，可采用非參數核密度估計法和極大似然估計分別求解各變量按季節分段后的邊緣分布uk=F1k(x),vk=F2k(y)和參數θk，最后采用歐氏距離[14]評估選取每段的最優Copula函數，組合生成最終的分段Copula函數。

分段Copula函數的采樣步驟如下詳述。

(1) 根據Ck(·)生成N維具有相關性的隨機向量Uk=[uk1,uk2,…,ukN]，Vk=[vk1,vk2,…,vkN]；

(2)根據原始數據求得其邊緣分布函數的逆函數xki=F1k-1(uki)，yki=F2k-1(vki)，i=1,2,…,N,進而獲得4N組滿足季節相關性和邊緣分布的隨機變量樣本{X’,Y’}，其中X’=[x11,…,x1N,x21,…,x2N,…,x41,…,x4N]T，Y’=[y11,…,y1N,y21,…,y2N,…,y41,…,y4N]T。

分段Copula模型能夠準確描繪風電相關性的季節變化特征，并且當已知分段Copula函數和某一風電場出力數據時，在條件概率下進行分層抽樣可獲取相鄰風電場的季度、年度出力。該模型可推廣至其他受氣候影響顯著的可再生能源領域。

2 負荷GMM

2.1 GMM

負荷的多樣性導致其分布難以用簡單的正態分布準確擬合，因此本文選用不受分布類型約束的GMM擬合[15]。一維隨機變量x的GMM表達式fG(x)為

(2)

2.2 基于改進K-means聚類的EM優化算法

通常須循環采用EM算法[8]求解不同n值下的GMM參數，并通過比較擬合度來選取最佳GMM，過程繁瑣。因此，為了提高建模效率，本文引進聚類有效性指標DB(-)[16]，提出一種基于改進K-means聚類的EM優化算法，主要步驟如下詳述。

(1)初始優化。考慮后續建模速度，設置高斯子成分個數n的取值范圍為2～10。在該范圍內采用K-means算法循環計算不同n值下的聚類結果及其DB(-)指標，選擇DB(-)最小值對應的n值作為最終GMM的子成分個數。各子成分的權重、期望和標準差初值則由該n值下的聚類結果決定。

(2)E-步。計算樣本容量為NL的原始數據xi(i=1,2,…,NL)由第k個子成分生成的概率p(i,k)：

(3)

(3)M-步。更新GMM各子成分的權重、期望和標準差：

(4)

(5)

(6)

(4)收斂性判斷。計算GMM的最大對數似然函數式L：

(7)

若式(7)收斂，計算結束；否則，重復迭代E-步和M-步。文獻[17]已證明DB(-)指標適合作為電力負荷曲線聚類的有效性指標。因此，上述改進的EM算法適用于負荷GMM的建立。

3 基于分段Copula函數和GMM的概率潮流計算方法

3.1 半不變量法概率潮流

傳統CM在輸入變量相互獨立的前提下，借助線性交流模型獲得輸出變量的概率分布特性。在確定性潮流所求基準點處對其進行泰勒級數展開，忽略2次及以上的高階項，得到線性化潮流方程：

(8)

式中：X、Z分別為節點狀態變量和支路狀態變量；X0、Z0分別為基準點的節點狀態變量和支路狀態變量；ΔW為節點注入功率變化量；J0為基準運行點的雅克比矩陣；G0為支路潮流對節點電壓的偏導；S0、T0為靈敏度矩陣。

為方便求得狀態變量的概率分布，在式(8)的基礎上將變量間復雜的卷積運算轉換為半不變量的計算：

(9)

傳統CM未計及輸入變量的相關性，并且當輸入變量波動范圍較大時，單點線性化的傳統CM會產生較大的誤差，而本文通過分段Copula函數和GMM分別獲得了具有相關性的風功率樣本和大范圍波動的負荷樣本，傳統CM不再適用，因此引入一種計及相關性的多段線性化處理方法。

3.2 潮流方程的多段線性化處理

根據輸入變量的波動范圍進行適當的分段，選取合適的基準點并在各段基準點處對潮流方程線性化。此處以實際負荷GMM為例進行說明，具體處理步驟如下詳述。

(1)假設k為某一實際負荷節點，依據其功率波動范圍合理均分為m段：Wk0～Wk1、Wk1～Wk2、…、Wk(m-1)～Wkm，其中Wki為功率波動區間的起止值，設各段功率期望值為Ek1，Ek2，…，Ekm；標準算例其他n-1個節點注入功率的期望分別為E1，E2，…，Ek-1，Ek+1，…，En。

(2)求系統m段的運行點At(t=1,2,…,m)。

At=(E1,…,Ek-1,Ekt,Ek+1,…,En)

(10)

式中Ekt為實際負荷節點的m段期望值。

此處除實際負荷節點外其他節點期望值不變。

(3)將潮流方程在At處線性化，得到：

ΔX=StΔW

(11)

式中St為確定性潮流在At處所得靈敏度矩陣。

當實際負荷注入功率在第t段期望值附近波動，其變化量ΔWk對狀態量變化量ΔXi具有一定的影響，影響系數S0ik根據St求得，潮流方程展開如下：

ΔXi=S0i1ΔW1+…+S0i(k-1)ΔWk-1+StikΔWk+

S0i(k+1)ΔWk+1+…+S0inΔWn

(12)

以上是以節點電壓為例分析，支路潮流與節點注入量的線性關系類似。對于功率波動范圍較大的實際負荷節點或者風電場，一次只考慮單個變量的注入功率變化，其余變量按期望值進行計算。

多段線性化處理破壞了狀態變量與輸入變量的各階半不變量的齊次性，不能簡單地根據式(9)計算，須進行修正[11-12]，修正方法如下詳述。

假設隨機變量Q與R之間存在如下的分段線性關系，本文即指狀態變量與輸入變量的分段線性化潮流方程：

Q=aR+b

(13)

R∈[r0,rm]。a、b取值：

(14)

(15)

(16)

根據式(16)，即可遞推依次得到系統多段線性化后隨機變量Q的k階半不變量κk：

(17)

3.3 隨機變量的相關性處理

由于本文風功率具有相關性，不滿足半不變量可加性的前提，因此在多段線性化CM計算前必須進行相關性處理,將具有相關性的輸入變量轉化為互不相關的變量。在式(8)的基礎上得到某一分段t下的表達式：

(18)

式中：ΔX(t)、ΔZ(t)分別為第t段的節點狀態變量變化量和支路狀態變量變化量；S0(t)、T0(t)分別為第t段的節點靈敏度矩陣；ΔW(t)為第t段的風電注入功率變化量；σ(t)為ΔW(t)的標準差；ΔW*(t)為ΔW(t)標準化后得到的風功率變化量。

(19)

(20)

根據式(18)—(20)，對某一分段t下靈敏度矩陣S0(t)進行修正：

(21)

式中S1(t)、T1(t)為第t段的靈敏度矩陣修正值。

對于多個風電和負荷波動因素，相關性處理后利用半不變量的可加性，將單個波動因素影響下的結果疊加即可求解最終狀態變量的各階半不變量[12]。本文引入C型Gram-Charlier級數展開擬合狀態變量的概率分布，有效避免概率密度函數(probability density function，PDF)出現負值，加快級數收斂。

3.4 算法流程

本文基于風功率的季節相關性模型和負荷GMM，采用多段線化求解PPF，具體流程如圖1所示。

圖1 基于分段Copula函數和GMM的PPF計算流程Fig.1 Flow chart of PPF calculation based on piecewise Copula and Gaussian mixture model

4 算例分析

4.1 算例說明

采用改進的IEEE 14節點系統驗證本文所提風電和負荷模型的準確性并說明不確定性對系統運行的影響。針對IEEE 14節點系統作如下改動：將某地區額定容量為20 MV·A的2個相鄰風電場A和B分別接入原系統節點13、14；以實際負荷節點K、M替代原系統節點7、13。系統的結構示意圖如附錄圖A1所示，根據2013年的實際風功率和負荷數據進行算例分析。

4.2 模型驗證

4.2.1風功率相關性模型

對風電場A和B的風功率進行季節性分段，按照1.2節所述最終選取的各段Copula函數及其參數見表1，其中各段具體參數見附錄表A1。

表1 各季節最佳Copula函數及其參數Table 1 The best Copula and its parameters for each season

注：季節1～4依次為春、夏、秋、冬。

單一Copula函數及其參數見附錄表A2。由附錄表A2可知，采用單一Copula函數建模時，t-Copula函數的擬合效果最好，對應的歐氏距離為2.051 7，秩相關系數為0.686 8，而分段Copula函數的歐氏距離之和僅為1.843 7；進一步，將風電場A的模擬數據和實際數據進行對比，結果如圖2所示。由圖2可知，分段Copula函數擬合原始數據的效果優于單一t-Copula函數。

圖2 風電場A出力概率密度曲線Fig.2 PDF of wind power output

由分段Copula函數的選取種類可知，僅秋季風功率的相關性符合對稱厚尾的t-Copula函數結構，其余季節風電出力存在著上尾獨立性，僅憑單一t-Copula函數無法描述。此外，由秩相關系數結果可知，春、秋、冬三季風電呈顯著相關，夏季風電呈低度相關，而單一t-Copula函數強化了夏季風電相關性，高估了2個風電場間的相互影響。

4.2.2負荷GMM

通過2.2節所述方法建立節點7、13的負荷有功GMM，其高斯子成分個數按照附表A3所示的聚類結果設為3和5，各子成分參數見表2。

表2 有功負荷的GMM參數Table 2 GMM parameters of the active power load

負荷有功功率模擬樣本和實際數據的PDF曲線如圖3所示。由圖3可知，2組數據總體擬合良好；經計算，模擬樣本與實際數據的期望值誤差小于0.03%，標準差誤差不超過1.2%；此外，基于改進K-means聚類的EM算法平均耗時為6.232 2 s，而傳統EM算法平均耗時29.941 9 s。因此，本文所建負荷GMM具有較高精度，且建模過程只需進行1次聚類檢驗和1次EM計算，提升了建模效率。

圖3 有功負荷概率密度曲線Fig.3 PDF of active load

4.3 算法性能測試

在上述模型的基礎上，采用本文所提算法對改進后的IEEE 14節點系統進行仿真分析。假設除負荷節點7、13外，其余節點負荷均服從正態分布，以原系統數據作為期望值，標準差為期望值的10%，本文暫不考慮節點負荷之間的相關性。

利用風功率分段Copula模型和負荷GMM分別生成樣本容量為30 000的風功率樣本和負荷樣本，并將風功率樣本分為2段，負荷樣本分為6段。

以30 000次MCSM所得的PPF結果作為基準，常規CM計算的PPF結果作為對比，采用狀態變量期望值和標準差的相對誤差指標[5]來衡量本文所提方法和常規CM方法的準確程度，結果見表3、4。

由表3、4可知，采用常規CM所得狀態變量數字特征的相對誤差指標的平均值最大為 22.010 5%，最大值最大為38.488 9%，而采用本文算法的相應結果分別為2.564 0%和8.498 4%，誤差明顯縮小。采用本文所提算法進行PPF計算平均耗時為15.18 s，而采用MCSM平均耗時1 687.36 s。由此說明在相同輸入條件下，采用本文所提算法計算PPF更加準確、快速。

表3 本文方法與MCSM的狀態變量數字特征的相對誤差Table 3 The relative error of the digital characteristic of the state variables between proposed method and MCSM

表4 常規CM與MCSM的狀態變量數字特征的相對誤差Table 4 The relative error of the digital characteristic of the state variables between normal CM and MCSM

為進一步說明風功率相關性模型、負荷GMM模型以及本文算法的適用性和準確性，設置5種運行場景進行分析比較，如下詳述。

場景一：考慮風功率分段Copula相關性模型和負荷GMM，采用MCSM進行PPF計算。

場景二：考慮風功率分段Copula相關性模型和負荷GMM，采用常規CM進行PPF計算。

場景三：考慮風功率分段Copula相關性模型和負荷GMM，采用本文所提算法進行PPF計算。

場景四：考慮風功率實際歷史數據和實際負荷數據，采用本文所提算法進行PPF計算。

場景五：考慮風功率單一t-Copula相關性模型和負荷正態分布的簡化模型，采用本文所提算法進行PPF計算。

以節點14電壓幅值和支路6—12有功功率的結果進行對比說明，結果如圖4、5所示。

根據圖4、5，對比場景一、場景二和場景三，采用相同風功率分段Copula相關性模型和負荷GMM，本文所提算法與MCSM所得曲線相似性較高，而常規CM與MCSM所得曲線偏差較大，說明在風功率和負荷波動范圍較大時，采用本文算法計算PPF結果更加準確；對比場景三、場景四和場景五，同樣采用本文所提算法，基于本文所述風功率及負荷模型的PPF結果，與實際數據計算所得PDF曲線幾乎重合，而場景五采用簡化模型計算的PDF曲線與場景三、場景四的偏差較大，再次驗證采用本文所述模型刻畫風功率相關性和負荷波動更加準確。

圖4 節點14電壓幅值概率密度曲線Fig.4 PDF of voltage magnitude at node 14

圖5 支路6—12有功功率概率密度曲線Fig.5 PDF of active power on branch 6-12

此外，采用常規CM的場景二和采用簡化模型的場景五中，節點電壓高低壓段和支路有功功率高值段的越界概率被低估，而支路有功功率低值段的越界概率被高估。模型和算法的誤差都會影響最終PPF計算的準確性，造成系統運行風險的誤判，從而導致運行和規劃人員決策失誤。

5 結 論

本文提出一種基于風功率分段Copula模型和負荷GMM的多段線性化PPF計算方法。理論和算例表明：

(1)分段Copula模型能夠清晰地描述風電出力在不同季節的相關結構，準確模擬相鄰風電場的季度、年度出力；

(2)采用改進K-means聚類的EM優化算法能夠快速準確地求解GMM的子成分個數及其參數，在保證模型精度的同時提高了建模效率；

(3)在獲得準確的風電和負荷數據后，采用Cholesky分解解決CM中的風功率相關性問題并對潮流方程進行多段線性化處理，有效降低了風功率與負荷的大范圍波動造成的PPF計算誤差，快速求解的同時，明顯提高了計算精度；

(4) 本文算法綜合考慮了輸入隨機變量的相關性和波動性，能夠快速準確地揭示系統運行中的薄弱環節和潛在風險，為電網運行和規劃提供有效的參考指標，在大規模風電等新能源并網的電力系統中具有良好的工程應用價值。