基于Adomian分解方法的暫態(tài)穩(wěn)定并行仿真研究*

劉云飛，江全元*，陳躍輝，張文磊，宋軍英

(1.浙江大學(xué)電氣工程學(xué)院，浙江杭州310027;2.湖南省電力公司，湖南長(zhǎng)沙410007)

0 引言

為滿(mǎn)足電源大規(guī)模集中投產(chǎn)和用電負(fù)荷增長(zhǎng)的需要，我國(guó)電網(wǎng)規(guī)模不斷擴(kuò)大。截至2011年底，除臺(tái)灣地區(qū)外，全國(guó)聯(lián)網(wǎng)格局已經(jīng)形成［1］，而且對(duì)電網(wǎng)的安全性和可靠性要求也在不斷提高。所以，尋找一種快速收斂、高度并行、適合大規(guī)模系統(tǒng)的暫態(tài)穩(wěn)定算法就顯得尤為重要。

國(guó)內(nèi)外很多學(xué)者在暫態(tài)穩(wěn)定算法方面開(kāi)展過(guò)相關(guān)研究［2-4］，其目的就是為了能夠快速求解大規(guī)模電力系統(tǒng)，實(shí)現(xiàn)(超)實(shí)時(shí)仿真。然而，隨著并行計(jì)算機(jī)的快速發(fā)展，并行算法已經(jīng)成為處理暫態(tài)穩(wěn)定問(wèn)題的主流方法。并行方式從原理上大致可分為3類(lèi):空間并行［5］、時(shí)間并行［6-7］和波形松弛［8-10］(Waveform Relaxation，WR)。但不論采用何種并行方式，都面臨對(duì)微分代數(shù)方程組(DAE)初值問(wèn)題的求解。交替求解和聯(lián)立求解是求解該類(lèi)問(wèn)題的主流方法，而由于聯(lián)立求解不存在交接誤差，在對(duì)仿真精度要求較高時(shí)被廣泛采用［11］。其基本過(guò)程為:先用隱式積分公式將微分方程組差分化，然后和代數(shù)方程組聯(lián)立，得到一個(gè)非線(xiàn)性方程組，而現(xiàn)有的研究往往采用牛頓法(Newton’s method)迭代求解該非線(xiàn)性方程組。其優(yōu)點(diǎn)是迭代函數(shù)簡(jiǎn)單，所以得到廣泛應(yīng)用。但是一些對(duì)收斂速度要求較高的場(chǎng)合，牛頓法往往不能滿(mǎn)足收斂性的要求。

近幾百年來(lái)，數(shù)學(xué)工作者們從來(lái)沒(méi)有停止對(duì)新型或高階收斂的迭代法的研究。上世紀(jì)80年代初，美國(guó)數(shù)學(xué)家G.Adomian［12］提出了一種解非線(xiàn)性方程(組)問(wèn)題，而后以他名字命名的Adomian分解方法。該方法的本質(zhì)是利用一個(gè)無(wú)窮級(jí)數(shù)去逼近問(wèn)題的精確解。通過(guò)構(gòu)造不同形式的Adomian多項(xiàng)式，可以得到不同的迭代格式和收斂精度。文獻(xiàn)［13］和文獻(xiàn)［14］分別提出了各自的分解格式，進(jìn)而得到了不同的迭代形式。Changbum Chu在文獻(xiàn)［15］中指出，利用Adomian級(jí)數(shù)技巧構(gòu)造的迭代，只要展開(kāi)的項(xiàng)數(shù)足夠多，就可以得到任意精度的解。文獻(xiàn)［16］在原Adomian分解法上加以修正，得到了一種收斂更快、精度更高的求解非線(xiàn)性方程的迭代算法，并可推廣到非線(xiàn)性系統(tǒng)的求解。

本研究采用波形松弛法作為并行框架，將Adomian分解方法配合非誠(chéng)實(shí)牛頓算法(Very dishonest Newton method，VDHN)應(yīng)用于各個(gè)子系統(tǒng)的求解。同時(shí)為了進(jìn)一步加快內(nèi)部收斂，本研究還使用預(yù)處理過(guò)程和初值猜測(cè)提高整體收斂性。此外，本研究采用OpenMP的共享內(nèi)存并行環(huán)境，以提高本研究的并行度。

1 Adomian分解方法

本節(jié)將對(duì)Adomian多項(xiàng)式和基于Adomian級(jí)數(shù)分解的迭代算法作簡(jiǎn)單介紹。

1.1 Adomian多項(xiàng)式及其計(jì)算

考慮非線(xiàn)性方程:

式中:f—連續(xù)可微函數(shù)。

將方程(1)改寫(xiě)成如下的形式:

式中:c—常數(shù)，N—非線(xiàn)性部分。

Adomian級(jí)數(shù)法首先把x分解為如下的無(wú)窮級(jí)數(shù)形式:

而非線(xiàn)性部分N分解為一個(gè)函數(shù)項(xiàng)級(jí)數(shù):

式中:Ai—依賴(lài)于x0，x1，x2，…，xi的Adomian多項(xiàng)式，其形式如下:

由式(2～4)可知:

從而級(jí)數(shù)形式解的每項(xiàng)xi與Adomian多項(xiàng)式之間的關(guān)系為:

解的部分和為:

且有:

1.2 非線(xiàn)性方程組的Adomian級(jí)數(shù)法

研究者得到了Adomian分解法的構(gòu)造，就可以通過(guò)該方法求解非線(xiàn)性方程式(1)。根據(jù)不同的分解方法，可以得到Adomian級(jí)數(shù)的不同形式，進(jìn)而得到相應(yīng)的迭代格式。其中Changbum Chu提出了如下分解:設(shè)α是方程式(1)的一個(gè)單根，γ是充分靠近α的一個(gè)初值。將f(x)在γ處一階Taylor展開(kāi)，可得到方程式(1)的耦合形式:

其中，余項(xiàng)為:

方程式(9)可以改寫(xiě)為式(2)的形式:

其中:

可以得到Adomian多項(xiàng)式的前幾項(xiàng)的具體表達(dá)式為:

由式(8)可知，Xm可以作為x的一個(gè)近似解，實(shí)際所需的展開(kāi)項(xiàng)很少時(shí)就可以達(dá)到所需要的精度。利用上述關(guān)系，取不同的項(xiàng)數(shù)來(lái)近似方程的解，就可以得到不同精度的各種迭代格式。

當(dāng)m=0時(shí):

從而得到二階收斂的Newton迭代序列:

當(dāng)m=1時(shí):

得到兩步三階收斂的迭代序列:

由于每一個(gè)子系統(tǒng)都需要求解一個(gè)相對(duì)較小規(guī)模的DAE系統(tǒng)，差分化后得一個(gè)非線(xiàn)性方程組。本研究采用式(15)的兩步迭代格式求解對(duì)各個(gè)子系統(tǒng)，該格式只需求解一次雅克比矩陣，與傳統(tǒng)牛頓迭代相比，僅僅增加了一次前推回代的過(guò)程。為了進(jìn)一步減少LU分解的次數(shù)，本研究采用VDHN算法，即在系統(tǒng)結(jié)構(gòu)不發(fā)生變化的情況下，迭代次數(shù)沒(méi)有超出設(shè)定閾值時(shí)，不更新雅克比矩陣。

2 暫態(tài)穩(wěn)定并行仿真流程

本研究采用波形松弛法作為本研究的并行方法，同時(shí)采用預(yù)處理和波形預(yù)測(cè)加速收斂。筆者采用商業(yè)軟件將大系統(tǒng)拆分成若干子系統(tǒng)，子系統(tǒng)內(nèi)部采用上節(jié)提到的基于Adomian級(jí)數(shù)分解的迭代算法求解，最后在共享內(nèi)存平臺(tái)上并行實(shí)現(xiàn)。

2.1 波形松弛方法及其加速方法

和許多大型電路系統(tǒng)的仿真模擬或多體動(dòng)力學(xué)模擬過(guò)程一樣，電力系統(tǒng)暫態(tài)穩(wěn)定問(wèn)題可以用一系列指標(biāo)1型、非線(xiàn)性、半隱式微分代數(shù)方程來(lái)描述，可以表示為:

式中:T＞0;x(t)∈Rn—狀態(tài)變量，用來(lái)描述發(fā)電機(jī)和勵(lì)磁調(diào)速系統(tǒng);V(t)∈Rm—代數(shù)變量，即節(jié)點(diǎn)電壓;x0—狀態(tài)變量初值。

用波形松弛法求解如式(16)所示的系統(tǒng)，可以將其分解成p個(gè)子系統(tǒng)，其過(guò)程可以表示為:

本研究采用雅克比(Jacobi)迭代求解，雅克比迭代是松弛算法中最基本而有效的迭代方法，而且是天然并行的算法，其過(guò)程可表示為:

本研究把整個(gè)時(shí)間段分為若干窗口，在每個(gè)窗口內(nèi)波形收斂后，再進(jìn)行下一窗口的迭代，以提高整體收斂性。窗口長(zhǎng)度的選取會(huì)影響到迭代的快慢和并行效率，但是由于最佳窗口長(zhǎng)度求解極為復(fù)雜［17］，工程上難以實(shí)現(xiàn)，故本研究選取固定的窗口長(zhǎng)度。

本研究同時(shí)采用文獻(xiàn)［18］中的預(yù)處理方法加速波形松弛法的收斂。其格式如下:

式中:P—預(yù)處理子，P=?g(x，V)/?V。

每次迭代之后增加如式(19)的一部求解，可以將每個(gè)子系統(tǒng)的電壓信息傳遞給其他分組，從而加速收斂。

由于波形松弛法在迭代開(kāi)始時(shí)每個(gè)子系統(tǒng)需要將其他子系統(tǒng)波形的“猜測(cè)值”作為初值，本研究采用3點(diǎn)Lagrange型多項(xiàng)式插值公式對(duì)所有變量在窗口內(nèi)的波形進(jìn)行預(yù)測(cè)，以達(dá)到減少迭代次數(shù)的目的，格式如下:

式中:yn—第n次迭代狀態(tài)變量和非狀態(tài)變量的值。

2.2 并行算法流程

綜上所述，本研究算法流程如圖1所示。

并行前，首先需要對(duì)系統(tǒng)進(jìn)行拆分，現(xiàn)有的系統(tǒng)分塊方法主要分為兩種:①采用商業(yè)軟件對(duì)系統(tǒng)進(jìn)行分塊;②按照地理分區(qū)對(duì)系統(tǒng)進(jìn)行分組。以上兩種方法都存在缺陷，對(duì)于空間并行方法，兩種分組方法無(wú)法在兼顧協(xié)調(diào)子網(wǎng)(關(guān)聯(lián)子網(wǎng))規(guī)模的大小的同時(shí)，將各個(gè)進(jìn)(線(xiàn))程的任務(wù)平均分配;而對(duì)于常規(guī)的波形松弛方法，兩種方法都不能保證良好的收斂性能。國(guó)內(nèi)外有

圖1 算法流程

部分學(xué)者針對(duì)電力系統(tǒng)分區(qū)策略開(kāi)展相關(guān)研究［19-20］。

由于本研究采用了預(yù)處理的加速方法，這樣每個(gè)節(jié)點(diǎn)的“信息”會(huì)在每次迭代后第一時(shí)間傳遞到其他分組，大幅縮短了信息的傳遞距離，提高了收斂性。故本研究采用分圖軟件METIS［21］對(duì)原系統(tǒng)進(jìn)行拆分，盡可能使各線(xiàn)程任務(wù)分配更加均勻，使各個(gè)子系統(tǒng)內(nèi)的變量數(shù)盡可能達(dá)到平衡。因此需要考慮每個(gè)節(jié)點(diǎn)的權(quán)重，即每個(gè)節(jié)點(diǎn)所含的狀態(tài)變量和電壓變量的個(gè)數(shù)。得到N個(gè)分組后，交給N個(gè)線(xiàn)程進(jìn)行并行求解。

在每次開(kāi)始新的窗口計(jì)算時(shí)，本研究先采用式(20)對(duì)波形進(jìn)行外推。各個(gè)子系統(tǒng)采用隱式梯形積分公式對(duì)狀態(tài)變量進(jìn)行差分化處理。每次求解之后，各個(gè)線(xiàn)程分別對(duì)各自的子系統(tǒng)進(jìn)行收斂判斷，并得到收斂標(biāo)志。在全部分組收斂后，系統(tǒng)開(kāi)始下一個(gè)窗口的計(jì)算。

2.3 并行環(huán)境的選擇

由于每次迭代，所有分組都需要交換數(shù)據(jù)、傳遞波形，如果采用分布式內(nèi)存的并行方式，就需要各個(gè)節(jié)點(diǎn)之間頻繁地通信，通信時(shí)間所占的比重會(huì)很大，會(huì)導(dǎo)致并行效率很低。

OpenMP［22］是一個(gè)共享存儲(chǔ)并行系統(tǒng)上的應(yīng)用編程接口，是基于線(xiàn)程的并行編程模型。OpenMP使用Fork-join并行執(zhí)行模型，該模型如圖2所示。由于OpenMP是基于共享內(nèi)存的編程模型，其更適合用于本研究這種細(xì)粒度、需要密集通信的并行環(huán)境。

圖2 OpenMP的Fork-join模型

3 算例測(cè)試

本研究的算法通過(guò)C++語(yǔ)言編程實(shí)現(xiàn)，采用的測(cè)試環(huán)境為Beowulf集群計(jì)算機(jī)，測(cè)試單節(jié)點(diǎn)為8路8核Intel Xeon CPU的SMP(對(duì)稱(chēng)多處理)結(jié)構(gòu)，同時(shí)采用Linux操作系統(tǒng)，并調(diào)用KLU［23］解法器求解線(xiàn)性方程組。分別采用Matpower［24］中的2 383節(jié)點(diǎn)算例與一個(gè)12 685節(jié)點(diǎn)算例對(duì)算法進(jìn)行測(cè)試，測(cè)試算例如表1所示。其中發(fā)電機(jī)采用6階模型，勵(lì)磁采用4階模型，調(diào)速器采用1階模型。仿真時(shí)間為20 s，仿真步長(zhǎng)0.02 s，窗口長(zhǎng)度0.1 s，收斂精度0.000 1。

表1 測(cè)試算例

本研究采用METIS對(duì)兩個(gè)算例進(jìn)行分組，最大分組數(shù)為16，分組結(jié)果如表2所示。為了表征分組的平均程度，定義如下變量:

當(dāng)變量Balance接近1時(shí)，表示分組越均勻。

表2 分組結(jié)果

為了驗(yàn)證算法精度，本研究采用經(jīng)典串行的VDHN算法作為比對(duì)算法，以Case2上萬(wàn)節(jié)點(diǎn)系統(tǒng)為例進(jìn)行比對(duì)。本研究選取了兩種故障，分別為0.1 s切除故障和0.5 s切除故障，其中某兩臺(tái)發(fā)電機(jī)的相對(duì)功角誤差如圖3所示，從圖3中可以看出，最大誤差為6.05×10-6rad，符合精度要求。

圖3 功角誤差曲線(xiàn)

為驗(yàn)證Adomian分解法的有效性，本研究仍采用VDHN算法求解子系統(tǒng)，作為比對(duì)，并行時(shí)間如表3所示。從表3可以看出，Adomian分解法可提高算法效率，最多可達(dá)39.17%，驗(yàn)證了本研究算法的有效性。

表3 兩種算法的求解時(shí)間比對(duì)(單位:s)

此外，從表3可以看出，通過(guò)本研究算法已實(shí)現(xiàn)了上萬(wàn)節(jié)點(diǎn)系統(tǒng)的超實(shí)時(shí)仿真。

本研究算法的并行度如圖4所示。其中并行加速比定義如下:

圖4 并行加速比

從圖4可以看出，本研究采用的OpenMP并行方式可以獲得較高的并行加速比。并且本研究算法對(duì)更大規(guī)模的系統(tǒng)有更好的適應(yīng)性，在16分組時(shí)可獲得10.65的并行度?，F(xiàn)有的空間并行算法中，最高并行度可達(dá)9.82［25］，現(xiàn)有波形松弛算法在相同加速比定義的情況下，最高加速比也低于10［26］(2 500節(jié)點(diǎn)系統(tǒng)，20進(jìn)程并行)。

4 結(jié)束語(yǔ)

本研究提出了一種基于Adomian分解的迭代算法求解大規(guī)模電力系統(tǒng)暫態(tài)穩(wěn)定問(wèn)題的并行算法。該算法采用波形松弛法作為本研究的并行框架，同時(shí)采用預(yù)處理和波形預(yù)測(cè)進(jìn)一步加快收斂過(guò)程。此外，本研究采用OpenMP的共享內(nèi)存并行環(huán)境，來(lái)提高本研究的并行度。

通過(guò)對(duì)兩個(gè)較大算例的測(cè)試，證明本研究算法在保證算法精度的同時(shí)，還獲得較高收斂速度和并行效率，最終實(shí)現(xiàn)了上萬(wàn)節(jié)點(diǎn)系統(tǒng)超實(shí)時(shí)仿真，充分說(shuō)明了本研究算法的有效性和對(duì)大規(guī)模系統(tǒng)的適應(yīng)性。

［1］劉振亞.中國(guó)電力與能源［M］.北京:中國(guó)電力出版社，2012:345.

［2］趙志奇，王建全.隱式精細(xì)積分算法在電力系統(tǒng)暫態(tài)穩(wěn)定分析中的應(yīng)用［J］.機(jī)電工程，2012，29(5):580-583.

［3］武同心，呂曉祥，王建全.Duhamel數(shù)值積分算法在電力系統(tǒng)暫態(tài)穩(wěn)定分析中的應(yīng)用［J］.機(jī)電工程，2013，30(6):741-745.

［4］FABOZZI D，CHIEH A S，HAUT B，et al.Accelerated and localized newton schemes for faster dynamic simulation of large power systems［J］.IEEE Transactions on Pouer Systens，2013，28(4):4936-4947.

［5］QUANYUAN J，HAN J.OpenMP-based parallel transient stability simulation for large-scale power systems［J］.Science China:Technological Sciences，2012(10):2837-2846.

［6］HONG C.Implementation of parallel-in-time newton method for transient stability analysis on a message passing multi-Computer［C］//Internation al Conferen ce on Power Syztem Teehnology，2002:1239-1243.

［7］洪潮.電力系統(tǒng)暫態(tài)穩(wěn)定計(jì)算的一種時(shí)間并行算法［J］.電網(wǎng)技術(shù)，2003(4):31-35.

［8］ILIC''-SPONG M，CROW M L，PAI M A.Transient stability simulation by waveform relaxation methods［J］.Power Systems，IEEE Transactions on，1987，2(4):943-949.

［9］薛巍，舒繼武，嚴(yán)劍峰，等.基于集群機(jī)的大規(guī)模電力系統(tǒng)暫態(tài)過(guò)程并行仿真［J］.中國(guó)電機(jī)工程學(xué)報(bào)，2003，23(8):38-43.

［10］林濟(jì)鏗，李楊春，羅萍萍，等.波形松弛法的電力系統(tǒng)暫態(tài)穩(wěn)定性并行仿真計(jì)算［J］.電工技術(shù)學(xué)報(bào)，2006，21(12):47-53，65.

［11］王錫凡，方萬(wàn)良，杜正春.現(xiàn)代電力系統(tǒng)分析［M］.北京:科學(xué)出版社，2003.

［12］ADOMIAN G，RACH R.On the solution of algebraic equations by the decomposition method［J］.Journal of Mathematical Analysis and Applications，1985，105(1):141-166.

［13］ABBASBANDY S.Improving Newton-Raphson method for nonlinear equations by modified Adomian decomposition method［J］.Applied Mathematics and Computation，2003，145(2-3):887-893.

［14］NOOR M A.New family of iterative methods for nonlinear equations［J］.Applied Mathematics and Computation，2007(190):553-558.

［15］CHUN C.Iterative methods improving newton's method by the decomposition method［J］.Computers＆Mathematics with Applications，2005，50(10-12):1559-1568.

［16］BABOLIAN E，BIAZAR J.Solution of nonlinear equations by modified adomian decomposition method［J］.Applied Mathematics and Computation，2002，132(1):167-172.

［17］蔣耀林.波形松弛方法［M］.北京:科學(xué)出版社，2009.

［18］PRUVOST F，LAURENT-GENGOUX P，MAGOULES F，et al.Accelerated waveform relaxation methods for power systems［C］//2011 International collferen on Electrical and Cohtrol Enginering，2011:2877-2888.

［19］舒繼武，薛巍，鄭緯民.一種電力系統(tǒng)暫態(tài)穩(wěn)定并行計(jì)算的優(yōu)化分區(qū)策略［J］.電力系統(tǒng)自動(dòng)化，2003(19):6-10.

［20］ZECEVIC A I，GACIC N.A partitioning algorithm for the parallel solution of differential-algebraic equations by waveform relaxation［J］.Circuits and Systems I:Fundamental Theory and Applications，IEEE Transactions on，1999，46(4):421-434.

［21］GEORGE K V K.METIS:A Software Package for Partitioning Unstructured Graphs，Partitioning Meshes，and Computing Fill-Reducing Orderings of Sparse Matrices Version4［M］.1998.

［22］BARBARA C G J A R.Using OpenMP:Portable Shared Memory Parallel Programming［M］.Cambridge，MAL:MIT Press，2007.

［23］DAVO T A，NATARAJAN E P.Algorithm 907:KLU a direct sparse solver for circuit simulation problems［J］.ACM Transactions on Mathematical Software，2010，3(37):1-17.

［24］ZIMMERMAN R D，MURILLO S X，NCHEZ C E，et al.MATPOWER:steady-state operations，planning，and analysis tools for power systems research and education［J］.Power Systems，IEEE Transactions on，2011，26(1):12-19.

［25］JIWU S，WEI X，WEIMIN Z.A parallel transient stability simulation for power systems［J］.Power Systems，IEEE Transactions on，2005，20(4):1709-1717.

［26］HOU L，BOSE A.Implementation of the waveform relaxation algorithm on a shared memory computer for the transient stability problem［J］.Power Systems，IEEE Transactions on，1997，12(3):1053-1060.