牟倩穎, 葉 華, 劉玉田

(電網(wǎng)智能化調(diào)度與控制教育部重點(diǎn)實(shí)驗(yàn)室(山東大學(xué)), 山東省濟(jì)南市 250061)

0 引言

隨著互聯(lián)電網(wǎng)規(guī)模的不斷增大,區(qū)域間低頻振蕩的問(wèn)題愈加突出[1-2]。基于廣域信息的阻尼控制器可以有效抑制區(qū)間低頻振蕩[3-5]。然而,廣域通信時(shí)滯是影響控制器性能的重要因素。在小干擾穩(wěn)定性分析[6]和廣域阻尼控制器設(shè)計(jì)[7]中需要考慮廣域通信時(shí)滯的影響[8-9]。

在頻域,時(shí)滯被轉(zhuǎn)化為指數(shù)項(xiàng),系統(tǒng)的特征方程成為超越方程。直接求解該方程得到系統(tǒng)的特征值變得異常困難[10]?，F(xiàn)有研究通常采用Padé近似有理多項(xiàng)式來(lái)逼近時(shí)滯環(huán)節(jié)[11]。然而,在大時(shí)滯情況下,Padé近似相位逼近誤差較大。近年來(lái),有學(xué)者提出了基于譜算子(包括無(wú)窮小生成元和解算子)離散化的時(shí)滯電力系統(tǒng)特征值計(jì)算方法,用于準(zhǔn)確求解系統(tǒng)的部分關(guān)鍵特征值。無(wú)窮小生成元離散化方法[9,12-14]能夠計(jì)算給定位移點(diǎn)附近的部分特征值。為了得到系統(tǒng)的全部機(jī)電振蕩模式,需要進(jìn)行多次特征值計(jì)算。解算子離散化(solution operator discretization,SOD)方法的優(yōu)勢(shì)在于通過(guò)一次計(jì)算便可得到系統(tǒng)最右側(cè)或阻尼比最小的部分特征值。文獻(xiàn)[15]提出基于解算子隱式龍格—庫(kù)塔離散化(SOD with implicit Runge-Kutta,SOD-IRK)的時(shí)滯電力系統(tǒng)特征值計(jì)算方法,但是該方法不能直接應(yīng)用于分析大規(guī)模時(shí)滯電力系統(tǒng)。并且,文獻(xiàn)[15]僅考慮了Radau-ⅡA方法求解常微分方程,并未研究其他隱式龍格—庫(kù)塔(IRK)方法的可用性和性能。

鑒于此,本文提出了一類基于SOD-IRK的大規(guī)模時(shí)滯電力系統(tǒng)關(guān)鍵特征值計(jì)算方法。一方面,采用多種IRK離散化方案離散解算子,得到高度稀疏和結(jié)構(gòu)化的解算子離散化矩陣;另一方面,通過(guò)充分利用系統(tǒng)增廣狀態(tài)矩陣和解算子離散化矩陣的稀疏特性,高效地計(jì)算大規(guī)模時(shí)滯電力系統(tǒng)阻尼比最小的部分關(guān)鍵特征值,以避免維數(shù)災(zāi)問(wèn)題。

1 問(wèn)題描述

時(shí)滯電力系統(tǒng)的狀態(tài)方程為[12]:

(1)

(2)

式(1)對(duì)應(yīng)的特征方程為超越方程,即

(3)

式中:λ和v分別為特征值和相應(yīng)的右特征向量。

(4)

(5)

式中:w∈Rl×1為輔助向量;In為單位矩陣。

為了避免直接求解式(3)和式(4),本文首先采用多種IRK方法對(duì)式(1)對(duì)應(yīng)的解算子進(jìn)行離散化,然后計(jì)算解算子離散化矩陣的部分特征值,進(jìn)而得到時(shí)滯電力系統(tǒng)阻尼比最小的部分關(guān)鍵特征值。

2 基于SOD-IRK的特征值計(jì)算方法

2.1 基本原理

設(shè)X:=C([-τmax,0],Rn×1)是由區(qū)間[-τmax,0]到n維實(shí)數(shù)空間Rn×1映射的連續(xù)函數(shù)構(gòu)成的巴拿赫空間,并賦有上確界范數(shù)。式(1)所示時(shí)滯電力系統(tǒng)的解算子T(h):X→X定義為將X上t時(shí)刻初值條件φ映射為t+h時(shí)刻系統(tǒng)狀態(tài)的線性算子,即

(T(h)φ)(t)=Δxh(t)=Δx(t+h)

(6)

式中:h為轉(zhuǎn)移步長(zhǎng),0≤h≤τmax。

在頻域,T(h)的特征值μ與時(shí)滯電力系統(tǒng)的特征值λ之間存在如下指數(shù)關(guān)系[16]:

μ=eλhμ∈σ(T(h)){0}

(7)

式中:σ(·)表示譜;表示集合差運(yùn)算。類似于Cayley變換[17],T(h)將s平面虛軸右側(cè)的特征值λ映射到z平面單位圓之外。

基于SOD-IRK的特征值計(jì)算的原理是,首先對(duì)T(h)進(jìn)行IRK離散化,得到其有限維離散化矩陣TNs,這是SOD-IRK方法的核心和關(guān)鍵;其次,計(jì)算TNs模值最大的部分特征值μ;最后,由式(7)得到時(shí)滯電力系統(tǒng)最右側(cè)的部分關(guān)鍵特征值λ。

T(h)的IRK離散化一般性步驟如下:首先,將時(shí)滯區(qū)間[-τmax,0]分為長(zhǎng)度為h=τmax/N的N個(gè)子區(qū)間;然后,利用p階s級(jí)IRK的s個(gè)橫坐標(biāo)c1,c2,…,cs對(duì)每個(gè)子區(qū)間進(jìn)行離散化,如圖1所示;最后,在集合ΩNs每一個(gè)離散點(diǎn)處,根據(jù)式(6)推導(dǎo)Δxh與φ之間表達(dá)式的系數(shù)矩陣,即得TNs。

(8)

圖1 離散點(diǎn)集合ΩNsFig.1 Discrete point set ΩNs

本文首先詳細(xì)說(shuō)明利用IRK的Radau-ⅡA方法對(duì)T(h)進(jìn)行離散化[18],然后對(duì)T(h)的其他IRK離散化方法進(jìn)行簡(jiǎn)要總結(jié)。

2.2 解算子的Radau-ⅡA離散化方案

1)當(dāng)t∈[-τmax,-h]時(shí),Δxh(t)為初始狀態(tài)φ。利用解算子的轉(zhuǎn)移特性,可以得到t=-jh+cqh(q=1,2,…,s;j=2,3,…,N)處Δx(t+h)的估計(jì)值為:

(9)

令Δx(0),Δx(1)∈RNsn×1分別為定義在區(qū)間[-τmax,0]和[h-τmax,h]上的系統(tǒng)的離散狀態(tài)變量,其表達(dá)式為:

(10)

于是,式(9)寫(xiě)成向量形式:

(11)

2)當(dāng)t∈[-h,0],t+h-τi<0,i=1,2,…,m時(shí),利用p階s級(jí)IRK的Radau-ⅡA法求解如下常微分方程:

(12)

從而,可推導(dǎo)得到Δx(1)與Δx(0)之間的顯式表達(dá)式為:

(13)

式中:RNs∈Rsn×sn,ΣNs∈Rsn×Nsn。

(14)

式中:A∈Rs×s為IRK方法Butcher表(A,b,c)的元素[19]。Radau-ⅡA方法的Butcher表具有如下特點(diǎn):①A可逆;②0

3)綜合式(11)和式(13),可得Δx(1)與Δx(0)之間的關(guān)系式為:

Δx(1)=TNsΔx(0)

(15)

(16)

(17)

式中:TNs∈RNsn×Nsn為解算子T(h)的IRK離散化矩陣。

2.3 解算子的其他IRK離散化方案

解算子的其他IRK離散化方案與Radau-ⅡA離散化方案的不同之處,可總結(jié)為以下兩點(diǎn)。

1)當(dāng)cs≠1時(shí),需要增加離散點(diǎn)cs+1=1以得到每個(gè)子區(qū)間最右端點(diǎn)處系統(tǒng)狀態(tài)的估計(jì)值。

2)當(dāng)c1=0時(shí),每個(gè)子區(qū)間左端點(diǎn)與相鄰子區(qū)間的右端點(diǎn)重合,即θj-1+c1h=θj+csh或θj-1+c1h=θj+cs+1h,其中j=2,3,…,N。因而,需要去掉這些重復(fù)的離散點(diǎn)。

以Lobatto-ⅢA/C離散化方案為例,其系數(shù)(A,b,c)有如下特點(diǎn):①0=c1

(18)

(19)

(20)

類似地,可處理Gauss-Legendre、單對(duì)角隱式龍格—庫(kù)塔法(SDIRK)、Radau-ⅠA和Lobatto-ⅢB等離散化方案(見(jiàn)附錄B)。此外,對(duì)于p階IRK方法,特征值計(jì)算的局部計(jì)算誤差為Ο(hp)[19]。

3 SOD-IRK方法的稀疏實(shí)現(xiàn)

3.1 旋轉(zhuǎn)—放大預(yù)處理

旋轉(zhuǎn)—放大預(yù)處理[20],一方面將時(shí)滯電力系統(tǒng)阻尼比小于ζ(=sinθ)的部分特征值λ乘以旋轉(zhuǎn)因子e-jθ,從而將它們變換到s平面虛軸右側(cè),進(jìn)而利用式(7)將它們映射到z平面單位圓之外;另一方面,將z平面的特征值進(jìn)行α次乘方,以增加它們之間相對(duì)距離。

預(yù)處理后,TNs變?yōu)門(mén)N′s:

(21)

式中:N′=τmax/(αh),其中 表示向上取整。

(22)

(23)

為了保證預(yù)處理后時(shí)滯為實(shí)數(shù),式(23)對(duì)時(shí)滯做了必要的近似,其對(duì)特征值的影響分析可參考文獻(xiàn)[20]。

在旋轉(zhuǎn)—放大預(yù)處理中,參數(shù)α取2或3,即可顯著增加T(h)特征值之間的相對(duì)距離,改善迭代特征值算法的收斂性。此外,α的取值還需要滿足:αh≤τmax。

參數(shù)θ的選取依據(jù)包括:待求部分關(guān)鍵特征值的分布情況和SOD-IRK方法的應(yīng)用目的。首先,希望得到的關(guān)鍵特征值應(yīng)當(dāng)位于等阻尼比線ζ的右側(cè)。其次,如果希望計(jì)算系統(tǒng)較多的關(guān)鍵特征值以設(shè)計(jì)控制器,則θ應(yīng)取較大的值,如θ≥5.74°;如果希望計(jì)算系統(tǒng)最右側(cè)的少量特征值用以快速、可靠地判斷系統(tǒng)的小干擾穩(wěn)定性,則θ應(yīng)取較小的值,如1.72°或2.87°。

3.2 稀疏特征值計(jì)算

采用隱式重啟動(dòng)Arnoldi算法來(lái)計(jì)算TN′s模值最大的r個(gè)特征值μ′。在計(jì)算過(guò)程中,計(jì)算量最大的操作是迭代生成Krylov向量。下面以Radau-ⅡA離散化方案為例,闡述由第j個(gè)Krylov向量qj∈CN′sn×1高效計(jì)算第j+1個(gè)Krylov向量qj+1的方法。

qj+1=TN′sqj

(24)

由于TN′s的次對(duì)角為單位矩陣,qj+1的后(N′-1)sn個(gè)元素可通過(guò)直接將qj的元素向后位移sn得到。前sn個(gè)元素qj+1(1:sn)的計(jì)算可以分解為式(25)和式(26)。

z=ΣN′sqj

(25)

(26)

式中:z∈Csn×1為中間向量。

令qj=vec(Q),其中Q∈Cn×N′s,vec(·)表示將矩陣按列重排為列向量的操作。利用Kronecker積的性質(zhì)[21],式(27)給出了式(25)的高效實(shí)現(xiàn)。

(27)

(28)

2)式(28)可改寫(xiě)為以下兩步:

(29)

在計(jì)算之前,僅對(duì)D0做一次稀疏三角分解。于是,式(29)的計(jì)算僅是一些稀疏矩陣與向量相乘和兩個(gè)三角方程的求解。

3.3 牛頓校驗(yàn)

(30)

4 算例分析

4.1 16機(jī)68節(jié)點(diǎn)系統(tǒng)

16機(jī)68節(jié)點(diǎn)算例系統(tǒng)的詳細(xì)數(shù)據(jù)參見(jiàn)文獻(xiàn)[22]?？紤]在G2和G5上分別安裝廣域阻尼控制器[20],反饋時(shí)滯和輸出時(shí)滯分別為τf1=150 ms,τo1=90.5 ms,τf2=70 ms和τo2=37.4 ms。所有測(cè)試均在8 GB RAM、3.4 GHz CPU臺(tái)式計(jì)算機(jī)上進(jìn)行。SOD-IRK方法和牛頓校驗(yàn)的收斂精度均為10-6。

首先,驗(yàn)證不進(jìn)行旋轉(zhuǎn)—放大預(yù)處理和牛頓校驗(yàn),采用各種IRK離散化方案時(shí),TNs逼近T(h)的準(zhǔn)確性。各參數(shù)分別為h=0.005 s,N=30,α=1和θ=0°。以2級(jí)Radau-ⅡA,Radau-ⅠA,Gauss-Legendre,SDIRK,Lobatto-ⅢA,Lobatto-ⅢB和Lobatto-ⅢC離散化方案為例,SOD-IRK方法分別能夠準(zhǔn)確估計(jì)T(h)模值大于0.744 7,0.861 7,0.727 0,0.777 3,0.937 8,0.903 3,0.938 2和0.903 6的特征值μ,分別對(duì)應(yīng)實(shí)部大于-59.41,-30.16,-65.45,-50.93,-19.95,-20.66,-12.97和-20.32的特征值λ。部分結(jié)果見(jiàn)附錄C。

通過(guò)對(duì)比分析可知,Radau-ⅡA,Radau-ⅠA,Gauss-Legendre和SDIRK離散化方案具有較好的特征值估計(jì)精度。此外,盡管各種IRK方法能夠準(zhǔn)確計(jì)算的特征值數(shù)目不同,這部分準(zhǔn)確特征值都足以分析系統(tǒng)的小干擾穩(wěn)定性。

接著,分析旋轉(zhuǎn)—放大預(yù)處理的影響。在應(yīng)用旋轉(zhuǎn)—放大預(yù)處理(α=2和θ=8.63°)后,各種SOD-IRK方法都具有快速計(jì)算系統(tǒng)全部(r=15)機(jī)電振蕩模式的能力,如圖2所示。

圖2 預(yù)處理后的系統(tǒng)特征值估計(jì)值及其準(zhǔn)確值λFig.2 Estimated value and exact value λ of eigenvalue after pre-processing

最后,對(duì)比分析各種SOD-IRK方法計(jì)算系統(tǒng)全部機(jī)電振蕩模式的效率。由表1可知,所有SOD-IRK方法的計(jì)算時(shí)間都在10 s以內(nèi),且單次迭代時(shí)間均小于基于解算子的偽譜離散化(SOD-PS)方法[20](其中切比雪夫多項(xiàng)式階數(shù)為3)。當(dāng)參數(shù)變化時(shí),隨著解算子離散化矩陣維數(shù)的增大,計(jì)算時(shí)間隨之增加。相比之下,前4種SOD-IRK方法的計(jì)算效率略低。

表1 16機(jī)68節(jié)點(diǎn)系統(tǒng)采用各種SOD方法的效率比較Table 1 Efficiency comparison among various SOD methods in 16-generator 68-bus system

4.2 華北—華中特高壓互聯(lián)電網(wǎng)

在某水平年下華北—華中特高壓互聯(lián)電網(wǎng)上進(jìn)一步驗(yàn)證各種SOD-IRK方法對(duì)大規(guī)模時(shí)滯電力系統(tǒng)的處理能力。系統(tǒng)含有33 028個(gè)節(jié)點(diǎn),2 405臺(tái)同步發(fā)電機(jī),16條高壓直流輸電線路?？紤]在蒙岱海1號(hào)和魯東海7號(hào)上分別裝設(shè)廣域阻尼控制器以抑制區(qū)域間低頻振蕩[13]。反饋和輸出時(shí)滯分別為τf1=120 ms,τo1=100 ms,τf2=100 ms和τo2=80 ms。系統(tǒng)狀態(tài)變量的維數(shù)為n=80 577,代數(shù)變量為l=162 718。

首先,利用各種SOD-IRK方法計(jì)算系統(tǒng)阻尼小于5%的r=20個(gè)特征值。圖3給出了2級(jí)Radau-ⅡA和Lobatto-ⅢA方法的計(jì)算結(jié)果。各參數(shù)分別設(shè)為h=0.006 s,N=20,α=2和θ=2.87°。圖中,部分對(duì)時(shí)滯靈敏度較大的特征值的估計(jì)誤差較大,需要通過(guò)牛頓校驗(yàn)予以消除。此外,系統(tǒng)仍然存在2對(duì)負(fù)阻尼區(qū)間振蕩模式和眾多失穩(wěn)和臨界失穩(wěn)局部模式。通過(guò)對(duì)比分析可知,各種SOD-IRK方法的精度基本相當(dāng)。

表2給出了不同SOD-IRK方法和SOD-PS方法計(jì)算20個(gè)特征值的IRA迭代次數(shù)和計(jì)算時(shí)間。通過(guò)對(duì)比分析可知:①SOD-IRK的單次迭代計(jì)算時(shí)間遠(yuǎn)小于SOD-PS;②Lobatto-ⅢA方法的計(jì)算效率最高(如表2中紅色數(shù)據(jù)所示)。

圖3 華北—華中特高壓電網(wǎng)特征值估計(jì)值及其準(zhǔn)確值λFig.3 Estimated value and exact value λ of eigenvalue in North China-Central China ultra-high voltage power grid

表2 華北—華中特高壓電網(wǎng)采用各種SOD方法的效率比較Table 2 Efficiency comparison among various SOD methods in North China-Central ultra-high voltage China power grid

5 結(jié)語(yǔ)

本文提出了一類基于SOD-IRK的大規(guī)模時(shí)滯電力系統(tǒng)特征值的高效計(jì)算方法。

2)Radau-ⅡA,Radau-ⅠA,Gauss-Legendre和SDIRK方法具有較好的計(jì)算精度,而Lobatto-ⅢA和Radau-ⅡA方法具有較高的計(jì)算效率。

本文提出的基于SOD-IRK的大規(guī)模時(shí)滯電力系統(tǒng)特征值計(jì)算方法中,系統(tǒng)所有的狀態(tài)變量都需要離散化,導(dǎo)致解算子離散化矩陣維數(shù)較大,特征值計(jì)算時(shí)間長(zhǎng)。因此,考慮只對(duì)與時(shí)滯相關(guān)的狀態(tài)變量進(jìn)行離散化以高效求解時(shí)滯電力系統(tǒng)的特征值是后續(xù)的研究方向。

本文受到山東大學(xué)青年學(xué)者未來(lái)計(jì)劃項(xiàng)目(2016WLJH06)支持,特此感謝!

附錄見(jiàn)本刊網(wǎng)絡(luò)版(http://www.aeps-info.com/aeps/ch/index.aspx)。