基于斯蒂芬遜迭代的剛性電力系統(tǒng)時(shí)域計(jì)算交替求解法

黃欽雄，張俊勃，劉 洋

（華南理工大學(xué)電力學(xué)院,廣東省廣州市 510641）

0 引言

電力系統(tǒng)時(shí)域計(jì)算的數(shù)學(xué)模型是一組非線性的微分代數(shù)方程組（differential algebra equations,DAEs）［1］,常用的求解方法為交替求解法［2-3］。在新型電力系統(tǒng)背景下,大量可再生能源與電力電子設(shè)備接入電網(wǎng),電力系統(tǒng)動態(tài)過程的時(shí)間常數(shù)差異增大,交替求解面臨剛性問題,導(dǎo)致計(jì)算過程的收斂性與計(jì)算效率變差［4-5］。因此,如何提高交替求解法收斂性,使算法能夠適應(yīng)系統(tǒng)動態(tài)過程的強(qiáng)剛性,是新型電力系統(tǒng)背景下時(shí)域計(jì)算的關(guān)鍵問題。

影響交替求解法收斂性的因素包括迭代初值、積分步長以及DAEs 的預(yù)處理格式等。在迭代初值設(shè)置方面,可通過顯式積分、人工神經(jīng)網(wǎng)絡(luò)等方法預(yù)估仿真變量初值,改善收斂性［6-7］。在積分步長方面,可基于方差和均值比最大值、局部截?cái)嗾`差等指標(biāo)進(jìn)行收斂性判別,實(shí)現(xiàn)自適應(yīng)步長調(diào)整［7-10］。在DAEs 預(yù)處理方面,已有研究從發(fā)電機(jī)注入電流表達(dá)式中分離出不同形式的虛擬導(dǎo)納,并將其并入網(wǎng)絡(luò)節(jié)點(diǎn)導(dǎo)納矩陣,以提高網(wǎng)絡(luò)代數(shù)方程求解收斂性［1,11-12］,或通過在差分方程中引入差分方程對應(yīng)的雅可比矩陣,以提高差分方程求解收斂性［13］。上述方法屬于改善計(jì)算收斂性的通用方法,并不能改善系統(tǒng)的剛性問題。此外,DAEs 的預(yù)處理方法還存在不能兼顧算法收斂性與計(jì)算效率的問題。

剛性問題主要有適應(yīng)系統(tǒng)剛性或弱化系統(tǒng)剛性兩種解決思路。前者包括:1）采用數(shù)值穩(wěn)定域與階數(shù)更高的數(shù)值積分方法,如龍格-庫塔法、Gear 法、泰勒級數(shù)法［14-16］等；2）縮小仿真步長或采用變步長方法以滿足系統(tǒng)中動態(tài)環(huán)節(jié)的最小時(shí)間常數(shù)限制［17-18］。后者主要將系統(tǒng)分為剛性與非剛性部分,不同部分采用不同的積分方法進(jìn)行求解,相關(guān)研究集中在如何對剛性與非剛性部分進(jìn)行分組以及如何選取不同組合算法［19-25］。各種數(shù)值積分或縮小步長算法能夠適應(yīng)系統(tǒng)的剛性,但也存在積分格式復(fù)雜和積分區(qū)間增加的問題,限制了計(jì)算效率。同時(shí),對系統(tǒng)進(jìn)行剛性與非剛性的區(qū)分能夠有效提高計(jì)算效率,但對不同仿真系統(tǒng)在計(jì)算之前都要進(jìn)行系統(tǒng)劃分,缺乏通用性,且劃分不當(dāng)時(shí)剛性部分仍存在迭代不收斂問題。

從不動點(diǎn)迭代理論出發(fā),上述提高收斂性以及解決剛性問題的思路均可以理解為對迭代過程中迭代函數(shù)的導(dǎo)數(shù)矩陣進(jìn)行優(yōu)化,但這些做法都不改變迭代過程的收斂階。若能對交替求解過程的迭代格式進(jìn)行改進(jìn),在DAEs 有解的情況下提高迭代過程的收斂階,擴(kuò)大各仿真時(shí)步的收斂范圍,則可能大幅提高迭代過程的收斂速度。

本文基于不動點(diǎn)迭代理論,分析交替求解法的迭代求解過程,在交替求解法流程中引入具有二階收斂速度的斯蒂芬遜迭代,使交替求解法能夠適應(yīng)新型電力系統(tǒng)強(qiáng)剛性導(dǎo)致的時(shí)域計(jì)算收斂性惡化。

1 交替求解法的收斂性分析

1.1 DAEs 的交替求解過程

電力系統(tǒng)機(jī)電暫態(tài)仿真的DAEs 可表示為:

式中:x和u分別為狀態(tài)變量和代數(shù)變量；Y為網(wǎng)絡(luò)節(jié)點(diǎn)導(dǎo)納矩陣；f為描述系統(tǒng)狀態(tài)變量動態(tài)元件的微分方程；i為網(wǎng)絡(luò)注入電流方程。

采用隱式梯形法將式（1）的微分方程離散化為差分方程后,計(jì)算n+1 時(shí)刻系統(tǒng)的狀態(tài)即求解如下非線性代數(shù)方程組:

式中:下標(biāo)n和n+1 分別表示仿真的第n和n+1 時(shí)刻；h為積分步長。

采用原始交替求解法對式（2）進(jìn)行迭代求解時(shí),迭代格式如下［3］:

式中:上標(biāo)中k和k+1 分別表示第k和k+1 次迭代。

式（3）中差分方程和代數(shù)方程的交替求解可在牛頓法下進(jìn)行理解。采用牛頓法進(jìn)行n+1 時(shí)刻第k次差分方程與代數(shù)方程求解,可構(gòu)造狀態(tài)變量和代數(shù)變量的殘差方程q、g如下:

根據(jù)殘差方程可分別構(gòu)造狀態(tài)變量x和代數(shù)變量u的修正方程。在差分方程和代數(shù)方程的交替求解中,分別將u和x認(rèn)為是已知量。因此,狀態(tài)變量x的修正量Δx只與x有關(guān),代數(shù)變量u的修正量Δu只與u有關(guān),即

其中,矩陣A、D為雅可比矩陣,其表達(dá)式為:

式中:I為單位矩陣。矩陣A、D的具體表達(dá)式與系統(tǒng)中動態(tài)元件模型有關(guān),其推導(dǎo)過程見附錄A。

根據(jù)第k次迭代的解x(k)n+1、u(k)n+1,得到第k+1 次迭代的解為:

式（4）、式（5）、式（7）即為牛頓法下交替求解法的基本方程。

對于原始交替求解法迭代格式,可認(rèn)為是在牛頓法下,對式（6）的雅可比矩陣進(jìn)行一定簡化。取A=I和D=Y,并將式（5）的修正方程代入式（7）后得到式（8）。

1.2 交替求解法的不動點(diǎn)迭代形式及新型電力系統(tǒng)背景下的收斂性惡化分析

將式（5）與式（7）聯(lián)立,可得單次迭代中差分方程與代數(shù)方程交替求解的迭代公式為:

由式（6）可知,式（9）等號右側(cè)的矩陣A、D分別為 關(guān) 于xn+1、un+1的 表 達(dá) 式,q(k)n+1、g(k)n+1也 分 別 為 關(guān)于xn+1、un+1的表達(dá)式。因此,式（9）具有標(biāo)準(zhǔn)的不動點(diǎn)迭代格式。狀態(tài)變量x與代數(shù)變量u的迭代函數(shù)Φx(xn+1)、Φu(un+1)分別為:

不動點(diǎn)迭代的收斂性定理指出,若迭代函數(shù)Φ(x)在定義域S內(nèi)有不動點(diǎn)x?,Φ(x)的各分量函數(shù)具有連續(xù)偏導(dǎo)數(shù),且Φ(x)的導(dǎo)數(shù)矩陣Φ′(x)的譜半徑滿足ρ(Φ′(x))<1,則對任意的x0∈S,迭代函數(shù)Φ(x)將收斂于不動點(diǎn)x?［26］。

因此,在不動點(diǎn)迭代方法下,對差分方程、網(wǎng)絡(luò)方程迭代求解的收斂性的討論就轉(zhuǎn)化為了對相應(yīng)迭代函數(shù)譜半徑大小的討論。對于式（10）的不動點(diǎn)迭代函數(shù),可得相應(yīng)的導(dǎo)數(shù)矩陣Φ′x(xn+1)、Φ′u(un+1)為:

根據(jù)式（13）可知,在交替求解過程中引入矩陣A、D后,相應(yīng)的不動點(diǎn)迭代函數(shù)導(dǎo)數(shù)矩陣中僅包含xn+1、un+1相關(guān)的高階項(xiàng),是相對稀疏的矩陣,故其具有較好的收斂性。這與文獻(xiàn)［13］中通過引入A矩陣來最大限度地提高狀態(tài)變量迭代收斂性的結(jié)論是一致的。

對于原始交替求解法,由于取A=I、D=Y,其相應(yīng)的導(dǎo)數(shù)矩陣變?yōu)?

從式（15）中可看出,原始交替求解法的迭代收斂性由差分方程以及代數(shù)方程的迭代收斂性共同決定。同時(shí),由于x、u在仿真過程中是時(shí)變的,導(dǎo)致導(dǎo)數(shù)矩陣譜半徑時(shí)變,從而仿真系統(tǒng)的迭代收斂性會隨仿真過程發(fā)生變化。

隨機(jī)擾動的出現(xiàn)使得代數(shù)方程的導(dǎo)數(shù)矩陣中出現(xiàn)了與擾動相關(guān)的項(xiàng),從而惡化了代數(shù)方程導(dǎo)數(shù)矩陣的譜分布。若計(jì)算過程中的收斂性主要取決于代數(shù)方程,則由于外部擾動不確定性的增強(qiáng),將使計(jì)算過程的收斂性變差。

此外,在新能源與電力電子模型中,各類動態(tài)過程中廣泛存在飽和、限幅等邏輯判斷環(huán)節(jié)。在邏輯環(huán)節(jié)處理層面,在每輪狀態(tài)變量與代數(shù)變量迭代計(jì)算結(jié)束后進(jìn)行邏輯判斷,根據(jù)限幅等環(huán)節(jié)的判斷結(jié)果更新本次迭代得到的變量結(jié)果。由于系統(tǒng)方程與邏輯環(huán)節(jié)的求解不是同步進(jìn)行的,實(shí)際計(jì)算需要在系統(tǒng)方程與邏輯環(huán)節(jié)間進(jìn)行更多次的迭代才能同時(shí)滿足系統(tǒng)方程與邏輯判斷環(huán)節(jié)的約束,導(dǎo)致系統(tǒng)迭代求解的收斂性變差。

根據(jù)式（13）,在交替求解法中引入完整的矩陣A、D或矩陣A、D的近似矩陣,有助于減小譜半徑,從而提高收斂性。一般來說,引入的A矩陣可分為以下3 種情況［27］:A=I/α,其中|α|<1 為松弛系數(shù)；A=diag(?qn+1/?xn+1)；A=?qn+1/?xn+1。引 入 的D矩陣根據(jù)注入電流對電壓偏導(dǎo)表達(dá)式的等效程度不同,可分為虛擬導(dǎo)納法、功角等效法［12］,相應(yīng)的D矩陣表達(dá)式分別為附錄A 中的D=Y+Yg和D=Y+Yg+Ytc,其中,Yg、Ytc分別為虛擬導(dǎo)納法和功角等效法分離出的定常導(dǎo)納矩陣。由于矩陣A、D與狀態(tài)變量、代數(shù)變量相關(guān),每次迭代都需要進(jìn)行更新,將產(chǎn)生巨大的計(jì)算開銷。實(shí)際計(jì)算中,通常采用非誠實(shí)牛頓（very dishonest Newton,VDHN）法,在不同時(shí)步間和迭代間傳遞矩陣A、D,以提高計(jì)算效率［9］。然而,采用VDHN 法對矩陣A、D進(jìn)行近似后,實(shí)際矩陣A、D與當(dāng)前迭代中采用的矩陣存在誤差ΔA、ΔD。對于剛性系統(tǒng),時(shí)間常數(shù)較小的動態(tài)環(huán)節(jié)對應(yīng)的狀態(tài)變量將在暫態(tài)過程中發(fā)生快速變化,使得ΔA、ΔD隨暫態(tài)過程的推進(jìn)變大,從而對譜半徑的影響逐漸增大,導(dǎo)致仿真收斂性變差。因此,VDHN 法雖然提高了時(shí)域計(jì)算的效率,但代價(jià)是犧牲了一定的收斂性,且對于剛性系統(tǒng)將導(dǎo)致顯著的收斂性問題。

1.3 提高交替求解法收斂性的討論

通過上述分析可知,交替求解法收斂性惡化的機(jī)理以及提高交替求解法收斂性的方法可以在迭代函數(shù)導(dǎo)數(shù)矩陣譜分布與譜半徑的視角下進(jìn)行理解。由此,可將提高交替求解法的收斂性方法分為3 類:

1）限制系統(tǒng)最大仿真步長,從而整體性壓縮導(dǎo)數(shù)矩陣的譜半徑,對于剛性系統(tǒng),需要將最大步長hmax縮小至10-4數(shù)量級或更?。?4］；

2）采用合理的初值預(yù)估,使迭代初值更接近真值,從而減小誤差；

3）優(yōu)化迭代格式,使得導(dǎo)數(shù)矩陣的元素相對稀疏,有效改善譜分布。

這3 類方法中,減小步長能夠有效改善收斂性,但增加了仿真時(shí)步數(shù)量,降低了計(jì)算效率；改善迭代初值能夠縮小初始值與收斂解的差,但如果初值預(yù)測策略選取過于復(fù)雜,同樣會增大計(jì)算量；優(yōu)化迭代格式能夠有效改善譜分布,但迭代格式的修正也在一定程度上引入了額外的計(jì)算量,降低了計(jì)算效率。此外,上述3 類方法在優(yōu)化譜分布時(shí),并不改變迭代過程的收斂階。若能通過改進(jìn)迭代格式提高交替求解過程的收斂階,則能夠在改善收斂性的同時(shí)兼顧計(jì)算效率。

2 基于斯蒂芬遜迭代的收斂性提升方法

2.1 斯蒂芬遜迭代過程

本節(jié)給出一種基于斯蒂芬遜迭代的交替求解收斂性提升方法,其本質(zhì)是通過構(gòu)造具有更高收斂階的迭代格式,改善交替求解的收斂性。斯蒂芬遜迭代的基本原理如下。

針對迭代函數(shù)Φ(x),設(shè)其不動點(diǎn)為x?,由x0出發(fā)構(gòu)造迭代公式:

式中:下標(biāo)i=0,1,2,…。

將式中的迭代函數(shù)Φ替換為Φx(xn+1)、Φu(un+1),即可得基于斯蒂芬遜迭代的差分方程與代數(shù)方程迭代格式。

已有研究證明,對于由不動點(diǎn)迭代函數(shù)Φ構(gòu)造的斯蒂芬遜迭代函數(shù)Ψ,若Φ存在不動點(diǎn)在x?的領(lǐng)域內(nèi)存在、連續(xù),且,則斯蒂芬遜迭代法至少具有二階收斂性［28］。

采用斯蒂芬遜迭代后,交替求解法的計(jì)算流程如下:

步驟1:初始化n+1 時(shí)刻狀態(tài)變量和代數(shù)變量的迭代初值x(0)n+1、u(0)n+1,迭代次數(shù)k=0。

步驟2:對于第k次迭代,若k=3i+1 或k=3i+2,則根據(jù)式（17）、式（18）依次計(jì)算第k次迭代的狀態(tài)變量和代數(shù)變量結(jié)果x(k)n+1、u(k)n+1；若k=3i+3,則根據(jù)式（19）通過斯蒂芬遜迭代求解第k次迭代的狀態(tài)變量和代數(shù)變量結(jié)果x(k)n+1、u(k)n+1。

步驟3:如果前后兩次迭代的狀態(tài)變量、代數(shù)變量最大迭代誤差max {||x(k+1)n+1-x(k)n+1||,||u(k+1)n+1-u(k)n+1||}小于允許值ε,則認(rèn)為本時(shí)步計(jì)算收斂,否則令k=k+1,返回步驟2,進(jìn)入下一輪迭代。

基于斯蒂芬遜迭代的時(shí)域計(jì)算交替求解流程如圖1 所示。由圖可知,算法在求解流程上與原始交替求解法的區(qū)別在于,在兩次常規(guī)狀態(tài)變量與代數(shù)變量迭代之后,插入一次標(biāo)量式的斯蒂芬遜迭代。

圖1 基于斯蒂芬遜迭代的交替求解法流程圖Fig.1 Flow chart of Steffensen iteration based alternating solving method

需要注意的是,在時(shí)域計(jì)算中,[x,u]為一組向量,相應(yīng)的迭代函數(shù)Φ也為向量函數(shù)。在上述流程中,當(dāng)k=3i+1 或k=3i+2 時(shí),是以向量的形式進(jìn)行迭代；當(dāng)k=3i+3 時(shí),是以標(biāo)量的形式進(jìn)行斯蒂芬遜迭代。因此,實(shí)際迭代過程比單變量迭代過程更加復(fù)雜,但標(biāo)量式斯蒂芬遜迭代的總體過程是不變的。

此外,由于限幅、死區(qū)等邏輯判斷環(huán)節(jié)的存在,實(shí)際DAEs 難以滿足待求方程的二階導(dǎo)數(shù)連續(xù)且存在的前提。實(shí)際求解時(shí),將求解過程分解為“不含邏輯判斷環(huán)節(jié)的方程求解”以及“根據(jù)邏輯判斷環(huán)節(jié)更新系統(tǒng)變量”兩個(gè)步驟。在“不含邏輯判斷環(huán)節(jié)的方程求解”中,由于不需要考慮限幅、死區(qū)等帶來的不連續(xù)及變量跳變問題,方程組仍滿足二階導(dǎo)數(shù)連續(xù)且存在的前提,但無法對系統(tǒng)方程與邏輯環(huán)節(jié)進(jìn)行同步求解,這也是算法的局限性所在。

2.2 算法的計(jì)算效率分析

在提高算法收斂性的同時(shí),還需要考慮算法是否會引入額外的計(jì)算量,從而導(dǎo)致時(shí)域計(jì)算的效率降低。因此,需要對比引入斯蒂芬遜迭代前后交替求解法的計(jì)算量差異。由于計(jì)算機(jī)乘除運(yùn)算耗時(shí)遠(yuǎn)大于加減運(yùn)算,算法計(jì)算量僅統(tǒng)計(jì)乘除運(yùn)算次數(shù)。

引入斯蒂芬遜迭代后,算法的計(jì)算量差異體現(xiàn)在以下兩方面:1）根據(jù)式（6）構(gòu)造A、D矩陣造成的計(jì)算量增加；2）各算法由于迭代公式的不同,導(dǎo)致進(jìn)行狀態(tài)變量與代數(shù)變量求解存在計(jì)算量差異。

對于方面1）,由于仿真過程中A、D矩陣各元素表達(dá)式差異較大,難以統(tǒng)計(jì)每個(gè)元素的計(jì)算量,只能對構(gòu)造A、D矩陣的整體計(jì)算量進(jìn)行評估。本文所提算法采用VDHN 法在各時(shí)步、各輪迭代間傳遞雅可比矩陣A,僅在收斂性較差時(shí)更新雅可比矩陣A。同時(shí),采用虛擬導(dǎo)納法,用D矩陣的常數(shù)部分Y+Yg近似D矩陣。因此,仿真過程中僅需要初始時(shí)刻以及拓?fù)渥兓臅r(shí)刻重新構(gòu)造D矩陣。若記仿真全過程A、D矩陣更新次數(shù)分別為nA、nD,單次更新A、D矩陣的計(jì)算量分別為OA、OD,則構(gòu)造A、D矩陣的計(jì)算量C為:

由于采用了VDHN 法,nA的值實(shí)際上較?。?7］,由于D矩陣的定常特性,nD的值也較小。因此,仿真全過程中構(gòu)造A、D矩陣引入的額外計(jì)算量并不是影響仿真過程計(jì)算效率的主要原因。

對于方面2）,設(shè)差分方程和網(wǎng)絡(luò)方程的階數(shù)分別為n、m。記單次迭代中常規(guī)狀態(tài)變量迭代與代數(shù)變量迭代的計(jì)算量為O1、O2。由式（19）可知,單次斯蒂芬遜迭代的計(jì)算包含1 次乘法操作和1 次除法操作。因此,對于n階差分方程和m階網(wǎng)絡(luò)方程,各需對n個(gè)狀態(tài)變量和m個(gè)代數(shù)變量進(jìn)行斯蒂芬遜迭代,計(jì)算量分別為2n和2m。假設(shè)某時(shí)步共需進(jìn)行K次迭代,則常規(guī)迭代與所提算法在求解差分方程與代數(shù)方程上的計(jì)算量對比如表1 所示。

表1 常規(guī)迭代與斯蒂芬遜迭代計(jì)算量差異Table 1 Difference in calculation overhead between regular iterations and Stephenson iterations

根據(jù)表1 可知,采用常規(guī)迭代格式與斯蒂芬遜迭代格式時(shí),計(jì)算量的差異與O1、O2的具體表達(dá)式相關(guān),而O1、O2的表達(dá)式則與實(shí)際計(jì)算中矩陣A與矩陣D的取值有關(guān)。需要作進(jìn)一步分析:

1）若采用式（3）的原始交替求解法,即取A=I、D=Y,對單個(gè)差分方程的求解涉及1 次乘法操作,故O1等于狀態(tài)變量的階數(shù)n；而代數(shù)方程的求解則涉及m階導(dǎo)納矩陣的求逆,故O2為kmm3,其中,m3為矩陣求逆的計(jì)算復(fù)雜度,而km為考慮導(dǎo)納矩陣稀疏度所取系數(shù)。此時(shí),常規(guī)迭代與斯蒂芬遜迭代中,差分方程計(jì)算量表達(dá)式均為n的一次表達(dá)式,但采用斯蒂芬遜迭代時(shí)一次項(xiàng)系數(shù)稍大；網(wǎng)絡(luò)方程計(jì)算量表達(dá)式均為m的三次表達(dá)式,但斯蒂芬遜迭代的三次項(xiàng)系數(shù)稍小。

2）若采用式（9）的迭代格式,引入完整矩陣A與近似矩陣D,則差分方程與網(wǎng)絡(luò)方程的求解計(jì)算量體現(xiàn)在對矩陣A、D的求逆上。因此,O1、O2分別為kAn3、kDm3,其中,kA、kD分別為考慮A、D矩陣稀疏度所取系數(shù)。此時(shí),差分方程與網(wǎng)絡(luò)方程計(jì)算量表達(dá)式均分別為n與m的三次表達(dá)式,但斯蒂芬遜迭代的三次項(xiàng)系數(shù)稍小。

從上述兩方面分析可知,由于采用了VDHN法,仿真過程中構(gòu)造A、D矩陣引入的計(jì)算量較??；對于迭代格式不同帶來的計(jì)算量差異,僅在常規(guī)迭代格式采用原始交替求解法時(shí),所提算法在差分方程求解的計(jì)算量上略高于常規(guī)迭代方法,而在其他迭代格式下所提算法均沒有引入額外的計(jì)算量。

此外,算法的計(jì)算效率是單次迭代計(jì)算量與迭代次數(shù)的綜合結(jié)果。即使原始交替求解法的單次迭代計(jì)算量低于斯蒂芬遜迭代,仍可能因?yàn)榈螖?shù)過高,而使總體計(jì)算效率低于基于斯蒂芬遜迭代的改進(jìn)交替求解法。

3 算例分析

3.1 算法正確性驗(yàn)證

在IEEE 標(biāo)準(zhǔn)16 機(jī)68 節(jié)點(diǎn)系統(tǒng)進(jìn)行算法的正確性驗(yàn)證。系統(tǒng)包含16 臺7 階同步發(fā)電機(jī),11 臺IEEE-DC4B 型 勵(lì) 磁 機(jī),1 臺IEEE-ST1A 型 勵(lì) 磁 機(jī),12 臺電力系統(tǒng)穩(wěn)定器,動態(tài)元件最小時(shí)間常數(shù)為0.01 s,雅可比矩陣A規(guī)模為227 階,網(wǎng)絡(luò)導(dǎo)納矩陣為68 階復(fù)數(shù)矩陣。測試電腦CPU 配置為2.3 GHz Xeon 5218,256 GB 內(nèi)存,操作系統(tǒng)為CentOS 7.5.1804。

案例測試中涉及的時(shí)域計(jì)算算法如表2 所示。表中:各算法在代數(shù)方程迭代格式中,D矩陣統(tǒng)一取D=Y+Yg。

表2 算例涉及的算法Table 2 Algorithms involved in case study

考慮到原始交替求解法,即表2 中算法1 存在最大步長限制,無法在0.010 s 的步長下取得收斂,在進(jìn)行正確性驗(yàn)證時(shí),算法1 采用0.005 s 固定步長進(jìn)行仿真,其余算法采用0.010 s 固定步長進(jìn)行仿真。擾動設(shè)置為1.0 s 時(shí)刻,220 kV 3 號母線發(fā)生三相短路故障,t=1.1 s 時(shí)故障清除。仿真時(shí)長為20 s,觀察變量為發(fā)電機(jī)G2、G4、G6、G8、G10、G12、和G16 的功角差。

G2 與G16 的功角差曲線如圖2 所示,其余發(fā)電機(jī)的功角差曲線見附錄B 圖B1。從圖中可以看出,采用不同算法進(jìn)行時(shí)域仿真得到的結(jié)果基本一致。

圖2 采用不同算法時(shí)的時(shí)域曲線Fig.2 Time-domain curves by using different algorithms

考慮到原始交替求解法的仿真固定步長小于其余4 種算法,后續(xù)算例中將不再對原始交替求解法進(jìn)行討論。

3.2 算法收斂性提升效果驗(yàn)證

為驗(yàn)證算法在含新能源場景下也能提高時(shí)域計(jì)算的收斂性與計(jì)算效率,設(shè)置68 節(jié)點(diǎn)系統(tǒng)的新能源場 景 為:將 系 統(tǒng) 中1 號、2 號、5 號、7 號、8 號 以 及10 號母線上的同步發(fā)電機(jī)替換為相同容量的雙饋風(fēng)機(jī),雙饋風(fēng)機(jī)模型及參數(shù)與MATLAB2018b/Simulink 中雙饋風(fēng)機(jī)相量模型相同。

針對無新能源場景和含新能源場景,分別對每個(gè)算法重復(fù)進(jìn)行10 次仿真,并對計(jì)算耗時(shí)求平均值。統(tǒng)計(jì)兩種場景下各算法的總迭代次數(shù)及耗時(shí)結(jié)果如表3 所示。

表3 采用不同算法時(shí)的總迭代次數(shù)與耗時(shí)對比Table 3 Comparison of total number of iterations and time consumption when using different algorithms

從表3 可以看出,不論是無新能源場景還是含新能源場景,在矩陣A格式相同的情況下,引入斯蒂芬遜迭代均可有效改善暫態(tài)時(shí)域仿真的迭代過程。在無新能源場景中:算法3 的仿真耗時(shí)相對算法2 減少70.5%；算法5 的仿真耗時(shí)相對算法4 減少63.8%；算法5 的仿真耗時(shí)相對算法2 減少70.8%。在含新能源場景中:算法3 的仿真耗時(shí)相對算法2減少40.2%；算法5 的仿真耗時(shí)相對算法4 減少38.4%；算法5 的仿真耗時(shí)相對算法2 減少40.7%。

結(jié)合2.3 節(jié)算法的計(jì)算效率分析可知,雖然算法3 的單次迭代計(jì)算量略高于算法2,但由于引入斯蒂芬遜迭代后計(jì)算過程的收斂性得到了較大改善,總的迭代次數(shù)遠(yuǎn)低于算法2。因此,在計(jì)算效率方面,算法3 相較算法2 取得了巨大的優(yōu)勢。同樣的分析也適用于算法4 與算法5。

此外,表3 中還給出了新能源場景下,固定步長為0.001 s 時(shí)算法2 至算法5 的性能對比。此時(shí),算法2 和算法4 無法在0.001 s 的固定步長下收斂,而算法3 和算法5 則仍能在0.001 s 的固定步長下收斂。這說明在新能源場景下,引入斯蒂芬遜迭代后,時(shí)域計(jì)算對新能源接入帶來的剛性有更好的適應(yīng)性。同時(shí),從迭代次數(shù)與仿真耗時(shí)來看,算法5 的收斂性優(yōu)于算法3,說明引入完整的矩陣A對于提高時(shí)域計(jì)算的收斂性是有益的。

圖3 給出了含新能源場景下20 s 仿真過程中,4 種算法每一時(shí)步的迭代次數(shù)比較。

圖3 采用不同算法時(shí)的迭代次數(shù)曲線Fig.3 Curves of iteration numbers when using different algorithms

從圖3 中可以看出,在仿真初期,4 種算法都可以在較低的迭代次數(shù)內(nèi)完成本時(shí)步的計(jì)算。當(dāng)1.0 s 發(fā)生短路故障后,4 種算法的迭代次數(shù)都出現(xiàn)明顯的提升,但算法3 和算法5 的迭代次數(shù)增長幅度明顯低于算法2 和算法4,且算法3 和算法5 在1.5 s左右穩(wěn)定于每時(shí)步進(jìn)行4 次迭代,而算法2 與算法4則在暫態(tài)后一直保持較高的迭代次數(shù)。

3.3 實(shí)際電網(wǎng)下的算法效果驗(yàn)證

在某2 970 節(jié)點(diǎn)區(qū)域電網(wǎng)驗(yàn)證本文所提算法對大型電力系統(tǒng)仿真的有效性。系統(tǒng)包含2 970 條母線、1 613 條線路、254 臺同步發(fā)電機(jī)與12 臺雙饋風(fēng)機(jī),包含勵(lì)磁、調(diào)速共16 種模型,負(fù)荷采用恒定阻抗+感應(yīng)電機(jī)模型,雅可比矩陣A規(guī)模為9 354 階,網(wǎng)絡(luò)導(dǎo)納矩陣規(guī)模為2 970 階,擾動設(shè)置為t=1.0 s某500 kV 母線三相短路故障,t=1.1 s 故障清除。

圖4 展示了采用本文算法以及BPA 進(jìn)行仿真時(shí),故障母線處發(fā)電機(jī)的功角差曲線。其他變量的對比見附錄B 圖B2?？梢钥闯?本文算法的仿真結(jié)果與BPA 的仿真結(jié)果具有一致性。

圖4 本文算法與BPA 時(shí)域曲線對比Fig.4 Comparison of time-domain curves by the proposed method in and BPA

采用表2 中的算法2、3、4 重復(fù)進(jìn)行10 次10 s 的時(shí)域仿真,統(tǒng)計(jì)各算法的總迭代次數(shù)以及耗時(shí)結(jié)果如表4 所示。

表4 新能源場景下采用不同算法時(shí)的總迭代次數(shù)與耗時(shí)對比Table 4 Comparison of total iteration numbers and time consumption when using different algorithms in renewable energy scenarios

從表4 可以看出,在大電網(wǎng)仿真中,算法2 與算法4 無法在較大步長條件下進(jìn)行仿真,需要將步長調(diào)整至10-6數(shù)量級,無法滿足機(jī)電暫態(tài)仿真的需要,而算法3 和算法5 則仍能以0.001 s 的固定步長進(jìn)行仿真。這說明引入斯蒂芬遜迭代對于時(shí)域計(jì)算收斂性提升及系統(tǒng)強(qiáng)剛性的適應(yīng)性。

4 結(jié)語

本文從不動點(diǎn)迭代的理論出發(fā),分析交替求解法的迭代求解過程,揭示了交替求解法不收斂的機(jī)理,提出了基于斯蒂芬遜迭代的電力系統(tǒng)時(shí)域計(jì)算交替求解法。通過算例分析得到以下結(jié)論:

1）對于剛性電力系統(tǒng),由于系統(tǒng)存在較小時(shí)間常數(shù)的動態(tài)環(huán)節(jié),傳統(tǒng)交替求解法及其改進(jìn)算法都將面臨微分方程強(qiáng)剛性帶來的收斂性問題。

2）引入斯蒂芬遜迭代能夠有效改善交替求解法對于系統(tǒng)剛性的適應(yīng)性,能夠采用比原始交替求解法更大的步長對系統(tǒng)進(jìn)行定步長仿真。

3）通過引入差分方程與代數(shù)方程的雅可比矩陣A、D,能夠較好地改善系統(tǒng)的收斂性,配合斯蒂芬遜迭代則能夠最大限度地兼顧交替求解過程的收斂性與效率。

4）在IEEE 標(biāo)準(zhǔn)68 節(jié)點(diǎn)系統(tǒng)中,同時(shí)引入矩陣A、D并采用斯蒂芬遜迭代進(jìn)行時(shí)域計(jì)算,在無新能源和有新能源接入的情況下,相比于原始交替求解法,仿真耗時(shí)分別減少70.8%和40.7%。在2 970 節(jié)點(diǎn)實(shí)際電網(wǎng)中,引入斯蒂芬遜迭代也能夠改善時(shí)域計(jì)算的收斂性,適應(yīng)系統(tǒng)的強(qiáng)剛性。

未來將進(jìn)一步研究非線性方程求解領(lǐng)域的其他高階迭代方法在時(shí)域計(jì)算中的適用性。

附錄見本刊網(wǎng)絡(luò)版（http：//www.aeps-info.com/aeps/ch/index.aspx），掃英文摘要后二維碼可以閱讀網(wǎng)絡(luò)全文。