曾　江， 蔡東陽， 黃德華

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

0　引言

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

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

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

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

1.1　諧波潮流方程

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

Ih=YhUh

(1)

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

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

(2)

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

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

1.2　半不變量

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

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

(3)

式中：gk為k階半不變量；ak為k階原點矩。

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

(4)

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

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

(5)

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

(6)

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

1.3　最大熵原理

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

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

(7)

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

(8)

(9)

(10)

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

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

2.1　諧波潮流方程線性化

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

(11)

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

U=f(I)

(12)

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

將式(12)在基準(zhǔn)運行點處進(jìn)行泰勒展開，并忽略2次及以上的高次項，可得

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

(13)

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

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

(14)

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

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

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

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

X=X0+ΔX

(15)

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

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

(16)

(17)

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

(18)

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

(19)

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

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

2.3　算法流程

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

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

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

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

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

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

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

7)計算諧波電壓幅值原點矩和相角原點矩。

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

3　算例分析

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4　結(jié)語

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

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