邵振國，吳瑾櫻，蘇文博

（福州大學 電氣工程與自動化學院，福建 福州 350116）

1 已有的諧波污染用戶模型

由于大量電力電子裝置及其他類型非線性負荷的不斷增加，電網(wǎng)中的諧波污染日益嚴重，降低了電網(wǎng)運行可靠性［1］。建立諧波污染用戶模型有助于分析用戶污染域［2］，設計治理方案。

在穩(wěn)定工作狀態(tài)下，污染用戶的諧波特性可統(tǒng)一表述成式（1）的形式［3］。

其中，Ih為用戶發(fā)出的h次諧波電流；U1和U2、…、Uh分別為用戶節(jié)點電壓中的基波和各次諧波分量；C為用戶控制參數(shù)。

式（1）考慮了不同次數(shù)諧波電壓對用戶諧波電流的影響，可稱之為耦合模型。理論上只要給出電壓波形和控制參數(shù)，通過供電電壓和負荷的伏安特性就能夠精確求解其注入系統(tǒng)的諧波電流值。但是在實際應用中往往不能得到式（1）的精確表達，并且也難以獲得參數(shù)C，因而工程中常常對該模型進行不同程度的簡化以得到一些實用化模型。文獻［4］提出了配電網(wǎng)建模的諾頓等效法，將諧波源負荷表示成阻抗、電流源的組合模型，用一個恒定的阻抗模擬負荷的線性部分，用一個諧波電流源模擬負荷的非線性部分。文獻［5］將不同次的諧波電壓、電流之間的耦合關系線性化，提出了基于交叉頻率導納矩陣的諧波源模型。文獻［6］將諧波源等效為恒流源和受控電流源，然后利用諧波監(jiān)測數(shù)據(jù)進行模型參數(shù)辨識。

文獻［7-8］利用時域狀態(tài)估計技術對通用模型進行參數(shù)辨識。文獻［9］利用模型中電壓、電流之間的物理關系，通過對電流的最佳平方進行逼近辨識出模型參數(shù)，可對各種類型的負荷進行諧波建模。文獻［10］利用現(xiàn)場測量數(shù)據(jù)來提高模型參數(shù)辨識的精度。這些參數(shù)辨識方法能夠準確地得到通用模型的各個參數(shù)，但諧波幅值和相角的測量值是模型參數(shù)辨識必不可少的條件。由于判斷諧波是否超標是以諧波幅值作為基準，同時為了節(jié)省存儲空間，當前的電能質(zhì)量監(jiān)測數(shù)據(jù)往往只保存諧波幅值，而不保存相角，這就限制了以上參數(shù)辨識方法的應用。此外，諧波污染用戶的工作狀態(tài)具有隨機性和時變性，模型參數(shù)不可能恒定不變。如果參數(shù)辨識樣本集包含用戶的多個工況數(shù)據(jù)，那么基于最佳估計思想的參數(shù)辨識結果并不能表征用戶運行特性的變化。

本文提出一種基于諧波幅值監(jiān)測數(shù)據(jù)的諧波用戶建模思路，同時考慮用戶模型參數(shù)的波動特點，建立諧波污染用戶的統(tǒng)計模型。該方法采用用戶的歷史監(jiān)測數(shù)據(jù)，不需要進行額外測量。計算結果為統(tǒng)計型模型參數(shù)，包含了用戶工況的變化特征，因而能夠更全面地評估用戶運行污染。

2 基于諧波幅值監(jiān)測數(shù)據(jù)的污染用戶建模

2.1 用戶污染模型

大電力用戶由變壓器和各種類型的用電設備構成，運行具有時變性和隨機性，在高次諧波情況下還存在多種耦合作用，因此無法精確模擬。在諧波分析中，通常采用圖1所示的負荷集總效應電路結構模擬污染用戶。

圖1 負荷諧波模型Fig.1 Model of load harmonic

圖1 中，Us、Zs分別為系統(tǒng)電壓和系統(tǒng)阻抗；Upcc、Ipcc分別為公共連接點（PCC）的電壓、電流；R、X 分別表示用戶負荷中電阻、電感元件部分的集總效應，一般電力用戶負荷為感性負荷，因此模型中不含有電容元件；諧波電流源Ich表示非線性負荷產(chǎn)生的諧波交叉耦合效應。該模型將諧波污染用戶負荷分為線性負荷部分和非線性負荷部分分別建模。在忽略系統(tǒng)背景諧波的情況下，可由圖2所示的基波等效電路和諧波等效電路分別求解。 其中，Upcc1、Upcch分別為PCC的基波電壓和諧波電壓監(jiān)測值，Ipcc1、Ipcch分別為PCC的基波電流和諧波電流監(jiān)測值。

圖2 諧波等效計算網(wǎng)絡Fig.2 Equivalent calculation networks for harmonic analysis

2.2 線性負荷模型參數(shù)

將用戶的線性負荷部分等效為并聯(lián)阻抗，則其基波阻抗與諧波阻抗的關系如式（3）所示。

其中，Zch為線性負荷諧波等值阻抗；Zc1為線性負荷基波等值阻抗。

基波阻抗值可以由用戶的基波電壓、電流和功率因數(shù)測量值計算得到。

根據(jù)式（2）可以推導出基波電抗的計算公式，如式（6）所示。

2.3 非線性負荷模型參數(shù)

根據(jù)圖2（b）所示電路，有如下關系式：

其中，Upcch、Ich分別為 Upcch、Ich的相量形式。

推導可得諧波電流源的幅值，如式（8）所示。

其中，Zsh為系統(tǒng)諧波等值阻抗，一般表示為串聯(lián)阻抗。系統(tǒng)諧波阻抗與系統(tǒng)基波阻抗?jié)M足如式（19）所示關系。

許多相關文獻已經(jīng)研究了辨識系統(tǒng)諧波阻抗參數(shù)的方法，如文獻［11］利用改進的偏最小二乘法求解系統(tǒng)阻抗，文獻［12］應用秩次回歸估計系統(tǒng)諧波阻抗。這些方法都能精確地估計出系統(tǒng)等效阻抗值。因此，本文在負荷諧波建模時假設系統(tǒng)阻抗已知。

由式（4）—（6）能夠從諧波污染用戶接入點的監(jiān)測數(shù)據(jù)計算得到用戶的線性負荷參數(shù)，由式（3）計算得到用戶的非線性負荷參數(shù)，根據(jù)式（8）可計算得到諧波污染用戶產(chǎn)生的諧波電流值。

3 基于拉丁超立方采樣的統(tǒng)計建模

電力系統(tǒng)中的負荷類型復雜，運行方式也處在不斷的變化之中，因而PCC的監(jiān)測數(shù)據(jù)包含大量的隨機波動成分，更適合采用統(tǒng)計型的模型參數(shù)描述用戶的污染特性。

利用蒙特卡羅模擬（MC）法，對監(jiān)測數(shù)據(jù)分組計算模型參數(shù)，進而進行概率統(tǒng)計分析，可得到統(tǒng)計型模型參數(shù)。根據(jù)GB／T19862—2005《電能質(zhì)量監(jiān)測設備通用要求》中的相關規(guī)定：諧波監(jiān)測的一個基波記錄周期為3s，存儲記錄周期為3min，即每3min產(chǎn)生一個保存記錄。因此諧波污染用戶接入點的監(jiān)測數(shù)據(jù)有成千上萬組，若對每組數(shù)據(jù)都進行計算，則存在計算量太大、占用內(nèi)存太多的缺點，并不實用。為了縮短計算時間，提高計算效率，本文引入基于分層采樣原理的拉丁超立方采樣（LHS）方法處理海量歷史監(jiān)測數(shù)據(jù)，該方法可以在保證計算精度的前提下大幅減小采樣數(shù)，提高計算效率。

3.1 拉丁超立方采樣方法

拉丁超立方采樣是M.D.Mckay等學者在1979年提出的一種分層采樣方法［13］，具有采樣效率高和穩(wěn)健性好等優(yōu)點。在進行相同規(guī)模的抽樣時，該方法可以獲得更高的精度。

拉丁超立方采樣方法分為采樣和排序2個步驟。采樣時把隨機變量累積分布函數(shù)的取值空間均分為N等份，從每個區(qū)間中隨機選取一個數(shù)（或選擇區(qū)間中間）作為累積分布函數(shù)的采樣值，通過求反函數(shù)得到隨機變量的采樣值，形成采樣矩陣S。

電力系統(tǒng)中輸入隨機變量的概率分布特性受眾多因素影響，難以建立概率分布函數(shù)的準確模型。一般只能得到大量離散的監(jiān)測數(shù)據(jù)。針對這種情況，本文利用離散拉丁超立方采樣法處理用戶接入點含有大量離散的歷史監(jiān)測數(shù)據(jù)而沒有其準確的概率分布函數(shù)的情形。

根據(jù)大數(shù)定律，隨機事件發(fā)生的頻率依概率收斂于該事件發(fā)生的概率。因此，離散拉丁超立方采樣并不計算具體的樣本值，而是計算在一定的樣本容量下某一狀態(tài)的數(shù)據(jù)。對于采樣規(guī)模N，輸入隨機變量X處于狀態(tài)xi的個數(shù)與其處于該狀態(tài)的概率成正比，即：

其中，xi為隨機變量X的第i個樣本值；n（xi）為隨機變量X處于狀態(tài)xi的個數(shù)；P（xi）為隨機變量X處于狀態(tài)xi的概率。

采樣過程具體實現(xiàn)如下。

將單個隨機變量的離散數(shù)據(jù)按從小到大的順序排序，并分成N等份，則每個子區(qū)間中包含n／N個樣本。 第 j個子區(qū)間（（j-1） ／N，j／N］的元素在數(shù)組中的分布范圍可由式（11）確定。由拉丁超立方采樣法的原理可知，第j個子區(qū)間對應的樣本可取區(qū)間內(nèi)任意的值，根據(jù)式（12）確定。

其中，a為子區(qū)間元素范圍下限；b為子區(qū)間元素范圍上限；k為第j個子區(qū)間采樣值在排序后的數(shù)字的位置；round為小數(shù)部分四舍五入取整函數(shù)；rand為隨機數(shù)生成函數(shù)。

所有隨機變量都完成采樣后，對采樣樣本進行排序，降低采樣值之間相關性對拉丁超立方采樣模擬精度的影響［14］。排序方法有多種，如Gram-Schmidt序列正交化方法［15］、Cholesky 分解法［16］等。 對于具有相關性的隨機變量，亦可通過排序方法［17-18］使得隨機變量采樣值間的相關性盡可能接近隨機變量間的相關性。對于離散的歷史監(jiān)測數(shù)據(jù)，各個變量之間的相關性可用其相關系數(shù)矩陣表示，結合Cholesky分解法進行排序能夠很方便地處理變量之間的相關性。

3.2 統(tǒng)計建模步驟

利用拉丁超立方采樣法進行諧波污染用戶統(tǒng)計建模的步驟如下：

a.輸入PCC的歷史監(jiān)測數(shù)據(jù)，包括電壓、電流幅值和功率因數(shù)；

b.利用上述單個隨機變量的采樣方法對離散數(shù)據(jù)中的各個變量進行采樣，得到初始樣本矩陣；

c.利用Cholesky分解法，根據(jù)相關系數(shù)矩陣對初始樣本矩陣進行排序，形成計算樣本；

d.對每一組樣本數(shù)據(jù)按式（4）—（6）和式（8）計算模型參數(shù)，并保存計算結果；

e.按步驟a—d完成對所有樣本數(shù)據(jù)的計算后，通過統(tǒng)計學方法求得負荷諧波模型中電阻、電抗、阻抗和諧波電流的期望值與方差等數(shù)字特征及其概率分布曲線。

4 仿真算例

為了驗證本文所提方法的正確性，通過仿真模擬了污染用戶的運行工況，此后采用MATLAB進行建模。仿真中系統(tǒng)側基波阻抗參數(shù)為：Rs1=0.03 Ω、Xs1=0.1 Ω；用戶側基波阻抗參數(shù)為：Rc1=5 Ω、Xc1=15 Ω，用戶5次諧波注入電流均值為14.14 A。假設諧波用戶負荷參數(shù)服從正態(tài)分布，其標準差為均值的10%。利用MATLAB的隨機數(shù)生成器得到負荷節(jié)點的運行樣本數(shù)據(jù)作為監(jiān)測數(shù)據(jù)。

利用拉丁超立方抽樣法對監(jiān)測數(shù)據(jù)進行抽樣，得到基波電壓、電流，諧波電壓、電流和功率因數(shù)的采樣數(shù)據(jù)，本文采樣次數(shù)取500次。對500組采樣數(shù)據(jù)分別進行確定性的負荷模型參數(shù)計算，得到500組結果，進行概率統(tǒng)計分析后可得到模型參數(shù)的期望值、方差和累積分布曲線等。將計算結果與實際的仿真參數(shù)進行對比，結果列在表1中。模型參數(shù)的概率分布曲線分別如圖3所示。根據(jù)圖中結果和表中數(shù)據(jù)對比可以看出，采用拉丁超立方抽樣的計算結果能夠準確地估計模型的統(tǒng)計特征參數(shù)，在大幅減小計算量的情況下，精度較高。

表1 模型參數(shù)計算值與真實值的對比Table1 Comparison between calculated and actual model parameters

圖3 模型參數(shù)的概率分布曲線Fig.3 Probability distribution curve for different model parameters

5 工程算例

實測數(shù)據(jù)來源于某變電站10 kV出線，每3min記錄一次，記錄數(shù)據(jù)為基波和各次諧波電壓、電流幅值，同時可以獲取運行時的功率因數(shù)。此處采用的是基波電壓、電流幅值與5次諧波電壓、電流幅值進行參數(shù)辨識。15d的6900組監(jiān)測數(shù)據(jù)波形如圖4所示，波形從上至下分別為基波電壓、電流幅值與5次諧波電壓、電流幅值。

圖4 電壓、電流幅值監(jiān)測數(shù)據(jù)Fig.4 Monitored measurements of voltage and current

分別采用蒙特卡羅模擬法和拉丁超立方采樣法進行參數(shù)辨識。其中，蒙特卡羅模擬法需要對每一組監(jiān)測數(shù)據(jù)進行一次計算，即需要6900次計算，而拉丁超立方采樣法在500次抽樣時就能達到較好的計算精度。對比可以看出，拉丁超立方采樣法降低計算量的效果十分明顯，而該方法所增加的只是對大量離散監(jiān)測數(shù)據(jù)的采樣過程，該過程額外增加的時間相比于節(jié)省的計算時間要少得多。由此可知，拉丁超立方采樣法具有減小蒙特卡羅法計算量的作用，并且監(jiān)測數(shù)據(jù)的規(guī)模越大，其效果越明顯。

蒙特卡羅模擬法和拉丁超立方采樣法進行統(tǒng)計建模的計算結果的期望值和標準差列在表2中，其概率分布曲線如圖5所示。

通過仿真對比可以看出，本文方法能從大量的歷史監(jiān)測數(shù)據(jù)中準確得到用戶模型參數(shù)的統(tǒng)計特征參數(shù)，從而快速地建立用戶的諧波污染統(tǒng)計模型，具有較好的工程應用價值。

表2 2種方法的計算結果Table 2 Calculated results by two methods

圖5 2種方法得到的模型參數(shù)的概率分布曲線對比Fig.5 Comparison of model parameter probability distribution between two methods

6 結論

針對電網(wǎng)中部分歷史監(jiān)測數(shù)據(jù)只保留用戶電壓、電流諧波幅值，而缺少諧波相角的情形，本文提出了一種基于諧波幅值的建模方法。該方法根據(jù)電壓、電流基波和諧波幅值以及功率因數(shù)計算負荷諧波模型參數(shù)，對實際工程較為實用。同時考慮用戶運行中具有的隨機波動特性，利用拉丁超立方采樣法由大量的歷史監(jiān)測數(shù)據(jù)計算模型參數(shù)的數(shù)字特征及其概率分布曲線。與傳統(tǒng)蒙特卡羅模擬法相比，計算量小且精度高，可用于實際的諧波污染用戶統(tǒng)計建模分析，具有很高的工程應用價值。仿真算例和工程實例證明了所提方法的有效性、準確性和快速性。

［1］林海雪.電力網(wǎng)中的諧波［M］.北京：中國電力出版社，1998：56-78.

［2］CHANG G，HATZIADONIU C，XU W，et al.Modeling devices with nonlinear voltage-current characteristics for harmonic studies［J］.IEEE Transactions on Power Delivery，2004，19 （4）：1802-1811.

［3］Task Force on Harmonics Simulation.Modeling and simulation of the propagation of harmonics in electric power networks，part I：concepts，models，and simulation techniques［J］.IEEE Transactions on Power Delivery，1996，11（1）：452-465.

［4］THUNBERG E，SODER L.A Norton approach to distribution network modeling for harmonic studies［J］.IEEE Transactions on Power Delivery，1999，14（1）：272-277.

［5］FAURI M.Harmonic modelling of non-linear load by means of crossed frequency admittance matrix［J］.IEEE Transactions on Power Systems，1997，12（4）：1632-1638.

［6］趙勇，張濤，李建華，等.一種新的諧波源簡化模型［J］.中國電機工程學報，2002，22（4）：46-51.ZHAO Yong，ZHANG Tao，LI Jianhua，et al.A new simplified harmonic source model for harmonic analysis and mitigation ［J］.Proceedings of the CSEE，2002，22（4）：46-51.

［7］吳篤貴，徐政. 電力負荷的諧波建模［J］. 電網(wǎng)技術，2004，28（3）：20-24.WU Dugui，XU Zheng.Harmonic modeling of electric loads［J］.Power System Technology，2004，28（3）：20-24.

［8］SOLIMAN S A，AL-KANDARI A M，EL-HAWARY M E.Time domain estimation techniques for harmonic load models［J］.Electric Machines and Power Systems，1997，25（8）：885-896.

［9］馮士剛，艾芊.一種諧波負荷建模的新方法［J］.電力系統(tǒng)自動化，2007，31（17）：26-30.FENG Shigang，AI Qian.A new approach for harmonic modeling ofelectricloads［J］.Automation ofElectric PowerSystems，2007，31（17）：26-30.

［10］RUDION K，GUO H，ABILDGAARD H，et al.Non-linear load modeling-requirements and preparation for measurement［C］∥2009 IEEE Power&Energy Society General Meeting.Calgary，AB，Canada：IEEE，2009：1-7.

［11］李麗，馬宏忠，姜寧，等.基于改進偏最小二乘法的系統(tǒng)諧波阻抗及諧波發(fā)射水平估算［J］.電力系統(tǒng)保護與控制，2011，39（1）：92-95.LI Li，MA Hongzhong，JIANG Ning，et al.Assessing harmonic impedance and the harmonic emission level based on improved partial least-squares regression method［J］.Power System Protection and Control，2011，39（1）：92-95.

［12］解紹鋒.基于秩次回歸的系統(tǒng)諧波阻抗分析與諧波評估方法［J］. 電力自動化設備，2010，30（11）：29-33.XIE Shaofeng.Harmonic impedance analysisand harmonic assessment based on rank regression［J］.Electric Power Automation Equipment，2010，30（11）：29-33.

［13］MCKAY M D.A comparison of three methods for selecting values of input variables in the analysis of output from a computer code［J］.Technometrics，2000，42（1）：55-61.

［14］蔡德福，石東源，陳金富.基于多項式正態(tài)變換和拉丁超立方采樣的概率潮流計算方法［J］. 中國電機工程學報，2013，33（13）：92-100.CAI Defu，SHI Dongyuan，CHEN Jinfu.Probabilistic load flow calculation method based on polynomial normal transformation and Latin hypercube sampling［J］.Proceedings of the CSEE，2013，33（13）：92-100.

［15］于晗，鐘志勇，黃杰波，等.采用拉丁超立方采樣的電力系統(tǒng)概率潮流計算方法［J］. 電力系統(tǒng)自動化，2009，33（21）：32-35.YU Han，CHUANG Chiyong，WONG Kitpo，et al.A probabilistic load flow calculation method with Latin hypercube sampling［J］.Automation of Electric Power Systems，2009，33（21）：32-35.

［16］YU Han，CHUANG Chiyong，WONG Kitpo.Probabilistic load flow evaluation with hybrid Latin hypercube sampling and Cholesky decomposition［J］.IEEE Transactions on Power Systems，2009，24（2）：661-667.

［17］陳雁，文勁宇，程時杰.考慮輸入變量相關性的概率潮流計算方法［J］. 中國電機工程學報，2011，31（22）：80-87.CHEN Yan，WEN Jinyu，CHENG Shijie.Probabilistic load flow analysis considering dependencies among input random variables［J］.Proceedings of the CSEE，2011，31（22）：80-87.

［18］徐瀟源，嚴正，馮冬涵，等.基于輸入變量秩相關系數(shù)的概率潮流計算方法［J］. 電力系統(tǒng)自動化，2014，38（12）：54-61.XU Xiaoyuan，YAN Zheng，F(xiàn)ENG Donghan，et al.Probabilistic load flow calculation based on rank correlation coefficient of input random variables［J］.Automation of Electric Power Systems，2014，38（12）：54-61.