基于RVM和全純嵌入法的考慮多風電場出力相關性的電-熱互聯(lián)概率能流計算方法

朱 溪，蘇晨博，劉教民

（華北電力大學 電氣與電子工程學院，北京102206）

0 引言

隨著電網中風電等可再生能源數量的不斷攀升，新能源的不確定性對IES運行狀態(tài)的影響愈發(fā)突出[1]，[2]。目前，主要采用概率能流（Probabilistic Energy Flow，PEF）研究不確定性對IES的影響[3]。PEF分析方法主要包括蒙特卡羅模擬法（Monte Carlo Method，MCM）[4]、半不變量法[5]和點估計法[6]。半不變量法根據輸入變量的半不變量結合級數展開求取輸出量的概率密度函數，該方法既能維持較高的運算效率，又能保證計算結果的精度[7]～[9]。然而傳統(tǒng)的半不變量計算方法都是基于NR法，當電-熱聯(lián)合系統(tǒng)包含大量新能源時，該算法容易出現迭代不收斂的情況。

為了解決上述問題，SD Rao提出了一種新的非迭代方法，稱為全純嵌入法（Holomorphic Embedding Method，HEM）[10]。HEM具有一個顯著的優(yōu)點：如果系統(tǒng)存在解，則該方法能夠保證最終找到一個可以運行的解；相反，如果系統(tǒng)沒有解，則計算出的有理近似值會產生數值振蕩[11]，[12]。

風機的輸出功率受到風速的影響，且同一地理區(qū)域內的風速存在非線性相關性[13]，因此在計算含風電場的電-熱聯(lián)合系統(tǒng)的概率能流時要充分考慮這種時空相關性[14]。根據Copula函數的特性和原理，多元變量的相關性模型可利用二元Copula的概率密度函數和邊緣概率密度函數相乘獲得[15]。但是特定的二元Copula很難準確地刻畫風速這類復雜變量的分布特性[16]～[18]。

為此，本文首先提出一種基于RVM的多維風場風速的聯(lián)合分布計算方法，該方法采用核函數能有效地避免參數估計帶來的二次誤差，從而提高計算精度，并且其權重系數矩陣為稀疏陣，能夠大幅度提升計算效率。其次，基于HEM方法提出一種考慮風電出力相關性的計算電-熱互聯(lián)系統(tǒng)概率能流的新算法。相較于傳統(tǒng)的NR法，該算法具有更好的收斂性和準確性。

1基于RVM的多維Copula函數模型

RVM以基函數為框架，以貝葉斯理論中的先驗分布和極大似然估計為基礎，將低維度的非線性問題轉化為高維空間的線性問題，并采用自動相關決策理論（Automatic Relevance Determination，ARD）來約束模型，進而獲得更加稀疏和精確的結果。假設給定訓練集{xi，ti}，其中，xi和ti分別表示輸入樣本矢量和相應的目標值；N為樣本總數，則RVM的數學模型可以表示為

式中：K（xi，xk）為核函數；序列ω=（ω1，ω2，…，ωn）T為每個核函數所對應的權重系數；ωik為對應于輸入矢量xi的第k個權重系數；ωi0為權重系數的偏差量（即常數項）；εi為與目標值之間的誤差。

權重系數ω的條件分布表達式如下：

式中：φ=[φ（x1），φ（x2），…，φ（xN）]為基函數序列，通過核函數構建，它們之間的關系可以表示為φ（xi）=[1，K（xi，x1），K（xi，x2），…，K（xi，xN）]T。

在本文中，輸入矢量是從不同的風速序列V=（v1，v2，…，vN）T中采樣所得。

根據統(tǒng)計學理論可知，極大似然法可以直接估計ω和σ2的值，但是很容易過擬合，導致擬合不準確。因此，采用標準正態(tài)分布的權重先驗概率分布來約束這兩個參數。根據貝葉斯理論可知，后驗概率分布等于加權先驗分布乘以極大似然估計[20]，則后驗概率分布p（ω|t，σ2，α）的表達式如下：

式中：A=diag（αi）；σ為樣本矢量xn的標準差；μ為樣本均值矢量；Σ為樣本協(xié)方差矩陣。

本文中的目標函數t（t1，t2，…，tn）T是從多維經驗Copula序列Cn中采樣獲得。根據R-Vine Copula理論[19]，多維經驗Copula函數可以由二維經驗Copula函數Cb、邊緣概率密度fi（vi）和條件分布fi|j（vi|vj）共同構建[21]。由RVM的定義可知，其回歸精度主要取決于核函數的種類。為了能夠更加準確地描述風速這類非線性變量的聯(lián)合分布，本文采用高斯核函數作為RVM的基函數[21]，因此，式（3）中的α和σ2可以由迭代獲得，分別為

式中：μi為后驗均值μ的第i個分量；∑ii為后驗協(xié)方差Σ中的第i個對角線分量。

將迭代得到的α和σ2代入式（3），可以得到權重系數ω的修正值。因此，基于RVM的多元Copula可以表示為

式中：ηi=［F1i（v1），F2i（v2），…，FNi（vN）］T為風速序列i的邊緣分布矢量。

利用計算出的多維風場風速的聯(lián)合分布計算得到相關系數矩陣ρ，結合Nataf變換[24]，獲得計及相關性的獨立風速分布序列，然后根據風速和風電功率的函數關系得到獨立的風功率序列。

綜上，RVM有兩個關鍵優(yōu)勢：①利用核函數為基本框架代替參數估計，能夠更加精確地刻畫風速這類非線性變量的聯(lián)合分布；②該模型中，每個權重參數ωi均有一個單獨的超參數αi，而不需要與其它權重參數共享超參數。在計算過程中，參數矩陣α中的大部分元素將增至無窮大，此時ω的后驗分布大量集中在零處，保證了權重矩陣的稀疏性，因此，該方法既能確保計算精度又能提高計算速度。

2 基于全純嵌入法的概率能流計算方法

2.1 電力網絡模型

電力系統(tǒng)采用經典交流系統(tǒng)模型，其直角坐標系下PQ和PV節(jié)點表達式為

式中：Si為節(jié)點i注入的功率，Si=Pi+jQi，Pi，Qi分別為有功功率和無功功率；Vj為節(jié)點j的電壓；Vi為節(jié)點i的電壓；V為平衡節(jié)點電壓幅值；Yij為導納矩陣。

HEM方法通過構造狀態(tài)量復數域的級數展開求解系統(tǒng)潮流問題，其表達式為[10]

式中：Yij，trans和Yi，shunt分別為導納矩陣Y的串聯(lián)支路部分和并聯(lián)支路部分；a*in為復冪級數系數ain的共軛；Vi（s）為電壓復冪級數。

為了求解上式中的級數序列，將Vi（s）的倒數表示為Wi（s）=bi0+bi1s+…+binsn。Vi（s）和Wi（s）的乘積可以表示為其系數的卷積，系數矩陣ain和bin的關系為

將式（9）代入式（8），得到節(jié)點功率的復冪級數展開式：

通過式（9）和（10）定義的遞推關系計算得到系數ain和bin，然后進一步計算PQ節(jié)點的電壓序列。PV節(jié)點的電壓幅值和有功功率均為已知量，而無功功率和節(jié)點電壓相角為未知量。假設節(jié)點i的無功功率為Qi，則其所對應的復冪級數展開式為

將式（9）帶入式（11），得到PV節(jié)點電壓和無功功率的復冪級數展開式：

將式（7）帶入式（12），并假設電壓復冪級數的常數項ai0=1.0，則其余項的系數可以表示為

式中：ain，re為電壓復冪級數的系數實部值；δn0和δn1為示性函數，其值為0或1。

結合式（9），（12）和（14），計算得到節(jié)點電壓復冪級數Vi（s）和Qi（s）的值。

2.2 熱力系統(tǒng)模型

在熱力系統(tǒng)模型中，水為主要的液體熱媒，在管道中實現熱能的傳遞與轉換。各管道中的熱水在傳輸過程中滿足流量連續(xù)性方程，根據能量守恒定律可得，熱水在熱力網絡中各個管道的壓頭損失的總和應該為零。基于此，則有：

式中：As為熱力網絡流量方向的關聯(lián)矩陣；m為各管道熱媒的流量；mq為各管道節(jié)點中流出的熱媒流量；Bh為管道的回路-支路關聯(lián)矩陣；hF為熱媒的壓頭損失向量；K為管道的阻力系數矩陣；mi為節(jié)點i流入的總流量；mj為管道j的流量。

針對熱力網絡中的負荷節(jié)點，基于節(jié)點供熱溫度Ts和節(jié)點輸出溫度To，可得熱力網絡的模型如下[22]：

式中：Φ為熱力網絡中流過各節(jié)點的熱功率；Cp為熱媒介的比熱容；Ta為熱力網絡外部的環(huán)境溫度；Tend為各條管道的末端溫度；λ為管道的熱傳導系數；L為各條管道的長度。

通過HEM方法，構造熱力網絡狀態(tài)量的復數域級數展開求解熱力網絡模型。首先通過熱力網回路-支路關聯(lián)矩陣求得hF的值，將式（15）展開為復數級數形式：

然后假設Ni（s）為Mi（s）的倒數，且其復數域級數展開式的系數有如下關系：

最后將式（18）和（19）代入式（17），得到熱力網絡各支路流量的復冪級數展開式。

因此，可以通過式（19）和（20）定義的遞推關系計算得到系數Min和Nin，進一步用式（16）計算出各節(jié)點的熱力功率。

2.3 電熱耦合環(huán)節(jié)

IES系統(tǒng)中存在著很多耦合環(huán)節(jié)，包括CHP機組、電鍋爐和熱鍋爐等，本文只考慮電-熱機組的耦合，其中電功率P和熱功率Φ的關系為

式中：cm為熱電比例系數，本文假設其為常數。

若節(jié)點k為電-熱耦合節(jié)點，則熱力網絡功率方程須修正為

將修正后得到的節(jié)點k的熱功率帶入式（8）和（11），對電力系統(tǒng)網絡的功率方程進行修正，由此可計算得到考慮耦合元件的電-熱網絡功率。其修正后的表達式如下：

2.4 基于全純嵌入法的半不變量模型

假設隨機輸入變量X的概率密度函數（Probability Density Function，PDF）和累積分布函數（Cumulative Distribution Function，CDF）分別為f（x）和F（x），則X的v階矩和v階中心矩的數學表達式為

式中：μ為X的平均值。

進而推出v階半不變量κv和αv的關系為

傳統(tǒng)的半不變量是基于牛頓-拉夫遜算法，為了將其更好地應用于全純嵌入方法，須要做進一步改進。假設Pˉ，Qˉ和hˉF分別為電-熱聯(lián)合網絡的有功功率、無功功率和節(jié)點壓頭損失向量的平均值，則聯(lián)合網絡能流的平均值模型為

當s=0時，式（26）的解與無負載和無分流元件的聯(lián)合網絡相對應，其解可作為復數域冪級數的常數項系數。若Padé[10]逼近在s=1時的解不發(fā)生數值振蕩，則系統(tǒng)收斂，此時系統(tǒng)狀態(tài)量的均值可以表示為各系數的和：

式中：G0和S0為由式（27）構成的系數矩陣。

式中：Zl0=Z0-G0S0W0。

因此，可以通過級數展開來計算電-熱網絡中狀態(tài)量的PDF。與Gram-Charlier級數和Edgeworth級數相比，Cornish-Fisher級數在處理非正態(tài)分布的隨機變量時具有更高的精度。本文采用Cornish-Fisher級數[23]計算熱力網絡狀態(tài)量的概率密度。根據Cornish-Fisher級數展開理論，可得其前五階的數學表達如下：

式中：Y（q）為輸出的分位數；τ（q）為標準正態(tài)分布的分位數，滿足τ（q）=Φ-1（q）。

3 算例分析

本文選取的IES系統(tǒng)如圖1所示，其中：電氣網絡節(jié)點數為33個；風電場為4個，各個風電場中風電機組的額定功率為80MW，其運行參數為中國張北地區(qū)采集的風場的實際歷史數據。本算例設定平衡節(jié)點的電壓標幺值為1 p.u.，功率因數為0.98。除此之外，系統(tǒng)還有13個熱力網絡節(jié)點，其中包含兩臺電鍋爐。

圖1 電-熱網絡結構圖Fig.1 Electrical-Heart structure network

熱-電互聯(lián)IES概率能流的計算過程如下。

①根據張北地區(qū)的實際采集數據，可以獲得一定時間間隔的風速序列。基于風速的歷史數據，利用核密度估計的方法計算得到各個風場中風速的邊緣分布；然后基于RVM算法和二元經驗Copula函數以及D-Vine結構計算四維風電場風速的聯(lián)合概率密度以及聯(lián)合分布。為了驗證所提方法的精確性，本算例采用K-S檢驗和歐式距離對比分析不同算法[最小二乘法（Ordinary Least Squares，OLS）、最大期望值估計法（Expectation-Maximization Algorithm，EM）和RVM算法]得出的聯(lián)合分布與經驗分布之間的差異。改變風速的采樣間隔，分別用每隔6 h，每隔1 h和每隔10min的風速序列作為采樣空間，用以驗證算例結果的普適性。本算例在采樣風速數據時，選取了5%的顯著水平，所得的d值和K值結果如表1所示。在K-S檢驗中，K值和歐式距離d越小，代表計算得到的各個風場的聯(lián)合概率分布與實際各個風場的聯(lián)合概率分布越接近。從表1中可以看出，所提算法擬合得到的聯(lián)合分布函數的d值和K值均為最小。結果表明，本文所提算法能夠更加精確地計算多維風電場風速的聯(lián)合分布。

表1 風電機群風速聯(lián)合分布對比Table 1 Comparison between differentmethods

續(xù)表1

注：CRVM（v1，v2，v3，v4）表示用RVM算法計算出的聯(lián)合分布函數；COLS（v1，v2，v3，v4）表示采用OLS算法計算出的聯(lián)合分布函數；CEM（v1，v2，v3，v4）表示用EM算法計算出的聯(lián)合分布函數。

②利用計算出的多維風場風速的聯(lián)合分布計算得到相關系數矩陣ρ，結合Nataf變換[24]得到計及相關性的獨立風速分布序列，帶入式（31）計算得到獨立風功率序列Pw，然后進一步采用所提出的基于HEM算法的半不變量計算概率能流。

式中：vwi為切入風速；vwo為切出風速；vr為額定風速；Pr為風電機組的額定功率。

本文計算時參考實際風電場，設vwi為4m/s，vwo為25m/s，vr為12m/s。因為模型中的輸入變量具有隨機不確定性，通過所提方法計算得到的能流結果為概率分布的形式。利用搜集到的網絡歷史數據作為參考值，計算輸出變量的期望值和標準差的相對誤差。

考慮到電-熱互聯(lián)IES概率能流模型的輸出變量數量較大，本文計算了每類輸出變量的相對誤差的平均值AVG以及最大值MAX，來量化和評估所提算法的精確性。

本文分別利用基于考慮風電相關性的NR法的半不變量與基于考慮風電相關性的HEM法的半不變量，計算電-熱網絡狀態(tài)量的概率密度并進行比較，結果如表2，3所示。結果表明，基于HEM法的半不變量計算得到的狀態(tài)量期望值和標準差的平均誤差小于基于NR法得到的平均誤差，驗證了所提算法的精確性；該方法將平均相對誤差限制在4%以內，最大相對誤差限制在10%以內；標準差的誤差大于期望值的誤差，這也符合半不變量法的特點。

表2 電力網絡概率能流誤差結果Table 2 The PEF results of power system

表3 熱力網絡概率能流誤差結果Table 3 The PEF results of heartnetwork

電-熱系統(tǒng)的節(jié)點電壓和管道流量期望誤差分別如圖2，3所示，隨著Padé逼近階數的增加，所提算法的精度也在提高。然而，Padé逼近的階數越高所需要的計算時間也越高，因此，在確保計算結果準確的前提下，應盡量減少其階數。

圖2 電-熱系統(tǒng)節(jié)點電壓期望誤差Fig.2 The error of node voltage in IES

圖3 電-熱系統(tǒng)管道流量期望誤差Fig.3 The error of pipeline flow in IES

4 結論

本文提出了一種基于RVM的多維風電場風速的聯(lián)合分布計算方法。在此基礎上，提出了一種基于HEM方法的新半不變量法用于電-熱聯(lián)合網絡的概率能流計算，該算法具有更好的收斂性和計算精度。經過采用修改的IEEE 33節(jié)點系統(tǒng)結合熱力網絡13節(jié)點系統(tǒng)進行算例分析，得出以下結論。

①在分析多能源系統(tǒng)的運行狀態(tài)時，應該考慮風速之間的相關性，其對系統(tǒng)能流具有重要影響。與傳統(tǒng)的參數估計方法相比，基于RVM的多元Copula函數以核函數為基礎，能夠有效地避免參數估計帶來的二次誤差，計算結果也更加接近真實風速的聯(lián)合分布。

②基于HEM和半不變量法的概率能流分析方法實現了對傳統(tǒng)PEF計算方法的改進，具有更準確的概率分布計算結果。與NR方法相比，它在冪級數收斂半徑以內或以外都能提供較好的收斂性，能夠保證更好的計算精度。