曹玉磊,李崇濤,何勁捷

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

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

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

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

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

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

1 電力系統(tǒng)特征分析基礎

下面將給出用于評估系統(tǒng)小干擾穩(wěn)定性的數(shù)學模型及用于特征值計算的Krylov-Schur算法。

1.1 電力系統(tǒng)特征分析模型

電力系統(tǒng)的動態(tài)性能通?？捎萌缦碌姆蔷€性微分代數(shù)方程來表示

(1)

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

(2)

(3)

根據(jù)A的不穩(wěn)定特征值便可判斷系統(tǒng)小干擾穩(wěn)定性?？紤]到矩陣A將失去系統(tǒng)的稀疏性,因此在特征值計算過程中將采用稀疏技術[23],并不直接計算A。

1.2 Krylov-Schur算法

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

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

(4)

根據(jù)文獻[15]可得到k階Krylov分解

(5)

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

(6)

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

(7)

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

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

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

(8)

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

(9)

(10)

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

(11)

(12)

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

2 多重譜變換方法

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

2.1 常用的譜變換方法

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

S=(A-sIn)-1

(13)

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

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

(14)

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

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

此外,凱萊變換的定義為

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

(15)

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

(16)

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

指數(shù)變換在不穩(wěn)定特征值計算過程中表現(xiàn)出了良好的譜特性,其定義為

(17)

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

μE=eλ和vC=v

(18)

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

2.2 存在的問題

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

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

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

2.3 多重譜變換

根據(jù)位移逆變換構造的多重位移逆變換如下

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

(19)

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

(20)

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

(21)

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

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

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

(22)

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

(23)

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

|eαλ|

(24)

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

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

(25)

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

(26)

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

3 算法流程

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

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

4 算例分析

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

4.1 不同譜變換譜性能對比

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

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

(a)Xingo6u系統(tǒng)

(b)Xingo3012系統(tǒng)

(a)Xingo6u系統(tǒng)

(b)Xingo3012系統(tǒng)

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

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

(a)Xingo6u系統(tǒng)

(a)Xingo6u系統(tǒng)

(b)Xingo3012系統(tǒng)

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

4.2 計算精度和效率測試

在得到多重凱萊矩陣的全部主導特征值及特征向量后,經(jīng)過逆變換便可得到Xingo6u和Xingo3012系統(tǒng)的不穩(wěn)定特征值λ及特征向量v。不穩(wěn)定特征值λ與誤差‖Av-λv‖見表1?？紤]到指數(shù)變換也可得到系統(tǒng)的全部不穩(wěn)定特征值,其結果也在表1中給出。

表1 系統(tǒng)的不穩(wěn)定特征值與誤差

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

(a)Xingo6u系統(tǒng)

(b)Xingo3012系統(tǒng)

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

表2 計算效率測試結果

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

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

5 結 論

文中提出了用以評估電力系統(tǒng)的小干擾穩(wěn)定性的多重譜變換方法,該方法的主要優(yōu)點如下。

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

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

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

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

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