電力系統小干擾穩定性評估的多重譜變換方法

曹玉磊,李崇濤,何勁捷

(西安交通大學電氣工程學院,710049,西安)

保證小干擾穩定是互聯電力系統安全穩定運行的前提[1-2]。目前特征分析方法被廣泛應用于評估大規模系統小干擾穩定性及優化阻尼控制器[3-5]。其中根據是否需要求解系統的全部特征值,相應的算法可以分為全特征值算法[6-8]和部分特征值算法[9-15]。由于全特征值算法需要較大的計算量和內存占用量,因此其往往適用于小規模系統。此外,在小干擾分析過程中通常無需關注占絕大多數的阻尼特性較好的特征值。因此在大規模系統特征分析過程中采用最多的是部分特征值算法。此外,控制器參數設置不合理甚至錯誤將會導致系統產生不穩定特征值,且漏掉其中任意的不穩定特征值都會對系統的小干擾穩定性產生誤判。

在現有小干擾穩定性分析軟件中,最常用的部分特征值算法主要是隱式重啟動Arnoldi算法[11-13]。其在PSD-SSAP、PEALS等國內外眾多小干擾商業軟件中均已采用該算法[16]。此外,與隱式重啟動Arnoldi方法相比,Krylov-Schur方法[14-15]具有相同的時間復雜度,但其在數值穩定性上表現更好,同樣該算法也被應用于Matlab等商業軟件中。這兩種算法的共同特點是均按模值遞減的順序收斂特征值,而判定系統小干擾穩定性的過程中所采用的是不穩定特征值,因此需要采用譜變換方法將不穩定特征值變換為具有最大模值的主導特征值。

目前在商業軟件中常用的譜變換方法主要包括位移逆變換和凱萊變換[17-19]。位移逆變換的特點是其將接近位移點的特征值變換為主導特征值。通過在虛軸或者弱阻尼比線上設定位移點,現有的商業軟件便可以得到幾乎全部弱阻尼特征值,但是這種方法可能會導致遺漏不穩定特征值。這是因為特征值的收斂性隨著遠離位移點而急速衰減,而右半平面上的特征值分布又無法預知。此外,凱萊變換可以將全部不穩定特征值變換成主導特征值,但經過映射后的不穩定特征值的主導性并不明顯且可能產生特征值聚集的情況,從而使得特征值計算困難[20]。這將導致對系統的小干擾穩定性產生誤判。

相應地,指數變換[20-22]在小干擾穩定性評估的過程中表現出了優異的譜性能。不穩定特征值經過指數變換后變得模值較大且遠離其他特征值,使得收斂性顯著提高,但該方法導致需要求解一個復雜的問題,即指數矩陣和向量的乘積。盡管可以利用Krylov子空間方法隱式地求解該乘積,但該乘積的計算始終非常耗時且截斷誤差始終存在,這將嚴重影響特征值計算的效率與精度。

本文提出了多重譜變換方法來評估電力系統的小干擾穩定性,該方法與指數變換類似,在不穩定特征值計算過程中同樣具有良好的譜特性,從而有利于特征值收斂;同時又避免了求解指數變換方法中引起的指數矩陣與向量的乘積,從而具有更高的計算效率和計算精度。最后采用來源于巴西互聯電網的Xingo6u和Xingo3012系統驗證了所提出方法的有效性。

1 電力系統特征分析基礎

下面將給出用于評估系統小干擾穩定性的數學模型及用于特征值計算的Krylov-Schur算法。

1.1 電力系統特征分析模型

電力系統的動態性能通常可用如下的非線性微分代數方程來表示

(1)

式中:x∈Rn×1和y∈Rm×1分別代表狀態和代數變量。然后在系統運行的平衡點處對式(1)進行線性化,便可得到用于小干擾穩定性評估的模型

(2)

(3)

根據A的不穩定特征值便可判斷系統小干擾穩定性。考慮到矩陣A將失去系統的稀疏性,因此在特征值計算過程中將采用稀疏技術[23],并不直接計算A。

1.2 Krylov-Schur算法

對于矩陣M∈Rn×n和非零向量u∈Rn,Krylov子空間Kk(M,u)可表示為

Kk(M,u)=span{u,Mu,…,Mk-1u}

(4)

根據文獻[15]可得到k階Krylov分解

(5)

在得到Bk的k個特征值對(μi,wi),i=1,…,k后,Ritz對(μi,Ukwi)被認為是矩陣M特征值對的近似解。相應的殘差ri可表示為

(6)

如果此時M的主導特征值并未滿足精度要求,就需要進一步采用重啟動技術。首先對Bk進行有序實Schur分解,得到

(7)

式中:Sk∈Rk×k為上三角矩陣;Qk∈Rk×k為變換矩陣。Sk的對角元素代表了Bk的特征值。此外可以通過對Qk的選取,使得Sk的對角元素重新排列。

接下來對k個Ritz值μi重新排列,得到

|μ1|≥…≥|μp|>|μp+1|≥…≥|μk|

(8)

其中不等式左側的p個Ritz值是在收縮過程中需要保留的。然后根據式(8)對矩陣Qk和Sk進行重新排列,得到

(9)

(10)

接下來根據式(10)可以得到如下p階Krylov分解

(11)

(12)

與式(5)相比,從式(12)中得到的Ritz對將是矩陣M特征值對的更優近似。注意到Krylov-Schur算法優先收斂到具有較大模值的特征值,因此需要采用譜變換方法將狀態矩陣A的最右特征值映射為矩陣M的主導特征值。

2 多重譜變換方法

下面將首先回顧常用的譜變換方法及其優缺點,并給出所提出的多重譜變換方法的定義及其與常用譜變換方法之間的聯系。

2.1 常用的譜變換方法

常用的譜變換方法包括位移逆變換和凱萊變換,不同譜變換示意圖如圖1所示。其中位移逆變換的定義為

S=(A-sIn)-1

(13)

式中:s表示位移點;In∈Rn×n為單位矩陣。位移逆矩陣S的特征值μS和特征向量vS與狀態矩陣A的特征值λ和特征向量v之間的對應關系滿足

μS=(λ-s)-1和vS=v

(14)

相應的譜變換示意圖如圖1(c)和圖1(f)所示。當λ平面上的特征值越來越遠離位移點s,對應的μS平面上的特征值主導性將迅速衰減。考慮到無法預知不穩定特征值的分布情況,因此位移逆變換在用于判斷小干擾穩定性時效果并不理想。

圖1 不同譜變換示意圖Fig.1 Schematic diagram of different spectral transforms

此外,凱萊變換的定義為

C=(A-s1In)(A-s2In)-1

(15)

式中:s1和s2表示位移點。凱萊矩陣C的特征值μC和特征向量vC與狀態矩陣A的特征值λ和特征向量v之間的對應關系滿足

(16)

從圖1(a)和圖1(f)中可以看出,λ右半平面上的全部不穩定特征值也均被變換成μC平面上模比1大的主導特征值。同樣地,當特征值λ遠離位移點s2,相應的特征值μC主導性也迅速衰減。同時該映射導致變換后產生特征值聚集的情況,從而使得特征值收斂性降低,求解變得困難。

指數變換在不穩定特征值計算過程中表現出了良好的譜特性,其定義為

(17)

相應地,指數矩陣eA的特征值μE和特征向量vE與狀態矩陣A的特征值λ和特征向量v之間滿足

μE=eλ和vC=v

(18)

從圖1(e)和圖1(f)中可以看出,當λ的實部增大,μE的主導性也隨之增大。這種特性使得不穩定特征值與其他特征值很好地分離,進而有利于收斂。

2.2 存在的問題

位移逆變換和凱萊變換在被用于特征值計算時,譜性能不佳,從而可能導致遺漏不穩定特征值,進而對電力系統的小干擾穩定性產生誤判。相較于位移逆變換和凱萊變換,指數變換在小干擾穩定性評估過程中具有更優異的譜特性,從而使得不穩定特征值易于收斂。

將指數變換作為Krylov-Schur算法的預處理方法時,在式(5)所示的子空間擴展過程中就不可避免的要計算指數矩陣和向量的乘積eAu。盡管文獻[24]中給出了多種求解eAu的方法,但是對于無論哪種方法來說,求解eAu都是非常耗費計算量的。此外,考慮到eA由無窮多項構成,計算eAu時截斷誤差的存在也將嚴重影響特征值計算精度。

下面分別根據位移逆變換和凱萊變換提出了多重譜變換方法,其在不穩定特征值計算過程中同樣具有類似于指數變換的良好譜特性,從而易于特征值收斂;同時又避免了求解指數變換方法中引起的指數矩陣與向量的乘積,具有更高的計算效率和計算精度。

2.3 多重譜變換

根據位移逆變換構造的多重位移逆變換如下

MS=[δ(A-δIn)-1]N

(19)

式中:位移系數δ=N/α,α為放大系數。相應地,多重位移逆矩陣MS的特征值μMS和特征向量vMS與狀態矩陣A的特征值λ和特征向量v之間滿足

(20)

可以看出當N=1時,多重位移逆變換便成為位移逆變換。此外,當N趨近于無窮大時,特征值μMS的主導性與指數變換后特征值eαλ的主導性一致。具體證明過程如下式

(21)

相應的譜變換示意圖如圖1(d)和圖1(f)所示,但注意到多重位移逆變換并不能嚴格保證λ右半平面上的全部不穩定特征值被變換成μMS平面上模比1大的主導特征值。

因此下面進一步根據凱萊變換構造了多重凱萊變換,其定義如下

MC=[(A+σIn)(A-σIn)-1]N

(22)

式中:位移系數σ=2N/α,α為放大系數。相應地,多重凱萊變換矩陣MC的特征值μMC和特征向量vMC與狀態矩陣A的特征值λ和特征向量v之間的一一對應關系滿足

(23)

可以看出當N=1時,多重凱萊變換就是凱萊變換。此外,當N趨近于無窮大時,特征值μMC的主導性也與指數變換后特征值eαλ的主導性一致。具體證明過程如下式

|eαλ|

(24)

從圖1(b)和圖1(f)可以看出,λ右半平面上的全部不穩定特征值均被多重凱萊變換變換成μMC平面上模比1大的主導特征值。考慮到多重位移逆變換并不能保證這一點,因此在以下分析中均以多重凱萊變換為例進行。此外,多重凱萊變換在具有類似于指數變換的譜特性的同時,又成功避免了求解eAu以及由此帶來的高計算量和計算誤差。

此外,在得到多重凱萊變換矩陣MC的特征值μMC和特征向量vMC后,考慮到λvMC=AvMC,相應的逆變換過程如下

(25)

需要說明的是,適當的增大系數α可以增強多重凱萊變換和指數變換的譜性能,從而更有利于特征值的收斂,但在指數變換中增大α將會增大求解eαAu計算量。下面給出簡要說明。eαAu為如下線性時不變系統在t=α時的解

(26)

其可以通過數值積分方法在區間[0,α]上積分得到。顯而易見,增大α將使得積分區間變長,進而增加求解eαAu的計算量。然而在多重凱萊變換中改變系數α并不會影響求解多重凱萊矩陣與向量乘積MCu的計算量。因此在多重凱萊變換中可以令α>1來增強其譜性能。

3 算法流程

在Krylov-Schur算法中采用多重凱萊變換方法作為預處理方法,也就是令M=MC。全部不穩定特征值被多重凱萊變換映射為主導特征值,從而易于收斂得到。相應的計算流程如圖2所示。其中,k表示Krylov分解的維數;ε表示收斂容差。此外重啟動過程中保留的Ritz值的數量p應該略大于模比1大的特征值數q。此外,相比于隱式重啟動Arnoldi算法,這里采用的Krylov-Schur算法除了具有數值穩定性高的優勢外,還可以保證子空間的擴展均在實數域進行,從而在盡可能地避免復數運算的同時進一步提高計算效率。

圖2 不穩定特征值計算的流程Fig.2 Procedures of unstable eigenvalues calculation

4 算例分析

在來源于巴西互聯電網的Xingo6u和Xingo3012系統上驗證所提出方法的有效性。其中巴西互聯電網[25]包含173臺同步機、5 056條支路和3 584條母線。Xingo6u和Xingo3012系統分別有2 940和3 012個狀態變量,17 798和17 932個代數變量,73 916和74 386個非零元素。相關測試在一臺裝有Intel i7處理器和16 GB內存的臺式機上采用Matlab編程實現。收斂容差設置為10-6。本文中所提出的多重凱萊變換中N設為10,放大系數α設為2。

4.1 不同譜變換譜性能對比

對比經過不同譜變換后特征值分布的差異。考慮到凱萊變換可以將全部不穩定特征值變換成主導特征值,下面將首先給出如圖3所示Xingo6u和Xingo3012系統經過凱萊變換(s1=-1,s2=1)后的特征值分布。盡管全部不穩定特征值被映射為主導特征值,但其并未與其他特征值很好地分離,反而相互聚集,導致特征值收斂困難。以Xingo6u系統為例,當Krylov子空間維數設為60,最大重啟動次數設為50時,仍無法收斂到全部不穩定特征值。這將會對系統的小干擾穩定性產生誤判。

測試系統經過位移逆變換(s=1)后的特征值分布如圖4所示。從圖4可以看出,位移逆變換并不能保證將全部不穩定特征值變換成模比1大的主導特征值。考慮到右半平面上的不穩定特征值分布情況又無法預知,因此位移逆變換在判斷小干擾穩定性時性能不佳,可能導致遺漏不穩定特征值的情況發生。

(a)Xingo6u系統

(b)Xingo3012系統

(a)Xingo6u系統

(b)Xingo3012系統

Xingo6u系統和Xingo3012系統經指數變換后的特征值分布如圖5所示。從圖5可以看出,作為對比,指數變換可以將全部不穩定特征值變換成主導特征值且與其他特征值較好地分離,從而有利于特征值的收斂。當Krylov子空間維數設為60,經過2次重啟動便可收斂到兩個測試系統全部不穩定特征值。但指數變換的缺陷在于其需要求解復雜的eAu。

對于所提出的多重凱萊變換來說,在具有類似于指數變換的譜性能的同時避免了直接求解eAu,從而顯著降低了計算的復雜度。經多重凱萊變換后的特征值分布如圖6所示。此外,本文通過引入放大系數α進一步提升了其譜性能。當α=2時,僅需要經過1次重啟動便可收斂到兩個測試系統的全部不穩定特征值。

(a)Xingo6u系統

(a)Xingo6u系統

(b)Xingo3012系統

綜上所述,與凱萊變換、位移逆變換和指數變換相比,通過引入放大系數α,本文提出的多重譜變換在譜主導性和譜稀疏性方面表現更好。這使得Krylov子空間方法更易收斂得到所有的不穩定特征值,從而避免對系統的小干擾穩定性產生誤判。

4.2 計算精度和效率測試

在得到多重凱萊矩陣的全部主導特征值及特征向量后,經過逆變換便可得到Xingo6u和Xingo3012系統的不穩定特征值λ及特征向量v。不穩定特征值λ與誤差‖Av-λv‖見表1。考慮到指數變換也可得到系統的全部不穩定特征值,其結果也在表1中給出。

表1 系統的不穩定特征值與誤差

從表1可以看出,相較于指數變換,多重凱萊變換的計算精度提高了約3個數量級。這是因為對于指數變換來說,求解eAu的過程中截斷誤差的存在嚴重影響了特征值計算的精度。多重凱萊變換避免了直接求解eAu,因此特征值求解精度顯著提高。此外,為了直觀起見,Xingo6u和Xingo3012系統的特征值分布如圖7所示。

(a)Xingo6u系統

(b)Xingo3012系統

進一步測試所提出的多重凱萊變換的計算效率。具體計算效率測試結果見表2,其中NIR表示重啟動次數。根據4.1節可知,位移逆變換和凱萊變換并不能保證收斂到全部不穩定特征值,因此并未對這兩種方法進行計算效率對比。

表2 計算效率測試結果

當Krylov子空間維數設為60時,多重凱萊變換經過1次重啟動便可得到Xingo6u和Xingo3012系統的全部不穩定特征值;當α=1時指數變換方法需要重啟動2次。此外當α進一步設置為2時,指數變換方法也僅需要重啟動1次。對于指數變換來說,增大α反而使得耗時增加。這是因為增大α盡管使得重啟動次數減少了,但求解e2Au也相比于求解eAu的計算量反而增大了約兩倍。

作為對比,多重凱萊變換的計算耗時減少了90%以上。這是因為一方面求解MCu的計算量相比于eAu大大降低;另一方面通過引入放大系數α增強了多重凱萊變換的譜性能,使重啟動次數相應減少。

5 結 論

文中提出了用以評估電力系統的小干擾穩定性的多重譜變換方法,該方法的主要優點如下。

(1)相較于位移逆變換和凱萊變換,多重譜變換方法由于具有類似于指數變換的良好譜特性,特征值易于收斂。

(2)多重譜變換方法又避免了求解指數變換方法引起的eAu以及由此帶來的高計算量和計算誤差,從而具有更高的計算效率和計算精度。

(3)通過引入放大系數α進一步提高了多重譜變換方法的收斂性,進而減少了特征值計算過程所采用的重啟動次數。

(4)根據所提出的多重譜變換方法,文中也揭示了指數變換與位移逆變換和凱萊變換的內在聯系。

(5)相較于隱式重啟動Arnoldi算法,文中通過采用Krylov-Schur算法使得所有子空間擴展均在實數域內進行,從而盡可能地避免了復數運算。