CN107918292B - 面向指数积分的电力电子电路暂态仿真gpu加速方法 - Google Patents

Abstract

本发明提出了一种面向指数积分的电力电子电路暂态仿真GPU加速方法。在开关状态频繁改变的电力电子电路的仿真过程中，对于指数积分算法，需要重新对

矩阵进行求解，这一部分消耗了大量的仿真时间。本发明提出的仿真加速策略在存储器中将曾经计算过的矩阵储存起来，当遇到相同开关状态时直接调用，以大幅节约

矩阵的计算时间；对考虑多重开关动作的线性插值法进行改进，使其仿真步长不变，更适用于矩阵存储；在仿真过程中，对于电气系统利用GPU并行计算资源进行高效求解，而将控制系统以及其他复杂的计算部分保留在CPU上计算。本发明所提方法可以有效地提升电力电子电路的仿真速度，为解决面向电力电子化的电力能源系统的新需求提供了新的思路。

Description

面向指数积分的电力电子电路暂态仿真GPU加速方法

技术领域

本发明涉及一种电力电子电路暂态仿真GPU加速方法。特别是涉及一种适用于电力系统 电磁暂态建模仿真应用的面向指数积分的电力电子电路暂态仿真GPU加速方法。

背景技术

近年来，随着可再生能源的大规模接入及交直流互联电网的运行规模不断增加，同时考 虑到电力电子化电力系统的发展趋势，电力能源系统正面临着前所未有的复杂动态特性的耦 合及其运行控制方面的挑战。而基于元件详细动态特性建模的电力系统电磁暂态仿真因其能 够准确刻画微秒级的系统快动态过程，正逐渐在新能源电力系统的分析、设计与运行等方面 获得更加广泛的重视与应用。同时，电磁暂态仿真在计算效率方面的矛盾也随着电网规模及 复杂性的增加和所关注动态过程时间尺度的延伸而愈发突出，传统的电磁暂态仿真方法已很 难满足要求。

目前，提高电磁暂态仿真的效率主要从以下两个维度考虑：一是通过仿真算法的改进， 如采用模型化简、高效的数值算法等；二是利用并行的硬件环境，如集群计算、FPGA(Field Programmable Gate Array，现场可编程门阵列)、GPU(Graphic Process Unit，中央处理器)等 并行计算资源提高计算效率。由于电磁暂态仿真的计算量随系统规模迅速增大，单纯依靠算 法的改进可能较难取得满意的加速效果。因此，并行计算资源及方法便成为解决电磁暂态仿 真效率问题的重要手段。其中，FPGA的灵活性较强，可以适应较为复杂的系统运算求解过 程，但与此同时也意味着更加繁复的系统设计方式。对于一个复杂的算法或元件模型，首先 要对目标进行分析，再使用高级硬件描述语言设计硬件模块，并对FPGA中大量的硬件电路 进行组合连接。而GPU的计算单元更丰富、峰值性能更强大，更适合处理规模较大但流程简 单的计算问题。此外，GPU的软件开发可以基于C/C++等成熟的编程语言，设计方式更加灵 活简单。近年来，图形处理器GPU的性能得到了飞速的发展，支持的运算也越来越复杂，可 编程性和功能都得到极大地扩展，已不仅仅局限于图形处理领域的应用。GPU由于拥有数量 众多的计算核心，在通用科学计算方面展示出了强大的计算潜力，其性能已超越同时期的 CPU。此外，CUDA(Compute Unified Device Architecture)的问世也推动了GPU在一般通用 计算中的应用。GPU提供的并行计算资源是当前及未来高性能计算的重要发展方向之一。

本发明采用的是CPU/GPU异构计算。程序的计算速度是由CPU与GPU共同决定的，对二者合理的分配任务以提高计算效率是非常必要的。根据CPU和GPU的结构特点，一般 将计算密度较高、逻辑分支简单的大规模并行任务交给GPU处理，而复杂逻辑和事物处理等串行计算由CPU解决，并在此基础上考虑尽可能减少二者之间的数据传递，从而最大限度的利用计算机的计算能力。程序根据编程者的设计在CPU与GPU上交替运行，并互相传输数据。

现代电力系统中，随着系统模型复杂程度的增加，动态过程的多时间尺度特性更加明显， 系统模型表现出较强的刚性。随着电磁暂态仿真的复杂程度不断提升，更多高效的数值算法 不断被提出并获得应用。

不同于EMTP类程序，指数积分方法主要基于状态分析框架，此时系统的仿真模型可以 用以下状态方程描述：

对于该非线性微分方程，一种典型的指数积分方法——指数欧拉法将其推导为如下形式：

其中指数欧拉公式在时间间隔[t_n,t_n+1]内，对式(2)中随时间变化的非线性部分g(t,x(t)) 近似处理为不变量g(t_n,x_n)，并将

族函数

引入积分公式。其中x_n∈R^N是精确值x(t_n) 的数值近似。h_n＝t_n+1-t_n为该步的仿真步长，

是

族函数的第一函数，定义如下：

其中，

族函数的递推关系如式(4)所示：

当前，随着可再生能源的广泛接入和系统电力电子化的演变趋势，电力电子装置已被广 泛地应用到电力系统中，如柔性交流输电技术(FACTS)、高压直流输电技术(HVDC)以及 各种可再生能源的分布式发电系统。针对电力电子器件的仿真已经成为了现阶段电磁暂态仿 真算法的核心问题之一，这不仅是传统积分算法面临的挑战，同时也对指数积分方法提出了 更高要求。在仿真过程中，电力电子开关的快速变化使得系统状态频繁发生改变。对于指数 积分算法，便需要在每次状态变化后重新生成状态矩阵，并对

族函数进行求解，这将会消耗 大量的仿真时间。另一方面，还需要对开关的确切动作时间进行确定，并重新计算开关动作 后的系统状态，此外还会遇到在一个仿真步长中开关多次变化的多重开关动作情况。因此对 于含大量电力电子变流装置的系统，提出高效的仿真策略具有重大的意义。

发明内容

本发明所要解决的技术问题是，提供一种能够提高仿真效率的面向指数积分的电力电子 电路暂态仿真GPU加速方法。

本发明所采用的技术方案是：一种面向指数积分的电力电子电路暂态仿真GPU加速方法， 包括如下步骤：

1)在状态分析框架下，建立待研究电力系统的电磁暂态仿真模型；

2)进行系统初始化，包括：设定仿真时间T，仿真步长ΔT＝h；设定当前时刻为仿真起 始时刻t_n＝t₀，依照仿真需要设置电气系统状态变量的仿真初值x_n＝x₀，并由y₀＝f(t₀,x₀)， 得到仿真起始时刻输出变量的值，写入输出文件；在图形处理器侧分配显存空间；

3)根据电磁暂态仿真模型中的矩阵A以及设定的仿真步长h在主机端计算

矩阵

4)将矩阵A、

矩阵

以及系统状态变量x数据由主机端内存传输到对应的设备端 显存上；

5)将当前状态下的

矩阵

存储在设备端的存储器上，用于当该开关状态再次出现 时直接调用；

6)利用图形处理器并行计算资源，使用指数积分并行算法对电气系统进行求解，得到电 气系统状态量x；

7)将得出的电气系统状态量x传回主机端，求解该时刻控制系统的状态量；

8)通过开关状态向量S判断开关状态在本时步内是否改变，若不变则转至步骤12)，否 则继续下一步；

9)开关状态改变时，采用适于矩阵存储的考虑多重开关动作的改进插值法确定开关动作 时刻，并计算该时刻电气系统状态量；

10)根据开关状态向量S判断该时刻的开关状态在之前的仿真过程中是否出现过，若是 则继续下一步，否则根据该时刻的开关状态对电磁暂态仿真模型中矩阵A的相应元素进行更 新后，返回步骤3)；

11)调用存储器中储存的与该时刻开关状态对应的

矩阵

返回步骤6)；

12)由y_n＝f(t_n,x_n)得到t_n时刻的输出向量并写入输出文件，更新仿真时刻t_n＝t_n+ΔT， 仿真向前推进一个步长；

13)比较当前时刻t_n与仿真时间T，判断是否已经抵达仿真结束时刻，若已经达到，则 仿真结束；若未达到，则返回步骤6)。

步骤1)所述的电力系统电磁暂态仿真模型，是采用如下非线性微分方程组描述：

其中，x为电气系统状态变量，矩阵A代表原微分方程组中的线性部分，含义是状态量之 间的相互耦合关系，包含了模型中的主要信息；而g(t,x)则代表系统状态量耦合的非线性部 分，y是通过系统状态变量x求出的输出量。

步骤6)所述的指数积分并行算法对电气系统进行求解，包括：

(1)在图形处理器上并行实现稀疏矩阵与向量的乘法运算dx←Ax_n；

(2)在图形处理器上并行完成dx←dx+g(t_n,x_n)部分的计算，即g(t_n,x_n)更新，和 g(t_n,x_n)与x的向量加法；

(3)在图形处理器上并行实现稠密矩阵与向量的乘法运算

(4)在图形处理器上并行计算向量加法x_n+1←x_n+dx

其中，x为电气系统状态变量，为矩阵，g(t,x)代表系统状态量耦合的非线性部分，h 为仿真步长。

步骤9)包括：

(1)在t时刻用指数积分算法仿真至t+ΔT，此时通过开关状态向量S的变化检测到开关 动作；

(2)由t与t+ΔT的状态变量线性插值得到电流的过零时刻t_a，t_a＝t+αΔT，

w为PWM变流器调制与载波信号的差，若t_a∈(t,t+ΔT)则求出该时刻电气系统 状态变量x(t_a)＝x(t)+α(x(t+ΔT)-x(t))，继续下一步，否则进入第(6)步；

(3)使用指数积分法仿真一个步长至t_a+ΔT；

(4)继续检测开关动作，若有动作则返回第(2)步，否则继续下一步：

(5)用线性插值法求出t+ΔT时刻电气系统状态变量，x(t+ΔT)＝x(t_a)+(1- α)(x(t_a+ΔT)-x(t_a))；

(6)使用指数积分法仿真至下一步t+2ΔT。

本发明的面向指数积分的电力电子电路暂态仿真GPU加速方法，是一种适合处理电磁暂 态模型刚性特征的数值积分方法，利用指数积分方法的矩阵运算特性，存储生成的矩阵，在 出现相同开关状态时直接调用；对考虑多重开关动作的插值算法针对指数积分和矩阵存储策 略进行相应改进；提出应用于电力电子器件的GPU并行仿真策略，图形处理器并行求解电气 系统，中央处理器串行求解控制系统及其他计算部分。对含电力电子变流装置的电力系统实 现了较好的加速效果。本发明的方法克服了系统刚性对算法数值稳定性的影响，使得仿真过 程可以采用较大的步长进行计算；同时，通过对系统的线性动态过程进行精确求解，避免了 传统数值积分方法所遇到的局部截断误差问题，使算法具有很高的数值精度。此外本发明的 方法尤其适于大规模高维非线性问题，且具有良好的数值性能，仿真计算速度得到明显提升。

附图说明

图1是本发明面向指数积分的电力电子电路暂态仿真GPU加速方法的流程图；

图2是永磁直驱风力发电系统仿真框图；

图3是适用于矩阵存储的考虑多重开关动作的改进插值法示意图；

图4是LCL滤波器电压的仿真结果图；

图5是LCL滤波器电压的仿真结果局部放大图；

图6是LCL滤波器电流的仿真结果图；

图7是LCL滤波器电流的仿真结果局部放大图。

具体实施方式

下面结合实施例和附图对本发明的面向指数积分的电力电子电路暂态仿真GPU加速方 法做出详细说明。

本发明的面向指数积分的电力电子电路暂态仿真GPU加速方法，基于指数积分算法通过 矩阵存储策略以及图形处理器并行计算，综合考虑多重开关动作，最终提高仿真效率，具有 良好的数值精度和稳定性。

如图1所示，本发明的面向指数积分的电力电子电路暂态仿真GPU加速方法，包括如下 步骤：

1)在状态分析框架下，建立待研究电力系统的电磁暂态仿真模型；

所述的电力系统电磁暂态仿真模型，是采用如下非线性微分方程组描述：

其中，x为电气系统状态变量，矩阵A代表原微分方程组中的线性部分，含义是状态量之 间的相互耦合关系，包含了模型中的主要信息；而g(t,x)则代表系统状态量耦合的非线性部 分，y是通过系统状态变量x求出的输出量。

2)进行系统初始化，包括：设定仿真时间T，仿真步长ΔT＝h；设定当前时刻为仿真起 始时刻t_n＝t₀，依照仿真需要设置电气系统状态变量的仿真初值x_n＝x₀，并由y₀＝f(t₀,x₀)， 得到仿真起始时刻输出变量的值，写入输出文件；在图形处理器侧分配显存空间；

3)根据电磁暂态仿真模型中的矩阵A以及设定的仿真步长h在主机端计算

矩阵

4)将矩阵A、

矩阵

以及系统状态变量x数据由主机端内存传输到对应的设备端 显存上；

5)将当前状态下的

矩阵

存储在设备端的存储器上，用于当该开关状态再次出现 时直接调用；在每一种开关状态首次出现时，将计算最耗时的

矩阵求解的结果存入设备端存 储器中，当该开关状态再次出现时直接调用，大大节约了计算时间。此外也减少了主机端与 设备端的数据传输；

6)利用图形处理器并行计算资源，使用指数积分并行算法对电气系统进行求解，得到电 气系统状态量x；

指数积分并行算法对电气系统的求解包含大量的矩阵向量乘法以及向量间加法运算，具 有高度的数据并行性。将该部分利用图形处理器的并行资源进行并行处理，而将其他复杂的 运算环节保留在中央处理器中串行计算。

所述的指数积分并行算法对电气系统进行求解，包括：

(1)在图形处理器上并行实现稀疏矩阵与向量的乘法运算dx←Ax_n；本发明实施例是 调用CUSPARSE库函数，在图形处理器上并行实现稀疏矩阵与向量的乘法运算dx←Ax_n。

(2)在图形处理器上并行完成dx←dx+g(t_n,x_n)部分的计算，即g(t_n,x_n)更新，和 g(t_n,x_n)与x的向量加法；

(3)在图形处理器上并行实现稠密矩阵与向量的乘法运算

本发明实 施例是调用CUBLAS库函数，在图形处理器上并行实现稠密矩阵与向量的乘法运算

(4)在图形处理器上并行计算向量加法x_n+1←x_n+dx，本发明实施例是调用CUBLAS库函数，在图形处理器上并行计算向量加法x_n+1←x_n+dx，

其中，x为电气系统状态变量，为矩阵，g(t,x)代表系统状态量耦合的非线性部分，h 为仿真步长。

7)将得出的电气系统状态量x传回主机端，求解该时刻控制系统的状态量；

8)通过开关状态向量S判断开关状态在本时步内是否改变，若不变则转至步骤12)，否 则继续下一步；

9)开关状态改变时，采用适于矩阵存储的考虑多重开关动作的改进插值法确定开关动作 时刻，并计算该时刻电气系统状态量；如图3所示，包括：

(1)在t时刻用指数积分算法仿真至t+ΔT，此时通过开关状态向量S的变化检测到开关 动作；

(2)由t与t+ΔT的状态变量线性插值得到电流的过零时刻t_a，t_a＝t+αΔT，

w为PWM变流器调制与载波信号的差，若t_a∈(t,t+ΔT)则求出该时刻电气系统 状态变量x(t_a)＝x(t)+α(x(t+ΔT)-x(t))，继续下一步，否则进入第(6)步；

(3)使用指数积分法仿真一个步长至t_a+ΔT；

(4)继续检测开关动作，若有动作则返回第(2)步，否则继续下一步：

(5)用线性插值法求出t+ΔT时刻电气系统状态变量，x(t+ΔT)＝x(t_a)+(1- α)(x(t_a+ΔT)-x(t_a))；

(6)使用指数积分法仿真至下一步t+2ΔT。

上述过程在仿真中不改变步长，只要通过开关向量判断出系统状态相同，那么由此状态 矩阵计算生成的

矩阵也完全相同，因此可以直接利用先前存储的矩阵，减少计算开销大的

函数运算，可以大幅提高仿真效率。

10)根据开关状态向量S判断该时刻的开关状态在之前的仿真过程中是否出现过，若是 则继续下一步，否则根据该时刻的开关状态对电磁暂态仿真模型中矩阵A的相应元素进行更 新后，返回步骤3)；

11)调用存储器中储存的与该时刻开关状态对应的

矩阵

返回步骤6)；

12)由y_n＝f(t_n,x_n)得到t_n时刻的输出向量并写入输出文件，更新仿真时刻t_n＝t_n+ΔT， 仿真向前推进一个步长；

13)比较当前时刻t_n与仿真时间T，判断是否已经抵达仿真结束时刻，若已经达到，则 仿真结束；若未达到，则返回步骤6)。

本发明实施例中，执行仿真计算的计算机硬件环境为Inter(R)Xeon E5-2623v3中央处理 器，内存容量32GB、Tesla K20C图形处理器，SM数为13，CUDA核心数为2496，CUDA 核心频率为706MHz，全局储存器容量4.67GB；软件环境为Linux14.04操作系统；使用CUDA 提供的面向稠密矩阵的CUBLAS库函数和面向稀疏矩阵的CUSPARSE库函数。

将理想电源直接接入5个永磁同步直驱风力发电系统，以此作为算例验证本发明的面向 指数积分的电力电子电路暂态仿真GPU加速方法。直驱风力发电系统采用永磁同步电机，采 用全功率双PWM变流器与电网相连作为并网结构，每台永磁同步电机包含6个电力电子开 关。永磁直驱风力发电系统仿真框图如图2所示。在dq0坐标系下，永磁同步电机电压方程 为：

磁链方程为：

电磁转矩方程为：

T_e＝i_sqψ_sd-i_sdψ_sq

将本发明的面向指数积分的电力电子电路暂态仿真GPU加速方法，与MATLAB版本的 指数积分一般性串行方法进行对比。LCL滤波器电压电流的仿真结果及其局部放大图如图4～ 图7所示，其中实线为MATLAB仿真结果，虚线为GPU仿真结果。从图中可以看出两者具有高度的一致性，说明文中所用方法可以有效地保证仿真精度。

表1给出了两种算法的计算时间对比。其中MATLAB版程序仿真耗时237.46s，而本发 明的面向指数积分的电力电子电路暂态仿真GPU加速方法仅需110.30s，可以达到2.15的加 速比，有效提升了计算效率。表2给出本发明的面向指数积分的电力电子电路暂态仿真GPU 加速方法，在0.2s的仿真过程里的矩阵存储情况。可以看出在电力电子元件的仿真中，系统 状态有着很高的重复率。对于指数积分算法，将曾经计算过的矩阵存储起来以便在遇到相同 开关状态时直接调用，这样可以有效地节约

矩阵的计算时间，从而大幅提升计算效率。

表1MATLAB与GPU程序计算时间

程序实现	计算时间(s)	加速比/GPU
			MATLAB	237.46	2.15
GPU	110.30	/

表2矩阵存储情况

风机台数N	存储矩阵数量	矩阵重复使用次数
			5	7035	50085

Claims (3)

1.一种面向指数积分的电力电子电路暂态仿真GPU加速方法，其特征在于，包括如下步骤：

1)在状态分析框架下，建立待研究电力系统的电磁暂态仿真模型；

2)进行系统初始化，包括：设定仿真时间T，仿真步长ΔT＝h；设定当前时刻为仿真起始时刻t_n＝t₀，依照仿真需要设置电气系统状态变量的仿真初值x_n＝x₀，并由y₀＝f(t₀，x₀)，得到仿真起始时刻输出变量的值，写入输出文件；在图形处理器侧分配显存空间；

3)根据电磁暂态仿真模型中的矩阵A以及设定的仿真步长h在主机端计算

矩阵

4)将矩阵A、

矩阵

以及系统状态变量x数据由主机端内存传输到对应的设备端显存上；

5)将当前状态下的

矩阵

存储在设备端的存储器上，用于当当前状态下的开关状态再次出现时直接调用；

6)利用图形处理器并行计算资源，使用指数积分并行算法对电气系统进行求解，得到电气系统状态量x；

7)将得出的电气系统状态量x传回主机端，求解该时刻控制系统的状态量；

8)通过开关状态向量S判断开关状态在本时步内是否改变，若不变则转至步骤12)，否则继续下一步；

9)开关状态改变时，采用适于矩阵存储的考虑多重开关动作的改进插值法确定开关动作时刻，并计算该时刻电气系统状态量；包括：

(1)在t时刻用指数积分算法仿真至t+ΔT，此时通过开关状态向量S的变化检测到开关动作；

(2)由t与t+ΔT的状态变量线性插值得到电流的过零时刻t_a，t_a＝t+αΔT，

w为PWM变流器调制与载波信号的差，若t_a∈(t,t+ΔT)则求出该时刻电气系统状态变量x(t_a)＝x(t)+α(x(t+ΔT)-x(t))，继续下一步，否则进入第(6)步；

(3)使用指数积分法仿真一个步长至t_a+ΔT；

(4)继续检测开关动作，若有动作则返回第(2)步，否则继续下一步：

(5)用线性插值法求出t+ΔT时刻电气系统状态变量，x(t+ΔT)＝x(t_a)+(1-α) ( x( t_a+ΔT) -x( t_a) ) ；

(6)使用指数积分法仿真至下一步t+2ΔT；

10)根据开关状态向量S判断该时刻的开关状态在之前的仿真过程中是否出现过，若是则继续下一步，否则根据该时刻的开关状态对电磁暂态仿真模型中矩阵A的相应元素进行更新后，返回步骤3)；

11)调用存储器中储存的与该时刻开关状态对应的

矩阵

返回步骤6)；

12)由y_n＝f(t_n,x_n)得到t_n时刻的输出向量并写入输出文件，更新仿真时刻t_n＝t_n+ΔT，仿真向前推进一个步长；

13)比较当前时刻t_n与仿真时间T，判断是否已经抵达仿真结束时刻，若已经达到，则仿真结束；若未达到，则返回步骤6)。

2.根据权利要求1所述的面向指数积分的电力电子电路暂态仿真GPU加速方法，其特征在于，步骤1)所述的电力系统电磁暂态仿真模型，是采用如下非线性微分方程组描述：

其中，x为电气系统状态变量，矩阵A代表原微分方程组中的线性部分，含义是状态量之间的相互耦合关系，包含了模型中的主要信息；而g(t,x)则代表系统状态量耦合的非线性部分，y是通过系统状态变量x求出的输出量。

3.根据权利要求1所述的面向指数积分的电力电子电路暂态仿真GPU加速方法，其特征在于，步骤6)所述的指数积分并行算法对电气系统进行求解，包括：

(1)在图形处理器上并行实现稀疏矩阵与向量的乘法运算dx←Ax_n；

(2)在图形处理器上并行完成dx←dx+g(t_n,x_n)部分的计算，即g(t_n,x_n)向量更新和g(t_n,x_n)与x的向量加法；

(3)在图形处理器上并行实现稠密矩阵与向量的乘法运算

(4)在图形处理器上并行计算向量加法x_n+1←x_n+dx

其中，x为电气系统状态变量，A为矩阵，g(t,x)代表系统状态量耦合的非线性部分，h为仿真步长。

Images