夏焰坤 ，唐文張 ，林欣懿

（1.西華大學(xué)電氣與電子信息學(xué)院，成都 610039；2.西南交通大學(xué)牽引動(dòng)力國家重點(diǎn)實(shí)驗(yàn)室，成都 610031）

由于當(dāng)前電力系統(tǒng)中存在電力電子設(shè)備的投切、高電壓直流輸電換流裝置的使用和電力機(jī)車的啟停等，導(dǎo)致電網(wǎng)中諧波源數(shù)量飛速增長(zhǎng)，諧波污染不可避免，電網(wǎng)運(yùn)行的穩(wěn)定性亦會(huì)因此降低[1-2]。為有效抑制諧波污染，國際上提出了獎(jiǎng)懲性方案，通過評(píng)判諧波發(fā)射水平并采用經(jīng)濟(jì)手段，對(duì)諧波負(fù)荷來源進(jìn)行懲罰或?qū)χC波受害者進(jìn)行補(bǔ)償[3]。

諧波責(zé)任分?jǐn)傄恢笔请娔苜|(zhì)量分析的重點(diǎn)問題[4-5]。諧波功率方向法和阻抗參數(shù)法是應(yīng)用于諧波責(zé)任劃分中最廣泛的兩大類方法。其中，諧波功率方向法又包括諧波有功功率方向法[6]、諧波無功功率方向法[7-8]及諧波視在功率識(shí)別法[9]等。這些方法均是以功率方向?yàn)榛A(chǔ)尋找諧波源的方法，因此諧波源的相位差對(duì)結(jié)果影響很大，特別是在大量新能源并網(wǎng)的今天，并不能保證此類方法始終適用。

阻抗參數(shù)法是當(dāng)今國內(nèi)外諧波水平研究中的熱門方向，其重點(diǎn)是如何有效計(jì)算出系統(tǒng)側(cè)諧波阻抗[10]。該方法由干預(yù)法和非干預(yù)法兩個(gè)大類組成。其中，干預(yù)法主要采用添加支路的方式對(duì)電路強(qiáng)加諧波電流，通過響應(yīng)的變化計(jì)算系統(tǒng)諧波阻抗。這種主動(dòng)注入諧波電流的方式在諧波阻抗計(jì)算的實(shí)驗(yàn)搭建和結(jié)果準(zhǔn)確性上具有顯著優(yōu)勢(shì)，但是此類方法在大系統(tǒng)中對(duì)系統(tǒng)運(yùn)行方式和穩(wěn)定性的影響也不可忽視，值得研究人員持續(xù)關(guān)注和進(jìn)一步研究。非干預(yù)法在近年來的研究過程中涌現(xiàn)出了一大批方法，例如波動(dòng)量法[11]、回歸法[12]、隨機(jī)獨(dú)立矢量協(xié)方差法[13]、盲源分離法[14]、支持向量機(jī)法[15]等。這些方法均在一定程度上解決了相應(yīng)的實(shí)際問題，例如文獻(xiàn)[16]提出了二元線性回歸法，將觀測(cè)值拆分為實(shí)部和虛部，通過線性回歸獲得系統(tǒng)諧波阻抗，原理簡(jiǎn)單，解決了諧波阻抗求解難的問題；文獻(xiàn)[17]提出了改進(jìn)的隨機(jī)獨(dú)立矢量協(xié)方差法，解決了觀測(cè)數(shù)據(jù)中相位缺失或者相位不能直接測(cè)得的問題；文獻(xiàn)[18]提出了支持向量機(jī)法，解決了觀測(cè)值樣本少、維數(shù)高、非線性及具有局部極小點(diǎn)等問題。背景諧波波動(dòng)性較大、系統(tǒng)側(cè)諧波與用戶側(cè)諧波具有關(guān)聯(lián)性等問題仍是諧波阻抗估計(jì)中的難點(diǎn)，故很難評(píng)判某些方法的泛化性能，需更進(jìn)一步研究。

為解決上述問題，本文提出一種優(yōu)化高斯過程回歸的諧波阻抗計(jì)算方法。高斯過程回歸[19-20]已普遍應(yīng)用于負(fù)荷預(yù)測(cè)、時(shí)間序列和計(jì)算機(jī)應(yīng)用等領(lǐng)域，受到大量專家學(xué)者的青睞。貝葉斯優(yōu)化算法[21]是一種全局優(yōu)化算法，作為一種優(yōu)秀的超參數(shù)調(diào)優(yōu)方式，將其運(yùn)用于高斯過程回歸中尋找最優(yōu)超參數(shù)。本文以公共連接點(diǎn)PCC（point of common coupling）處諧波電壓和諧波電流為基礎(chǔ)，通過貝葉斯優(yōu)化高斯過程回歸BO-GPR（Bayesian optimized Gaussian process regression）獲得系統(tǒng)側(cè)諧波阻抗，通過仿真分析與實(shí)例分析驗(yàn)證了該方法的有效性、準(zhǔn)確性及泛化性能。

1 基本原理

1.1 諧波阻抗估計(jì)原理

戴維南等效電路和諾頓等效電路是電力系統(tǒng)諧波分析中兩種常用的等效電路，如圖1所示。其中，Us為系統(tǒng)側(cè)等效諧波電壓源；Zs為系統(tǒng)側(cè)諧波阻抗；Is為系統(tǒng)側(cè)等效諧波電流源；Ic為用戶側(cè)等效諧波電流源；Zc為用戶側(cè)諧波阻抗。

圖1 等效電路Fig.1 Equivalent circuits

對(duì)諧波等效電路分析可得以下方程組：

式中，Upcc、Ipcc分別為PCC處的諧波電壓和諧波電流。通過優(yōu)化高斯過程回歸獲得Upcc、Ipcc與Zs的內(nèi)在聯(lián)系，即可估計(jì)系統(tǒng)側(cè)諧波阻抗。

1.2 高斯過程回歸

高斯過程回歸在處理低維、小樣本和非線性問題中具有顯著優(yōu)勢(shì)，作為一個(gè)隨機(jī)分布，其任意隨機(jī)變量的有限子集均服從聯(lián)合高斯分布。

式中：f()為映射函數(shù)；εi為服從高斯分布N～(0 ,σ2)的高斯白噪聲。

在數(shù)據(jù)集T的有限集合中，一個(gè)聯(lián)合高斯分布由[f(x1),f(x2),…,f(xi)]組成，其高斯過程模型可完全由均值函數(shù)與協(xié)方差函數(shù)決定，即

式中：GP()為高斯分布；X、X′為任意隨機(jī)變量，X,X′∈Rk；m()為均值函數(shù)；k(X,X′)為協(xié)方差函數(shù)，即核函數(shù)；E()為期望。高斯過程回歸常見的核函數(shù)包括平方指數(shù)核、Matern32核、指數(shù)核、有理二次核等。本文選用平方指數(shù)核作為核函數(shù)，具體形式為

式中：σf為標(biāo)準(zhǔn)偏差；σl為特征長(zhǎng)度標(biāo)尺。

當(dāng)式（3）中m(X)=0時(shí)，觀測(cè)值Y的先驗(yàn)分布可表示為

式中：K(X,X)為n維對(duì)稱陣；為噪聲協(xié)方差矩陣，其中σn為方差，In為單位陣；kij為向量Xi和Xj的相似程度，kij越大表示兩個(gè)變量越相似。

給定測(cè)試樣本X*，則其預(yù)測(cè)值Y*與訓(xùn)練樣本的觀測(cè)值Y組成的聯(lián)合先驗(yàn)分布為

式中：N(0,·)為高斯分布；X為訓(xùn)練樣本輸入；X*為測(cè)試樣本輸入；K(X,X)為訓(xùn)練樣本輸入的協(xié)方差矩陣；K(X,X*)為訓(xùn)練樣本與測(cè)試樣本間的協(xié)方差矩陣；K(X*,X)=K(X,X*)T；K(X*,X*)為測(cè)試樣本輸入的協(xié)方差矩陣。

通過邊緣化上述聯(lián)合分布，可得到預(yù)測(cè)值Y*的后驗(yàn)分布為

式中：*為預(yù)測(cè)均值；σ2(Y*)為預(yù)測(cè)方差。

一般情況下，高斯過程回歸中的超參數(shù)集合θ=[σf,σl,σn]可采用極大似然估計(jì)求解最優(yōu)值，其對(duì)數(shù)似然函數(shù)L(θ)為

式中，n為矩陣維度。通過共軛梯度法或牛頓法等求解極大似然函數(shù)即可獲得超參數(shù)集合。

1.3 貝葉斯優(yōu)化

最優(yōu)化超參數(shù)本質(zhì)上就是求取一個(gè)函數(shù)的極值，若函數(shù)形式未知或函數(shù)為非凸函數(shù)，則梯度優(yōu)化等方法無法解決。貝葉斯優(yōu)化恰能解決這類問題，貝葉斯優(yōu)化算法是一種基于模型的序貫優(yōu)化方法，其優(yōu)勢(shì)在于只需要不斷采樣來推測(cè)函數(shù)的極值，經(jīng)過一次評(píng)估之后再進(jìn)行下一次評(píng)估，僅需少數(shù)次目標(biāo)函數(shù)評(píng)估即可獲得最佳參數(shù)，是一種有效的全局優(yōu)化方法。故貝葉斯優(yōu)化在實(shí)際問題中可能比其他優(yōu)化方法具有更好的效果。

貝葉斯優(yōu)化可在一定范圍內(nèi)求解目標(biāo)函數(shù)的最優(yōu)值，即

式中：x*為最優(yōu)值；A為變量x的搜索空間；f(x)為由概率代理模型擬合的目標(biāo)函數(shù)。

貝葉斯優(yōu)化中有兩個(gè)關(guān)鍵部分：①使用概率代理模型代替復(fù)雜的目標(biāo)函數(shù)；②通過代理模型構(gòu)造采集函數(shù)。

1.3.1 概率代理模型

貝葉斯優(yōu)化中概率代理模型的參數(shù)更新依據(jù)貝葉斯定理，可表示為

式中：p(f|D1:i)為f的后驗(yàn)概率分布；p(D1:i|f)為y的似然分布；p(f)為f的先驗(yàn)概率分布；p(D1:i)為邊緣似然分布；D1:i為觀測(cè)樣本集，D1:i={(x1,y1),(x2,y2),…,(xi,yi)}。

概率代理模型可分為參數(shù)模型和非參數(shù)模型。常見的參數(shù)模型有線性模型、廣義線性模型、貝塔-伯努利模型。相比于參數(shù)固定的參數(shù)模型，非參數(shù)模型更具靈活性和可擴(kuò)展性，在貝葉斯優(yōu)化中不易出現(xiàn)“過擬合”，非參數(shù)模型主要包括高斯過程、隨機(jī)森林、深度神經(jīng)網(wǎng)絡(luò)等。

1.3.2 采集函數(shù)

采集函數(shù)是一種根據(jù)后驗(yàn)概率分布主動(dòng)選擇下一個(gè)樣本點(diǎn)進(jìn)行優(yōu)化的策略，包括基于提升策略的概率改善 PI（probability of improvement）、期望改善EI（expected improvement）及基于置信邊界策略的置信上限UCB（upper confidence bound）等。

本文中采用基于EI策略的采集函數(shù)，可表示為

式中：V*為最優(yōu)值；μi(x)為均值；?(M)為高斯分布累計(jì)密度函數(shù)；M=[V*-μi(x)]/σi(x)；σi(x)為標(biāo)準(zhǔn)差；α(x;D1:i)為目標(biāo)函數(shù)。

1.4 BO-GPR計(jì)算諧波阻抗

將BO-GPR用于諧波阻抗計(jì)算，計(jì)算流程如圖2所示。具體步驟如下。

圖2 BO-GPR流程Fig.2 Flow chart of BO-GPR

步驟1通過PCC處測(cè)量數(shù)據(jù)獲得一定數(shù)量的樣本，選擇部分樣本為訓(xùn)練樣本，輸入為諧波電壓Upcc與諧波電流Ipcc組成的二維向量，輸出為諧波阻抗Zs。再選擇一定數(shù)量的樣本作為測(cè)試樣本，組成結(jié)構(gòu)與訓(xùn)練樣本一致。

步驟2根據(jù)訓(xùn)練樣本進(jìn)行高斯過程回歸獲得模型。

步驟3構(gòu)造類似式（12）的代理模型及式（13）和式（14）的采集函數(shù)，使用貝葉斯優(yōu)化獲得噪聲標(biāo)準(zhǔn)差σn的最優(yōu)值。將σn的最優(yōu)值代入高斯過程回歸模型中的式（8）和式（9）得到最終的回歸模型。

步驟4對(duì)比BO-GPR與其他方法的計(jì)算結(jié)果，評(píng)價(jià)本文方法的魯棒性、準(zhǔn)確性及泛化性能。

2 仿真分析

仿真分析中，分別采用4種方法對(duì)諧波模型進(jìn)行阻抗估計(jì)。其中，方法1為二元線性回歸法；方法2為支持向量機(jī)法；方法3為獨(dú)立矢量協(xié)方差法；方法4為本文方法（BO-GPR法）。

2.1 背景諧波為非高斯分布的仿真

根據(jù)圖1中等效電路在Matlab軟件中建立背景諧波為非高斯分布的諧波分析模型（仿真1），系統(tǒng)頻率為50 Hz，具體參數(shù)設(shè)置如下。

（1）系統(tǒng)側(cè)諧波電壓源為 100 V∠53.13°。

（2）用戶側(cè)諧波電流源為Ic=(4.73+j4.74)A，并在實(shí)部添加0.42 A的擾動(dòng)電流，虛部添加0.41 A的擾動(dòng)電流。

（3）諧波阻抗中系統(tǒng)側(cè)諧波阻抗為Zs=(6+j25)Ω，并在實(shí)部添加0.1 Ω的擾動(dòng)阻抗，虛部添加0.3 Ω的擾動(dòng)阻抗；用戶側(cè)諧波阻抗為Zc=(6+j25)Ω，并在實(shí)部添加0.1 Ω的擾動(dòng)阻抗，虛部添加0.3 Ω的擾動(dòng)阻抗。

表1 阻抗均值Tab.1 Mean impedance

由表1結(jié)果可知，除方法3的計(jì)算阻抗誤差較大外，其余3種方法均獲得了很好的計(jì)算結(jié)果，但本文方法（方法4）相對(duì)而言最精確，其計(jì)算阻抗與真實(shí)阻抗完全一致，故當(dāng)背景諧波為非高斯分布時(shí)，本文方法在諧波阻抗計(jì)算中具有高準(zhǔn)確性。

2.2 背景諧波為高斯分布的仿真

為探究本文方法在不同背景諧波下的效果，根據(jù)圖1中諾頓等效電路建立背景諧波為高斯分布的仿真模型，系統(tǒng)頻率為50 Hz，具體參數(shù)設(shè)置如下。

（1）系統(tǒng)側(cè)諧波阻抗Zs服從高斯分布Zs～N(2+j6,0.1+j0.3)，Ω。

（2）用戶側(cè)諧波阻抗Zc服從高斯分布Zc～N(20+j160,0.5+j5)，Ω。

（3）用戶側(cè)諧波電流源Ic服從高斯分布Ic～N(6+j8,0.3+j0.4)，A。

（4）系統(tǒng)側(cè)諧波電流源Is為Ic的k倍，其中k從0.2～1.0每隔0.2取1個(gè)值。

表2 阻抗均值實(shí)部Tab.2 Real part of mean impedance

表3 阻抗均值虛部Tab.3 Imaginary part of mean impedance

對(duì)表2和表3的結(jié)果縱向比較，方法1的誤差隨k的增大而增大，甚至超過了100%，說明該方法在背景諧波非線性程度較高時(shí)效果較差；方法2的誤差隨k的變化也有一定的波動(dòng)，誤差波動(dòng)范圍在7%以內(nèi)，相對(duì)較穩(wěn)定；方法3的誤差隨k的變化趨勢(shì)類似于方法1，其誤差相對(duì)于方法1更小，但仍然具有較大的誤差，獨(dú)立矢量協(xié)方差法只能在Is與Ipcc具有弱相關(guān)性或者近似獨(dú)立時(shí)實(shí)現(xiàn)，而該仿真下這項(xiàng)前提條件并不成立，故出現(xiàn)較大誤差；本文方法（方法4）是4種方法中誤差最小的，且方法4的誤差幾乎不隨k值的變化而變化，具有顯著的穩(wěn)定性。

再結(jié)合表1中的數(shù)據(jù)進(jìn)行對(duì)比可知，方法1在不同的背景諧波下效果差異明顯；方法2在背景諧波特性變化時(shí)效果較好，也具備較高的穩(wěn)定性；方法3在兩種背景諧波下均出現(xiàn)較大的誤差，計(jì)算結(jié)果在源數(shù)據(jù)相關(guān)性較大時(shí)失去參考意義；本文方法（方法4）在兩種背景諧波下均具有極小的誤差，很大程度地抑制了背景諧波變化以及諧波電流源呈一定相關(guān)性帶來的誤差，體現(xiàn)出較好的泛化性能。

3 實(shí)例分析

實(shí)例分析數(shù)據(jù)來自某變電所帶牽引負(fù)荷的110 kV母線。使用電能分析儀測(cè)量PCC處的電壓和諧波電流，進(jìn)行快速傅里葉變換得到諧波電壓與諧波電流，本文以3次諧波為例。在PCC處每隔3 s測(cè)量1次，在10 h內(nèi)共得到12 000個(gè)樣本，分析基波電壓及電流，去除樣本中的空載點(diǎn)、輕載點(diǎn)及離群點(diǎn)，選擇800個(gè)負(fù)載樣本作為研究對(duì)象，得到PCC處帶負(fù)載3次諧波電壓和諧波電流如圖3所示。

圖3 帶負(fù)載時(shí)的3次諧波電壓與電流Fig.3 Third harmonic voltage and harmonic current with load

實(shí)例分析中參考阻抗為（26.45+j42.77）Ω，值得注意的是，此前的研究中常采用二維圖像展示結(jié)果，考慮到輸入數(shù)據(jù)是二維矢量，輸出為一維矢量，使用三維圖像更能凸顯各方法之間的差異。在方法4的計(jì)算過程中，將篩選后的諧波電壓與諧波電流進(jìn)行優(yōu)化高斯過程回歸，并對(duì)結(jié)果進(jìn)行插值，得到3次諧波阻抗如圖4所示。使用均值、平均誤差、均值最大偏差、標(biāo)準(zhǔn)差、均方根誤差及運(yùn)算時(shí)間等指標(biāo)評(píng)價(jià)其結(jié)果，本文方法得到的指標(biāo)參數(shù)如表4所示。

圖4 通過BO-GPR得到的諧波阻抗Fig.4 Harmonic impedance estimated by BO-GPR

前3種方法實(shí)部和虛部的估計(jì)阻抗如圖5～圖7所示。前3種方法各項(xiàng)指標(biāo)參數(shù)見表5～表7。

圖5 通過二元線性回歸得到的諧波阻抗Fig.5 Harmonic impedance estimated by bivariate linear regression

圖6 通過SVM得到的諧波阻抗Fig.6 Harmonic impedance estimated by SVM

圖7 通過獨(dú)立矢量協(xié)方差得到的諧波阻抗Fig.7 Harmonic impedance estimated by random vector covariance

表5 二元線性回歸的評(píng)價(jià)指標(biāo)Tab.5 Evaluation indicators for bivariate linear regression

表6 SVM的評(píng)價(jià)指標(biāo)Tab.6 Evaluation indicators for SVM

表7 獨(dú)立矢量協(xié)方差的評(píng)價(jià)指標(biāo)Tab.7 Evaluation indicators for random vector covariance

對(duì)比圖4～圖7及表4～表7可知，方法1計(jì)算方式簡(jiǎn)單，運(yùn)行時(shí)間短，雖然均值接近參考阻抗，但實(shí)測(cè)諧波數(shù)據(jù)不平穩(wěn)且線性程度低，導(dǎo)致其他指標(biāo)均不理想；方法2在運(yùn)算中可能陷入局部最優(yōu)，使得其結(jié)果出現(xiàn)偏差，并且顯示出一定的波動(dòng)性；方法3在4種方法中效果最差，在系統(tǒng)側(cè)與用戶側(cè)諧波發(fā)射有一定聯(lián)系的情況下，因不滿足方法3的前提條件而出現(xiàn)不適用性；本文方法（方法4）各項(xiàng)指標(biāo)均優(yōu)于其他3種方法，具有一定有效性和實(shí)用性，從圖4可以看出，BO-GPR的結(jié)果可視化后得到的諧波阻抗實(shí)部、虛部在諧波電壓與諧波電流的正負(fù)半軸上具有一定的對(duì)稱性，諧波阻抗具有相對(duì)較小的波動(dòng)性，唯一的不足是計(jì)算時(shí)間較長(zhǎng)，不具實(shí)時(shí)性。

4 結(jié)語

本文將高斯過程回歸與貝葉斯優(yōu)化相結(jié)合并用于計(jì)算諧波阻抗中，提出了一種BO-GPR的諧波阻抗計(jì)算方法，旨在解決背景諧波波動(dòng)較大、系統(tǒng)側(cè)諧波與用戶側(cè)諧波有一定關(guān)聯(lián)程度情況下的諧波阻抗計(jì)算問題。仿真分析驗(yàn)證了本文方法準(zhǔn)確性和穩(wěn)定性較好，實(shí)例分析驗(yàn)證本文方法具有一定的有效性和實(shí)用性。本文方法具有顯著的優(yōu)勢(shì)和不足，其優(yōu)勢(shì)在于：①可以降低非線性背景諧波對(duì)結(jié)果的影響；②能夠避免陷入局部最優(yōu)；③實(shí)例分析中系統(tǒng)側(cè)諧波與用戶側(cè)諧波或有一定關(guān)聯(lián)程度，本文方法對(duì)PCC處諧波電壓和電流的相關(guān)性沒有要求；④在背景諧波特性不同的情況下可得到較好結(jié)果，具有較好的泛化性能。但是，本文方法計(jì)算用時(shí)相對(duì)較長(zhǎng)，在系統(tǒng)穩(wěn)定運(yùn)行時(shí)可根據(jù)已生成模型計(jì)算阻抗，而若系統(tǒng)出現(xiàn)故障或拓?fù)浣Y(jié)構(gòu)發(fā)生改變，此法難以快速反應(yīng)。如何在保證結(jié)果準(zhǔn)確的情況下增強(qiáng)方法的時(shí)效性將是進(jìn)一步研究的方向。