周書宇，蔡國偉，楊德友，王麗馨

（東北電力大學電氣工程學院，吉林省吉林市 132012）

0 引言

隨著現(xiàn)代電力系統(tǒng)的快速發(fā)展，大型區(qū)域電網(wǎng)互聯(lián)使得電力系統(tǒng)小干擾穩(wěn)定問題日益突出，這給電力系統(tǒng)的安全穩(wěn)定運行帶來了巨大挑戰(zhàn)。準確評估小干擾穩(wěn)定特征參數(shù)對系統(tǒng)穩(wěn)定運行狀態(tài)下的低頻振蕩預(yù)警、及時采取有效調(diào)控措施抑制低頻振蕩現(xiàn)象具有重要意義［1-4］。

近年來，同步相量測量單元（PMU）逐步成為監(jiān)測電力系統(tǒng)動態(tài)運行的主要工具，將PMU 采集到的量測數(shù)據(jù)按照數(shù)據(jù)性質(zhì)可分為暫態(tài)振蕩數(shù)據(jù)和隨機響應(yīng)數(shù)據(jù)。暫態(tài)振蕩數(shù)據(jù)指系統(tǒng)受到大擾動沖擊后的自由振蕩響應(yīng)。隨機響應(yīng)數(shù)據(jù)由負荷隨機波動等自然擾動所激發(fā)，外在表現(xiàn)為類噪聲數(shù)據(jù)，這種數(shù)據(jù)在系統(tǒng)中時刻存在，且包含系統(tǒng)小干擾穩(wěn)定特征參數(shù)，利用這類數(shù)據(jù)辨識特征參數(shù)更符合系統(tǒng)振蕩預(yù)警要求，避免了傳統(tǒng)方法只能在大擾動特殊工況下進行特征參數(shù)提取，針對這類數(shù)據(jù)的辨識受到國內(nèi)外研究人員的廣泛關(guān)注［5-14］。

基于隨機響應(yīng)數(shù)據(jù)的辨識方法大致可分為預(yù)處理方法和直接處理方法。預(yù)處理方法先通過預(yù)處理技術(shù)將隨機響應(yīng)信號處理為自由振蕩信號，再利用針對暫態(tài)振蕩信號的方法進行參數(shù)辨識。文獻［8］利用隨機量減技術(shù)（RDT）從系統(tǒng)隨機響應(yīng)信號中獲取自由振蕩信號，再通過Prony 方法辨識振蕩模式。文獻［9］將自然激勵技術(shù)（NExT）與總體最小二乘-旋轉(zhuǎn)不變技術(shù)的信號參數(shù)估計（TLSESPRIT）方法相結(jié)合，通過滑窗的方式獲取模式參數(shù)。 文獻［10］引入NExT 和多參考點復(fù)指數(shù)（PRCE）法在獲得振蕩頻率和阻尼比的基礎(chǔ)上得到模式模態(tài)，同時討論了信噪比對辨識結(jié)果的影響。預(yù)處理方法中預(yù)處理技術(shù)和參數(shù)辨識方法是2 個獨立算法，兩者在數(shù)據(jù)運算過程中誤差的疊加可能導(dǎo)致難以保證辨識精度。直接處理方法則可避免算法之間的累加誤差，直接利用隨機響應(yīng)數(shù)據(jù)進行參數(shù)辨識。文獻［11］將自回歸滑動平均模型（ARMA）法應(yīng)用于特征參數(shù)辨識，但這一辨識過程中存在模型定階問題。文獻［12］提出一種貝葉斯方法在頻域中辨識模態(tài)參數(shù)，大大簡化了識別過程。文獻［13］利用隨機子空間（SSI）算法獲得系統(tǒng)狀態(tài)矩陣，通過特征分解求得模式特征信息。文獻［14］則在SSI 算法基礎(chǔ)上引入遞歸過程來降低LQ 分解計算的復(fù)雜度，提高了算法的在線辨識能力。

傳統(tǒng)的動態(tài)模式分解（DMD）法能從全局角度提取振蕩特征參數(shù)［15］，但由于隨機響應(yīng)數(shù)據(jù)中小干擾穩(wěn)定特征參數(shù)較弱，DMD 法難以對系統(tǒng)隨機響應(yīng)數(shù)據(jù)進行有效辨識，且DMD 法只能選擇在固定的低維正交空間下構(gòu)建系統(tǒng)近似狀態(tài)矩陣，加大了系統(tǒng)低維模型的誤差。為了克服傳統(tǒng)DMD 法的不足，本文提出了一種子空間最優(yōu)模式分解（subspace optimal dynamic mode decomposition，Sub-OpMD）算法，該方法將基于正交投影的矩陣線性變換方法以及共軛梯度優(yōu)化算法引入DMD 法，用于在隨機響應(yīng)數(shù)據(jù)中辨識系統(tǒng)小干擾穩(wěn)定特征參數(shù)。利用IEEE 16 機68 節(jié)點系統(tǒng)和實際電網(wǎng)PMU 數(shù)據(jù)進行仿真驗證。相較于SSI 算法，本文所提方法辨識準確度更高，充分滿足了現(xiàn)代大型區(qū)域互聯(lián)系統(tǒng)對安全穩(wěn)定評估的精度要求。

1 系統(tǒng)隨機響應(yīng)信號模型

在電力系統(tǒng)小干擾穩(wěn)定分析中，運行點處每臺發(fā)電機的振蕩行為可表示為［12］:

式中:M、D、K分別為慣性時間常數(shù)、阻尼系數(shù)和同步力矩系數(shù)；x(t)為發(fā)電機在平衡點附近的功角偏差；F(t)為系統(tǒng)中的自然激勵，主要表現(xiàn)為負荷的隨機波動。

在自然激勵條件下式（1）的解x(t)可表示為:

式中:ψi為模式i的模態(tài)；ηi(t)為模式i滿足運動解耦方程的標量模態(tài)。

式中:ξi為模式i的阻尼比；ωi=2πfi，其中fi為模式i的振蕩頻率；pi(t) 為模式i的激勵，pi(t)=(t)/(ψiTMψi)。

在實際電力系統(tǒng)中由于受到負荷投切、新能源出力不穩(wěn)定等自然激勵的影響，PMU 量測的系統(tǒng)隨機響應(yīng)數(shù)據(jù)外在表現(xiàn)為混亂的類噪聲信號。通過觀察式（2）和式（3）可知，在這種自然激勵條件下，系統(tǒng)看似雜亂的隨機響應(yīng)數(shù)據(jù)中實則包含了豐富的系統(tǒng)動態(tài)信息，利用可靠的隨機數(shù)據(jù)辨識方法可以有效地從中提取系統(tǒng)小干擾特征信息（振蕩頻率、阻尼比、模態(tài)）。

2 小干擾穩(wěn)定評估方法

2.1 DMD 法

DMD 法作為一種高維系統(tǒng)降階算法已廣泛應(yīng)用于流體力學等研究領(lǐng)域［16］，在2015 年首次由Barocio 等學者引入電力系統(tǒng)振蕩模式辨識［15］。下面簡單介紹一下DMD 法原理。

在研究電力系統(tǒng)小干擾穩(wěn)定分析時，可以將電力系統(tǒng)用狀態(tài)空間模型表示為［7］:

式中:A為系統(tǒng)高維映射矩陣。對A進行特征值計算可得到系統(tǒng)特征信息。

由于PMU 采集的信號均為離散量，則可以將式（4）離散化表示為xk+1=Axk，其中k為離散化后的一個采樣點，設(shè)置連續(xù)狀態(tài)變量的采樣時間間隔為Δt，則存在xk=x(kΔt)。利用m個所采集的數(shù)據(jù)序列構(gòu)造2 個輸入矩陣X1=[x1，x2，…，xm]和X2=[x2，x3，…，xm+1]。根據(jù)式（4）的離散化表示可得出X1和X2之間存在一個映射矩陣A，使得:

將高維狀態(tài)矩陣A投影到固定的低維正交空間上，得到A的低維近似狀態(tài)矩陣。

式中:Uq為低維正交空間，可以通過對X1進行奇異值分解（singular value decomposition，SVD）得到，q為子矩陣X1的秩。

式中:Uq和Vq分別為X1的左奇異矩陣和右奇異矩陣；Sq為X1的奇異值矩陣。

通過求解Frobenius 范數(shù)的最小值優(yōu)化問題以計算低維近似矩陣A～。

則系統(tǒng)低維近似狀態(tài)矩陣A～ 為:

DMD 法僅僅適用于從大擾動數(shù)據(jù)中辨識小干擾特征參數(shù)，對外在表現(xiàn)為類噪聲的隨機響應(yīng)數(shù)據(jù)適應(yīng)性較差，無法有效辨識。此外，DMD 法利用SVD 求取固定維數(shù)的低維正交空間，而低維正交空間的維數(shù)選擇不當會導(dǎo)致A～ 無法獲取最佳的系統(tǒng)動態(tài)信息。

2.2 Sub-OpMD 算法

本文提出的Sub-OpMD 算法在DMD 法的基礎(chǔ)上引入隨機子空間算法中基于正交投影的矩陣線性變換［17］以及共軛梯度算法［18-19］，來彌補DMD 法的不足。

構(gòu)造數(shù)據(jù)矩陣X3=[x3，x4，…，xm+2]和X4=[x4，x5，…，xm+3]，并形成2 個新的數(shù)據(jù)矩陣Xp和Xf。

將Xf的行向量投影到Xp的行空間上，構(gòu)造正交投影矩陣O。

對正交投影矩陣O進行SVD 分析:

式中:U和V分別為O的左奇異矩陣和右奇異矩陣；S為O的奇異值矩陣。

正交投影變換能有效壓縮數(shù)據(jù)，使數(shù)據(jù)降維，并減少數(shù)據(jù)降維過程中的信息丟失，SVD 可以在數(shù)據(jù)降維的同時提取數(shù)據(jù)所含有的特征，這2 種矩陣線性變換的應(yīng)用使算法可以有效辨識隨機響應(yīng)數(shù)據(jù)。本文定義U1和U2分別為U的前n行子矩陣和后n行子矩陣，U1和U2之間同樣存在一個包含系統(tǒng)動態(tài)信息的高維映射矩陣，滿足U2U-11→。

引入正交矩陣L，與式（8）類似，則系統(tǒng)低維近似矩陣構(gòu)建問題可表示為:

式（13）可進一步改寫為:

式中:Z=UT2L(LTU2UT2L)-1LTU2。

則式（13）的雙變量（L，Aˉ*）優(yōu)化問題可以轉(zhuǎn)化為單變量（L）優(yōu)化問題:

共軛梯度算法首先計算初始梯度G0與初始搜索方向H0。

式中:L0為正交矩陣L的初始矩陣。

在梯度方向G上將L參數(shù)化表示為:

Lk+1(t)?LkQΓQT+PTQTt∈[0，1]（17）式中:Lk為第k次迭代下的正交矩陣；P和Q分別為第k次迭代下的搜索方向矩陣Hk的左奇異矩陣和右奇異矩陣，即Hk=PNQT，其中N為Hk的奇異值矩陣；Γ和T分別為矩陣Nt對角線元素的余弦值和正弦值構(gòu)成的矩陣。

迭代更新梯度G和搜索方向H:

其中，修正項Δk+1為:

式中:· 表示求內(nèi)積；τ(Gk)=Gk-[LkQTk+P(I-Γk)]PTGk，其中Γk和Tk分別為矩陣Ntk對角線元素的余弦值和正弦值構(gòu)成的矩陣，tk為第k次迭代優(yōu)化t的最小值；τ(Hk)=-(LkQTk+PΓk)NQT。

利用式（18）和式（19）迭代更新得到L*，直至g(Lk+1)和g(Lk)之間的差值在誤差允許范圍內(nèi)，即g(Lk+1)-g(Lk)＜ε時，停止迭代。利用迭代優(yōu)化得到L*，可進一步得到最優(yōu)低維近似狀態(tài)矩陣。

利用共軛梯度算法可以實現(xiàn)式（15）單變量優(yōu)化問題的求解，目的是計算最佳維度的低維正交空間L*，并在此正交空間上找到系統(tǒng)低維近似模型Aˉ*與L*的最優(yōu)組合［18-19］。所構(gòu)建的最優(yōu)系統(tǒng)模型Aˉ*可以更好地描述系統(tǒng)動態(tài)特征。

3 系統(tǒng)特征參數(shù)辨識

小干擾穩(wěn)定性分析過程中特征參數(shù)可以確定某振蕩模式的頻率、振蕩幅值衰減速度（阻尼比）和狀態(tài)變量之間的相對振蕩關(guān)系（模態(tài)）。通過對Sub-OpMD 算法得到的最優(yōu)近似狀態(tài)矩陣Aˉ*進行特征值分解，可以得到上述振蕩特征參數(shù)。

式中:Λ=diag(λ1，λ2，…，λi，…，λl)，其中λi為系統(tǒng)模式i的特征值，l為最優(yōu)近似狀態(tài)矩陣Aˉ*的維數(shù)；Φ=[φ1，φ2，…，φi，…，φl]，其中φi為系統(tǒng)模式i的右特征向量。

計算系統(tǒng)模式i的振蕩頻率fi和阻尼比ξi［20］。

模式i的模態(tài)ψi表示為［18］:

利用Sub-OpMD 算法計算小干擾特征參數(shù)的詳細過程如圖1 所示。

圖1 Sub-OpMD 算法辨識小干擾特征參數(shù)過程Fig.1 Process of Sub-OpMD algorithm for identifying small signal characteristic parameters

4 算例分析

4.1 IEEE 16 機68 節(jié)點系統(tǒng)

本文采用IEEE 16 機68 節(jié)點系統(tǒng)作為仿真算例，該系統(tǒng)可劃分為5 個區(qū)域，詳細接線圖見附錄A圖A1。通過將系統(tǒng)在工作點附近線性化并進行特征分解可以得到15 種振蕩模式，其中4 種區(qū)域間振蕩模式是本文研究的重點。為了全面評價Sub-OpMD 算法的性能，將基于模型的小信號穩(wěn)定分析（SSSA）所得振蕩頻率、阻尼比結(jié)果作為參考值，即模式1 下振蕩頻率為0.375 2 Hz，阻尼比為11.35%；模式2 下振蕩頻率為0.504 7 Hz，阻尼比為7.25%；模式3 下振蕩頻率為0.554 9 Hz，阻尼比為2.80%；模式4 下振蕩頻率為0.749 1 Hz，阻尼比為6.99%。

為了模擬系統(tǒng)中負荷隨機波動所激發(fā)的系統(tǒng)隨機響應(yīng)，本文在系統(tǒng)所有負荷中加入基準值為2%的隨機擾動，以16 臺發(fā)電機的角頻率、功角、發(fā)電機出口功率作為輸入信號，在OptiPlex 7080 工作站（Intel Core i7-10700 CPU，2.90 GHz）上驗證所提算法的有效性。

利用Sub-OpMD 算法對10 min 內(nèi)系統(tǒng)隨機響應(yīng)信號進行辨識，采用滑動數(shù)據(jù)窗方法，分析窗長為2 min，滑動窗長為0.2 s，共得到2 400 個窗口的辨識結(jié)果。4 種模式辨識結(jié)果波形圖見附錄A 圖A2。由圖A2 可以看出，系統(tǒng)振蕩頻率和阻尼比在均值附近波動，這種波動具有一定隨機性，且阻尼比的波動程度明顯大于振蕩頻率。

對10 min 內(nèi)系統(tǒng)隨機響應(yīng)信號辨識結(jié)果進行統(tǒng)計分析，將統(tǒng)計結(jié)果與SSI 算法和SSSA 進行比較，比較結(jié)果如表1 所示。從表1 可知，Sub-OpMD算法得到的4 種區(qū)域間振蕩模式的頻率和阻尼比均值更接近SSSA 得到的參考結(jié)果，同時標準差的統(tǒng)計結(jié)果小于SSI 算法得出的結(jié)果，說明本文提出的方法在振蕩特征參數(shù)辨識上具有較高的可靠性，這種現(xiàn)象在強阻尼模式下更為明顯。

表1 區(qū)域間小干擾特征參數(shù)辨識結(jié)果對比Table 1 Comparison of identification results for inter-area small signal characteristic parameters

本文驗證了Sub-OpMD 算法在不同幅值負荷波動下參數(shù)辨識的魯棒性。由于強阻尼模式表示結(jié)果波動較大，驗證效果更加顯著，故以阻尼比較高的模式1 為例，在原有負荷波動2%的基礎(chǔ)上，加入0.5%、1.0%和3.0%的負荷波動下Sub-OpMD 算法以及SSI 算法的辨識結(jié)果如表2 所示。從表2 可以看出，Sub-OpMD 算法對這4 種不同幅值負荷波動所激勵的隨機響應(yīng)數(shù)據(jù)均具有魯棒性，且辨識結(jié)果精度優(yōu)于SSI 算法。此外，根據(jù)大量仿真實驗，對于負荷波動0.5%以下的信號，本文所提算法的辨識誤差明顯增大，甚至不能有效辨識小干擾特征參數(shù)。

表2 不同負荷波動小干擾特征參數(shù)辨識結(jié)果對比Table 2 Comparison of identification results for small signal characteristic parameters with different load fluctuations

本文針對上文所采用的10 min 數(shù)據(jù)還驗證了分析窗長和滑動窗長的選擇對辨識結(jié)果的影響。首先研究分析窗長對辨識結(jié)果的影響，固定滑動窗長為0.2 s，并分別對分析窗長為0.5、1.0、2.0、3.0 min的數(shù)據(jù)進行辨識。由附錄A 圖A2 可知，阻尼比波動程度明顯大于振蕩頻率，故此次驗證以模式1 的阻尼比辨識結(jié)果為例，驗證所得結(jié)果如圖2 所示。從圖2 可以看出，隨著分析窗長長度的增加阻尼比辨識結(jié)果均值更接近參考值，且標準差更小。而對于滑動窗長而言，固定分析窗長為2 min，分別調(diào)整滑動窗長為0.2、0.3、0.4、0.8 s。經(jīng)過仿真驗證，對于相同一段10 min 數(shù)據(jù)的參數(shù)辨識，不同滑動窗長的選擇對辨識結(jié)果精度影響不大。

圖2 不同分析窗長下阻尼比辨識結(jié)果Fig.2 Identification results of damping ratio with different analysis window lengths

在辨識系統(tǒng)振蕩頻率和阻尼比的同時，還可以得到振蕩模式的模態(tài)，即發(fā)電機組之間的振蕩關(guān)系，它反映了機群發(fā)電機轉(zhuǎn)子之間的非同步搖擺。Sub-OpMD 算法辨識得到了4 個區(qū)域間模式下發(fā)電機角頻率之間的振蕩關(guān)系，統(tǒng)計結(jié)果為每個分析窗口辨識結(jié)果的平均值，詳細模態(tài)圖見附錄A 圖A3。模式1中區(qū)域A、B 發(fā)電機群（G1 至G13）與區(qū)域C、D、E 發(fā)電機群（G14 至G16）相互振蕩；模式2 中區(qū)域C 發(fā)電機（G14）與區(qū)域E 發(fā)電機（G16）相互振蕩；模式3 中區(qū)域A 發(fā)電機群（G1 至G9）與區(qū)域B 發(fā)電機群（G10至G13）相互振蕩；模式4 中區(qū)域D 發(fā)電機（G15）與區(qū)域E 發(fā)電機（G16）相互振蕩。以上Sub-OpMD 算法模態(tài)辨識結(jié)果與SSSA 計算得到的發(fā)電機振蕩分組一致。

4.2 實際系統(tǒng)算例

以某實際系統(tǒng)中PMU 實測數(shù)據(jù)為例，進一步驗證Sub-OpMD 算法在振蕩模式辨識的有效性。選取該系統(tǒng)中4 個PMU 采集的母線頻率、母線電壓幅值、電壓相角數(shù)據(jù)作為輸入數(shù)據(jù)，采集到的10 min內(nèi)電壓幅值波形如圖3（a）所示，對這些采集到的數(shù)據(jù)進行功率譜分析，得到功率譜密度波形如圖3（b）所示，可以看到功率譜密度在0.319 4 Hz 處有明顯尖峰，表示該振蕩模式在這4 個PMU 測量信號中具有很強的可觀性［21-22］，本文重點對這一可觀性較強的模式進行分析。

圖3 實測數(shù)據(jù)波形和功率密度譜Fig.3 Waveform and power density spectrum of measured data

利用Sub-OpMD 算法對10 min 內(nèi)實測數(shù)據(jù)采用滑動數(shù)據(jù)窗方法進行辨識，分析窗長為2 min，滑動窗長為0.5 s，共得到1 000 個窗口的辨識結(jié)果。時域辨識結(jié)果以及相應(yīng)的頻次直方圖如圖4 所示，從圖4 可以看出，辨識結(jié)果呈正態(tài)分布，且多集中于均值附近，這在振蕩頻率上體現(xiàn)得更為明顯。

圖4 實測數(shù)據(jù)振蕩頻率和阻尼比頻次Fig.4 Oscillation frequency and damping ratio frequency of measured data

對實測數(shù)據(jù)頻率和阻尼比的辨識結(jié)果進行統(tǒng)計，如表3 所示。從表3 可以看出，所提方法辨識出的模式振蕩頻率均值為0.332 2 Hz，阻尼比均值為2.90%，相較于SSI 算法更接近功率譜分析得出的真實值（0.319 4 Hz），且標準差更小。該振蕩模式在振蕩關(guān)系上表現(xiàn)為PMU4 相對于PMU1 至PMU3振蕩。

表3 實測數(shù)據(jù)辨識結(jié)果對比Table 3 Comparison of identification results of measured data

5 結(jié)語

本文在分析系統(tǒng)隨機響應(yīng)信號特性的基礎(chǔ)上，將基于正交投影的矩陣線性變換與共軛梯度算法引入DMD 法，提出了隨機數(shù)據(jù)驅(qū)動的Sub-OpMD 算法?；贗EEE 16 機68 節(jié)點系統(tǒng)和某實際測量數(shù)據(jù)的辨識結(jié)果表明，Sub-OpMD 算法可用于多通道數(shù)據(jù)，在小幅隨機激勵的環(huán)境下能從全局角度一次性獲得全部小干擾特征參數(shù)，且辨識結(jié)果可以較準確地反映系統(tǒng)運行狀態(tài)，有利于系統(tǒng)振蕩模式的在線監(jiān)測。通過與SSI 算法進行對比驗證，表明Sub-OpMD 算法的辨識結(jié)果均值更靠近SSSA 的結(jié)果，且標準差更小。

考慮實時評估系統(tǒng)振蕩模式參數(shù)的要求，縮短輸入數(shù)據(jù)矩陣線性變換時間，提高Sub-OpMD 實時應(yīng)用能力將是后續(xù)研究的重點。

附錄見本刊網(wǎng)絡(luò)版（http：//www.aeps-info.com/aeps/ch/index.aspx），掃英文摘要后二維碼可以閱讀網(wǎng)絡(luò)全文。