朱 溪，蘇晨博，劉教民

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

0 引言

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

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

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

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

1基于RVM的多維Copula函數(shù)模型

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

式中：K（xi，xk）為核函數(shù)；序列ω=（ω1，ω2，…，ωn）T為每個(gè)核函數(shù)所對(duì)應(yīng)的權(quán)重系數(shù)；ωik為對(duì)應(yīng)于輸入矢量xi的第k個(gè)權(quán)重系數(shù)；ωi0為權(quán)重系數(shù)的偏差量（即常數(shù)項(xiàng)）；εi為與目標(biāo)值之間的誤差。

權(quán)重系數(shù)ω的條件分布表達(dá)式如下：

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

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

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

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

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

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

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

式中：ηi=［F1i（v1），F(xiàn)2i（v2），…，F(xiàn)Ni（vN）］T為風(fēng)速序列i的邊緣分布矢量。

利用計(jì)算出的多維風(fēng)場(chǎng)風(fēng)速的聯(lián)合分布計(jì)算得到相關(guān)系數(shù)矩陣ρ，結(jié)合Nataf變換[24]，獲得計(jì)及相關(guān)性的獨(dú)立風(fēng)速分布序列，然后根據(jù)風(fēng)速和風(fēng)電功率的函數(shù)關(guān)系得到獨(dú)立的風(fēng)功率序列。

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

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

2.1 電力網(wǎng)絡(luò)模型

電力系統(tǒng)采用經(jīng)典交流系統(tǒng)模型，其直角坐標(biāo)系下PQ和PV節(jié)點(diǎn)表達(dá)式為

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

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

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

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

將式（9）代入式（8），得到節(jié)點(diǎn)功率的復(fù)冪級(jí)數(shù)展開式：

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

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

將式（7）帶入式（12），并假設(shè)電壓復(fù)冪級(jí)數(shù)的常數(shù)項(xiàng)ai0=1.0，則其余項(xiàng)的系數(shù)可以表示為

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

結(jié)合式（9），（12）和（14），計(jì)算得到節(jié)點(diǎn)電壓復(fù)冪級(jí)數(shù)Vi（s）和Qi（s）的值。

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

在熱力系統(tǒng)模型中，水為主要的液體熱媒，在管道中實(shí)現(xiàn)熱能的傳遞與轉(zhuǎn)換。各管道中的熱水在傳輸過程中滿足流量連續(xù)性方程，根據(jù)能量守恒定律可得，熱水在熱力網(wǎng)絡(luò)中各個(gè)管道的壓頭損失的總和應(yīng)該為零?；诖?，則有：

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

針對(duì)熱力網(wǎng)絡(luò)中的負(fù)荷節(jié)點(diǎn)，基于節(jié)點(diǎn)供熱溫度Ts和節(jié)點(diǎn)輸出溫度To，可得熱力網(wǎng)絡(luò)的模型如下[22]：

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

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

然后假設(shè)Ni（s）為Mi（s）的倒數(shù)，且其復(fù)數(shù)域級(jí)數(shù)展開式的系數(shù)有如下關(guān)系：

最后將式（18）和（19）代入式（17），得到熱力網(wǎng)絡(luò)各支路流量的復(fù)冪級(jí)數(shù)展開式。

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

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

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

式中：cm為熱電比例系數(shù)，本文假設(shè)其為常數(shù)。

若節(jié)點(diǎn)k為電-熱耦合節(jié)點(diǎn)，則熱力網(wǎng)絡(luò)功率方程須修正為

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

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

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

式中：μ為X的平均值。

進(jìn)而推出v階半不變量κv和αv的關(guān)系為

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

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

式中：G0和S0為由式（27）構(gòu)成的系數(shù)矩陣。

式中：Zl0=Z0-G0S0W0。

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

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

3 算例分析

本文選取的IES系統(tǒng)如圖1所示，其中：電氣網(wǎng)絡(luò)節(jié)點(diǎn)數(shù)為33個(gè)；風(fēng)電場(chǎng)為4個(gè)，各個(gè)風(fēng)電場(chǎng)中風(fēng)電機(jī)組的額定功率為80MW，其運(yùn)行參數(shù)為中國(guó)張北地區(qū)采集的風(fēng)場(chǎng)的實(shí)際歷史數(shù)據(jù)。本算例設(shè)定平衡節(jié)點(diǎn)的電壓標(biāo)幺值為1 p.u.，功率因數(shù)為0.98。除此之外，系統(tǒng)還有13個(gè)熱力網(wǎng)絡(luò)節(jié)點(diǎn)，其中包含兩臺(tái)電鍋爐。

圖1 電-熱網(wǎng)絡(luò)結(jié)構(gòu)圖Fig.1 Electrical-Heart structure network

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

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

表1 風(fēng)電機(jī)群風(fēng)速聯(lián)合分布對(duì)比Table 1 Comparison between differentmethods

續(xù)表1

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

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

式中：vwi為切入風(fēng)速；vwo為切出風(fēng)速；vr為額定風(fēng)速；Pr為風(fēng)電機(jī)組的額定功率。

本文計(jì)算時(shí)參考實(shí)際風(fēng)電場(chǎng)，設(shè)vwi為4m/s，vwo為25m/s，vr為12m/s。因?yàn)槟Ｐ椭械妮斎胱兞烤哂须S機(jī)不確定性，通過所提方法計(jì)算得到的能流結(jié)果為概率分布的形式。利用搜集到的網(wǎng)絡(luò)歷史數(shù)據(jù)作為參考值，計(jì)算輸出變量的期望值和標(biāo)準(zhǔn)差的相對(duì)誤差。

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

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

表2 電力網(wǎng)絡(luò)概率能流誤差結(jié)果Table 2 The PEF results of power system

表3 熱力網(wǎng)絡(luò)概率能流誤差結(jié)果Table 3 The PEF results of heartnetwork

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

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

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

4 結(jié)論

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

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

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