江友華，劉子瑜，張 煜，楊興武，吳衛(wèi)民

（1.上海電力大學，上海 200120；2.國網(wǎng)山東省電力公司棗莊供電公司，山東 棗莊 277000；3.上海海事大學，上海 200120）

0 引言

隨著智能電網(wǎng)、主動配電網(wǎng)等概念的提出和逐步實現(xiàn)，我國電網(wǎng)架構(gòu)和運行方式都較以往發(fā)生了巨大變化，大量新能源和電力電子設(shè)備的接入，使電網(wǎng)的電能質(zhì)量問題，尤其是諧波污染問題日益嚴重[1-4]。針對諧波污染的獎懲機制建立問題，國內(nèi)外已經(jīng)進行了廣泛而深刻的探討，而獎懲機制建立的前提是諧波責任的合理區(qū)分[5]。

諧波責任區(qū)分的關(guān)鍵問題是諧波阻抗的有效估計。目前，諧波阻抗的估計大多基于線性回歸法[6-10]。一般的線性回歸法在背景諧波電壓變化的情況下，有較大誤差。針對這一問題已進行了很多探討：文獻[11]利用聚類的方法，篩選背景諧波電壓穩(wěn)定的有效數(shù)據(jù)段，以期降低背景諧波電壓波動的影響，但系統(tǒng)中存在未知諧波源時，聚類效果變差；文獻[12]將傳統(tǒng)線性回歸模型回歸常數(shù)項修改為一個關(guān)于時間的函數(shù)，利用泰勒公式近似后，采用核估計的方法找到其表達式，據(jù)此估計諧波阻抗，該方法僅在背景諧波電壓波動不大的情況下有較理想的效果；文獻[13]采用波形匹配的方法篩選背景諧波電壓相對穩(wěn)定的數(shù)據(jù)，對諧波監(jiān)測數(shù)據(jù)進行分類處理，分別進行回歸分析。

由以上分析可知：現(xiàn)有的研究往往需要對數(shù)據(jù)做人為處理或假設(shè)，可能引入誤差；對諧波監(jiān)測數(shù)據(jù)分類標準不同，方法各異，實際應(yīng)用時常有局限性，普適性并不理想。針對以上問題，提出一種普適性諧波阻抗估計方法。考慮背景諧波電壓難以直接測量的特點，在一般線性模型的基礎(chǔ)上將背景諧波電壓視為隱變量，用GMM（高斯混合模型）建模，并指出GMM 參數(shù)在實際工程背景下的意義。將GMM 參數(shù)、線性模型參數(shù)建模為DPM（狄利克雷過程混合模型），推導(dǎo)基于諧波監(jiān)測數(shù)據(jù)的各參數(shù)的后驗分布，利用MCMC（馬爾科夫鏈-蒙特卡洛）采樣方法從各后驗分布中抽取樣本，進行貝葉斯估計，求解諧波阻抗和背景諧波電壓工況數(shù)，并評估諧波責任。

1 含隱變量的多諧波源系統(tǒng)模型

1.1 含未知諧波源的多諧波源系統(tǒng)

近年來隨著光伏、風力機等新能源的并網(wǎng)，電力電子裝置的大量應(yīng)用，系統(tǒng)中單一諧波源的諧波發(fā)射水平可能并不顯著，但大量此類諧波源的共同作用不能忽略。顯然，要逐一監(jiān)測這些諧波源是不現(xiàn)實的，多諧波源系統(tǒng)中存在未知諧波源。本文將多諧波源系統(tǒng)中的諧波源劃分為關(guān)注諧波源和背景諧波源，背景諧波源進一步劃分為可監(jiān)測諧波源和未知諧波源，其關(guān)系如圖1 所示。

這種含有未知諧波源的多諧波源系統(tǒng)的等效電路如圖2 所示，其中，電流和阻抗Z 的下標表示相應(yīng)的諧波源，c1—ck為關(guān)注諧波源，ck+1—cq為可監(jiān)測諧波源，cq+1—cn為未知諧波源。

圖1 諧波源的分類

圖2 含未知諧波源的多諧波源系統(tǒng)

對于某次諧波，PCC（公共耦合點）處的諧波電壓可表示為：

式中：R 為諧波責任；I 為諧波電流幅值；Upcc為PCC 處諧波電壓幅值；為諧波阻抗的模值；θ為諧波源在PCC 處產(chǎn)生的諧波電壓與的夾角。本文假定系統(tǒng)中各諧波源不存在耦合或已經(jīng)過解耦處理，將各饋線電流視為諧波源的發(fā)出電流。按以下原則評估諧波責任，若R 為正，則對應(yīng)諧波源發(fā)出諧波，應(yīng)按比例承擔相應(yīng)責任；若R 為負，則對應(yīng)濾波源吸收諧波，不應(yīng)受到懲罰或應(yīng)得獎勵。

實際工程中，諧波電壓和諧波電流的相角通常難以準確測量，參考文獻[13]的方法，以為基準，將其旋轉(zhuǎn)至正實軸，其余相量旋轉(zhuǎn)同樣的角度，此時取式（1）的實部則有：

1.2 背景諧波電壓GMM 及其參數(shù)的意義

背景諧波電壓常常包含未知諧波源，因此難以直接測量背景諧波電壓數(shù)據(jù)，故將背景諧波電壓視為隱變量。用線性回歸法估計諧波阻抗時，式（4）可寫作式（5）：

式中：Upcc，i為PCC 處諧波電壓，即回歸模型的響應(yīng)變量（i 表示第i 個數(shù)據(jù)，下同）；Ii為諧波源的諧波電流；Uo，i為背景諧波電壓，是隱變量；εi為隨機誤差，εi～N（0，σε2），服從均值為0、方差為σε2的高斯分布，滿足協(xié)方差Cor（εi，εj）=0，i≠j，均值E（εiIi）=0，E（εiUo，i）=0。

背景諧波電壓的分布往往為非標準高斯分布，呈現(xiàn)峰值偏移、多峰、長尾等特征[15]，僅采用核估計法難以準確建模背景諧波電壓分布。高斯混合分布常用于近似表示復(fù)雜的非標準分布，故將背景諧波電壓建模為GMM，即Uo，i～N（μi，σi2），且εi～N（0，σε2），由于高斯分布的線性性質(zhì)，Uo，i+εi服從高斯混合分布，則Upcc，i-ZIi服從高斯混合分布，即：

式中：f（g）為高斯分布的密度函數(shù)，其參數(shù)為Θc；πc為權(quán)重參數(shù)，表示數(shù)據(jù)來自第c 個高斯分布的概率；混合模型有K 個不同的高斯分布，K≤N，N 為數(shù)據(jù)容量。對于每個數(shù)據(jù)i，其對應(yīng)參數(shù)Θc，i由μc，i和組成。則有：

式（7）即為背景諧波電壓的GMM。

在實際的工程背景下，背景諧波電壓GMM參數(shù)K，c，πc和Θc都具有實際意義：參數(shù)K 在GMM 中代表模型中混合的單個高斯分布的個數(shù)，在實際背景中代表背景諧波電壓工況的個數(shù)，若諧波監(jiān)測數(shù)據(jù)屬于同一種背景諧波電壓工況，則其背景諧波電壓服從同一個高斯分布，在一定范圍內(nèi)保持穩(wěn)定；參數(shù)c 在GMM 中表示數(shù)據(jù)來自第c 個高斯分布，在實際中表示諧波監(jiān)測數(shù)據(jù)屬于第c 個背景諧波電壓工況；參數(shù)πc在GMM 中代表數(shù)據(jù)來自第c 個高斯分布的概率，在實際中表示諧波監(jiān)測數(shù)據(jù)屬于第c 個背景諧波電壓工況的概率；參數(shù)Θc在GMM 中是第c 個高斯分布的參數(shù)，在實際中是第c 個背景諧波電壓工況對應(yīng)分布的參數(shù)。

1.3 基于狄利克雷過程混合模型的參數(shù)建模

在含未知諧波源的背景諧波電壓GMM 中，混合分布的元素個數(shù)K 難以直接確定?；诓恢付ū尘爸C波電壓服從分布類型的前提，即在“更弱、更一般”的情況下，對參數(shù)Θ 引入一個先驗G，G 是一個狄利克雷過程混合分布。Θ 的共軛先驗，即G 的基分布G0如下：

完整的DPM 可以表示為：

式（10）—（14）中：Discrete（g）為離散分布；INVGamma（g）為逆Gamma 分布，a，b 為其參數(shù)；Dir（g）為狄利克雷分布，α 為其參數(shù)。

基于DPM，可通過從分布中抽樣，用參數(shù)估計的方式確定K 及其他參數(shù)的值。對DPM 抽樣需要已知參數(shù)的后驗分布由貝葉斯定理，各參數(shù)的后驗分布可統(tǒng)一表示為：

2 基于MCMC 采樣的諧波責任區(qū)分

2.1 模型參數(shù)的后驗分布

諧波阻抗的計算和諧波責任區(qū)分應(yīng)建立在對前文模型進行非參數(shù)貝葉斯估計的基礎(chǔ)上，需要估計的參數(shù)由狄利克雷過程參數(shù){c，Θ，α}和線性模型參數(shù){Z，Uo}兩部分組成。

在貝葉斯理論中，參數(shù)估計需要基于參數(shù)的后驗分布進行。目前DPM 的后驗分布已有成熟的研究成果[16]，由于其推導(dǎo)較為復(fù)雜繁瑣，本文將根據(jù)各參數(shù)的先驗分布，在已知諧波監(jiān)測數(shù)據(jù)的前提下直接給出各參數(shù)的后驗分布。

2.1.1 參數(shù)c

參數(shù)c 和參數(shù)Θ 存在一一對應(yīng)的關(guān)系，即ci對應(yīng)Θi，若ci=cj，則Θi=Θj。首先定義：

其中#（g）為滿足條件的元素個數(shù)。由狄利克雷過程定義，ci的先驗為：

定義F（Upcc，i，Θc）為似然函數(shù)，b 為歸一化參數(shù)，則ci的后驗為：

2.1.2 參數(shù)Θ

背景諧波電壓工況對應(yīng)分布的參數(shù)Θ 由μ和σ2組成。其先驗已經(jīng)在前文給出，如式（8）、式（9）所示。

參數(shù)Θ 的后驗分布為：

2.1.3 參數(shù)Z

線性模型參數(shù)Z，即諧波阻抗的先驗分布可以通過一般的線性回歸分析進行粗估，將粗估結(jié)果設(shè)為諧波阻抗的先驗分布：

在獲得諧波監(jiān)測數(shù)據(jù)后，可以得到諧波阻抗的后驗分布：

2.1.4 參數(shù)Uo

Uo的后驗分布為：

2.1.5 參數(shù)α

參數(shù)α 為狄利克雷過程的集中參數(shù)，在實際模型中是背景諧波電壓工況的個數(shù)的初始值。由式（14）和式（15）可推導(dǎo)α 的后驗分布為：

式中：Γ（g）為Gamma 函數(shù)；α′為α 所有可能的取值。

式（26）即為參數(shù)α 的后驗分布。

2.2 MCMC 采樣和諧波責任區(qū)分

在推導(dǎo)出各模型參數(shù)的后驗分布后，即可使用MCMC 采樣方法對模型各參數(shù)進行抽樣。在各模型參數(shù)抽樣完畢后，依據(jù)線性模型參數(shù)Z，即諧波阻抗的抽樣結(jié)果，計算其均值作為諧波阻抗的估計值。進而將諧波阻抗Z，即cosθ 帶入諧波責任定義式（2）中，可求出關(guān)注諧波源在PCC 處的諧波責任。MCMC 采樣及諧波責任評估流程如圖3 所示。

圖3 MCMC 采樣及諧波責任評估流程

3 案例分析

為驗證文章方法的有效性，利用Simulink 平臺，對IEEE 14 節(jié)點測試系統(tǒng)和實際案例中測得的諧波發(fā)射水平數(shù)據(jù)進行測試。IEEE 14 節(jié)點測試系統(tǒng)網(wǎng)絡(luò)結(jié)構(gòu)如圖4 所示。

圖4 IEEE 14 節(jié)點系統(tǒng)

將節(jié)點4 視為PCC，在節(jié)點2、節(jié)點9、節(jié)點12 處分別接入諧波源HS1，HS2，HS3，作為關(guān)注諧波源，諧波源諧波發(fā)射水平分別參照文獻[17]中案例22、案例23、案例25 中5 次諧波電流數(shù)據(jù)設(shè)置；在節(jié)點10 處接入HS4 作為濾波源，按案例29[17]的濾波設(shè)備數(shù)據(jù)設(shè)置。為模擬圖2 中含有未知諧波源的多諧波源系統(tǒng)，在節(jié)點11 接入諧波源HSB，作為背景諧波源。為模擬背景諧波電壓的不同工況，將案例21[17]中5 次諧波發(fā)射水平的100%，80%，50%，30%視為4 種不同工況，并加上方差為0.5 的高斯隨機波動，將4 種工況的背景諧波電壓數(shù)據(jù)等數(shù)量混合作為背景諧波電壓。采集3 600 個數(shù)據(jù)點，包括諧波源HS1，HS2，HS3，HS4 的諧波電流和PCC 處的諧波電壓，組成諧波測量數(shù)據(jù)集。

根據(jù)采集到的諧波數(shù)據(jù)，使用MCMC 采樣方法對模型參數(shù)進行抽樣，依據(jù)抽樣結(jié)果估計諧波阻抗。進行MCMC 抽樣前，需要設(shè)置模型參數(shù)的初值，目前還沒有一般性的初值設(shè)置方法。對于諧波阻抗初值，先采用一般線性回歸法進行粗估，將粗估結(jié)果作為MCMC 采樣的初值，其他參數(shù)給出表1 所示的初值參考。

表1 參數(shù)初值參考值

文獻[18]指出，一條MCMC 采樣鏈需要經(jīng)過多次迭代后才能趨于穩(wěn)定，其抽樣結(jié)果才具有用于貝葉斯推斷的意義。因此，對于不同的MCMC采樣鏈，需設(shè)置不同的burn-in 值來剔除MCMC采樣穩(wěn)定前的數(shù)據(jù)，burn-in 值需要通過具體測試才能確定。通過測試，1 500 次迭代后，MCMC采樣鏈會趨于穩(wěn)定，因此將burn-in 值設(shè)定為1 500，共迭代2 000 次，取最后500 次作為推斷依據(jù)數(shù)據(jù)。

采用非參數(shù)貝葉斯估計法和一般線性回歸法分別重復(fù)進行20 次諧波阻抗估計試驗。以理論計算諧波阻抗為真實值，將非參數(shù)貝葉斯估計法和一般線性回歸法估計諧波源HS1，HS2，HS3諧波阻抗的相對誤差進行對比，結(jié)果如圖5—7所示。

圖5 HS1 諧波阻抗估計誤差

圖6 HS2 諧波阻抗估計誤差

圖7 HS3 諧波阻抗估計誤差

由圖5 可知，對于諧波源HS1 的諧波阻抗估計，非參數(shù)貝葉斯估計法的估計誤差大部分在5%以內(nèi)，少數(shù)誤差較大的結(jié)果不超過10%，一般線性回歸法的估計誤差則在16%～25%。由圖6 可知，對于諧波源HS2 的諧波阻抗估計，非參數(shù)貝葉斯估計法的估計誤差在6%以內(nèi)，一般線性回歸法的估計誤差大部分在6%～12%。由圖7 可知，對于諧波源HS3 的諧波阻抗估計，非參數(shù)貝葉斯估計法的估計誤差在6%以內(nèi)，一般線性回歸法的估計誤差大部分在6%～13%?？梢钥闯?，非參數(shù)貝葉斯估計法較一般線性回歸法在估計諧波阻抗的精度方面更具優(yōu)勢。

取20 次重復(fù)試驗計算所得諧波責任的平均值作為最終結(jié)果。以理論計算諧波責任為真實值，將非參數(shù)貝葉斯估計法和一般線性回歸法估計諧波源HS1，HS2，HS3，HS4 的諧波責任進行對比，結(jié)果見表2。

由表2 可知，在評估諧波責任方面，非參數(shù)貝葉斯估計法比一般線性回歸法精度更高。此外，從諧波責任真實值可以看出，諧波源HS1，HS2，HS3，HS4 的諧波責任總和沒有達到100%，其余部分由背景諧波電壓貢獻，非參數(shù)貝葉斯估計法的評估結(jié)果較一般線性回歸法更能反映這一事實。

表2 諧波責任評估數(shù)據(jù)

傳統(tǒng)線性回歸法在處理背景諧波電壓波動的情形時，需要借助聚類算法分別計算各段背景諧波電壓穩(wěn)定時的諧波責任，然后再進行綜合評估，聚類簇數(shù)K 的確定常常需要人為指定。非參數(shù)貝葉斯估計法在估計諧波阻抗的同時，可以通過迭代給出背景諧波電壓工況的個數(shù)K，即線性回歸法的聚類簇數(shù)，而不需要人為指定。重復(fù)20次非參數(shù)貝葉斯估計法試驗，所得的背景諧波電壓諧波工況數(shù)K 的頻數(shù)分布如圖8 所示。

圖8 背景諧波電壓工況數(shù)頻數(shù)分布

可以看出，K=4 的頻數(shù)為10，K=3 和K=5 的頻數(shù)分別為4 和3，非參數(shù)貝葉斯估計法對諧波工況個數(shù)的估計較為準確。

4 結(jié)語

線性回歸法在背景諧波電壓波動時，對諧波阻抗的估計常常有較大誤差，從而導(dǎo)致諧波責任評估失準；運用聚類-線性回歸法估計諧波阻抗時，在聚類簇個數(shù)的設(shè)定上又存在主觀性。本文針對以上兩點提出使用非參數(shù)貝葉斯估計法估計諧波阻抗，評估諧波責任，較線性回歸法有以下優(yōu)勢：

（1）非參數(shù)貝葉斯估計法估計諧波阻抗的精度高于一般線性回歸法。一般線性回歸法在背景諧波電壓數(shù)值波動較大時，誤差較大。實驗表明非參數(shù)貝葉斯估計法在非單一背景諧波電壓工況的情況下，仍能保持較為理想的評估精度。

（2）非參數(shù)貝葉斯估計法在“更弱、更一般”的假設(shè)前提下，直接利用諧波監(jiān)測數(shù)據(jù)估計諧波阻抗，評估諧波責任，無需人為干預(yù)，無需人工確定聚類簇數(shù)，是一種更具普適性的方法。非參數(shù)貝葉斯估計法還可以給出背景諧波電壓工況數(shù)，為背景諧波電壓的分析與研究提供了一種可行的方法。

此外，與傳統(tǒng)方法相比，非參數(shù)貝葉斯估計法含有多次迭代步驟，在運行速度方面并未表現(xiàn)出明顯優(yōu)勢，因此在對運行速度要求較高的場合，非參數(shù)貝葉斯估計法仍需改進或結(jié)合其他方法使用，以取得理想效果。