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

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

0 引言

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

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

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

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

1 Adomian分解方法

本節(jié)將對Adomian多項式和基于Adomian級數(shù)分解的迭代算法作簡單介紹。

1.1 Adomian多項式及其計算

考慮非線性方程:

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

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

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

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

而非線性部分N分解為一個函數(shù)項級數(shù):

式中:Ai—依賴于x0，x1，x2，…，xi的Adomian多項式，其形式如下:

由式(2～4)可知:

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

解的部分和為:

且有:

1.2 非線性方程組的Adomian級數(shù)法

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

其中，余項為:

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

其中:

可以得到Adomian多項式的前幾項的具體表達式為:

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

當(dāng)m=0時:

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

當(dāng)m=1時:

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2.2 并行算法流程

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

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

圖1 算法流程

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

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

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

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

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

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

圖2 OpenMP的Fork-join模型

3 算例測試

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

表1 測試算例

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

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

表2 分組結(jié)果

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

圖3 功角誤差曲線

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

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

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

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

圖4 并行加速比

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

4 結(jié)束語

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

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

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

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

［3］武同心，呂曉祥，王建全.Duhamel數(shù)值積分算法在電力系統(tǒng)暫態(tài)穩(wěn)定分析中的應(yīng)用［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］.電網(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］薛巍，舒繼武，嚴劍峰，等.基于集群機的大規(guī)模電力系統(tǒng)暫態(tài)過程并行仿真［J］.中國電機工程學(xué)報，2003，23(8):38-43.

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

［11］王錫凡，方萬良，杜正春.現(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)定并行計算的優(yōu)化分區(qū)策略［J］.電力系統(tǒ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.