電磁暫態(tài)仿真程序中高斯噪聲的注入方法研究

周博英，唐可翾，沈沉

（清華大學(xué) 電機(jī)工程與應(yīng)用電子技術(shù)系，北京 100084）

含高比例可再生能源發(fā)電和高比例電力電子設(shè)備（“雙高”）正成為電力系統(tǒng)發(fā)展的重要趨勢(shì)和關(guān)鍵特征?！半p高”電力系統(tǒng)的動(dòng)態(tài)特性中，隨機(jī)性和非線性特性進(jìn)一步凸顯，其中隨機(jī)性越來(lái)越成為影響電力系統(tǒng)安全和穩(wěn)定運(yùn)行的重要因素。傳統(tǒng)電力系統(tǒng)穩(wěn)定性往往使用確定性的時(shí)域仿真法對(duì)系統(tǒng)的微分方程進(jìn)行數(shù)值分析，并基于概率論采用蒙特卡洛等分析方法對(duì)系統(tǒng)狀態(tài)進(jìn)行重復(fù)模擬及統(tǒng)計(jì)分析。這種基于給定的系統(tǒng)模型和數(shù)據(jù)的方法很難全面刻畫系統(tǒng)中的隨機(jī)因素。

隨機(jī)響應(yīng)分析側(cè)重于研究系統(tǒng)受動(dòng)態(tài)時(shí)變隨機(jī)干擾的情況，通過(guò)注入隨機(jī)噪聲并對(duì)系統(tǒng)隨機(jī)響應(yīng)統(tǒng)計(jì)特征進(jìn)行分析，可在電磁暫態(tài)仿真中實(shí)現(xiàn)電力系統(tǒng)隨機(jī)穩(wěn)定性動(dòng)態(tài)分析。電力系統(tǒng)隨機(jī)因素大多滿足維納過(guò)程（布朗運(yùn)動(dòng)），可通過(guò)高斯型噪聲進(jìn)行建?！，F(xiàn)有的研究很少關(guān)注如何在電磁暫態(tài)仿真中隨仿真步長(zhǎng)變化生成隨機(jī)噪聲，且缺少對(duì)噪聲注入后隨機(jī)響應(yīng)仿真結(jié)果正確性的驗(yàn)證。本文擬在電磁暫態(tài)仿真中實(shí)現(xiàn)對(duì)高斯型隨機(jī)噪聲的正確仿真，為電力系統(tǒng)基于隨機(jī)響應(yīng)的動(dòng)態(tài)分析提供支撐。

本文使用的電磁暫態(tài)仿真平臺(tái)為CloudPSS 平臺(tái)［1］，由清華大學(xué)和清華四川能源互聯(lián)網(wǎng)研究院共同研發(fā)，是一套面向新型電力系統(tǒng)、多時(shí)間尺度建模的高性能電磁暫態(tài)仿真云平臺(tái)。

1 電力系統(tǒng)高斯白噪聲建模

在電力系統(tǒng)研究中，通常假設(shè)隨機(jī)過(guò)程為高斯分布的白噪聲。對(duì)于受到隨機(jī)噪聲干擾的系統(tǒng)，如果注入的隨機(jī)噪聲干擾滿足維納過(guò)程，則在動(dòng)力學(xué)分析時(shí)可以通過(guò)高斯白噪聲模型進(jìn)行模擬［2］。白噪聲是一種最簡(jiǎn)單的隨機(jī)過(guò)程，是由一系列不相關(guān)的隨機(jī)變量組成的理想化隨機(jī)過(guò)程。

電力系統(tǒng)的隨機(jī)因素大多滿足維納過(guò)程，連續(xù)隨機(jī)擾動(dòng)因素大部分可以用高斯白噪聲隨機(jī)擾動(dòng)進(jìn)行描述［3］，例如系統(tǒng)諧波的處理、負(fù)荷的隨機(jī)波動(dòng)、風(fēng)電功率的隨機(jī)變化、機(jī)械功率的隨機(jī)扭振等，均可視為高斯白噪聲擾動(dòng)［4］。

高斯白噪聲的標(biāo)準(zhǔn)矢量隨機(jī)微分方程如式（1）所示［5］：

其中，X(t)為n維矢量隨機(jī)變量，W(t)表示n維維納過(guò)程，W(t)的形式導(dǎo)數(shù)定義為W(t)/dt=η(t)，為高斯過(guò)程。本文高斯白噪聲統(tǒng)一用符號(hào)η表示。

線性系統(tǒng)高斯白噪聲隨機(jī)擾動(dòng)可視作系統(tǒng)輸入，向系統(tǒng)的狀態(tài)方程和輸出方程進(jìn)行注入。注入擾動(dòng)后線性系統(tǒng)的隨機(jī)微分方程如式（2）所示：

其中，x為狀態(tài)變量，y為輸出變量，η為服從標(biāo)準(zhǔn)正態(tài)分布的高斯白噪聲，a，b，c，d為對(duì)應(yīng)參數(shù)。

電力系統(tǒng)通常由高維非線性微分代數(shù)方程組描述，注入噪聲后的微分代數(shù)方程如式（3）所示：

發(fā)電側(cè)的高斯型隨機(jī)擾動(dòng)主要指原動(dòng)機(jī)隨機(jī)變化導(dǎo)致的發(fā)電機(jī)輸入機(jī)械功率的隨機(jī)波動(dòng)［6］。以單機(jī)無(wú)窮大系統(tǒng)為例［7］，加入隨機(jī)擾動(dòng)項(xiàng)后的系統(tǒng)方程如式（4）所示：

用電側(cè)負(fù)荷也具有較強(qiáng)的不確定性和波動(dòng)性。在研究隨機(jī)擾動(dòng)對(duì)電力系統(tǒng)電壓穩(wěn)定的影響時(shí)，需要研究電力系統(tǒng)無(wú)功負(fù)荷的隨機(jī)波動(dòng)。在無(wú)功負(fù)荷中添加高斯白噪聲后，其定義如式（5）所示［8］：

此項(xiàng)無(wú)功負(fù)荷的波動(dòng)會(huì)反映到系統(tǒng)電壓幅值和相角的方程中，導(dǎo)致系統(tǒng)電壓出現(xiàn)隨機(jī)波動(dòng)。

電力系統(tǒng)在日常運(yùn)行過(guò)程中，母線電壓、相角和線路功率等信號(hào)均存在類似噪聲的微小幅值波動(dòng)。這類信號(hào)定義為類噪聲信號(hào)。若將負(fù)載變化看做完全隨機(jī)的過(guò)程，則類噪聲信號(hào)可看做白噪聲信號(hào)激勵(lì)電力系統(tǒng)產(chǎn)生的響應(yīng)。通過(guò)對(duì)類噪聲信號(hào)進(jìn)行隨機(jī)系統(tǒng)建模分析，有望估計(jì)出系統(tǒng)的動(dòng)態(tài)模型［9］。

2 含高斯白噪聲微分方程一般求解

為探究電磁暫態(tài)仿真過(guò)程中高斯白噪聲的注入方法及其正確性，首先向只含一個(gè)狀態(tài)變量的一階RL電路中注入高斯白噪聲隨機(jī)擾動(dòng)。含噪聲的系統(tǒng)微分方程如式（6）所示：

在電磁暫態(tài)仿真中，動(dòng)態(tài)元件經(jīng)過(guò)諾頓等效可視為諾頓等效電導(dǎo)與諾頓等效歷史電流源并聯(lián)形式［10］，電路方程如式（7）所示：

其中，諾頓等效電導(dǎo)Geff=0.5ΔtL?1，諾頓等效歷史電流源

根據(jù)隨機(jī)微分方程求解法則，當(dāng)積分步長(zhǎng)Δt較小時(shí)，(Wt?Wt?Δt)可以用半階無(wú)窮小的高斯過(guò)程來(lái)代替。(Wt?Wt?Δt)為獨(dú)立變量，且與當(dāng)前時(shí)刻電路狀態(tài)無(wú)關(guān)。得到遞推公式如式（8）所示：

其中，σN(0，1)為添加的隨機(jī)擾動(dòng)項(xiàng)。

為了最小程度地改變仿真拓?fù)洌S機(jī)擾動(dòng)項(xiàng)可以視作電流源iweff并聯(lián)于原電感元件，即可以將含有隨機(jī)擾動(dòng)的電流源拆分成2 個(gè)元件，得到含高斯白噪聲的RL電路諾頓等效電路圖如圖1所示。

圖1 含高斯白噪聲的RL電路諾頓等效電路圖Fig.1 Norton equivalent circuit diagram of RL circuit with Gaussian white noise

該電路的約束方程如式（9）所示：

原本的電感支路iL0(t)記錄了原電感元件暫態(tài)電流，而iweff(t)記錄了每一時(shí)刻高斯白噪聲注入后隨機(jī)擾動(dòng)等效電流源對(duì)系統(tǒng)的電流貢獻(xiàn)。按照該等效電路圖的思路，在CloudPSS平臺(tái)上對(duì)含高斯白噪聲的RL電路進(jìn)行建模與仿真，仿真電路拓?fù)浣Y(jié)構(gòu)如圖2所示。

圖2 CloudPSS平臺(tái)仿真含高斯白噪聲的RL電路Fig.2 CloudPSS simulation of RL circuit with Gaussian white noise

根據(jù)隨機(jī)分析的數(shù)學(xué)知識(shí)，一階RL簡(jiǎn)單電路的隨機(jī)微分方程解析解如式（10）所示：

一階隨機(jī)微分方程的狀態(tài)變量中注入高斯白噪聲，iL(t)解析解也為高斯過(guò)程，服從正態(tài)分布，即

計(jì)算該高斯過(guò)程的方差及自協(xié)方差函數(shù)理論值，得到一階RL電路系統(tǒng)統(tǒng)計(jì)特征解析解如式（12）所示：

將仿真結(jié)果的方差和自協(xié)方差函數(shù)與理論值進(jìn)行對(duì)比，結(jié)果如圖3 所示，驗(yàn)證了該仿真方法在對(duì)含高斯白噪聲的RL電路進(jìn)行仿真時(shí)的正確性。

圖3 RL電路統(tǒng)計(jì)特征仿真值與理論值對(duì)比Fig.3 Comparison between simulation value and theoretical value of RL circuit statistical characteristics

3 高斯白噪聲注入微分代數(shù)方程所引發(fā)的問(wèn)題

沿用上節(jié)中的處理方法，在單機(jī)單負(fù)荷系統(tǒng)的負(fù)荷側(cè)有功功率中注入高斯白噪聲，即在這個(gè)電路中將高斯白噪聲注入代數(shù)變量中進(jìn)行仿真。

在圖4 所采用的單機(jī)單負(fù)荷電路中，E=30 kV，R+jX=(0.1+1.6×10?6j)Ω，負(fù)荷數(shù)據(jù)為Pset=65 MW，Qset=20 MVar。在電力系統(tǒng)由微分代數(shù)方程描述的情況下，同樣用式（8）的半階無(wú)窮小方法對(duì)高斯白噪聲進(jìn)行處理，對(duì)代數(shù)變量中注入隨機(jī)擾動(dòng)后的系統(tǒng)隨機(jī)響應(yīng)進(jìn)行仿真。

圖4 含高斯白噪聲的單機(jī)單負(fù)荷電路拓?fù)鋱DFig.4 Single machine single load circuit topology with Gaussian white noise

仿真預(yù)期的正確結(jié)果為變量的方差不隨仿真步長(zhǎng)的變化而變化，但在CloudPSS仿真平臺(tái)上測(cè)試后發(fā)現(xiàn)（仿真結(jié)果如圖5 所示），無(wú)論是否進(jìn)行半階無(wú)窮小處理，系統(tǒng)無(wú)功功率的方差均會(huì)隨著仿真步長(zhǎng)的變化而變化。

圖5 無(wú)功功率方差隨仿真步長(zhǎng)變化圖Fig.5 Reactive power variance changes with simulation step

仿真步長(zhǎng)變化導(dǎo)致噪聲響應(yīng)的性質(zhì)發(fā)生變化是不符合物理規(guī)律的。這一現(xiàn)象說(shuō)明，第2 節(jié)中針對(duì)隨機(jī)微分方程的求解方法并不適用于隨機(jī)微分代數(shù)方程的求解。實(shí)際上，當(dāng)系統(tǒng)由微分代數(shù)方程描述時(shí)，無(wú)論是否采用式（8）的半階無(wú)窮小處理方法，均會(huì)出現(xiàn)噪聲響應(yīng)統(tǒng)計(jì)特征隨步長(zhǎng)變化而變化的問(wèn)題。因此，必須研究新的噪聲處理方法以求解微分方程與代數(shù)方程共同約束的系統(tǒng)。

4 高斯有色噪聲的構(gòu)造

為了解決含高斯白噪聲的微分代數(shù)方程求解問(wèn)題，嘗試將注入的高斯白噪聲替換為高斯有色噪聲。事實(shí)上，現(xiàn)實(shí)中白噪聲過(guò)程并不能很好地描述自然界的各種擾動(dòng)［11］，理想的白噪聲在物理上是很難實(shí)現(xiàn)的，并且實(shí)際工程中測(cè)量的數(shù)據(jù)所包含的噪聲往往是有色噪聲。使用PMU（Phase Measurement Unit）對(duì)實(shí)際電力系統(tǒng)進(jìn)行數(shù)據(jù)采集時(shí)，會(huì)在理想高斯白噪聲的基礎(chǔ)上混雜一部分高斯有色噪聲［12］，這部分高斯有色噪聲是PMU 防混疊濾波器產(chǎn)生的［13］。同時(shí)，高斯白噪聲經(jīng)過(guò)低通濾波處理也將形成高度相關(guān)的高斯有色噪聲［14］。

高斯白噪聲的功率譜密度呈全頻域均勻分布，其自相關(guān)函數(shù)為δ沖激函數(shù)，方差為無(wú)窮大。而高斯有色噪聲的功率譜密度和自相關(guān)函數(shù)則如圖6所示［15］，其方差為有限值。不同于高斯白噪聲，有色噪聲序列中不同時(shí)刻的數(shù)據(jù)具有相關(guān)性。有色噪聲可視為白噪聲序列經(jīng)過(guò)線性環(huán)節(jié)的輸出，是理想白噪聲的一種近似。

圖6 高斯白噪聲與有色噪聲的功率譜密度及自相關(guān)函數(shù)Fig.6 Curves of power spectral density and autocorrelation function of Gaussian white noise and colored noise

考慮到電磁暫態(tài)仿真計(jì)算過(guò)程的單向遞推特性，高斯有色噪聲的構(gòu)造采取由高斯白噪聲插值的方法。圖7 展示了如何利用插值方法由高斯白噪聲構(gòu)造高斯有色噪聲。(dT)?1為高斯白噪聲注入的頻率，為保證噪聲的性質(zhì)相同，dT是不隨仿真步長(zhǎng)變化而變化的，即在不同的仿真步長(zhǎng)下，高斯白噪聲注入頻率相同。dt為實(shí)際仿真步長(zhǎng)，則(dT/dt?1)即為構(gòu)造有色噪聲時(shí)每個(gè)dT時(shí)段中插值點(diǎn)的數(shù)量。隨著仿真步長(zhǎng)的減小，插入數(shù)據(jù)的密度增加。

圖7 插值構(gòu)造高斯有色噪聲示意圖Fig.7 Diagram of Gaussian colored noise construction by interpolation

由這種方法構(gòu)造出的有色噪聲序列相關(guān)性與白噪聲略有差異，不同dT時(shí)段的數(shù)據(jù)之間沒有相關(guān)性，但每段dT中插入的數(shù)據(jù)點(diǎn)之間具有相關(guān)性。

為了方便在電磁暫態(tài)仿真中實(shí)現(xiàn)，選取階梯插值與分段線性插值2 種插值方法。經(jīng)插值后產(chǎn)生的高斯有色噪聲序列如圖8所示。

圖8 不同插值方法產(chǎn)生的高斯有色噪聲Fig.8 Gauss colored noise produced by different interpolation methods

圖9 不同插值點(diǎn)數(shù)的噪聲序列自相關(guān)函數(shù)Fig.9 Autocorrelation function of noise sequence with different interpolation points

圖10 不同插值點(diǎn)數(shù)的噪聲序列方差Fig.10 Variance of noise sequence with different interpolation points

圖11 不同仿真步長(zhǎng)下輸入噪聲的自相關(guān)函數(shù)Fig.11 Autocorrelation function of input noise sequence with different simulation steps

圖12 不同插值點(diǎn)數(shù)的輸入噪聲方差Fig.12 Variance of input noise sequence with different interpolation points

圖13 不同插值點(diǎn)數(shù)的系統(tǒng)無(wú)功功率響應(yīng)方差Fig.13 Variance of systemreactive power response with different interpolation points

圖8中，高斯白噪聲的注入頻率為5 × 105Hz，2 個(gè)相鄰的白噪聲數(shù)據(jù)中插值199 個(gè)數(shù)據(jù)點(diǎn)，即生成的高斯有色噪聲序列頻率為1×108Hz。

對(duì)不同插值方法產(chǎn)生的序列方差進(jìn)行理論分析。由于階梯插值所插入的數(shù)據(jù)點(diǎn)均與原白噪聲相同，因此所產(chǎn)生的有色噪聲均值和方差應(yīng)當(dāng)均與原白噪聲相同。而對(duì)于線性插值產(chǎn)生的有色噪聲，首先考慮隨機(jī)變量X與Y相互獨(dú)立且均服從正態(tài)分布N（0，1）時(shí)，隨機(jī)變量aX+(1?a)Y的方差為=a2+(1?a)2。類似的，在服從正態(tài)分布N（0，1）的高斯白噪聲序列中的兩點(diǎn)(xi，yi)與(xi+1，yi+1)之間進(jìn)行線性插值時(shí)，兩端的數(shù)據(jù)點(diǎn)可視為相互獨(dú)立的隨機(jī)變量，中間的插值數(shù)據(jù)點(diǎn)[kyi+(1?k)yi+1]，k∈(0，1)服從正態(tài)分布，其均值為0，方差為=k2+(1?k)2。因此，當(dāng)插值點(diǎn)數(shù)為m時(shí)，每段dT中的全部數(shù)據(jù)點(diǎn)為[j(m+1)?1yi+(m+1?j)(m+1)?1yi+1]，j=0，1，2，3，…，m。經(jīng)過(guò)必要的近似處理（如忽略樣本方差與整體方差的差別及忽略首末端數(shù)據(jù)的影響等）之后，可以計(jì)算出線性插值有色噪聲序列的方差如式（13）所示：

計(jì)算結(jié)果表明，對(duì)于線性插值后的有色噪聲，當(dāng)插值點(diǎn)數(shù)較大時(shí)，有色噪聲序列方差比原高斯白噪聲序列方差有所減小，并隨插值點(diǎn)數(shù)的增加趨近于固定值的2/3。

為驗(yàn)證噪聲插值后統(tǒng)計(jì)特征的變化規(guī)律，在Matlab中對(duì)不同仿真步長(zhǎng)（插值點(diǎn)數(shù)）下的有色噪聲統(tǒng)計(jì)特征進(jìn)行數(shù)值驗(yàn)證，結(jié)果如圖9-10所示。

在數(shù)值實(shí)驗(yàn)中，為避免插值前噪聲序列的偶然隨機(jī)性對(duì)方差結(jié)果產(chǎn)生持續(xù)的影響，在不同的插值點(diǎn)數(shù)下需要重新生成正態(tài)分布的隨機(jī)噪聲序列。

Matlab 實(shí)驗(yàn)結(jié)果表明，插值前后噪聲序列的統(tǒng)計(jì)特征變化符合預(yù)期，這為進(jìn)一步在電磁暫態(tài)仿真中將有色噪聲作為系統(tǒng)輸入提供了正確性保證。

對(duì)于階梯插值，同一噪聲序列插值所得到的新序列與原白噪聲序列的自相關(guān)函數(shù)趨勢(shì)及方差大小均基本相同，且不隨插值點(diǎn)數(shù)量的變化而變化。

對(duì)于線性插值，在插值點(diǎn)數(shù)較大（≥20）時(shí)，插值后序列的統(tǒng)計(jì)特征基本不隨插值點(diǎn)數(shù)的變化而變化，且方差大致穩(wěn)定在原序列方差的2/3；在插值點(diǎn)數(shù)較?。?20）時(shí)，插值點(diǎn)數(shù)的變化對(duì)于插值后序列的方差及自相關(guān)函數(shù)的峰值有著較大的影響，插值后序列的方差可由式（13）進(jìn)行定量計(jì)算。在實(shí)際仿真中使用線性插值方法得到有色噪聲序列時(shí)，若想要保持插值后序列的統(tǒng)計(jì)特征穩(wěn)定（不隨仿真步長(zhǎng)變化而變化），則仿真頻率應(yīng)比噪聲頻率高20 倍以上；若想要保持插值后序列的方差大小與原白噪聲序列基本相同，則可以基于式（13）的定量計(jì)算結(jié)果通過(guò)其他方式補(bǔ)償。

5 高斯有色噪聲的注入

使用2 種插值方法分別構(gòu)造有色噪聲，并在CloudPSS 平臺(tái)上分別作為隨機(jī)干擾注入到圖4 所示的單機(jī)單負(fù)荷系統(tǒng)的負(fù)荷有功功率中。電路參數(shù)與上文相同，隨機(jī)數(shù)據(jù)的注入頻率為1 × 104Hz。采用不同仿真步長(zhǎng)（插值點(diǎn)數(shù)）進(jìn)行仿真，每次仿真均重新生成高斯白噪聲序列，插值后輸入系統(tǒng)的有色噪聲序列和系統(tǒng)響應(yīng)（負(fù)荷無(wú)功功率）的統(tǒng)計(jì)特征如圖11-13所示。

由圖11-13 可見，仿真結(jié)果的理論計(jì)算結(jié)果和數(shù)值實(shí)驗(yàn)結(jié)果基本相符。

對(duì)于階梯插值，插值點(diǎn)數(shù)（仿真步長(zhǎng)）不同時(shí)，插值后噪聲序列的自相關(guān)函數(shù)形態(tài)及方差大小基本不變，系統(tǒng)無(wú)功功率響應(yīng)的方差在不同的插值點(diǎn)數(shù)下也基本保持不變。對(duì)于線性插值，在插值點(diǎn)數(shù)較大時(shí)，所得到的有色噪聲方差趨近于原高斯白噪聲方差的2/3，噪聲的性質(zhì)以及注入噪聲后系統(tǒng)響應(yīng)的統(tǒng)計(jì)特征在不同的仿真步長(zhǎng)下均基本保持不變；在插值點(diǎn)數(shù)較小時(shí)，仿真步長(zhǎng)的變化對(duì)噪聲序列及系統(tǒng)響應(yīng)的統(tǒng)計(jì)特征則有較大影響。仿真結(jié)果驗(yàn)證了通過(guò)插值方法生成高斯有色噪聲并注入電磁暫態(tài)仿真的可行性與合理性。

6 結(jié)論

電力系統(tǒng)中的隨機(jī)擾動(dòng)多為高斯型噪聲。在電磁暫態(tài)仿真中注入高斯型噪聲隨機(jī)擾動(dòng)，需要保證系統(tǒng)的隨機(jī)響應(yīng)統(tǒng)計(jì)特征不隨仿真步長(zhǎng)變化而變化。在RL基礎(chǔ)電路中可以將高斯白噪聲隨機(jī)擾動(dòng)等效為并聯(lián)半階無(wú)窮小電流源的形式加入系統(tǒng)，噪聲注入系統(tǒng)后的隨機(jī)響應(yīng)統(tǒng)計(jì)特征與理論值相符。但該方法無(wú)法應(yīng)用于求解含高斯白噪聲的微分-代數(shù)系統(tǒng)。采用階梯插值以及插值點(diǎn)數(shù)較大的線性插值方法所產(chǎn)生的有色噪聲統(tǒng)計(jì)特征在一定程度上與插值前的白噪聲相似，通過(guò)將注入的白噪聲替換為經(jīng)插值產(chǎn)生的有色噪聲，能夠規(guī)避隨機(jī)微分代數(shù)方程求解中所遇到的問(wèn)題。在電磁暫態(tài)仿真中，通過(guò)向系統(tǒng)的微分代數(shù)方程組中注入本文方法所產(chǎn)生的隨機(jī)有色噪聲，可以實(shí)現(xiàn)對(duì)系統(tǒng)隨機(jī)響應(yīng)統(tǒng)計(jì)特征的正確分析。