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

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

0 引言

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

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

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

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

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

1.1 數(shù)學模型

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

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

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

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

需要說明的是,狀態(tài)矩陣A失去了稀疏特性。在實際計算中,關(guān)于狀態(tài)矩陣A的所有運算都采用稀疏技術(shù)［13］,而不直接計算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ū)域。除此之外,另一個區(qū)域為位于虛軸右側(cè)的不穩(wěn)定區(qū)域。在實際電力系統(tǒng)中,控制器參數(shù)調(diào)節(jié)不恰當可能會導(dǎo)致電力系統(tǒng)出現(xiàn)不穩(wěn)定特征值。而遺漏其中任何一個不穩(wěn)定特征值都將導(dǎo)致對電力系統(tǒng)的小干擾穩(wěn)定性產(chǎn)生錯誤的判斷。因此,本文提出了一種雙層Arnoldi 算法來求解電力系統(tǒng)的全部不穩(wěn)定特征值。

2 外層隱式重啟動Arnoldi 算法

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

2.1 基本的Arnoldi 方法

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

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

2.2 指數(shù)變換

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

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

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

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

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

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

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

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

2.3 重啟動策略

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

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

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

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

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

2.3.2 重啟動策略

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

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

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

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

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

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

3.3 求解eAv 的算法流程

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

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

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

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

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

4 算法實現(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的計算利用內(nèi)層Arnoldi 算法實現(xiàn)。

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

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

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

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

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

5 算例分析

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

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

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

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

表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é)果與精確解的對比如圖4 所示?？梢钥闯?指數(shù)映射將所有不穩(wěn)定特征值映射為模值較大的主導(dǎo)特征值,且不同主導(dǎo)特征值之間的間距較遠。因此,利用所提雙層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為位移點。在本測試中,采用位移點為s=1+j10 的復(fù)位移-逆變換以及位移點為s1=-1和s2=1 的實Cayley 變換。需要說明的是,對于復(fù)位移-逆變換來說,只需要關(guān)注具有非負虛部的不穩(wěn)定特征值即可,即圖5（a）中的μ1、μ2～μ6、μ12、μ14。

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

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

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

5.2 效率測試

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

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

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

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

6 結(jié)語

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

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

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

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

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