王啟国，徐 晋，汪可友，周建其，樊 涛

（1.电力传输与功率变换控制教育部重点实验室（上海交通大学），上海市 200240；2.国网浙江省电力有限公司嘉兴供电公司，浙江省 嘉兴市 314000；3.国家电网有限公司数字化工作部，北京市 100032）

0 引言

随着“碳达峰·碳中和”目标的提出，电力系统正在向以新能源为主体、大规模交直流互联的复杂网络转变［1-2］。电力电子设备占比的提高和多控制环节之间交互作用引起的多模态振荡问题也逐渐凸显［3-4］。这种全局多模态振荡不仅涉及中高频段，也涉及新能源场站内部的交互作用，设备级的简化模型［5-6］忽略了开关器件的高频切换过程，场站级的聚合模型［7-9］与详细模型的阻抗特性存在差异［10］。因此，有必要以满足精确模拟开关高频动作过程的小步长进行精细化建模和仿真。

系统详细建模会使得仿真规模急剧增长，为提升大规模系统的仿真效率，一系列并行加速算法被提出。从数学的角度看，目前广泛使用的节点分析法（nodal analysis method，NAM）主要受限于节点电压方程的求解效率。文献［11］提出了并行Crout 分解算法，但未对大规模系统进行效率测试。文献［12］提出一种节点映射结构，将大规模系统转换为多个小型子系统进行求解，但其仿真效率与系统规模存在函数关系。文献［13］提出一种利用图形处理器（graphics processor unit，GPU）加速的LU 分解方法实现了循环的并行化求解。针对该方法的数据依赖性，文献［14］提出了改进的稀疏LU 分解方法。上述方法极大地依赖于系统节点导纳矩阵的稀疏程度，其在电力电子设备高渗透下强耦合、强非线性电力系统中的适用性还有待进一步研究。

从物理的角度看，复杂电力系统可以划分为多个独立求解的子系统，由多个处理器进行并行计算。长传输线解耦［15］方法利用了长传输线的自然解耦特征实现并行仿真，但其分网要求较为严格。在此基础上，文献［16］利用全隐式积分和内插值实现了多速率仿真；文献［17-18］提出了一种多区域戴维南等值（multi-area Thevenin equivalent，MATE）的分网方法，但需要利用集中参数元件来进行分网；文献［19］提出了分网多速率电磁暂态并行仿真方法，在分网的基础上支持2 种及2 种以上仿真步长；文献［20］提出了节点分裂法，提高了分网并行的灵活性；文献［21］进一步提出了基于节点分裂接口的多速率仿真方法，提高了系统仿真精度。上述分网并行求解方法，存在子系统规模均匀分配的问题，需要找到子系统数目和子系统规模之间的效率最优点，才能更好地发挥细粒度并行硬件的优势。

延迟插入法（latency insertion method，LIM）是大规模集成电路中提高设计效率的一种仿真方法［22］，其具有细粒度并行优势。目前，已在集成电路中的温度场分析［23］、强耦合传输线［24］、三相耦合传输线［25］以及电力电子系统［26］中得到应用研究。但受限于LIM 的模型要求［22］，现有研究未对含多样化设备的电力系统进行统一化建模和仿真。

本文在LIM 基础上，针对新能源电力系统中无法直接表示为LIM 拓扑形式的设备提出了支持受控源接入的细粒度建模方法。该方法采用的基本并行单元为三相节点、三相支路，从而可以充分发挥GPU 的计算资源优势，使得新能源电力系统的仿真效率不会随着系统规模的增长而明显变化。

1 LIM 原理

LIM 是大规模集成电路设计中常用的一种快速瞬变仿真方法［22］，要求电路由支路拓扑结构和节点拓扑结构组成，并要求每条支路包含一个电感、每个节点包含一个电容，并在各自的拓扑中产生一个延迟［25］。两种拓扑结构如图1 所示。图中：Ca、Ga、Ia分别为节点a的对地电容、电导和电流源值；Ia1、Ia2、Ia3、Ia4为流出节点a的相邻支路的电流；Lab、Rab、Vab分别为支路ab上的电感、电阻和电压源值。

图1 LIM 的2 种基本拓扑结构Fig.1 Two basic topologies of LIM

采用蛙跳法［27］对电压、电流方程进行离散化，因为蛙跳法求解须错开半个仿真步长，所以得到的结果［25］分别如式（1）和式（2）所示。

式中：n为电流求解的当前仿真步；n-1/2 为电压求解的当前仿真步；Va，n-1/2和Va，n+1/2分别为节点a当前仿真步和下一仿真步的电压；Vb，n+1/2为节点b下一仿真步的电压；h为仿真步长；Ma为与节点a相邻的支路数量；o为与节点a相邻支路的编号；Iao，n为与节点a相邻支路o当前仿真步的电流；Iab，n和Vab，n+1/2分别为支路ab当前仿真步的电流和下一仿真步的电压源值；Ia，n为节点a当前仿真步的电流源值。下文变量下标的仿真步长含义与此类似。

在含多个节点、多条支路的网络中，式（1）和式（2）可以整理为相应的矩阵形式。依据式（1）和式（2）可知，节点电压的更新只与上一仿真步的本地节点电压以及相邻支路的电流有关。因此，各个节点电压的更新互相独立。另外，支路电流的更新也只与上一仿真时步的本地支路电流以及两端的节点电压有关，故各支路电流的更新也互相独立。综上，LIM 具有良好的并行优势。

2 细粒度建模方法

通过将LIM 基本拓扑结构［20］中的支路电压源和节点电流源拓展为受控电压源和受控电流源，以输电网络为上层LIM 建模框架，将电力系统中无法直接表示为LIM 形式的设备以受控源形式接入上层LIM 框架中。新能源电力系统的细粒度建模框架如图2 所示。图中：GT 代表同步发电机及其出口变压器。

图2 新能源电力系统的细粒度建模框架Fig.2 Fine-grained modeling framework of power system with renewable energy

2.1 接口建模

以新能源出口滤波电路中的电容作为解耦元件，可以将系统分为常规动态过程的传统交流电网部分和快速动态过程的电力电子部分。如图2 中接口模型所示，通过获取网侧变换器交流侧的支路电流和变换器相连的网侧支路电流，对接口电容上的电压进行更新，然后将更新的接口电压反馈至两侧子系统。因此，接口模型涉及的是单个电容的电压更新，其差分化模型如式（3）所示。

式中：VCf，n-1/2为接口电容的电压向量；Cif为接口电容值组成的对角矩阵；Igsc，n为网侧变换器交流侧的支路电流向量；Ib，n为与变换器相连的网侧支路电流组成的向量。

网络接口的引入使得快速动态系统与常规动态系统之间的求解相互独立，为2 类子系统的并行求解奠定了基础。

2.2 常规动态系统建模

对常规动态系统来说，其主要组成元件包括同步发电机、变压器、输电网络和负荷，各主要元件建模如下。

1）对于输电网络来说，各条线路采用π 型电路模型，并建模为图1 所示的2 种基本拓扑组合形式［25］。

2）对于负荷来说，采用恒阻抗模型可以等效为电阻-电感（RL）串联形式，对应图2 中负荷模型所示的对地支路，其接地点不参与节点电压的计算，在求解过程中直接置零。

3）对于输电网络中的变压器来说，忽略变压器的励磁回路，用变压器的漏抗串联一个无损耗的理想变压器来模拟变压器，进而基于LIM 支路单元的形式进行求解。

4）对于同步发电机来说，其建模过程如图2 中GT 模型所示，将同步发电机及其出口变压器等效为受控电流源形式，通过其数学模型接入上层LIM建模网络，列写同步发电机电压和磁链方程、电感电压方程，并进行差分化，推导得到的GT 支路输出电流表达式如式（4）所示，具体推导过程见附录A。

式中：IGT，n为由GT 支路的三相电流和同步发电机转子电流组成的向量；VGT，n+1/2为由GT 支路所在节点的三相电压和同步发电机励磁电压组成的向量；A1和A2为GT 支路电流更新的2 个系数矩阵。

2.3 快速动态系统建模

以直驱风电机组并网系统为例，快速动态系统包括直驱永磁同步电机、机侧变换器和控制系统以及网侧变换器和控制系统。

1）直驱式永磁同步电机及其串联电感以GT 支路的形式建模。

2）变换器主电路的建模是通过直流侧电容将机侧变换器和网侧变换器进行划分，并根据开关元件的输入信号将其等效为受控源的形式接入上层LIM 框架中，如图2 中变换器模型所示。

其主电路求解过程与输电线路类似，但相间耦合的存在使其求解略有差异。以网侧变换器为例进行推导，得到的表达式如式（5）所示，具体推导过程见附录B。机侧变换器模型的推导过程与网侧变换器类似。

式中：B1和B2为转换矩阵；A3为系数矩阵，具体定义见附录B；Vgsc，n+1/2=[Vab1，n+1/2Vbc1，n+1/2]T，为网侧变换器交流侧等效受控电压源值组成的向量，其中，Vab1，n+1/2和Vbc1，n+1/2分别为直流电压在网侧变换器交流侧的a、b 两相和b、c 两相之间形成的等效受控电压源值。

另外，以网侧变换器为例，列写其解耦两侧受控源的更新方程如式（6）所示，机侧变换器与其类似。

式中：B3、B4、B5、B6为转换矩阵，用来求解相间等效受控电压源值，其中，Vdc，n-1/2为直流电压；Idcgsc，n+1为网侧变换器交流侧在直流侧的等效受控电流源值；S为由网侧变换器开关控制信号组成的矩阵，S=[SaSbSc]，其中，Sa、Sb、Sc分别为a 相、b 相、c 相上桥臂开关控制信号，当a 相上桥臂开关导通时，Sa为1，下桥臂开关导通时，Sa为0，Sb、Sc定义与此类似。

3）机侧变换器采用附录C 图C1 所示的外环功率控制、内环电流控制方案。

4）网侧变换器采用附录C 图C2 所示的基于电网电压定向的矢量控制技术。

对于直接表示为LIM 基本拓扑的设备均以系统拓扑中存在的电容为解耦元件，通过其底层数学模型接入上层LIM 建模框架。从仿真求解角度看，各类变量求解相互独立，符合GPU 细粒度并行求解的要求。为了满足解耦的有效性，需要选取合适的解耦元件参数和仿真步长。

3 多速率并行仿真算法设计

3.1 网络接口数据交互方法

网络接口的引入使得不同子系统的求解相互独立，可以依据其时间尺度特征差异进行多速率求解。因此，采用插值的方法引入常规动态系统估计值，协助其进行求解。

为更好说明网络接口的数据交互方法，设接口的状态变量为Xf1，仿真步长为h；快速动态系统的状态变量为Xf2，仿真步长为h；常规动态系统的状态变量为Xs，仿真步长为mh，即2 类系统的步长比为m，l为0～m之间的整数，通过插值过程可以得到在n+lh时步Xs的估计值，如式（7）所示。

具体来说，在非同步时刻，只需对快速动态系统进行求解，但需要考虑常规动态系统参数对其计算过程的影响。首先，利用Xs，n和Xf2，n以及Xf1，n-h/2来更新Xf1，n+h/2；进而，利用Xf1，n+h/2来更新Xf2，n+h；然后，利用插值公式得到常规动态系统n+h时步的状态变量估计值，与Xf2，n+h和Xf1，n+h/2共同求解Xf1，n+3h/2。以此类推，直到n+mh的同步时刻，此时，在更新快速动态系统参数的同时也需要对常规动态系统参数进行更新，即利用Xf1，n+mh-h/2更新Xs，n+mh和Xf2，n+mh，后续又进入到非同步时刻，如前所述，进行下一个大仿真步长的更新计算。以此循环，直到仿真结束。

3.2 细粒度并行与多速率仿真算法流程

基于所提的细粒度建模方法设计新能源电力系统的细粒度并行与多速率电磁暂态仿真算法，算法流程如图3 所示。图中：nmax为最大仿真步数，w为自然数。其中，中央处理器（central processing unit，CPU）主要完成模型初始化和系数矩阵预存储，GPU 主要完成快速/常规动态系统、网络接口的变量求解。

图3 细粒度并行与多速率电磁暂态仿真算法流程图Fig.3 Flow chart of fine-grained parallel and multi-rate electromagnetic transient simulation algorithm

在GPU 并行求解时，针对不同模型设计不同核函数与之对应，每个核函数求解所需线程块数量及各线程块所含线程的数量依据其求解的复杂程度决定。例如，对于普通线路的节点电压（或支路电流）的求解来说，由于节点（或支路）之间不存在耦合关系，依据所需求解子系统的节点数（或支路数）分配线程即可；对于控制器输出的更新来说，控制器之间求解相互独立，依据系统控制器的数量来分配线程的数量即可；对于GT 支路电流的更新来说，由于各GT 支路的求解相互独立，每条GT 支路分配一个线程块进行计算即可，其中，每个线程块含有的线程数用来并行执行矩阵的乘法和加法运算。对于矩阵乘法来说，每个线程负责一对元素的乘积计算；对于矩阵加法来说，每个线程负责同位置元素的加法计算。

4 步长选择

4.1 快速动态系统的步长选择

新能源电力系统的多速率仿真效率与快速动态系统的仿真步长密切相关。文献［28］利用李雅普诺夫函数给出了由电导-电感-电容（GLC）所组成电路的稳定性判据；文献［29］从波的角度利用电报方程给出了纯电阻-电感-电容（RLC）网络的数值稳定性判据，以上2 种方案均未考虑网络中存在受控源的情况。文献［30］给出了分区LIM 的一般性稳定性判据，考虑了存在受控源接入的情况，但受控源并没有反映到具体设备的数学模型中。上述稳定性判据均只是针对简单电路进行研究，并没有考虑到复杂电力系统含大量无法被LIM 基本拓扑表达的设备。

本文结合LIM 基本拓扑的结构特征，首先，依据精确模拟开关动作的基本条件初步确定快速动态系统的仿真步长；然后，利用混合数值稳定性判据优化选择快速动态系统的步长。

满足精确模拟开关动作的仿真步长约为开关周期的1%，所以仿真步长初选条件如式（8）所示。

式中：fs为开关频率；hinit为初选仿真步长。

利用混合数值稳定性判据对快速动态系统的步长进行优化选择。首先，对含受控源接入的支路，基于局部离散状态空间模型进行数值稳定性判断。以风电机组变换器为例，根据前文中式（1）至式（3）、式（5）、式（6），联立其机侧、网侧、直流侧、接口及其并网点的电压和电流方程，可以得到其局部系统的方程如式（9）所示。

式中：Vn，n+1/2为新能源并网处节点电压组成的向量；L和R分别为线路电感值、电阻值组成的电感矩阵和电阻矩阵；Cdc为直流电容值；Vpmsg，n+1/2为直驱式永磁同步电机出口电压向量；A4为机侧变换器电流求解的系数矩阵；Imsc，n和Vmsc，n+1/2分别为机侧变换器支路电流和机侧变换器交流侧等效受控电压源值组成的向量，其定义与网侧变换器类似；Q为机侧变换器开关控制信号，其定义与S类似。

对式（9）进行整理，可以得到如下的状态空间表达式：

式中：Xsys，n为全系统状态变量；J为状态矩阵；K为输入矩阵。

推导过程以及状态矩阵J和输入矩阵K的定义如附录D 所示。若要该局部离散时间系统稳定，则状态矩阵J的所有特征值须全部位于单位圆内，即

式中：ρ(J)为状态矩阵J的谱半径。

对不含受控源的节点和支路，通过李雅普诺夫能量函数判断其数值稳定性条件，保障系统整体的互联稳定性，其步长选择与电容电感值相关［30］，有

式中：i=1，2，…，Nn，其中，Nn为系统节点数；p=1，2，…，Nb，i，其中，Nb，i为与节点i相连的支路数；Ci为节点i的对地电容值；Li，p为与节点i相连的支路p上的电感值。

式（12）和式（13）共同构成系统的混合数值稳定性判据。式（12）可依据快速动态系统的模型一致性并行求解多个受控源接入点的局部离散状态矩阵的特征值，式（13）可依据系统拓扑信息并行求解各节点电容及其相邻支路电感对步长的约束条件。

4.2 最优步长比的选择

当常规动态系统的局部截断误差达到快速动态系统可接受的误差上限时，系统最优步长比就可以确定。采用梯形法进行离散化时，快速和常规动态系统的局部截断误差分别如式（14）和式（15）所示。

在多速率仿真情况下，快速动态系统在n+mh时步的误差为：

快速和常规动态系统的三阶导数可以依据各系统中4 个不同时步的参数变量进行计算，如式（17）所示。

步长比的选取应该使常规动态系统的误差与快速动态系统的误差具有相同的阶数［31］，由此可以得到最优步长比m的选择公式，如式（18）所示。

从步长选择条件来看，LIM 仿真步长的选取与拓扑所含电容、电感的值有关，其步长通常小于常规电磁暂态仿真步长，更适用于需要以小步长进行详细建模和仿真的大规模可再生能源电力系统。

5 算例分析

5.1 精度验证

为了说明所提方法的准确性，以附录E 中所示改进39 节点系统为仿真算例，分别设计了基于节点分析法的测试程序以及基于所提方法的测试程序。其中，开关频率为5 kHz。硬件平台参数见附录F。

首先，依据混合数值稳定性判据确定快速动态系统的仿真步长，分析中发现，电容参数对系统数值稳定性的影响大，故以接口电容和仿真步长为未知变量求解状态矩阵J的特征值，得到主导特征值（变化最显著的特征值）的参数空间如图4（a）所示，其中，稳定边界为谱半径为1 的平面与主导特征值参数空间的交点线。根据接口电容值可得到对应的仿真步长，此处50 μF 的接口电容选用2 μs 步长进行仿真。以该步长执行阶段化仿真，基于历史仿真步的变量值确定常规动态系统的仿真步长，依据2 类系统局部截断误差可以确定常规动态系统的仿真步长为8 μs，详细的局部截断误差比较见附录G。

图4 仿真结果Fig.4 Simulation results

考虑节点2 在2 s 时发生三相接地短路故障，2.05 s 时清除，所提多速率方法与节点分析法的仿真结果如图4（b）和图4（c）所示。可以发现，所提算法无论在稳态过程中还是在暂态过程中，都与节点分析法的仿真结果基本一致。

5.2 效率验证

5.2.1 系统规模对效率的影响

为了验证所提仿真方法在效率上的优越性，考虑对不同规模的系统进行仿真效率测试。系统规模的变化主要体现在2 个方面，分别是接入电网的可再生能源数量、网侧系统的节点数和机组数。以含多台风电机组接入的改进39 节点级联系统为例，测试系统的仿真效率，系统拓扑见附录E。相关效率测试结果如图5 所示。为保证模型的一致性，节点分析法也采用开关函数对电力电子设备进行建模，对新能源设备进行详细建模时所能接受的仿真步长在10 μs 左右［32］。因此，节点分析法以10 μs 步长进行仿真，所提方法以2 μs 步长进行单速率仿真，以10 μs 为一个仿真时段，对比2 种方法在不同平台上的计算效率。

图5 本文方法与不同平台下节点分析法的效率对比Fig.5 Efficiency comparison among proposed method and nodal analysis method with different platforms

可以发现，相比于节点分析法，所提方法有效解决了系统规模增长所带来的计算效率显著降低的问题，其仿真效率不会随着系统规模的增大而显著降低，但由于所提方法存在步长约束，其在大规模系统仿真中更具有优势。需要说明的是，仿真效率在系统节点个数为468 时会有明显的变化，因为此时仿真所需的最大线程数超出硬件可用的最大线程数，GPU 会自动串行执行新增的计算任务。

5.2.2 仿真步长对效率的影响

为了验证所提仿真方法采用多速率求解方式在效率上的优越性，考虑对含10 个可再生能源接入的改进39 节点系统进行仿真效率测试。主要考虑对快速和常规动态系统采用不同的仿真步长进行仿真，系统在不同仿真步长组合下的效率测试结果如表1 所示。可以发现，多速率求解方式的引入可以进一步提升系统的仿真效率。

表1 不同仿真步长组合下单步求解时间Table 1 Single-step solution time with different combinations of simulation step sizes

6 结语

本文提出了一种针对新能源电力系统的细粒度并行与多速率电磁暂态仿真方法。一方面，通过对快速和常规动态系统进行细粒度建模，并结合GPU的并行计算网格，解决了系统的均衡划分问题，提升了系统求解的并行度，降低了求解复杂度，实现了各子系统的全局细粒度求解；另一方面，所提方法能将系统中不同时间尺度的设备进行解耦，解决了多速率求解方法存在的分网问题，并基于混合数值稳定性判据及局部截断误差的方法确定了各子系统的仿真步长。在此基础上，利用插值进行了多速率仿真，实现了系统仿真效率的进一步提升。但所提方法的仿真步长相对于常规电磁暂态仿真步长来说较小，会导致仿真计算量增加。因此，在小规模、不关注开关级动态过程的电力系统仿真中是不必要的。它更适用于需要进行详细建模和小步长仿真的大规模新能源场站的振荡分析场合。未来如何提升其仿真效率是可以研究的重要方向。

本文研究得到国网浙江省电力有限公司嘉兴供电公司项目(5211JX230004)的资助，谨此致谢！

附录见本刊网络版（http：//www.aeps-info.com/aeps/ch/index.aspx），扫英文摘要后二维码可以阅读网络全文。