曾　江， 蔡?hào)|陽(yáng)， 黃德華

(華南理工大學(xué)電力學(xué)院， 廣東省廣州市 510641)

0　引言

隨著大量非線性負(fù)載的接入，電網(wǎng)中的諧波污染日益嚴(yán)重，對(duì)電網(wǎng)中諧波傳播與分布的研究工作具有重要意義，諧波潮流計(jì)算是研究諧波分布特性的有效手段[1-3]。傳統(tǒng)諧波潮流計(jì)算旨在求解特定工況下電網(wǎng)各節(jié)點(diǎn)基波、諧波電壓以及支路潮流的確定值。但實(shí)際電網(wǎng)中存在著諸多不確定因素，如負(fù)荷變化、系統(tǒng)隨機(jī)故障、網(wǎng)絡(luò)結(jié)構(gòu)變化等。此外，新能源發(fā)電技術(shù)的快速發(fā)展也加劇了電力系統(tǒng)的不確定性[4-5]，確定性諧波潮流已不能滿足工程應(yīng)用的需求。因此，考慮系統(tǒng)的不確定性，進(jìn)行概率諧波潮流計(jì)算[6-7]，對(duì)電網(wǎng)諧波管理、薄弱環(huán)節(jié)分析具有重要意義。

目前，概率諧波潮流算法主要有解析法、模擬法和擬合法三種。解析法[8-9]指構(gòu)建待求隨機(jī)變量的概率密度解析式進(jìn)行求解，概率諧波潮流的求解實(shí)質(zhì)上是多個(gè)隨機(jī)變量的疊加，從理論上可以用卷積法進(jìn)行求解，但卷積涉及大量積分計(jì)算，當(dāng)隨機(jī)變量數(shù)量增多時(shí)，運(yùn)算將變得十分困難，因此不適用于大規(guī)模網(wǎng)絡(luò)的概率諧波潮流計(jì)算。模擬法[10-11]主要包含蒙特卡洛法，其思路是根據(jù)網(wǎng)絡(luò)注入量的概率密度函數(shù)進(jìn)行反復(fù)抽樣計(jì)算，得到多組諧波潮流的結(jié)果，最后進(jìn)行統(tǒng)計(jì)分析得到節(jié)點(diǎn)狀態(tài)量的概率密度。為了保證計(jì)算結(jié)果的精度，可能需要進(jìn)行上萬(wàn)甚至幾十萬(wàn)次的運(yùn)算，對(duì)運(yùn)算時(shí)間和內(nèi)存空間消耗巨大，且難以進(jìn)行靈敏度分析。但由于蒙特卡洛法具有編程簡(jiǎn)單、準(zhǔn)確度可靠等優(yōu)點(diǎn)，通常作為對(duì)照實(shí)驗(yàn)方法使用。擬合法包括點(diǎn)估計(jì)法[12]、仿射法[13]、級(jí)數(shù)擬合法[14]等，該方法不直接求出隨機(jī)變量的概率密度，而是通過(guò)求取隨機(jī)變量的數(shù)字特征，通過(guò)數(shù)值擬合方法，獲得概率密度函數(shù)，具有運(yùn)算量小、精度可靠等特點(diǎn)。其中仿射法常用于改善區(qū)間分析的保守性，能夠快速計(jì)算獲得諧波的分布區(qū)間，但無(wú)法分析區(qū)間內(nèi)的概率分布；點(diǎn)估計(jì)法通常需要構(gòu)造2m或2m+1個(gè)場(chǎng)景(m為隨機(jī)變量個(gè)數(shù))進(jìn)行計(jì)算，當(dāng)隨機(jī)變量數(shù)量增多時(shí)，計(jì)算量會(huì)明顯增大，有一定的局限性。

本文提出一種基于半不變量及最大熵的概率諧波潮流算法，根據(jù)諧波電流的監(jiān)測(cè)樣本數(shù)據(jù)計(jì)算電網(wǎng)諧波電壓的概率密度函數(shù)，屬于擬合法的范疇。在計(jì)算多個(gè)隨機(jī)諧波電流共同作用產(chǎn)生的諧波電壓概率密度函數(shù)時(shí)，為了避免卷積運(yùn)算，引入半不變量，半不變量能夠與樣本數(shù)據(jù)的高階矩互相轉(zhuǎn)換，易于求取，利用半不變量的可加性與齊次性簡(jiǎn)化概率求解過(guò)程，且計(jì)算量不會(huì)隨著隨機(jī)變量的增多而明顯增加。根據(jù)最大熵原理能夠在給定的約束條件下確定一種主觀成分最少、最接近實(shí)際情況的概率密度。以待求諧波的高階矩作為約束條件構(gòu)建最大熵模型，可以快速準(zhǔn)確地?cái)M合出諧波概率密度函數(shù)，進(jìn)而求取諧波95%概率值等特征參數(shù)。將本文算法在4節(jié)點(diǎn)系統(tǒng)上與解析法比較，以及在IEEE 57節(jié)點(diǎn)系統(tǒng)上與蒙特卡洛法比較，驗(yàn)證本文算法的有效性。

1　諧波潮流方程及數(shù)學(xué)基礎(chǔ)

1.1　諧波潮流方程

在諧波潮流計(jì)算中，通常將諧波源等值為恒流源模型或諾頓等效模型[15]，本文采用恒流源模型進(jìn)行等值。假設(shè)諧波潮流計(jì)算中的諧波次數(shù)為h，計(jì)算電網(wǎng)各元件諧波參數(shù)，生成諧波導(dǎo)納矩陣Yh。根據(jù)諧波注入電流情況和諧波導(dǎo)納矩陣，可以列出諧波潮流方程：

Ih=YhUh

(1)

式中：Ih為h次諧波注入電流；Yh為h次諧波導(dǎo)納矩陣；Uh為節(jié)點(diǎn)的h次諧波電壓。

當(dāng)諧波源注入電流和網(wǎng)絡(luò)參數(shù)已知時(shí)，可運(yùn)用稀疏技術(shù)對(duì)式(1)進(jìn)行求解，得到電網(wǎng)各節(jié)點(diǎn)諧波電壓。若對(duì)諧波導(dǎo)納矩陣求逆可得諧波阻抗矩陣，相應(yīng)地，式(1)可改寫(xiě)為：

(2)

式中：Zh為h次諧波阻抗矩陣。

使用計(jì)算機(jī)求解諧波潮流時(shí)，采用大規(guī)模電網(wǎng)的計(jì)算框架，諧波導(dǎo)納矩陣可根據(jù)電網(wǎng)拓?fù)浣Y(jié)構(gòu)和基波參數(shù)直觀地求得，采用稀疏矩陣技術(shù)儲(chǔ)存數(shù)據(jù)和求解諧波潮流方程，有利于節(jié)省內(nèi)存和提高計(jì)算速度，與基波潮流的成熟算法類(lèi)似，以上方法適用于各種規(guī)模的電網(wǎng)計(jì)算。

1.2　半不變量

當(dāng)考慮諧波源注入電流的隨機(jī)性時(shí)，式(2)則無(wú)法滿足計(jì)算要求。此時(shí)將諧波電流視為隨機(jī)變量進(jìn)行概率諧波潮流計(jì)算，并引入半不變量簡(jiǎn)化計(jì)算過(guò)程。

半不變量是隨機(jī)變量很重要的一種數(shù)字特征，又稱(chēng)累積量，半不變量序列可唯一確定隨機(jī)變量的分布規(guī)律[16]。根據(jù)定義，半不變量可通過(guò)對(duì)隨機(jī)變量的特征函數(shù)取自然對(duì)數(shù)，再進(jìn)行泰勒展開(kāi)求取而得[17]，計(jì)算過(guò)程復(fù)雜，在實(shí)際運(yùn)用中，通常將半不變量與高階原點(diǎn)矩進(jìn)行互相轉(zhuǎn)換。對(duì)隨機(jī)變量X，其半不變量與高階原點(diǎn)矩的關(guān)系如下：

(3)

式中：gk為k階半不變量；ak為k階原點(diǎn)矩。

利用隨機(jī)變量X的樣本數(shù)據(jù)可直接求得ak，有

(4)

式中：xi為隨機(jī)變量X的第i個(gè)可能取值；pi表示取值為xi時(shí)的概率。半不變量具有可加性和齊次性兩個(gè)重要性質(zhì)，在隨機(jī)分析的簡(jiǎn)化計(jì)算中起到關(guān)鍵作用。

可加性：若X(t)為n個(gè)互相獨(dú)立的隨機(jī)變量之和，即X(t)=X(1)+X(2)+…+X(n)，則X(t)的k階半不變量為：

(5)

齊次性：若隨機(jī)變量X(t)與隨機(jī)變量X(1)呈線性關(guān)系X(t)=aX(1)+b，則X(t)的k階半不變量為：

(6)

可以看出，半不變量與高階矩之間的轉(zhuǎn)換、隨機(jī)變量的半不變量之間的轉(zhuǎn)換均可輕易實(shí)現(xiàn)，因此半不變量在隨機(jī)變量分析中具有獨(dú)特優(yōu)勢(shì)，能夠大大簡(jiǎn)化運(yùn)算過(guò)程。

1.3　最大熵原理

本文采用最大熵原理，實(shí)現(xiàn)高階矩對(duì)概率密度函數(shù)的擬合。最大熵又稱(chēng)最大信息熵，其擬合隨機(jī)變量概率密度的中心思想是：在滿足給定的約束條件下，使得信息熵取最大值所對(duì)應(yīng)的概率密度分布是最接近實(shí)際情況的分布。自提出以來(lái)，最大信息熵原理已在通信、交通、氣象等領(lǐng)域獲得成功應(yīng)用，在電力系統(tǒng)中也越來(lái)越受到關(guān)注[18-20]。

假設(shè)X是一個(gè)離散型隨機(jī)變量，取值為xi時(shí)對(duì)應(yīng)的概率為pX(xi)，則X的信息熵為：

(7)

根據(jù)最大熵原理的思想，當(dāng)信息熵取得最大值時(shí)，得到的隨機(jī)變量概率密度函數(shù)最接近事實(shí)，含主觀誤差最小。以隨機(jī)變量X的k階原點(diǎn)矩ak和歸一化條件作為約束，對(duì)信息熵建立最大值規(guī)劃模型：

(8)

(9)

(10)

諧波電流、電壓屬于連續(xù)型隨機(jī)變量，但諧波監(jiān)測(cè)樣本數(shù)據(jù)是一組離散數(shù)列，因此對(duì)諧波概率密度函數(shù)的擬合也轉(zhuǎn)換為對(duì)離散值概率分布的擬合。

2　基于半不變量及最大熵的概率諧波潮流

2.1　諧波潮流方程線性化

省略諧波次數(shù)h，將式(2)的諧波潮流方程用相量的形式表示：

(11)

假設(shè)電網(wǎng)結(jié)構(gòu)不發(fā)生改變，即上式中Zij不變，諧波源諧波電流Ij為節(jié)點(diǎn)注入量，節(jié)點(diǎn)諧波電壓Ui是節(jié)點(diǎn)狀態(tài)量，可以寫(xiě)出方程：

U=f(I)

(12)

假設(shè)諧波電流的監(jiān)測(cè)樣本中電流幅值數(shù)據(jù)為{I}，選取幅值期望值作為諧波潮流基準(zhǔn)運(yùn)行點(diǎn)，即I0=E(I)，代入式(11)得到基準(zhǔn)運(yùn)行點(diǎn)節(jié)點(diǎn)諧波電壓為U0。

將式(12)在基準(zhǔn)運(yùn)行點(diǎn)處進(jìn)行泰勒展開(kāi)，并忽略2次及以上的高次項(xiàng)，可得

U=U0+ΔU=f(I0)+f′(I0)ΔI

(13)

式中：ΔI為諧波源諧波電流誤差量；ΔU為節(jié)點(diǎn)諧波電壓誤差量；f′(I0)為ΔU對(duì)ΔI的靈敏度系數(shù)。

將上式中的諧波電流電壓寫(xiě)成幅值和相角的形式，靈敏度系數(shù)用矩陣表示為：

(14)

式中：R的元素Rij=?Ui/?Ij;S的元素Sij=?Ui/?θIj;J的元素Jij=?θUi/?Ij;T的元素Tij=?θUi/?θIj。

上述參數(shù)均是基于諧波潮流基準(zhǔn)運(yùn)行點(diǎn)進(jìn)行的計(jì)算，結(jié)果為常數(shù)，可見(jiàn)式(14)是一組線性方程組。

2.2　隨機(jī)變量數(shù)字特征求取

采用線性化諧波方程時(shí)，諧波源諧波電流的幅值和相角與節(jié)點(diǎn)諧波電壓的幅值和相角可記為基準(zhǔn)量和誤差量之和，為了統(tǒng)一分析，將這4個(gè)隨機(jī)變量表示為：

X=X0+ΔX

(15)

式中：X為狀態(tài)量；X0為基準(zhǔn)量；ΔX為誤差量。

由于X0是一個(gè)常數(shù)列，由式(3)和式(4)可得k階原點(diǎn)矩和半不變量為：

(16)

(17)

根據(jù)半不變量的可加性，可得X與ΔX的k階半不變量之間的關(guān)系如下：

(18)

根據(jù)半不變量的可加性和齊次性，利用式(14)的線性關(guān)系，可計(jì)算得到諧波電壓幅值和相角誤差量的半不變量為：

(19)

式中：gΔUk為節(jié)點(diǎn)諧波電壓幅值誤差量的半不變量；為節(jié)點(diǎn)諧波電壓相角誤差量的半不變量;gΔIk為節(jié)點(diǎn)諧波電流幅值誤差量的半不變量；為節(jié)點(diǎn)諧波電流相角誤差量的半不變量。

得到諧波電壓半不變量后，根據(jù)式(3)轉(zhuǎn)化為原點(diǎn)矩，利用最大熵原理進(jìn)行擬合獲得概率分布。

2.3　算法流程

基于半不變量和最大熵的線性化概率諧波潮流算法步驟如下。

1)讀取電網(wǎng)基礎(chǔ)信息，確定待求諧波潮流諧波次數(shù)，列出諧波潮流方程。

2)讀取諧波源電流樣本數(shù)據(jù)，確定諧波潮流基準(zhǔn)運(yùn)行點(diǎn)，計(jì)算該基準(zhǔn)運(yùn)行點(diǎn)下的諧波潮流。

3)在基準(zhǔn)運(yùn)行點(diǎn)處，將諧波潮流方程線性化，計(jì)算諧波電壓誤差量對(duì)電流誤差量的靈敏度矩陣。

4)處理諧波源電流樣本數(shù)據(jù)，計(jì)算諧波電流幅值誤差量的半不變量和相角誤差量的半不變量。

5)利用半不變量的可加性和齊次性，根據(jù)式(19)計(jì)算待求節(jié)點(diǎn)諧波電壓幅值誤差量的半不變量和相角誤差量的半不變量。

6)與基準(zhǔn)值進(jìn)行疊加，根據(jù)式(18)計(jì)算諧波電壓幅值半不變量和相角半不變量。

7)計(jì)算諧波電壓幅值原點(diǎn)矩和相角原點(diǎn)矩。

8)建立最大熵模型，使用合適的數(shù)值分析方法求解，擬合出諧波電壓的概率分布。

3　算例分析

3.1　4節(jié)點(diǎn)系統(tǒng)算例

為了說(shuō)明本文方法的適用性，對(duì)一個(gè)簡(jiǎn)單4節(jié)點(diǎn)系統(tǒng)進(jìn)行分析，從而便于采用卷積算法計(jì)算準(zhǔn)確的概率分布曲線。該系統(tǒng)包含4個(gè)節(jié)點(diǎn)、3條線路、2臺(tái)發(fā)電機(jī)和1臺(tái)變壓器，網(wǎng)絡(luò)結(jié)構(gòu)見(jiàn)附錄A圖A1?；ㄏ碌木W(wǎng)絡(luò)參數(shù)可查閱文獻(xiàn)[21]。

節(jié)點(diǎn)1和3均帶有非線性負(fù)載，兩者互相獨(dú)立。為了使卷積計(jì)算簡(jiǎn)單，只研究幅值的概率密度。假設(shè)節(jié)點(diǎn)1的5次諧波電流幅值服從正態(tài)分布I1～N(0.08,0.022)，節(jié)點(diǎn)3的5次諧波電流幅值服從正態(tài)分布I3～N(0.06,0.012)。計(jì)算網(wǎng)絡(luò)中各節(jié)點(diǎn)5次諧波電壓的概率密度。

首先，采用卷積法進(jìn)行求解，獲得準(zhǔn)確的概率曲線。由于只研究幅值的概率密度，諧波方程中的復(fù)數(shù)運(yùn)算變?yōu)槟Ｖ颠\(yùn)算，根據(jù)正態(tài)分布的卷積特性，可直接計(jì)算各節(jié)點(diǎn)諧波電壓的概率密度函數(shù)。

利用本文方法計(jì)算節(jié)點(diǎn)諧波電壓的概率密度，與卷積法進(jìn)行對(duì)比，兩種方法得到的概率密度曲線基本一致，具體可參考附錄A圖A2，證明了本文方法的有效性。當(dāng)網(wǎng)絡(luò)規(guī)模以及諧波源數(shù)量進(jìn)一步增大時(shí)，使用直接法時(shí)卷積計(jì)算量和復(fù)雜程度大大增加，不利于實(shí)際應(yīng)用，此時(shí)運(yùn)用本文的方法將具有明顯優(yōu)勢(shì)。

3.2　IEEE 57節(jié)點(diǎn)系統(tǒng)算例

采用IEEE 57節(jié)點(diǎn)系統(tǒng)作為研究對(duì)象。該系統(tǒng)包含57個(gè)節(jié)點(diǎn)、7臺(tái)發(fā)電機(jī)、63條輸電線路以及17臺(tái)變壓器，拓?fù)浣Y(jié)構(gòu)如圖1所示。節(jié)點(diǎn)19，20，26，45，54，55帶有非線性負(fù)載，產(chǎn)生并向電網(wǎng)中注入5次諧波。電網(wǎng)中發(fā)電機(jī)、變壓器等設(shè)備及其他節(jié)點(diǎn)產(chǎn)生的諧波忽略不計(jì)。

圖1IEEE 57節(jié)點(diǎn)電力系統(tǒng)
Fig.1IEEE 57-bus power system

以上6個(gè)諧波注入點(diǎn)的諧波樣本取自廣東某電網(wǎng)諧波監(jiān)測(cè)點(diǎn)的監(jiān)測(cè)數(shù)據(jù)。每個(gè)節(jié)點(diǎn)的諧波樣本包含約16 000個(gè)5次諧波電流幅值數(shù)據(jù)，其概率密度詳見(jiàn)附錄B圖B1。

運(yùn)用蒙特卡洛法進(jìn)行104次重復(fù)計(jì)算，結(jié)果與本文方法進(jìn)行對(duì)比。由于該電力系統(tǒng)中節(jié)點(diǎn)較多，受篇幅所限，每個(gè)節(jié)點(diǎn)的計(jì)算結(jié)果無(wú)法一一列出，本文隨機(jī)選取節(jié)點(diǎn)4，10，15，25，38，53進(jìn)行對(duì)比分析。

節(jié)點(diǎn)5次諧波電壓概率密度計(jì)算結(jié)果對(duì)比如圖2所示。從圖中可以看出，本文方法與蒙特卡洛法得到的節(jié)點(diǎn)諧波電壓概率密度曲線基本吻合，從而證明本文方法能夠有效評(píng)估諧波電壓的概率密度。

基于本文方法和蒙特卡洛法評(píng)估的諧波電壓幅值統(tǒng)計(jì)數(shù)據(jù)如表1所示，兩種方法得到的諧波電壓期望值和標(biāo)準(zhǔn)差相差較小，具有近似的平均值和離散度。

圖2　諧波電壓幅值概率密度曲線Fig.2　Probability density curves of harmonic voltage amplitude

表1諧波電壓幅值統(tǒng)計(jì)數(shù)據(jù)
Table1Statisticaldataofharmonicvoltageamplitude

為了進(jìn)一步評(píng)估本文方法的準(zhǔn)確性，以蒙特卡洛法的結(jié)果作為對(duì)照組，計(jì)算節(jié)點(diǎn)諧波電壓概率密度的平均誤差和最大誤差，用來(lái)評(píng)估概率密度曲線吻合程度，如表2所示，除了節(jié)點(diǎn)10以外，其他節(jié)點(diǎn)的概率密度評(píng)估值誤差都很小，與蒙特卡洛法評(píng)估結(jié)果有較好的吻合程度。

表2　諧波電壓概率密度誤差分析Table 2　Analysis on probability density error of harmonic voltage

表3展示了兩種方法計(jì)算而得的節(jié)點(diǎn)諧波電壓幅值的95%概率值，本文方法得到的結(jié)果與蒙特卡洛法的結(jié)果基本相同，誤差在5%以內(nèi)。

表3　諧波電壓幅值95%概率值Table 3　95% probability values of harmonic voltage amplitude

蒙特卡洛法因?yàn)楹?jiǎn)單準(zhǔn)確的優(yōu)點(diǎn)被廣泛應(yīng)用于隨機(jī)分析方法的有效性驗(yàn)證中，但其計(jì)算量巨大，且隨著網(wǎng)絡(luò)規(guī)模的擴(kuò)大和隨機(jī)變量的增加，計(jì)算量呈指數(shù)性增長(zhǎng)。對(duì)比蒙特卡洛法和本文方法在同一平臺(tái)上的計(jì)算耗時(shí)，排除網(wǎng)絡(luò)拓?fù)浜蛥?shù)讀入、諧波導(dǎo)納矩陣生成所需的時(shí)間，蒙特卡洛法耗時(shí)138.28 s，而本文方法耗時(shí)1.53 s，在計(jì)算速度方面具有明顯的優(yōu)勢(shì)。

綜上所述，本文方法得到的結(jié)果具有較高的準(zhǔn)確度，能夠合理地評(píng)估節(jié)點(diǎn)諧波電壓的概率分布。評(píng)估得到的諧波電壓期望值、標(biāo)準(zhǔn)差、95%概率值均與蒙特卡洛法的計(jì)算結(jié)果基本相同，證明本文方法能對(duì)節(jié)點(diǎn)諧波的分布預(yù)期、離散度做出較準(zhǔn)確的判斷，對(duì)分布規(guī)律的評(píng)估可能會(huì)出現(xiàn)局部誤差，原因在于諧波方程線性化的線性誤差、最大熵?cái)M合的計(jì)算誤差以及蒙特卡洛法本身的誤差。本文方法較其他評(píng)估方法的優(yōu)點(diǎn)主要在于能夠避免大量重復(fù)模擬運(yùn)算，運(yùn)算時(shí)間和占用內(nèi)存空間少，對(duì)于電網(wǎng)結(jié)構(gòu)或諧波源的變化，僅需修改相應(yīng)參數(shù)即可重新計(jì)算，且計(jì)算過(guò)程中無(wú)需對(duì)隨機(jī)變量的分布規(guī)律做任何假設(shè)，能夠避免主觀判斷誤差，具有客觀準(zhǔn)確性。

4　結(jié)語(yǔ)

由于電力系統(tǒng)中的諧波分布具有明顯的隨機(jī)性，確定性諧波潮流計(jì)算無(wú)法滿足實(shí)際應(yīng)用的需求，本文結(jié)合數(shù)據(jù)半不變量、線性化諧波潮流方程以及最大熵規(guī)劃模型，提出一種概率諧波計(jì)算方法。該方法根據(jù)諧波源電流的概率分布，計(jì)算并擬合出電網(wǎng)各個(gè)節(jié)點(diǎn)的諧波電壓的概率分布曲線，進(jìn)而能夠?yàn)殡娋W(wǎng)薄弱環(huán)節(jié)和潛在風(fēng)險(xiǎn)的評(píng)估、諧波管理等工作提供指導(dǎo)意見(jiàn)，具有重要的實(shí)際應(yīng)用價(jià)值。本文方法能夠避免復(fù)雜的卷積運(yùn)算和大量的蒙特卡洛模擬計(jì)算，具有較高的準(zhǔn)確性。如何減少線性化和曲線擬合等環(huán)節(jié)的誤差，仍需進(jìn)一步研究。

附錄見(jiàn)本刊網(wǎng)絡(luò)版(http://www.aeps-info.com/aeps/ch/index.aspx)。