基于指數(shù)變換的電力系統(tǒng)不穩(wěn)定特征值計(jì)算方法

陳亦平，曹玉磊，趙利剛，肖 亮，李崇濤

（1. 中國(guó)南方電網(wǎng)電力調(diào)度控制中心,廣東省廣州市 510663；2. 中國(guó)南方電網(wǎng)科學(xué)研究院有限責(zé)任公司,廣東省廣州市 510663）

0 引言

特征分析法是電力系統(tǒng)小干擾穩(wěn)定性分析過程中的一個(gè)有效方法［1］。其中,QR/BR 等全特征值計(jì)算方法具有數(shù)值穩(wěn)定性好、收斂速度快且不漏根的優(yōu)勢(shì)［2］。然而,隨著現(xiàn)代電力系統(tǒng)規(guī)模的不斷擴(kuò)大,系統(tǒng)狀態(tài)變量的數(shù)量可以達(dá)到幾千階甚至數(shù)萬階,全特征值計(jì)算方法已經(jīng)不能滿足大規(guī)模系統(tǒng)特征分析的需要。在小干擾穩(wěn)定分析過程中,通常只有部分關(guān)鍵特征值是所關(guān)心的。因此,部分特征值計(jì)算方法成為進(jìn)行大規(guī)模電力系統(tǒng)小干擾穩(wěn)定性分析的切實(shí)可行方法。其中,電力系統(tǒng)小干擾穩(wěn)定分析中關(guān)鍵特征值通常可以分為兩類［3］,即弱阻尼特征值和不穩(wěn)定特征值。

目前,部分特征值計(jì)算方法包括序貫法、子空間法［4］。其中,利用子空間法可以得到系統(tǒng)的一組模值最大的特征值。子空間方法中最具有代表性的是基于Krylov 子空間的隱式重啟動(dòng)Arnoldi（implicitly restarted Arnoldi,IRA）方法［5-7］、基于非Krylov 子空間的Jacobi-Davidson（JD）方法［8］和子空間加速瑞利商迭代方法［9］。

IRA 方法具有算法適應(yīng)性強(qiáng)、收斂速度快的優(yōu)點(diǎn),是求解電力系統(tǒng)部分特征值最成功的方法之一,已經(jīng)被國(guó)內(nèi)外多種電力系統(tǒng)特征值計(jì)算程序和小干擾分析程序采用,如PEALS、PSAT 和PSASP等［10］。該方法的主要特征是可以優(yōu)先收斂到模值最大的特征值,在實(shí)際應(yīng)用中需要配合位移-逆變換或者Cayley 變換將關(guān)鍵特征值映射為模值最大的主導(dǎo)特征值［11］。因此,通過在弱阻尼區(qū)域設(shè)置一系列的位移點(diǎn),便可以得到幾乎所有的弱阻尼特征值。但位移-逆變換或者Cayley 變換在求解不穩(wěn)定特征值時(shí)性能并不好［12］。這是因?yàn)樘卣髦档闹鲗?dǎo)性隨著離位移點(diǎn)的距離增加而減弱,而不穩(wěn)定模態(tài)在右半平面上的分布無法提前估計(jì)。

作為以上特征值計(jì)算方法的補(bǔ)充,本文提出一種基于指數(shù)變換的雙層Arnoldi 算法來求解大規(guī)模電力系統(tǒng)的不穩(wěn)定特征值。首先,該方法采用了一種指數(shù)變換的譜變換方式。由于指數(shù)變換在求解不穩(wěn)定特征值時(shí)具有良好的譜主導(dǎo)性和譜稀疏性,因此所提方法可以不遺漏地求解電力系統(tǒng)的全部不穩(wěn)定特征值。此外,本文采用了內(nèi)層Arnoldi 算法隱式地求解指數(shù)矩陣與向量的乘積,避免了顯式計(jì)算指數(shù)矩陣,使得所提方法在保證計(jì)算精度的前提下具有較高的計(jì)算效率,從而可以滿足在大規(guī)模實(shí)際系統(tǒng)中應(yīng)用的需要。

1 小干擾穩(wěn)定性分析基礎(chǔ)

1.1 數(shù)學(xué)模型

電力系統(tǒng)的小干擾穩(wěn)定性分析的數(shù)學(xué)模型通?？梢杂萌缦挛⒎执鷶?shù)方程來描述:

式中:A?∈Rn×n、B?∈Rn×m、C?∈Rm×n、D?∈Rm×m為稀疏矩陣,其中n為狀態(tài)變量Δx的個(gè)數(shù),m為代數(shù)變量Δy的個(gè)數(shù)。

消去式（1）中代數(shù)變量Δy,得到:

式中:A∈Rn×n為系統(tǒng)的狀態(tài)矩陣,其具體形式如式（3）所示。

需要說明的是,狀態(tài)矩陣A失去了稀疏特性。在實(shí)際計(jì)算中,關(guān)于狀態(tài)矩陣A的所有運(yùn)算都采用稀疏技術(shù)［13］,而不直接計(jì)算A。在本文中,為了敘述方便,算法中仍然采用A來描述。

1.2 關(guān)鍵特征值分類

電力系統(tǒng)小干擾穩(wěn)定分析中關(guān)鍵特征值所在區(qū)域通?？梢苑譃槿踝枘釁^(qū)域和不穩(wěn)定區(qū)域［3］,如圖1所示。

圖1 關(guān)鍵特征值所在區(qū)域分類Fig.1 Region classification of critical eigenvalues

從圖1 可以看出,弱阻尼區(qū)域是由關(guān)鍵阻尼比ζc、虛軸以及最大振蕩頻率fm所決定的窄三角區(qū)域。除此之外,另一個(gè)區(qū)域?yàn)槲挥谔撦S右側(cè)的不穩(wěn)定區(qū)域。在實(shí)際電力系統(tǒng)中,控制器參數(shù)調(diào)節(jié)不恰當(dāng)可能會(huì)導(dǎo)致電力系統(tǒng)出現(xiàn)不穩(wěn)定特征值。而遺漏其中任何一個(gè)不穩(wěn)定特征值都將導(dǎo)致對(duì)電力系統(tǒng)的小干擾穩(wěn)定性產(chǎn)生錯(cuò)誤的判斷。因此,本文提出了一種雙層Arnoldi 算法來求解電力系統(tǒng)的全部不穩(wěn)定特征值。

2 外層隱式重啟動(dòng)Arnoldi 算法

本章首先回顧了基本的Arnoldi 方法；然后,給出了指數(shù)變換的定義；最后,基于指數(shù)變換給出了計(jì)算電力系統(tǒng)不穩(wěn)定特征值的外層Arnoldi 算法和相應(yīng)的隱式重啟動(dòng)策略。

2.1 基本的Arnoldi 方法

給定一個(gè)矩陣M∈Rn×n和一個(gè)初始非零向量v∈Rn,則k階的Krylov 子空間可以表示為:

矩陣M關(guān)于向量v的k階Arnoldi 分解具有如下形式:

2.2 指數(shù)變換

指數(shù)變換的定義如下:

式中:In∈Rn×n為單位矩陣。假設(shè)A是可對(duì)角化的,那么可以得到狀態(tài)矩陣A和指數(shù)矩陣eA的特征值對(duì)之間的對(duì)應(yīng)關(guān)系。如果(λ,u)是A的一個(gè)特征值對(duì),那么(μ,u)將是eA的一個(gè)特征值對(duì),且滿足:

相反,如果(μ,u)是eA的一個(gè)特征值對(duì),那么(λ,u)是A的一個(gè)特征值對(duì),且滿足:

式中:N為參數(shù),是一個(gè)具體的正整數(shù)。

從譜變換的角度來看,指數(shù)變換將狀態(tài)矩陣A的不穩(wěn)定特征值映射為指數(shù)矩陣eA的模值大于1 的主導(dǎo)特征值,如圖2 所示。Re(λ2)

圖2 指數(shù)變換示意圖Fig.2 Schematic diagram of exponential transformation

作為對(duì)比,Cayley 變換的示意圖如圖3 所示?？梢钥闯?Cayley 變換也可以將不穩(wěn)定特征值映射為模值大于1 的特征值,但這種變換降低了映射后特征值之間的距離,使得不穩(wěn)定特征值收斂較為困難［12］。

圖3 Cayley 變換示意圖Fig.3 Schematic diagram of Cayley transformation

2.3 重啟動(dòng)策略

Krylov 子空間方法的特征是其可以優(yōu)先收斂到模值最大的特征值。在外層Arnoldi 算法中,令

也就是對(duì)指數(shù)矩陣eA進(jìn)行Arnoldi 分解,然后便可以得到其主導(dǎo)特征值。逆變換后便可以得到狀態(tài)矩陣A的不穩(wěn)定特征值。

2.3.1 逆變換求解參數(shù)N

假設(shè)(μi,ui)是M的一個(gè)收斂的特征對(duì),那么需要計(jì)算式（9）中的N值來得到對(duì)應(yīng)的特征值λi。考慮到特征值收斂時(shí)誤差ri的值已經(jīng)非常小,根據(jù)式（9）可以得到:

式中:Round(·)為取整函數(shù)；H 表示共軛轉(zhuǎn)置。得到N后,λi=lnμi+j2πN,λi和ui分別為狀態(tài)矩陣A的特征值和特征向量。

2.3.2 重啟動(dòng)策略

當(dāng)?shù)玫骄仃嘙關(guān)于式（5）的Arnoldi 分解后,如果此時(shí)主導(dǎo)特征值并沒有全部收斂,就需要采用重啟動(dòng)策略來提高特征值的收斂性。對(duì)Hk的特征值按模值遞減排序,即

此外,在利用外層隱式重啟動(dòng)Arnoldi 算法求解狀態(tài)矩陣A的不穩(wěn)定特征值的過程中,需要計(jì)算指數(shù)矩陣eA與向量v的乘積?？紤]到直接求解大規(guī)模指數(shù)矩陣eA是不現(xiàn)實(shí)的,下面采用內(nèi)層Arnoldi 算法來隱式求解eAv。

3 內(nèi)層采用加速收斂技術(shù)的Arnoldi 算法

3.1 內(nèi)層Arnoldi 算法求解eAv

本節(jié)利用Arnoldi 算法隱式地求解eAv的值［14］。給定狀態(tài)矩陣A和向量v（其中||v||2=1）,則k階的Krylov 子空間可以表示為:

3.2 加速收斂技術(shù)

實(shí)際的電力系統(tǒng)通常為剛性系統(tǒng),也就是同時(shí)存在模值很大和模值很小的特征值,這會(huì)導(dǎo)致內(nèi)層Arnoldi 算法的收斂速度很慢［16］。因此,為了提高Arnoldi 算法的收斂速度,本節(jié)采用基于位移逆矩陣S的加速收斂技術(shù)［15-16］。位移逆矩陣S定義為:

式中:τ為位移系數(shù),通常取值為很小的正實(shí)數(shù),如0.02。

相應(yīng)的Arnoldi 分解具有如下形式:

如果誤差rk小于收斂容差εinner,那么VkeGke1可以被認(rèn)為是eAv的一個(gè)很好的近似解。需要說明的是,eGk是一個(gè)低階指數(shù)矩陣,可通過文獻(xiàn)［17］中的算法高效求解。

3.3 求解eAv 的算法流程

利用內(nèi)層Arnoldi 算法求解eAv流程如下:

步 驟1:給 定A、v、εinner、kinner,其 中kinner為 內(nèi) 層Arnoldi 算法的最大擴(kuò)展維數(shù),令M=S。

步驟2:令k=1,2,…,kinner,根據(jù)式（25）得到k階Arnoldi 分 解,根 據(jù) 式（31）計(jì) 算 誤 差rk,若rk<εinner,結(jié)束迭代。

步驟3:根據(jù)式（29）計(jì)算得到eAv。

在步驟2 的Arnoldi 分解過程中,需要計(jì)算Sv的值,為了減小計(jì)算量,這里采用稀疏技術(shù)而不直接計(jì)算狀態(tài)矩陣A［13］。此外,考慮到Arnoldi 分解的階數(shù)k較小時(shí)誤差rk通常較大,為了減小計(jì)算量,可以當(dāng)k較大（k>40）時(shí)再計(jì)算誤差rk。

4 算法實(shí)現(xiàn)

利用雙層Arnoldi 算法求解電力系統(tǒng)不穩(wěn)定特征值的算法流程如下:

步 驟1:給 定A、v、kouter、εouter,其 中kouter為 外 層Arnoldi 算 法 的 維 數(shù),εouter需 要 設(shè) 定 的 比εinner小。令v←v/||v||2。

步驟2:形成矩陣eA的kouter階Arnoldi 分解,其中eAv的計(jì)算利用內(nèi)層Arnoldi 算法實(shí)現(xiàn)。

步驟3:計(jì)算得到Hk的所有特征值對(duì)(μi,zi),i=1,2,…,kouter,然后對(duì)特征值μi按式（13）進(jìn)行排序。

步驟4:令q為模值大于1 的特征值μi的數(shù)量,根據(jù)式（6）判斷前q個(gè)特征值μi是否收斂。如果收斂,則結(jié)束迭代,轉(zhuǎn)向步驟6。

步驟5:令前p（p>q）個(gè)特征值為期望特征值,根據(jù)式（16）得到前p階Arnoldi 關(guān)系,然后擴(kuò)展到kouter階Arnoldi 分解,轉(zhuǎn)向步驟3。

步驟6:根據(jù)式（9）采用逆變換,得到狀態(tài)矩陣A的所有不穩(wěn)定特征值。

考慮到實(shí)際系統(tǒng)中由于控制器參數(shù)設(shè)置不合理等因素產(chǎn)生的不穩(wěn)定特征值數(shù)目較少,kouter設(shè)置為60 就已足夠。重啟動(dòng)過程中期望特征值的個(gè)數(shù)p需要設(shè)定的比q略大。

5 算例分析

本文利用Xing6u 和NF2016 這2 個(gè)實(shí)際系統(tǒng)來驗(yàn)證所提方法的有效性。其中,Xing6u 是開源的系統(tǒng)［9］,NF2016 是 以 中 國(guó) 某 區(qū) 域 電 網(wǎng)2016 年 運(yùn) 行 方式為基礎(chǔ)的實(shí)際系統(tǒng)。Xing6u 系統(tǒng)有2 940 個(gè)狀態(tài)變量、17 798 個(gè)代數(shù)變量、73 946 個(gè)非零元素。NF2016 系 統(tǒng) 有40 819 個(gè) 狀 態(tài) 變 量、71 656 個(gè) 代 數(shù) 變量、486 605 個(gè)非零元素。本文所進(jìn)行的測(cè)試都是在一臺(tái)裝有Intel i7 CPU、內(nèi)存為16 GB 的臺(tái)式機(jī)上進(jìn)行的,并且所有測(cè)試都用C++代碼編寫。

5.1 不穩(wěn)定特征值計(jì)算

首先,在Xing6u 和NF2016 這2 個(gè)系統(tǒng)上測(cè)試內(nèi)層Arnoldi 算法的有效性。任取初始向量v,隨著內(nèi)層Arnoldi 算法維數(shù)的擴(kuò)展,根據(jù)式（31）得到計(jì)算eAv時(shí)的誤差。當(dāng)eAv的計(jì)算精度為10-8時(shí),Arnoldi分解的維數(shù)通常小于60,即內(nèi)層Arnoldi 算法可以快速且準(zhǔn)確地求解eAv。因此,收斂容差εinner可以設(shè)置為10-8。對(duì)于極少數(shù)初始向量v,當(dāng)擴(kuò)展維度小于80 時(shí),內(nèi)層Arnoldi 算法可能收斂不到10-8。此時(shí),為了提高計(jì)算效率,可以選取最小誤差時(shí)對(duì)應(yīng)的維度來求解eAv。此外,為了減小計(jì)算量,可以當(dāng)內(nèi)層Arnoldi 算法的維數(shù)k較大（k>40）時(shí)再計(jì)算誤差rk。

接下來,利用雙層Arnoldi 算法計(jì)算Xing6u 系統(tǒng)的不穩(wěn)定特征值。外層Arnoldi 算法的維數(shù)kouter設(shè)置為60,εouter設(shè)置為10-6。經(jīng)過2 次重啟動(dòng)后,算法2 的迭代終止。指數(shù)矩陣eA模值大于1 的主導(dǎo)特征值μ1～μ15如表1 所示。經(jīng)過逆變換后,狀態(tài)矩陣A的不穩(wěn)定特征值如表2 所示。由表2 可以看出,該系統(tǒng)的15 個(gè)不穩(wěn)定特征值λ1～λ15均被得到,包括2.94、0.08±j4.32、0.11±j4.95 以 及5 對(duì) 重 特 征 值1.01±j8.08。與QR 方法得到精確解的對(duì)比顯示,所提方法具有較高的計(jì)算精度。

表1 指數(shù)矩陣的主導(dǎo)特征值Table 1 Dominant eigenvalues of exponential matrix

表2 狀態(tài)矩陣A 的不穩(wěn)定特征值Table 2 Unstable eigenvalues of state matrix A

所得結(jié)果與精確解的對(duì)比如圖4 所示。可以看出,指數(shù)映射將所有不穩(wěn)定特征值映射為模值較大的主導(dǎo)特征值,且不同主導(dǎo)特征值之間的間距較遠(yuǎn)。因此,利用所提雙層Arnoldi 算法非常容易得到系統(tǒng)所有的不穩(wěn)定特征值。

圖4 Xing6u 系統(tǒng)指數(shù)矩陣和狀態(tài)矩陣的特征值分布Fig.4 Distribution of eigenvalues for exponential matrix and state matrix in Xing6u system

最后,為了直觀地展示不穩(wěn)定特征值經(jīng)過位移-逆變換和Cayley 變換后的分布情況,Xing6u 系統(tǒng)的特征值經(jīng)過相應(yīng)變換后的結(jié)果如圖5 所示。

圖5 Xing6u 系統(tǒng)經(jīng)位移-逆變換和Cayley 變換后的特征值分布Fig.5 Distribution of eigenvalues after shift-invert transformation and Cayley transformation in Xing6u system

位移-逆變換和Cayley 變換的形式分別為A→(A-sI)-1和A→(A-s1I)(A-s2I)-1,其中,s、s1和s2為位移點(diǎn)。在本測(cè)試中,采用位移點(diǎn)為s=1+j10 的復(fù)位移-逆變換以及位移點(diǎn)為s1=-1和s2=1 的實(shí)Cayley 變換。需要說明的是,對(duì)于復(fù)位移-逆變換來說,只需要關(guān)注具有非負(fù)虛部的不穩(wěn)定特征值即可,即圖5（a）中的μ1、μ2～μ6、μ12、μ14。

從圖5（a）可以看出,在位移-逆變換譜上,靠近位移點(diǎn)s=1+j10 的5 個(gè)重特征值成為模值最大的主導(dǎo)特征值。然而,遠(yuǎn)離該位移點(diǎn)的特征值μ1、μ12、μ14的模值非常小。與位移-逆變換不同,Cayley 變換可以將所有的不穩(wěn)定特征值映射為模值大于1 的主導(dǎo)特征值。但從圖5（b）可以看出,大多數(shù)不穩(wěn)定特征值尤其是5 對(duì)重特征值與其他特征值聚集在一起,使得IRA 方法收斂的難度較大。

對(duì)Cayley 變換來說,當(dāng)位移點(diǎn)選擇合適時(shí)也可以得到系統(tǒng)的所有不穩(wěn)定特征值。但由于事先并不能獲取系統(tǒng)不穩(wěn)定特征值的數(shù)量以及位置分布情況,需要不斷改變位移點(diǎn)的位置,直到可以收斂到系統(tǒng)的全部不穩(wěn)定特征值。依次選取位移點(diǎn)(s1,s2)為(-0.5,0.5)、 (-1,1)、 (-1.5,1.5)、 (-2,2)、 (-2.5,2.5)和(-3,3),最大重啟動(dòng)次數(shù)設(shè)置為50,子空間維數(shù)設(shè)置為60,收斂容差設(shè)置為10-6。當(dāng)位移點(diǎn)(s1,s2)選擇為(-2.5,2.5)和(-3,3)時(shí),分別經(jīng)過43次和40 次重啟動(dòng)收斂到所有的不穩(wěn)定特征值。而對(duì)于其他4 組位移點(diǎn),當(dāng)達(dá)到最大重啟動(dòng)次數(shù)時(shí),仍存在遺漏不穩(wěn)定特征值的情況。

與圖4 相比,在求解不穩(wěn)定特征值時(shí),指數(shù)變換在譜主導(dǎo)性和譜稀疏性方面都具有較大的優(yōu)勢(shì)。這是因?yàn)橹笖?shù)變換可以將所有不穩(wěn)定特征值映射為主導(dǎo)特征值,且主導(dǎo)性隨著特征值實(shí)部的增大而增大。當(dāng)子空間維數(shù)設(shè)置為60 時(shí),經(jīng)過2 次重啟動(dòng),所提雙層Arnoldi 方法就可以收斂到Xing6u 系統(tǒng)所有的不穩(wěn)定特征值。

5.2 效率測(cè)試

在本節(jié)中,利用Xing6u 和NF2016 這2 個(gè)實(shí)際系統(tǒng)來測(cè)試所提方法的計(jì)算效率。需要說明的是,NF2016 系統(tǒng)是正常運(yùn)行的系統(tǒng),并不包含不穩(wěn)定特征值。因此,在測(cè)試過程中通過修改控制器參數(shù)使其變?yōu)椴环€(wěn)定系統(tǒng),并包含3 對(duì)不穩(wěn)定特征值,分別 為 0.447 8±j9.596 3、0.315 6±j10.081 7 和0.134 7±j13.542 2。當(dāng)維數(shù)kouter設(shè)置為60、收斂容差εouter設(shè)置為10-6時(shí),利用所提雙層Arnoldi 算法,未經(jīng)過重啟動(dòng)過程就可以收斂到全部的特征值。求解Xing6u 和NF2016 系統(tǒng)不穩(wěn)定特征值的計(jì)算耗時(shí)如表3 所示。其中,NUnstable為不穩(wěn)定特征值的個(gè)數(shù),NIR為重啟動(dòng)的次數(shù)。

表3 Xing6u 和NF2016 系統(tǒng)的計(jì)算耗時(shí)Table 3 Calculation time of Xing6u and NF2016 systems

由表3 可知,雙層Arnoldi 算法求解Xing6u 系統(tǒng)不穩(wěn)定特征值的計(jì)算耗時(shí)為4.2 s,而QR 方法的計(jì)算耗時(shí)為13 s。對(duì)于更大規(guī)模的NF2016 系統(tǒng)來說,計(jì)算耗時(shí)為15.7 s。此外,由于NF2016 系統(tǒng)規(guī)模較大,基于本文所采用的測(cè)試環(huán)境無法用QR 方法求解全部特征值。可以看出,所提雙層Arnoldi 算法可以適用于較大規(guī)模的系統(tǒng)。

此外,當(dāng)位移點(diǎn)選擇合適時(shí),Cayley 變換也可以較為快速地得到系統(tǒng)的所有不穩(wěn)定特征值。例如,當(dāng)位移點(diǎn)選擇為s1=-3 和s2=3 時(shí),Cayley 變換求解Xing6u 系統(tǒng)全部不穩(wěn)定特征值的耗時(shí)為1.3 s。但對(duì)任意一個(gè)系統(tǒng)來說,由于事先不知道不穩(wěn)定特征值的分布,并不能保證一開始選取的位移點(diǎn)就是合適的位移點(diǎn)。如果用于試錯(cuò)的位移點(diǎn)數(shù)目較多,Cayley 變換的計(jì)算耗時(shí)則可能會(huì)超過所提方法。另一方面,對(duì)于大規(guī)模系統(tǒng)來說,所提方法的計(jì)算效率已經(jīng)可以滿足實(shí)際應(yīng)用的需要?？紤]到所采用的方法在譜主導(dǎo)性和譜稀疏性上具有較大的優(yōu)勢(shì),其可以作為現(xiàn)有特征值計(jì)算方法的補(bǔ)充,以確保在電力系統(tǒng)部分特征值計(jì)算過程中能夠不遺漏不穩(wěn)定特征值。

6 結(jié)語

本文提出了一種基于指數(shù)變換的雙層Arnoldi算法來求解大規(guī)模電力系統(tǒng)的不穩(wěn)定特征值。該方法的主要優(yōu)勢(shì)在于:

1）可以將電力系統(tǒng)所有的不穩(wěn)定特征值映射為主導(dǎo)特征值,并且特征值的主導(dǎo)性隨實(shí)部的增大而增大。

2）結(jié)合指數(shù)變換的譜主導(dǎo)性和譜稀疏性的優(yōu)勢(shì),外層Arnoldi 算法易于收斂到所有不穩(wěn)定特征值,且內(nèi)層Arnoldi 算法可以快速地計(jì)算其中所需的指數(shù)矩陣和向量的乘積。因此,所提方法具有較高的計(jì)算效率。

3）利用Xing6u 和NF2016 系統(tǒng)的數(shù)值測(cè)試表明,所提方法具有較高的計(jì)算精度。

此外,雙層Arnoldi 方法是針對(duì)求解不穩(wěn)定特征值提出的,而指數(shù)變換的應(yīng)用使得該方法不適用于求解弱阻尼特征值。這是因?yàn)橹笖?shù)變換將所有的弱阻尼特征值映射為非主導(dǎo)特征值?？紤]到位移-逆變換或者Cayley 變換在求解弱阻尼特征值時(shí)是非常高效的,可以針對(duì)不同的區(qū)域采用不同的譜變換方法。例如,在弱阻尼區(qū)域采用位移-逆變換或者Cayley 變換,而在不穩(wěn)定區(qū)域采用指數(shù)變換。綜上,所提方法可以作為現(xiàn)有特征值計(jì)算方法的補(bǔ)充,以確保在大規(guī)模電力系統(tǒng)部分特征值計(jì)算過程中能夠不遺漏不穩(wěn)定特征值。需要說明的是,研究更高效的指數(shù)矩陣與向量乘積計(jì)算方法對(duì)于所提方法計(jì)算效率的提升具有顯著意義,這也是需進(jìn)一步探索的內(nèi)容。