王啟國，徐 晉，汪可友，周建其，樊 濤

（1.電力傳輸與功率變換控制教育部重點實驗室（上海交通大學(xué)），上海市 200240；2.國網(wǎng)浙江省電力有限公司嘉興供電公司，浙江省 嘉興市 314000；3.國家電網(wǎng)有限公司數(shù)字化工作部，北京市 100032）

0 引言

隨著“碳達峰·碳中和”目標(biāo)的提出，電力系統(tǒng)正在向以新能源為主體、大規(guī)模交直流互聯(lián)的復(fù)雜網(wǎng)絡(luò)轉(zhuǎn)變［1-2］。電力電子設(shè)備占比的提高和多控制環(huán)節(jié)之間交互作用引起的多模態(tài)振蕩問題也逐漸凸顯［3-4］。這種全局多模態(tài)振蕩不僅涉及中高頻段，也涉及新能源場站內(nèi)部的交互作用，設(shè)備級的簡化模型［5-6］忽略了開關(guān)器件的高頻切換過程，場站級的聚合模型［7-9］與詳細模型的阻抗特性存在差異［10］。因此，有必要以滿足精確模擬開關(guān)高頻動作過程的小步長進行精細化建模和仿真。

系統(tǒng)詳細建模會使得仿真規(guī)模急劇增長，為提升大規(guī)模系統(tǒng)的仿真效率，一系列并行加速算法被提出。從數(shù)學(xué)的角度看，目前廣泛使用的節(jié)點分析法（nodal analysis method，NAM）主要受限于節(jié)點電壓方程的求解效率。文獻［11］提出了并行Crout 分解算法，但未對大規(guī)模系統(tǒng)進行效率測試。文獻［12］提出一種節(jié)點映射結(jié)構(gòu)，將大規(guī)模系統(tǒng)轉(zhuǎn)換為多個小型子系統(tǒng)進行求解，但其仿真效率與系統(tǒng)規(guī)模存在函數(shù)關(guān)系。文獻［13］提出一種利用圖形處理器（graphics processor unit，GPU）加速的LU 分解方法實現(xiàn)了循環(huán)的并行化求解。針對該方法的數(shù)據(jù)依賴性，文獻［14］提出了改進的稀疏LU 分解方法。上述方法極大地依賴于系統(tǒng)節(jié)點導(dǎo)納矩陣的稀疏程度，其在電力電子設(shè)備高滲透下強耦合、強非線性電力系統(tǒng)中的適用性還有待進一步研究。

從物理的角度看，復(fù)雜電力系統(tǒng)可以劃分為多個獨立求解的子系統(tǒng)，由多個處理器進行并行計算。長傳輸線解耦［15］方法利用了長傳輸線的自然解耦特征實現(xiàn)并行仿真，但其分網(wǎng)要求較為嚴(yán)格。在此基礎(chǔ)上，文獻［16］利用全隱式積分和內(nèi)插值實現(xiàn)了多速率仿真；文獻［17-18］提出了一種多區(qū)域戴維南等值（multi-area Thevenin equivalent，MATE）的分網(wǎng)方法，但需要利用集中參數(shù)元件來進行分網(wǎng)；文獻［19］提出了分網(wǎng)多速率電磁暫態(tài)并行仿真方法，在分網(wǎng)的基礎(chǔ)上支持2 種及2 種以上仿真步長；文獻［20］提出了節(jié)點分裂法，提高了分網(wǎng)并行的靈活性；文獻［21］進一步提出了基于節(jié)點分裂接口的多速率仿真方法，提高了系統(tǒng)仿真精度。上述分網(wǎng)并行求解方法，存在子系統(tǒng)規(guī)模均勻分配的問題，需要找到子系統(tǒng)數(shù)目和子系統(tǒng)規(guī)模之間的效率最優(yōu)點，才能更好地發(fā)揮細粒度并行硬件的優(yōu)勢。

延遲插入法（latency insertion method，LIM）是大規(guī)模集成電路中提高設(shè)計效率的一種仿真方法［22］，其具有細粒度并行優(yōu)勢。目前，已在集成電路中的溫度場分析［23］、強耦合傳輸線［24］、三相耦合傳輸線［25］以及電力電子系統(tǒng)［26］中得到應(yīng)用研究。但受限于LIM 的模型要求［22］，現(xiàn)有研究未對含多樣化設(shè)備的電力系統(tǒng)進行統(tǒng)一化建模和仿真。

本文在LIM 基礎(chǔ)上，針對新能源電力系統(tǒng)中無法直接表示為LIM 拓撲形式的設(shè)備提出了支持受控源接入的細粒度建模方法。該方法采用的基本并行單元為三相節(jié)點、三相支路，從而可以充分發(fā)揮GPU 的計算資源優(yōu)勢，使得新能源電力系統(tǒng)的仿真效率不會隨著系統(tǒng)規(guī)模的增長而明顯變化。

1 LIM 原理

LIM 是大規(guī)模集成電路設(shè)計中常用的一種快速瞬變仿真方法［22］，要求電路由支路拓撲結(jié)構(gòu)和節(jié)點拓撲結(jié)構(gòu)組成，并要求每條支路包含一個電感、每個節(jié)點包含一個電容，并在各自的拓撲中產(chǎn)生一個延遲［25］。兩種拓撲結(jié)構(gòu)如圖1 所示。圖中：Ca、Ga、Ia分別為節(jié)點a的對地電容、電導(dǎo)和電流源值；Ia1、Ia2、Ia3、Ia4為流出節(jié)點a的相鄰支路的電流；Lab、Rab、Vab分別為支路ab上的電感、電阻和電壓源值。

圖1 LIM 的2 種基本拓撲結(jié)構(gòu)Fig.1 Two basic topologies of LIM

采用蛙跳法［27］對電壓、電流方程進行離散化，因為蛙跳法求解須錯開半個仿真步長，所以得到的結(jié)果［25］分別如式（1）和式（2）所示。

式中：n為電流求解的當(dāng)前仿真步；n-1/2 為電壓求解的當(dāng)前仿真步；Va，n-1/2和Va，n+1/2分別為節(jié)點a當(dāng)前仿真步和下一仿真步的電壓；Vb，n+1/2為節(jié)點b下一仿真步的電壓；h為仿真步長；Ma為與節(jié)點a相鄰的支路數(shù)量；o為與節(jié)點a相鄰支路的編號；Iao，n為與節(jié)點a相鄰支路o當(dāng)前仿真步的電流；Iab，n和Vab，n+1/2分別為支路ab當(dāng)前仿真步的電流和下一仿真步的電壓源值；Ia，n為節(jié)點a當(dāng)前仿真步的電流源值。下文變量下標(biāo)的仿真步長含義與此類似。

在含多個節(jié)點、多條支路的網(wǎng)絡(luò)中，式（1）和式（2）可以整理為相應(yīng)的矩陣形式。依據(jù)式（1）和式（2）可知，節(jié)點電壓的更新只與上一仿真步的本地節(jié)點電壓以及相鄰支路的電流有關(guān)。因此，各個節(jié)點電壓的更新互相獨立。另外，支路電流的更新也只與上一仿真時步的本地支路電流以及兩端的節(jié)點電壓有關(guān)，故各支路電流的更新也互相獨立。綜上，LIM 具有良好的并行優(yōu)勢。

2 細粒度建模方法

通過將LIM 基本拓撲結(jié)構(gòu)［20］中的支路電壓源和節(jié)點電流源拓展為受控電壓源和受控電流源，以輸電網(wǎng)絡(luò)為上層LIM 建?？蚣埽瑢㈦娏ο到y(tǒng)中無法直接表示為LIM 形式的設(shè)備以受控源形式接入上層LIM 框架中。新能源電力系統(tǒng)的細粒度建?？蚣苋鐖D2 所示。圖中：GT 代表同步發(fā)電機及其出口變壓器。

圖2 新能源電力系統(tǒng)的細粒度建?？蚣蹻ig.2 Fine-grained modeling framework of power system with renewable energy

2.1 接口建模

以新能源出口濾波電路中的電容作為解耦元件，可以將系統(tǒng)分為常規(guī)動態(tài)過程的傳統(tǒng)交流電網(wǎng)部分和快速動態(tài)過程的電力電子部分。如圖2 中接口模型所示，通過獲取網(wǎng)側(cè)變換器交流側(cè)的支路電流和變換器相連的網(wǎng)側(cè)支路電流，對接口電容上的電壓進行更新，然后將更新的接口電壓反饋至兩側(cè)子系統(tǒng)。因此，接口模型涉及的是單個電容的電壓更新，其差分化模型如式（3）所示。

式中：VCf，n-1/2為接口電容的電壓向量；Cif為接口電容值組成的對角矩陣；Igsc，n為網(wǎng)側(cè)變換器交流側(cè)的支路電流向量；Ib，n為與變換器相連的網(wǎng)側(cè)支路電流組成的向量。

網(wǎng)絡(luò)接口的引入使得快速動態(tài)系統(tǒng)與常規(guī)動態(tài)系統(tǒng)之間的求解相互獨立，為2 類子系統(tǒng)的并行求解奠定了基礎(chǔ)。

2.2 常規(guī)動態(tài)系統(tǒng)建模

對常規(guī)動態(tài)系統(tǒng)來說，其主要組成元件包括同步發(fā)電機、變壓器、輸電網(wǎng)絡(luò)和負荷，各主要元件建模如下。

1）對于輸電網(wǎng)絡(luò)來說，各條線路采用π 型電路模型，并建模為圖1 所示的2 種基本拓撲組合形式［25］。

2）對于負荷來說，采用恒阻抗模型可以等效為電阻-電感（RL）串聯(lián)形式，對應(yīng)圖2 中負荷模型所示的對地支路，其接地點不參與節(jié)點電壓的計算，在求解過程中直接置零。

3）對于輸電網(wǎng)絡(luò)中的變壓器來說，忽略變壓器的勵磁回路，用變壓器的漏抗串聯(lián)一個無損耗的理想變壓器來模擬變壓器，進而基于LIM 支路單元的形式進行求解。

4）對于同步發(fā)電機來說，其建模過程如圖2 中GT 模型所示，將同步發(fā)電機及其出口變壓器等效為受控電流源形式，通過其數(shù)學(xué)模型接入上層LIM建模網(wǎng)絡(luò)，列寫同步發(fā)電機電壓和磁鏈方程、電感電壓方程，并進行差分化，推導(dǎo)得到的GT 支路輸出電流表達式如式（4）所示，具體推導(dǎo)過程見附錄A。

式中：IGT，n為由GT 支路的三相電流和同步發(fā)電機轉(zhuǎn)子電流組成的向量；VGT，n+1/2為由GT 支路所在節(jié)點的三相電壓和同步發(fā)電機勵磁電壓組成的向量；A1和A2為GT 支路電流更新的2 個系數(shù)矩陣。

2.3 快速動態(tài)系統(tǒng)建模

以直驅(qū)風(fēng)電機組并網(wǎng)系統(tǒng)為例，快速動態(tài)系統(tǒng)包括直驅(qū)永磁同步電機、機側(cè)變換器和控制系統(tǒng)以及網(wǎng)側(cè)變換器和控制系統(tǒng)。

1）直驅(qū)式永磁同步電機及其串聯(lián)電感以GT 支路的形式建模。

2）變換器主電路的建模是通過直流側(cè)電容將機側(cè)變換器和網(wǎng)側(cè)變換器進行劃分，并根據(jù)開關(guān)元件的輸入信號將其等效為受控源的形式接入上層LIM 框架中，如圖2 中變換器模型所示。

其主電路求解過程與輸電線路類似，但相間耦合的存在使其求解略有差異。以網(wǎng)側(cè)變換器為例進行推導(dǎo)，得到的表達式如式（5）所示，具體推導(dǎo)過程見附錄B。機側(cè)變換器模型的推導(dǎo)過程與網(wǎng)側(cè)變換器類似。

式中：B1和B2為轉(zhuǎn)換矩陣；A3為系數(shù)矩陣，具體定義見附錄B；Vgsc，n+1/2=[Vab1，n+1/2Vbc1，n+1/2]T，為網(wǎng)側(cè)變換器交流側(cè)等效受控電壓源值組成的向量，其中，Vab1，n+1/2和Vbc1，n+1/2分別為直流電壓在網(wǎng)側(cè)變換器交流側(cè)的a、b 兩相和b、c 兩相之間形成的等效受控電壓源值。

另外，以網(wǎng)側(cè)變換器為例，列寫其解耦兩側(cè)受控源的更新方程如式（6）所示，機側(cè)變換器與其類似。

式中：B3、B4、B5、B6為轉(zhuǎn)換矩陣，用來求解相間等效受控電壓源值，其中，Vdc，n-1/2為直流電壓；Idcgsc，n+1為網(wǎng)側(cè)變換器交流側(cè)在直流側(cè)的等效受控電流源值；S為由網(wǎng)側(cè)變換器開關(guān)控制信號組成的矩陣，S=[SaSbSc]，其中，Sa、Sb、Sc分別為a 相、b 相、c 相上橋臂開關(guān)控制信號，當(dāng)a 相上橋臂開關(guān)導(dǎo)通時，Sa為1，下橋臂開關(guān)導(dǎo)通時，Sa為0，Sb、Sc定義與此類似。

3）機側(cè)變換器采用附錄C 圖C1 所示的外環(huán)功率控制、內(nèi)環(huán)電流控制方案。

4）網(wǎng)側(cè)變換器采用附錄C 圖C2 所示的基于電網(wǎng)電壓定向的矢量控制技術(shù)。

對于直接表示為LIM 基本拓撲的設(shè)備均以系統(tǒng)拓撲中存在的電容為解耦元件，通過其底層數(shù)學(xué)模型接入上層LIM 建?？蚣?。從仿真求解角度看，各類變量求解相互獨立，符合GPU 細粒度并行求解的要求。為了滿足解耦的有效性，需要選取合適的解耦元件參數(shù)和仿真步長。

3 多速率并行仿真算法設(shè)計

3.1 網(wǎng)絡(luò)接口數(shù)據(jù)交互方法

網(wǎng)絡(luò)接口的引入使得不同子系統(tǒng)的求解相互獨立，可以依據(jù)其時間尺度特征差異進行多速率求解。因此，采用插值的方法引入常規(guī)動態(tài)系統(tǒng)估計值，協(xié)助其進行求解。

為更好說明網(wǎng)絡(luò)接口的數(shù)據(jù)交互方法，設(shè)接口的狀態(tài)變量為Xf1，仿真步長為h；快速動態(tài)系統(tǒng)的狀態(tài)變量為Xf2，仿真步長為h；常規(guī)動態(tài)系統(tǒng)的狀態(tài)變量為Xs，仿真步長為mh，即2 類系統(tǒng)的步長比為m，l為0～m之間的整數(shù)，通過插值過程可以得到在n+lh時步Xs的估計值，如式（7）所示。

具體來說，在非同步時刻，只需對快速動態(tài)系統(tǒng)進行求解，但需要考慮常規(guī)動態(tài)系統(tǒng)參數(shù)對其計算過程的影響。首先，利用Xs，n和Xf2，n以及Xf1，n-h/2來更新Xf1，n+h/2；進而，利用Xf1，n+h/2來更新Xf2，n+h；然后，利用插值公式得到常規(guī)動態(tài)系統(tǒng)n+h時步的狀態(tài)變量估計值，與Xf2，n+h和Xf1，n+h/2共同求解Xf1，n+3h/2。以此類推，直到n+mh的同步時刻，此時，在更新快速動態(tài)系統(tǒng)參數(shù)的同時也需要對常規(guī)動態(tài)系統(tǒng)參數(shù)進行更新，即利用Xf1，n+mh-h/2更新Xs，n+mh和Xf2，n+mh，后續(xù)又進入到非同步時刻，如前所述，進行下一個大仿真步長的更新計算。以此循環(huán)，直到仿真結(jié)束。

3.2 細粒度并行與多速率仿真算法流程

基于所提的細粒度建模方法設(shè)計新能源電力系統(tǒng)的細粒度并行與多速率電磁暫態(tài)仿真算法，算法流程如圖3 所示。圖中：nmax為最大仿真步數(shù)，w為自然數(shù)。其中，中央處理器（central processing unit，CPU）主要完成模型初始化和系數(shù)矩陣預(yù)存儲，GPU 主要完成快速/常規(guī)動態(tài)系統(tǒng)、網(wǎng)絡(luò)接口的變量求解。

圖3 細粒度并行與多速率電磁暫態(tài)仿真算法流程圖Fig.3 Flow chart of fine-grained parallel and multi-rate electromagnetic transient simulation algorithm

在GPU 并行求解時，針對不同模型設(shè)計不同核函數(shù)與之對應(yīng)，每個核函數(shù)求解所需線程塊數(shù)量及各線程塊所含線程的數(shù)量依據(jù)其求解的復(fù)雜程度決定。例如，對于普通線路的節(jié)點電壓（或支路電流）的求解來說，由于節(jié)點（或支路）之間不存在耦合關(guān)系，依據(jù)所需求解子系統(tǒng)的節(jié)點數(shù)（或支路數(shù)）分配線程即可；對于控制器輸出的更新來說，控制器之間求解相互獨立，依據(jù)系統(tǒng)控制器的數(shù)量來分配線程的數(shù)量即可；對于GT 支路電流的更新來說，由于各GT 支路的求解相互獨立，每條GT 支路分配一個線程塊進行計算即可，其中，每個線程塊含有的線程數(shù)用來并行執(zhí)行矩陣的乘法和加法運算。對于矩陣乘法來說，每個線程負責(zé)一對元素的乘積計算；對于矩陣加法來說，每個線程負責(zé)同位置元素的加法計算。

4 步長選擇

4.1 快速動態(tài)系統(tǒng)的步長選擇

新能源電力系統(tǒng)的多速率仿真效率與快速動態(tài)系統(tǒng)的仿真步長密切相關(guān)。文獻［28］利用李雅普諾夫函數(shù)給出了由電導(dǎo)-電感-電容（GLC）所組成電路的穩(wěn)定性判據(jù)；文獻［29］從波的角度利用電報方程給出了純電阻-電感-電容（RLC）網(wǎng)絡(luò)的數(shù)值穩(wěn)定性判據(jù)，以上2 種方案均未考慮網(wǎng)絡(luò)中存在受控源的情況。文獻［30］給出了分區(qū)LIM 的一般性穩(wěn)定性判據(jù)，考慮了存在受控源接入的情況，但受控源并沒有反映到具體設(shè)備的數(shù)學(xué)模型中。上述穩(wěn)定性判據(jù)均只是針對簡單電路進行研究，并沒有考慮到復(fù)雜電力系統(tǒng)含大量無法被LIM 基本拓撲表達的設(shè)備。

本文結(jié)合LIM 基本拓撲的結(jié)構(gòu)特征，首先，依據(jù)精確模擬開關(guān)動作的基本條件初步確定快速動態(tài)系統(tǒng)的仿真步長；然后，利用混合數(shù)值穩(wěn)定性判據(jù)優(yōu)化選擇快速動態(tài)系統(tǒng)的步長。

滿足精確模擬開關(guān)動作的仿真步長約為開關(guān)周期的1%，所以仿真步長初選條件如式（8）所示。

式中：fs為開關(guān)頻率；hinit為初選仿真步長。

利用混合數(shù)值穩(wěn)定性判據(jù)對快速動態(tài)系統(tǒng)的步長進行優(yōu)化選擇。首先，對含受控源接入的支路，基于局部離散狀態(tài)空間模型進行數(shù)值穩(wěn)定性判斷。以風(fēng)電機組變換器為例，根據(jù)前文中式（1）至式（3）、式（5）、式（6），聯(lián)立其機側(cè)、網(wǎng)側(cè)、直流側(cè)、接口及其并網(wǎng)點的電壓和電流方程，可以得到其局部系統(tǒng)的方程如式（9）所示。

式中：Vn，n+1/2為新能源并網(wǎng)處節(jié)點電壓組成的向量；L和R分別為線路電感值、電阻值組成的電感矩陣和電阻矩陣；Cdc為直流電容值；Vpmsg，n+1/2為直驅(qū)式永磁同步電機出口電壓向量；A4為機側(cè)變換器電流求解的系數(shù)矩陣；Imsc，n和Vmsc，n+1/2分別為機側(cè)變換器支路電流和機側(cè)變換器交流側(cè)等效受控電壓源值組成的向量，其定義與網(wǎng)側(cè)變換器類似；Q為機側(cè)變換器開關(guān)控制信號，其定義與S類似。

對式（9）進行整理，可以得到如下的狀態(tài)空間表達式：

式中：Xsys，n為全系統(tǒng)狀態(tài)變量；J為狀態(tài)矩陣；K為輸入矩陣。

推導(dǎo)過程以及狀態(tài)矩陣J和輸入矩陣K的定義如附錄D 所示。若要該局部離散時間系統(tǒng)穩(wěn)定，則狀態(tài)矩陣J的所有特征值須全部位于單位圓內(nèi)，即

式中：ρ(J)為狀態(tài)矩陣J的譜半徑。

對不含受控源的節(jié)點和支路，通過李雅普諾夫能量函數(shù)判斷其數(shù)值穩(wěn)定性條件，保障系統(tǒng)整體的互聯(lián)穩(wěn)定性，其步長選擇與電容電感值相關(guān)［30］，有

式中：i=1，2，…，Nn，其中，Nn為系統(tǒng)節(jié)點數(shù)；p=1，2，…，Nb，i，其中，Nb，i為與節(jié)點i相連的支路數(shù)；Ci為節(jié)點i的對地電容值；Li，p為與節(jié)點i相連的支路p上的電感值。

式（12）和式（13）共同構(gòu)成系統(tǒng)的混合數(shù)值穩(wěn)定性判據(jù)。式（12）可依據(jù)快速動態(tài)系統(tǒng)的模型一致性并行求解多個受控源接入點的局部離散狀態(tài)矩陣的特征值，式（13）可依據(jù)系統(tǒng)拓撲信息并行求解各節(jié)點電容及其相鄰支路電感對步長的約束條件。

4.2 最優(yōu)步長比的選擇

當(dāng)常規(guī)動態(tài)系統(tǒng)的局部截斷誤差達到快速動態(tài)系統(tǒng)可接受的誤差上限時，系統(tǒng)最優(yōu)步長比就可以確定。采用梯形法進行離散化時，快速和常規(guī)動態(tài)系統(tǒng)的局部截斷誤差分別如式（14）和式（15）所示。

在多速率仿真情況下，快速動態(tài)系統(tǒng)在n+mh時步的誤差為：

快速和常規(guī)動態(tài)系統(tǒng)的三階導(dǎo)數(shù)可以依據(jù)各系統(tǒng)中4 個不同時步的參數(shù)變量進行計算，如式（17）所示。

步長比的選取應(yīng)該使常規(guī)動態(tài)系統(tǒng)的誤差與快速動態(tài)系統(tǒng)的誤差具有相同的階數(shù)［31］，由此可以得到最優(yōu)步長比m的選擇公式，如式（18）所示。

從步長選擇條件來看，LIM 仿真步長的選取與拓撲所含電容、電感的值有關(guān)，其步長通常小于常規(guī)電磁暫態(tài)仿真步長，更適用于需要以小步長進行詳細建模和仿真的大規(guī)?？稍偕茉措娏ο到y(tǒng)。

5 算例分析

5.1 精度驗證

為了說明所提方法的準(zhǔn)確性，以附錄E 中所示改進39 節(jié)點系統(tǒng)為仿真算例，分別設(shè)計了基于節(jié)點分析法的測試程序以及基于所提方法的測試程序。其中，開關(guān)頻率為5 kHz。硬件平臺參數(shù)見附錄F。

首先，依據(jù)混合數(shù)值穩(wěn)定性判據(jù)確定快速動態(tài)系統(tǒng)的仿真步長，分析中發(fā)現(xiàn)，電容參數(shù)對系統(tǒng)數(shù)值穩(wěn)定性的影響大，故以接口電容和仿真步長為未知變量求解狀態(tài)矩陣J的特征值，得到主導(dǎo)特征值（變化最顯著的特征值）的參數(shù)空間如圖4（a）所示，其中，穩(wěn)定邊界為譜半徑為1 的平面與主導(dǎo)特征值參數(shù)空間的交點線。根據(jù)接口電容值可得到對應(yīng)的仿真步長，此處50 μF 的接口電容選用2 μs 步長進行仿真。以該步長執(zhí)行階段化仿真，基于歷史仿真步的變量值確定常規(guī)動態(tài)系統(tǒng)的仿真步長，依據(jù)2 類系統(tǒng)局部截斷誤差可以確定常規(guī)動態(tài)系統(tǒng)的仿真步長為8 μs，詳細的局部截斷誤差比較見附錄G。

圖4 仿真結(jié)果Fig.4 Simulation results

考慮節(jié)點2 在2 s 時發(fā)生三相接地短路故障，2.05 s 時清除，所提多速率方法與節(jié)點分析法的仿真結(jié)果如圖4（b）和圖4（c）所示?？梢园l(fā)現(xiàn)，所提算法無論在穩(wěn)態(tài)過程中還是在暫態(tài)過程中，都與節(jié)點分析法的仿真結(jié)果基本一致。

5.2 效率驗證

5.2.1 系統(tǒng)規(guī)模對效率的影響

為了驗證所提仿真方法在效率上的優(yōu)越性，考慮對不同規(guī)模的系統(tǒng)進行仿真效率測試。系統(tǒng)規(guī)模的變化主要體現(xiàn)在2 個方面，分別是接入電網(wǎng)的可再生能源數(shù)量、網(wǎng)側(cè)系統(tǒng)的節(jié)點數(shù)和機組數(shù)。以含多臺風(fēng)電機組接入的改進39 節(jié)點級聯(lián)系統(tǒng)為例，測試系統(tǒng)的仿真效率，系統(tǒng)拓撲見附錄E。相關(guān)效率測試結(jié)果如圖5 所示。為保證模型的一致性，節(jié)點分析法也采用開關(guān)函數(shù)對電力電子設(shè)備進行建模，對新能源設(shè)備進行詳細建模時所能接受的仿真步長在10 μs 左右［32］。因此，節(jié)點分析法以10 μs 步長進行仿真，所提方法以2 μs 步長進行單速率仿真，以10 μs 為一個仿真時段，對比2 種方法在不同平臺上的計算效率。

圖5 本文方法與不同平臺下節(jié)點分析法的效率對比Fig.5 Efficiency comparison among proposed method and nodal analysis method with different platforms

可以發(fā)現(xiàn)，相比于節(jié)點分析法，所提方法有效解決了系統(tǒng)規(guī)模增長所帶來的計算效率顯著降低的問題，其仿真效率不會隨著系統(tǒng)規(guī)模的增大而顯著降低，但由于所提方法存在步長約束，其在大規(guī)模系統(tǒng)仿真中更具有優(yōu)勢。需要說明的是，仿真效率在系統(tǒng)節(jié)點個數(shù)為468 時會有明顯的變化，因為此時仿真所需的最大線程數(shù)超出硬件可用的最大線程數(shù)，GPU 會自動串行執(zhí)行新增的計算任務(wù)。

5.2.2 仿真步長對效率的影響

為了驗證所提仿真方法采用多速率求解方式在效率上的優(yōu)越性，考慮對含10 個可再生能源接入的改進39 節(jié)點系統(tǒng)進行仿真效率測試。主要考慮對快速和常規(guī)動態(tài)系統(tǒng)采用不同的仿真步長進行仿真，系統(tǒng)在不同仿真步長組合下的效率測試結(jié)果如表1 所示?？梢园l(fā)現(xiàn)，多速率求解方式的引入可以進一步提升系統(tǒng)的仿真效率。

表1 不同仿真步長組合下單步求解時間Table 1 Single-step solution time with different combinations of simulation step sizes

6 結(jié)語

本文提出了一種針對新能源電力系統(tǒng)的細粒度并行與多速率電磁暫態(tài)仿真方法。一方面，通過對快速和常規(guī)動態(tài)系統(tǒng)進行細粒度建模，并結(jié)合GPU的并行計算網(wǎng)格，解決了系統(tǒng)的均衡劃分問題，提升了系統(tǒng)求解的并行度，降低了求解復(fù)雜度，實現(xiàn)了各子系統(tǒng)的全局細粒度求解；另一方面，所提方法能將系統(tǒng)中不同時間尺度的設(shè)備進行解耦，解決了多速率求解方法存在的分網(wǎng)問題，并基于混合數(shù)值穩(wěn)定性判據(jù)及局部截斷誤差的方法確定了各子系統(tǒng)的仿真步長。在此基礎(chǔ)上，利用插值進行了多速率仿真，實現(xiàn)了系統(tǒng)仿真效率的進一步提升。但所提方法的仿真步長相對于常規(guī)電磁暫態(tài)仿真步長來說較小，會導(dǎo)致仿真計算量增加。因此，在小規(guī)模、不關(guān)注開關(guān)級動態(tài)過程的電力系統(tǒng)仿真中是不必要的。它更適用于需要進行詳細建模和小步長仿真的大規(guī)模新能源場站的振蕩分析場合。未來如何提升其仿真效率是可以研究的重要方向。

本文研究得到國網(wǎng)浙江省電力有限公司嘉興供電公司項目(5211JX230004)的資助，謹(jǐn)此致謝！

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