新能源電力系統(tǒng)細(xì)粒度并行與多速率電磁暫態(tài)仿真

王啟國(guó)，徐 晉，汪可友，周建其，樊 濤

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

0 引言

隨著“碳達(dá)峰·碳中和”目標(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)振蕩不僅涉及中高頻段，也涉及新能源場(chǎng)站內(nèi)部的交互作用，設(shè)備級(jí)的簡(jiǎn)化模型［5-6］忽略了開關(guān)器件的高頻切換過程，場(chǎng)站級(jí)的聚合模型［7-9］與詳細(xì)模型的阻抗特性存在差異［10］。因此，有必要以滿足精確模擬開關(guān)高頻動(dòng)作過程的小步長(zhǎng)進(jìn)行精細(xì)化建模和仿真。

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

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

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

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

1 LIM 原理

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

圖1 LIM 的2 種基本拓?fù)浣Y(jié)構(gòu)Fig.1 Two basic topologies of LIM

采用蛙跳法［27］對(duì)電壓、電流方程進(jìn)行離散化，因?yàn)橥芴ㄇ蠼忭氬e(cuò)開半個(gè)仿真步長(zhǎng)，所以得到的結(jié)果［25］分別如式（1）和式（2）所示。

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

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

2 細(xì)粒度建模方法

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

圖2 新能源電力系統(tǒng)的細(xì)粒度建模框架Fig.2 Fine-grained modeling framework of power system with renewable energy

2.1 接口建模

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

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

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

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

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

1）對(duì)于輸電網(wǎng)絡(luò)來(lái)說(shuō)，各條線路采用π 型電路模型，并建模為圖1 所示的2 種基本拓?fù)浣M合形式［25］。

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

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

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

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

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

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

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

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

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

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

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

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

對(duì)于直接表示為L(zhǎng)IM 基本拓?fù)涞脑O(shè)備均以系統(tǒng)拓?fù)渲写嬖诘碾娙轂榻怦钤ㄟ^其底層數(shù)學(xué)模型接入上層LIM 建?？蚣?。從仿真求解角度看，各類變量求解相互獨(dú)立，符合GPU 細(xì)粒度并行求解的要求。為了滿足解耦的有效性，需要選取合適的解耦元件參數(shù)和仿真步長(zhǎng)。

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

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

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

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

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

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

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

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

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

4 步長(zhǎng)選擇

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

5 算例分析

5.1 精度驗(yàn)證

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

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

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

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

5.2 效率驗(yàn)證

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

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

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

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

5.2.2 仿真步長(zhǎng)對(duì)效率的影響

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

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

6 結(jié)語(yǔ)

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

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

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