曾　江， 蔡東陽， 黃德華

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

0　引言

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

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

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

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

1.1　諧波潮流方程

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

Ih=YhUh

(1)

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

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

(2)

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

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

1.2　半不變量

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

半不變量是隨機變量很重要的一種數(shù)字特征，又稱累積量，半不變量序列可唯一確定隨機變量的分布規(guī)律[16]。根據(jù)定義，半不變量可通過對隨機變量的特征函數(shù)取自然對數(shù)，再進行泰勒展開求取而得[17]，計算過程復雜，在實際運用中，通常將半不變量與高階原點矩進行互相轉(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ù)的擬合。最大熵又稱最大信息熵，其擬合隨機變量概率密度的中心思想是：在滿足給定的約束條件下，使得信息熵取最大值所對應的概率密度分布是最接近實際情況的分布。自提出以來，最大信息熵原理已在通信、交通、氣象等領域獲得成功應用，在電力系統(tǒng)中也越來越受到關(guān)注[18-20]。

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

(7)

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

(8)

(9)

(10)

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

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

2.1　諧波潮流方程線性化

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

(11)

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

U=f(I)

(12)

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

將式(12)在基準運行點處進行泰勒展開，并忽略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ù)均是基于諧波潮流基準運行點進行的計算，結(jié)果為常數(shù)，可見式(14)是一組線性方程組。

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

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

X=X0+ΔX

(15)

式中：X為狀態(tài)量；X0為基準量；Δ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)化為原點矩，利用最大熵原理進行擬合獲得概率分布。

2.3　算法流程

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

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

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

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

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

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

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

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

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

3　算例分析

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

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

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

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

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

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

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

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

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

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

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

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

為了進一步評估本文方法的準確性，以蒙特卡洛法的結(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

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

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

4　結(jié)語

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

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