基于優化高斯過程回歸的諧波阻抗估計方法

夏焰坤 ，唐文張 ，林欣懿

（1.西華大學電氣與電子信息學院，成都 610039；2.西南交通大學牽引動力國家重點實驗室，成都 610031）

由于當前電力系統中存在電力電子設備的投切、高電壓直流輸電換流裝置的使用和電力機車的啟停等，導致電網中諧波源數量飛速增長，諧波污染不可避免，電網運行的穩定性亦會因此降低[1-2]。為有效抑制諧波污染，國際上提出了獎懲性方案，通過評判諧波發射水平并采用經濟手段，對諧波負荷來源進行懲罰或對諧波受害者進行補償[3]。

諧波責任分攤一直是電能質量分析的重點問題[4-5]。諧波功率方向法和阻抗參數法是應用于諧波責任劃分中最廣泛的兩大類方法。其中，諧波功率方向法又包括諧波有功功率方向法[6]、諧波無功功率方向法[7-8]及諧波視在功率識別法[9]等。這些方法均是以功率方向為基礎尋找諧波源的方法，因此諧波源的相位差對結果影響很大，特別是在大量新能源并網的今天，并不能保證此類方法始終適用。

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

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

1 基本原理

1.1 諧波阻抗估計原理

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

圖1 等效電路Fig.1 Equivalent circuits

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

式中，Upcc、Ipcc分別為PCC處的諧波電壓和諧波電流。通過優化高斯過程回歸獲得Upcc、Ipcc與Zs的內在聯系，即可估計系統側諧波阻抗。

1.2 高斯過程回歸

高斯過程回歸在處理低維、小樣本和非線性問題中具有顯著優勢，作為一個隨機分布，其任意隨機變量的有限子集均服從聯合高斯分布。

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

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

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

式中：σf為標準偏差；σl為特征長度標尺。

當式（3）中m(X)=0時，觀測值Y的先驗分布可表示為

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

給定測試樣本X*，則其預測值Y*與訓練樣本的觀測值Y組成的聯合先驗分布為

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

通過邊緣化上述聯合分布，可得到預測值Y*的后驗分布為

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

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

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

1.3 貝葉斯優化

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

貝葉斯優化可在一定范圍內求解目標函數的最優值，即

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

貝葉斯優化中有兩個關鍵部分：①使用概率代理模型代替復雜的目標函數；②通過代理模型構造采集函數。

1.3.1 概率代理模型

貝葉斯優化中概率代理模型的參數更新依據貝葉斯定理，可表示為

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

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

1.3.2 采集函數

采集函數是一種根據后驗概率分布主動選擇下一個樣本點進行優化的策略，包括基于提升策略的概率改善 PI（probability of improvement）、期望改善EI（expected improvement）及基于置信邊界策略的置信上限UCB（upper confidence bound）等。

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

式中：V*為最優值；μi(x)為均值；?(M)為高斯分布累計密度函數；M=[V*-μi(x)]/σi(x)；σi(x)為標準差；α(x;D1:i)為目標函數。

1.4 BO-GPR計算諧波阻抗

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

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

步驟1通過PCC處測量數據獲得一定數量的樣本，選擇部分樣本為訓練樣本，輸入為諧波電壓Upcc與諧波電流Ipcc組成的二維向量，輸出為諧波阻抗Zs。再選擇一定數量的樣本作為測試樣本，組成結構與訓練樣本一致。

步驟2根據訓練樣本進行高斯過程回歸獲得模型。

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

步驟4對比BO-GPR與其他方法的計算結果，評價本文方法的魯棒性、準確性及泛化性能。

2 仿真分析

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

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

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

（1）系統側諧波電壓源為 100 V∠53.13°。

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

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

表1 阻抗均值Tab.1 Mean impedance

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

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

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

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

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

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

（4）系統側諧波電流源Is為Ic的k倍，其中k從0.2～1.0每隔0.2取1個值。

表2 阻抗均值實部Tab.2 Real part of mean impedance

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

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

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

3 實例分析

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

圖3 帶負載時的3次諧波電壓與電流Fig.3 Third harmonic voltage and harmonic current with load

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

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

前3種方法實部和虛部的估計阻抗如圖5～圖7所示。前3種方法各項指標參數見表5～表7。

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

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

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

表5 二元線性回歸的評價指標Tab.5 Evaluation indicators for bivariate linear regression

表6 SVM的評價指標Tab.6 Evaluation indicators for SVM

表7 獨立矢量協方差的評價指標Tab.7 Evaluation indicators for random vector covariance

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

4 結語

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