CN105574809B - 基于矩阵指数的电磁暂态仿真图形处理器并行计算方法 - Google Patents

Abstract

一种基于矩阵指数的电磁暂态仿真图形处理器并行计算方法，在状态分析框架下建立待研究电力系统的整体电磁暂态仿真模型，结合矩阵指数算法的数据并行性以及图形处理器并行计算的性能优势，实现电力系统电磁暂态快速仿真。本发明保留了矩阵指数积分方法良好的数值精度和刚性处理能力，对电力系统元件的非线性环节具有一般性的建模和仿真能力，利用了矩阵指数积分算法的数据高度并行性特点，实现了矩阵指数积分算法的在大规模电力系统电磁暂态仿真领域的高效性。本发明在状态分析框架下，基于矩阵指数运算实现了一般电力系统电磁暂态模型的仿真图形处理器并行计算，提高了矩阵指数电磁暂态仿真方法的计算速度。

Description

基于矩阵指数的电磁暂态仿真图形处理器并行计算方法

技术领域

本发明涉及一种电磁暂态仿真图形处理器并行计算方法。特别是涉及一种适用于电力系统电磁暂态建模的基于矩阵指数的电磁暂态仿真图形处理器并行计算方法。

背景技术

电力系统电磁暂态仿真主要反映系统中电场与磁场的相互影响产生的电气量的变化过程，得到从工频到几十kHz频谱范围内的三相电压电流瞬时值波形，需要采用详细的动态建模和微秒级的仿真步长才能准确刻画，使得系统的仿真规模扩大和计算量增大。与此同时，随着互联电网规模的不断扩大，分布式电源的大规模接入，为了快速准确估测系统运行状态，满足对电网进行更快的安全控制，对实时仿真及快速计算能力提出了严峻挑战，提高电磁暂态仿真程序计算速度成为当务之急。

提高电磁暂态仿真计算速度主要有两种途径，算法设计的改进，如降维近似、不同的积分算法等；并行的硬件环境，如集群计算、图形处理器(Graphic Process Unit，GPU)并行计算等。单纯地采用算法的改进到目前已进入瓶颈，难以取得突破性进展，为此，并行计算成为了解决此类问题的有效途径之一。近年来，图形处理器发展迅猛。图形处理器由于拥有数量众多的计算核心，其计算性能早已超过同时期的中央处理器(Central ProcessingUnit，CPU)，在通用科学计算领域展示出强大的计算潜力。与此同时，CUDA(ComputeUnified Device Architecture，统一计算设备架构)的问世及快速发展，使得图形处理器具备了通用计算的能力，迅速成为开始广泛应用的一种并行计算工具。有学者预测，基于图形处理器的并行计算代表着未来高性能计算的发展趋势。

统一计算设备架构编程模型将中央处理器作为主机端，图形处理器看作设备端。主机端和设备端都有独立的内存，分别称为主机端的内存和设备端的显存。数据可在主机端与设备端之间传递。主机端在中央处理器上执行，设备端在图形处理器上执行。运行在图形处理器上的并行计算函数称为内核函数。内核函数以线程网格的形式组织，每个线程网格由若干个线程块组成，而每个线程块又由若干个线程组成。统一计算设备架构还为矩阵运算提供了面向稠密矩阵的CUBLAS函数库和面向稀疏矩阵的CUSPARSE函数库，两者均是“基础代数数学函数集”接口在图形处理器上的移植，包含了诸如矩阵相乘、矩阵与向量相加相乘、向量内积等基本运算函数，为实现基于图形处理器的电力系统仿真奠定了良好的基础。

电力系统电磁暂态仿真本质上可归结为对动力学系统时域响应的求取，它包括系统本身的数学模型和与之相适应的数值算法。

当前，电力系统电磁暂态仿真基本框架可分为两类，包括节点分析法(NodalAnalysis)和状态变量分析法(State-VariableAnalysis)。基于节点分析框架的电磁暂态仿真方法可概括为先采用某种数值积分方法(通常为梯形积分法)将系统中动态元件的特性方程差分化，得到等效的计算电导与历史项电流源并联形式的诺顿等效电路，此时联立整个电气系统的元件特性方程形成节点电导矩阵，如式(1)所示，对其求解即可得到系统中各节点电压的瞬时值。

Gu＝i (1)

式(1)所示的节点电导矩阵为线性方程组，可使用各种成熟的线性稀疏矩阵算法库进行求解。节点分析法广泛应用于EMTP、PSCAD/EMTDC等专业的电力系统电磁暂态仿真程序中，工程上也称基于节点分析框架的电磁暂态仿真工具为EMTP类程序。节点分析法的主要优势体现在程序实现难度和仿真计算效率方面，但由于式(1)的节点电导方程本身已将数值积分方法与系统模型融为一体，导致EMTP类程序在求解算法选择方面缺乏灵活性与开放性，同时式(1)已不能给出系统本身的特征信息。

状态变量分析框架与节点分析法不同，状态变量分析法属于一般性建模方法(general purpose modeling)，不仅适于电路与电力系统仿真，同样也适于其它工程领域的动力学系统的建模与仿真。Matlab/SimPowerSystems软件是状态变量分析框架下暂态仿真程序的典型代表。与节点分析框架相比，状态方程在模型的计算求解方面具有高度的开放性和灵活性，可方便地选择与问题相适应的数值积分方法，同时能够提供关于系统各种特征的丰富信息(如系统的特征值)，进而能够从全局角度了解系统的动态特性，为各种快速、准确、高效的仿真算法的开发与测试工作提供了便利条件。

应用状态变量分析的基础是形成式(2)所示标准形式的状态-输出方程，此时系统中的电源作为输入u。

在电力系统仿真领域，式(2)可以由改进节点法Modified Nodal Analysis(MNA)通过KCL、KVL等约束关系以及元件伏安特性进行构造而得到MNA模型，再经过一定的正规化处理(regularization)转化而来。MNA模型是形如式(3)的状态-输出方程。

也可以采用一般支路等方法，如Automated State Model Generator(ASMG)方法直接构造得到。基于这些方法得到的电力系统模型能够很容易的与本发明所采用的状态变量分析框架下的电磁暂态仿真程序进行接口。

在数值算法方面，传统数值积分方法可分为显式和隐式两类，不同积分方法所具有的数值稳定性和数值精度各不相同。一般来说，隐式方法处理仿真模型中刚性特征的能力较强。电力系统由于动态过程时间尺度差异较大，系统模型表现出一定刚性，这使得主流电磁暂态软件EMTP类程序采用隐式方法以保证数值稳定性。从计算开销方面来看，隐式方法在每一时步内需求解线性方程组，极大限制了其仿真在大规模系统的能力。与之相对的，传统显式方法无需迭代，在每一时步内的运算量较小，但其有限的数值稳定域使得仿真步长受到约束，综合来看对刚性系统的仿真性能不佳。对于现代电力系统来说，系统中既存在微秒级的电力电子开关动态过程，又存在同步机组的励磁、调速等秒级的机电动态过程，时间尺度差异极大，系统刚性特征十分显著。充分利用状态方程框架在数值算法选择方面的灵活性，结合电力系统电磁暂态仿真的应用场景与特殊需求，发展合适的数值积分方法，是提高电力系统电磁暂态快速仿真计算和应用前景的重要前提。

矩阵指数积分方法(Exponential Integrator)是近年来从应用数学领域发端的一种数值积分方法。它使用矩阵指数算子e^hA精确描述动态系统的线性变化规律，可以准确求解形如

的线性动态系统，并具有刚性处理能力强、计算过程数据并行性高等特点。随着现代电力系统中分布式电源和电力电子装置的大量接入，元件模型复杂，数量众多，网络规模庞大，结构各异。充分利用矩阵指数积分方法良好的数据并行性，基于图形处理器平台进行并行计算成为进行大规模复杂结构电力系统仿真的有效途径。

发明内容

本发明所要解决的技术问题是，提供一种基于矩阵指数的电磁暂态仿真图形处理器并行计算方法。其结合了矩阵指数对于线性动态系统的精确仿真能力和刚性处理能力，利用了其数据高度并行性的特点，实现了电磁暂态仿真的高效性。

本发明所采用的技术方案是：一种基于矩阵指数的电磁暂态仿真图形处理器并行计算方法，在状态分析框架下建立待研究电力系统的整体电磁暂态仿真模型，结合矩阵指数算法的数据并行性以及图形处理器并行计算的性能优势，实现电力系统电磁暂态快速仿真；具体包括如下步骤：

1)在状态分析框架下，建立待研究电力系统整体电磁暂态仿真模型，表示为

其中，x是包含当前时刻系统内所有储能元件与记忆元件状态的状态向量，A为稀疏矩阵，t为时间，Ax和f(x,t)分别表示系统动态特性中的线性和非线性部分，y是仿真使用者所要求的仿真输出向量，输出函数g(x,t)根据仿真研究关注点不同，由使用者任意指定；

2)设定仿真时间T，仿真步长Δt，设定当前时刻为仿真起始时刻t_n＝t₀，依照仿真需要设置仿真初值x_n＝x₀，并由y₀＝g(x₀,t₀)，得到仿真起始时刻输出变量的值y₀，写入输出文件；

3)分别在中央处理器上分配内存空间、图形处理器上分配显存空间，在所述内存空间和显存空间内用于存放数据：稀疏矩阵A、稠密矩阵初始状态变量x₀、中间计算状态向量、常数项以及常数向量，输出数据t_n+1时刻的状态变量x_n+1、输出向量y，在中央处理器上预先处理非线性部分f(x,t)元素内常数数据间的运算，将所需数据由中央处理器内存空间传输到对应的图形处理器显存上，采用完全基于图形处理器的并行计算；

4)从统一计算设备架构中调用CUSPARSE库函数，在图形处理器上并行计算当前时刻t_n的稀疏矩阵向量乘法dx₁←Ax_n，将计算结果数据存放在图形处理器显存空间中对应的区域；

5)调用自定义内核函数，在图形处理器上并行计算当前时刻t_n的更新向量f(x_n,t_n)并实现叠加dx₂←dx₁+f(x_n,t_n)，将计算结果数据存放在图形处理器显存空间中对应的区域；

6)从统一计算设备架构中调用CUBLAS库函数，在图形处理器上并行计算当前时刻t_n的稠密向量乘法将计算结果数据存放在图形处理器显存空间中对应的区域；

7)从统一计算设备架构中调用CUBLAS库函数，在图形处理器上并行计算当前时刻t_n的向量加法x_n+1←x_n+dx₃，得到t_n+1时刻的状态变量x_n+1，将计算结果数据存放在图形处理器显存空间中对应的区域；

8)将结果向量x_n+1拷贝到中央处理器内存空间，由y_n+1＝g(x_n+1,t_n+1)得到t_n+1时刻输出向量的值并写入输出数据文件，更新当前时刻为下一时刻t_n+1＝t_n+Δt，仿真向前推进一个步长；

9)比较当前时刻t_n与仿真时间T，判断仿真是否结束，是则仿真结束；否则返回到步骤4)继续进行计算，依此循环迭代，直到仿真结束。

步骤3)中是以压缩稀疏行格式存储稀疏矩阵A，包含矩阵中的非零元素、非零元素所在行的序号、非零元素所在的列索引，以列为主方式存储稠密矩阵

步骤4)中当调用CUSPARSE库函数时，是调用CUSPARSE库函数中的csrmv格式函数进行dx₁←Ax_n计算。

步骤5)具体包括如下步骤：

(1)根据硬件配置参数，设置与硬件配置参数以及数据规模相对应的内核执行参数，线程数选为32的倍数；

(2)在每一个线程内计算当前步的f(x_n,t_n)内每一元素，为了实现指令优化，使用在中央处理器预先计算的f(x,t)所需的常量；

(3)将显存空间中存放当前步的中间计算向量dx₁与f(x_n,t_n)叠加，实现dx₂←dx₁+f(t_n,x_n)，将计算结果数据存放在图形处理器显存空间中对应的区域。

步骤6)中当调用CUBLAS库函数时，是调用CUBLAS库函数中的gemv格式函数进行计算。

步骤7)中当调用CUBLAS库函数时，是调用CUBLAS库函数中的axpy格式函数进行x_n+1←x_n+dx₃计算。

本发明的基于矩阵指数的电磁暂态仿真图形处理器并行计算方法，保留了矩阵指数积分方法良好的数值精度和刚性处理能力，对电力系统元件的非线性环节具有一般性的建模和仿真能力，利用了矩阵指数积分算法的数据高度并行性特点，实现了矩阵指数积分算法的在大规模电力系统电磁暂态仿真领域的高效性。本发明采用图形处理器并行计算方法进行基于矩阵指数算法的电磁暂态仿真的求解，有较高的计算效率。本发明在状态分析框架下，基于矩阵指数运算实现了一般电力系统电磁暂态模型的仿真图形处理器并行计算，提高了矩阵指数电磁暂态仿真方法的计算速度。

附图说明

图1是测试风电场系统原理图；

图2是本发明基于矩阵指数的电磁暂态仿真图形处理器并行计算方法的流程图；

图3a是系统编号为1风机单元的定子绕组q轴磁变；

图3b是图3a中所示subplot-(a)的局部放大图；

图3c是图3a各版本仿真结果之间的绝对误差；

图4a是系统编号为1风机单元的定子绕组d轴磁变；

图4b是图4a中所示subplot-(b)的局部放大图；

图4c是图4a各版本仿真结果之间的绝对误差；

图5a是系统编号为1风机单元的转子绕组q轴磁变；

图5b是图5a中所示subplot-(c)的局部放大图；

图5c是图5a各版本仿真结果之间的绝对误差；

图6a是系统编号为1风机单元的转子绕组d轴磁变；

图6b是图6a中所示subplot-(d)的局部放大图；

图6c是图6a各版本仿真结果之间的绝对误差；

图7a是系统编号为1风机单元的转速；

图7b是图7a中所示subplot-(e)的局部放大图；

图7c是图7a各版本仿真结果之间的绝对误差。

具体实施方式

下面结合实施例和附图对本发明的基于矩阵指数的电磁暂态仿真图形处理器并行计算方法做出详细说明。

如图2所示，本发明的基于矩阵指数的电磁暂态仿真图形处理器并行计算方法，在状态分析框架下建立待研究电力系统的整体电磁暂态仿真模型，结合矩阵指数算法的数据并行性以及图形处理器并行计算的性能优势，实现电力系统电磁暂态快速仿真；具体包括如下步骤：

1)在状态分析框架下，建立待研究电力系统整体电磁暂态仿真模型，表示为

其中，x是包含当前时刻系统内所有储能元件与记忆元件状态的状态向量，A为稀疏矩阵，t为时间，Ax和f(x,t)分别表示系统动态特性中的线性和非线性部分，Ax和f(x,t)由各元件子系统的电磁暂态仿真模型和各元件子系统之间互联关系确定。y是仿真使用者所要求的仿真输出向量，输出函数g(x,t)根据仿真研究关注点不同，由使用者任意指定；

2)设定仿真时间T，仿真步长Δt，设定当前时刻为仿真起始时刻t_n＝t₀，依照仿真需要设置仿真初值x_n＝x₀，并由y₀＝g(x₀,t₀)，得到仿真起始时刻输出变量的值y₀，写入输出文件；

3)分别在中央处理器上分配内存空间、图形处理器上分配显存空间，在所述内存空间和显存空间内用于存放数据：稀疏矩阵A、稠密矩阵初始状态变量x₀、中间计算状态向量、常数项以及常数向量，输出数据t_n+1时刻的状态变量x_n+1、输出向量y，在中央处理器上预先处理非线性部分f(x,t)元素内常数数据间的运算，将所需数据由中央处理器内存空间传输到对应的图形处理器显存上，采用完全基于图形处理器的并行计算；本发明是以压缩稀疏行(CSR)格式存储稀疏矩阵A，包含矩阵中的非零元素、非零元素所在行的序号、非零元素所在的列索引，以列为主方式存储稠密矩阵

4)从统一计算设备架构(Compute Unified Device Architecture)中调用CUSPARSE库函数，在图形处理器上并行计算当前时刻t_n的稀疏矩阵向量乘法dx₁←Ax_n，将计算结果数据存放在图形处理器显存空间中对应的区域；所述的调用CUSPARSE库函数，是调用CUSPARSE库函数中的csrmv格式函数进行dx₁←Ax_n计算。

5)调用自定义内核函数，在图形处理器上并行计算当前时刻t_n的更新向量f(x_n,t_n)并实现叠加dx₂←dx₁+f(t_n,x_n)，将计算结果数据存放在图形处理器显存空间中对应的区域；具体包括如下步骤：

(1)根据硬件配置参数，设置与硬件配置参数以及数据规模相对应的内核执行参数，线程数选为32的倍数；

(2)在每一个线程内计算当前步的f(x_n,t_n)内每一元素，为了实现指令优化，使用在中央处理器预先计算的f(x,t)所需的常量；

(3)将显存空间中存放当前步的中间计算向量dx₁与f(x_n,t_n)叠加，实现dx₂←dx₁+f(t_n,x_n)，将计算结果数据存放在图形处理器显存空间中对应的区域。

6)从统一计算设备架构(Compute Unified Device Architecture)中调用CUBLAS库函数，在图形处理器上并行计算当前时刻t_n的稠密向量乘法将计算结果数据存放在图形处理器显存空间中对应的区域；所述的调用CUBLAS库函数，是调用CUBLAS库函数中的gemv格式函数进行计算。

7)从统一计算设备架构(Compute Unified Device Architecture)中调用CUBLAS库函数，在图形处理器上并行计算当前时刻t_n的向量加法x_n+1←x_n+dx₃，得到t_n+1时刻的状态变量x_n+1，将计算结果数据存放在图形处理器显存空间中对应的区域；所述的调用CUBLAS库函数，是调用CUBLAS库函数中的axpy格式函数进行x_n+1←x_n+dx₃计算。

8)将结果向量x_n+1拷贝到中央处理器内存空间，由y_n+1＝g(x_n+1,t_n+1)得到t_n+1时刻输出向量的值并写入输出数据文件，更新当前时刻为下一时刻t_n+1＝t_n+Δt，仿真向前推进一个步长；

9)比较当前时刻t_n与仿真时间T，判断仿真是否结束，是则仿真结束；否则返回到步骤4)继续进行计算，依此循环迭代，直到仿真结束。

下面，以测试风电场系统作为算例进行分析，系统结构原理图如图1所示。此风电场包含17台风机单元，并配有本地的功率因数补偿并联电容器，接入风电场馈线网络。

本发明提出的算法流程图如图2所示，现详细说明如下：

1)在状态分析框架下，分别建立待研究电力系统各元件子系统的电磁暂态仿真模型。

A.异步发电机的风力发电机组单元建模；

B.功率因数校正并联电容器的建模；

C.风电场中其他无源元件的建模

其中，x_net是电网络的状态变量，包含网络中独立的电感电流和电容电压，A_net是电网络子系统的状态矩阵；u_net是来自并网非线性元件的输入电流，B_net是对应的输入矩阵；y_net是电网络与并网非线性元件接口处的端口电压(相-相)，C_net是对应的输出矩阵。

D.由元件子系统模型生成研究系统的整体状态方程模型。模型形式为：

式中其中，A和f(x,t)分别是系统动态特性中的线性和非线性部分，如图3a、图3b、图3c所示，x包含子系统的所有状态变量。

2)启动仿真程序，设定仿真时间T，仿真步长Δt。设定当前时刻为仿真起始时刻t_n＝t₀，依照仿真需要设置仿真初值x_n＝x₀，并由y₀＝g(x₀,t₀)，得到仿真起始时刻的输出变量，写入输出文件。分别在主机端和设备端分配内存、显存空间，用于存放输入输出数据以及中间计算变量，将所需数据文件由主机端内存传输到设备端显存上，设置合适的内核函数的执行参数。本实例从初始零时刻以零状态启动，设定仿真时间3s，仿真步长50us；

3)调用CUSPARSE库的CSRMV格式函数在图形处理器平台上进行dx₁←Ax_n计算，以压缩稀疏行(CSR)格式存储稀疏矩阵A，包含矩阵中的非零元素、非零元素所在行的序号、非零元素所在的列索引作为函数参数输入，结果数据写到显存中对应的区域；

4)调用内核函数进行dx₂←dx₁+f(t_n,x_n)运算，定义内核函数具体的实现步骤为：

(1)根据硬件配置参数，设置合适的内核执行参数，选择使用较小线程数的多个线程块，线程数选为32的倍数；

(2)在每一个线程内计算当前步的f(x_n,t_n)内每一元素，并将显存空间中存放上一步的中间计算向量dx与f(x_n,t_n)叠加，实现dx₂←dx₁+f(t_n,x_n)，结果数据存放在对应显存空间；

5)调用CUBLAS库GEMV格式函数进行计算，稠密矩阵以按1索引以列为主方式存储，将结果写到显存中对应的区域；

6)调用CUBLAS库AXPY格式函数实现向量加法x_n+1←x_n+dx₃，得到下一时刻的状态变量x_n+1，结果数据存放在对应的显存空间；

7)将显存结果x_n+1输出到内存中，由y_n+1＝g(x_n+1,t_n+1)得到t_n+1时刻输出向量的值并写入输出数据文件。更新当前时刻为下一时刻t_n+1＝t_n+Δt，仿真向前推进一个步长；

8)比较当前时刻t_n与仿真时间T，判断是否已经抵达仿真结束时刻。若已经达到，则仿真结束；若未达到，则回到步骤3)继续进行计算。依此循环迭代，直到仿真结束。

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

分别将本发明所提出的GPU版本基于矩阵指数的电磁暂态并行计算仿真方法与商业仿真软件MATLAB/SimPowerSystems仿真结果、CPU版本采用GOTOBLAS2库且只在单个图形处理器上运行的仿真结果进行比较。以MATLAB版本仿真结果作为验证正确性的比对结果。

图4a、图4b、图4c、图5a、图5b、图5c、图6a、图6b、图6c、图7a、图7b、图7c为测试风电场系统中编号为1的风机单元的状态变量结果，其中红色线代表MATLAB版本仿真结果、绿色线代表CPU版本仿真结果、蓝色代表GPU版本仿真结果。从局部放大图中可以看出，各版本风机单元的磁链变化仿真结果相差很小。从绝对误差图中可以看出，GPU版与CPU版磁链变化仿真结果基本相同，绝对误差数量级在1e-14左右，达到双精度浮点数运算的舍入误差级别，可忽略不计；GPU版与MATLAB版本、CPU版本与MATLAB版本仿真结果最大绝对误差数量级在1e-2左右；GPU版与CPU版转速的仿真最大误差数量级在1e-15左右；GPU版与MATLAB版本、CPU版本与MATLAB版本仿真结果最大绝对误差数量级在1e-6左右，表明GPU版本、CPU版本仿真程序具有很高的准确性。

表3.不同版本仿真程序的计算时间

不同版本电磁暂态仿真的时间如表3所示。从表3中可以看出，在此测试风电场系统规模下，GPU版本的仿真计算时间为3.11s，程序效率是MATLAB版本、CPU版本的2倍左右，这表明GPU版本仿真程序具有高效性，GPU对基于矩阵指数算法的电磁暂态仿真具有较好的加速效果，对测试例这一较大规模的电力系统具备超实时仿真的条件。大量研究已表明，随着系统规模的不断增加，GPU的加速效果会明显的增加。在此387节点规模的风电场系统下加速比为2倍左右，在更大规模的系统下，有望达到更高的加速效率，对快速实现大规模电力系统的电磁暂态仿真具有重要的意义。

以上算例测试结果证明，本发明提出的基于矩阵指数电磁暂态仿真的GPU并行计算方法的准确性与高效性，对提高电磁暂态仿真程序计算速度提供了一种很好的方法。

Claims (6)

1.一种基于矩阵指数的电磁暂态仿真图形处理器并行计算方法，其特征在于，在状态分析框架下建立待研究电力系统的整体电磁暂态仿真模型，结合矩阵指数算法的数据并行性以及图形处理器并行计算的性能优势，实现电力系统电磁暂态快速仿真；具体包括如下步骤：

1)在状态分析框架下，建立待研究电力系统整体电磁暂态仿真模型，表示为

其中，x是包含当前时刻系统内所有储能元件与记忆元件状态的状态向量，A为稀疏矩阵，t为时间，Ax和f(x,t)分别表示系统动态特性中的线性和非线性部分，y是仿真使用者所要求的仿真输出向量，输出函数g(x,t)根据仿真研究关注点不同，由使用者任意指定；

2)设定仿真时间T，仿真步长Δt，设定当前时刻为仿真起始时刻t_n＝t₀，依照仿真需要设置仿真初值x_n＝x₀，并由y₀＝g(x₀,t₀)，得到仿真起始时刻输出变量的值y₀，写入输出文件；

3)分别在中央处理器上分配内存空间、图形处理器上分配显存空间，在所述内存空间和显存空间内用于存放数据：稀疏矩阵A、稠密矩阵初始状态变量x₀、中间计算状态向量、常数项以及常数向量，输出数据t_n+1时刻的状态变量x_n+1、输出向量y，在中央处理器上预先处理非线性部分f(x,t)元素内常数数据间的运算，将所需数据由中央处理器内存空间传输到对应的图形处理器显存上，采用完全基于图形处理器的并行计算；

4)从统一计算设备架构中调用CUSPARSE库函数，在图形处理器上并行计算当前时刻t_n的稀疏矩阵向量乘法dx₁←Ax_n，将计算结果数据存放在图形处理器显存空间中对应的区域；

5)调用自定义内核函数，在图形处理器上并行计算当前时刻t_n的更新向量f(x_n,t_n)并实现叠加dx₂←dx₁+f(x_n,t_n)，将计算结果数据存放在图形处理器显存空间中对应的区域；

6)从统一计算设备架构中调用CUBLAS库函数，在图形处理器上并行计算当前时刻t_n的稠密向量乘法将计算结果数据存放在图形处理器显存空间中对应的区域；

7)从统一计算设备架构中调用CUBLAS库函数，在图形处理器上并行计算当前时刻t_n的向量加法x_n+1←x_n+dx₃，得到t_n+1时刻的状态变量x_n+1，将计算结果数据存放在图形处理器显存空间中对应的区域；

8)将结果向量x_n+1拷贝到中央处理器内存空间，由y_n+1＝g(x_n+1,t_n+1)得到t_n+1时刻输出向量的值并写入输出数据文件，更新当前时刻为下一时刻t_n+1＝t_n+Δt，仿真向前推进一个步长；

9)比较当前时刻t_n与仿真时间T，判断仿真是否结束，是则仿真结束；否则返回到步骤4)继续进行计算，依此循环迭代，直到仿真结束。

2.权利要求1所述的一种基于矩阵指数的电磁暂态仿真图形处理器并行计算方法，其特征在于，步骤3)中是以压缩稀疏行格式存储稀疏矩阵A，包含矩阵中的非零元素、非零元素所在行的序号、非零元素所在的列索引，以列为主方式存储稠密矩阵

3.权利要求1所述的一种基于矩阵指数的电磁暂态仿真图形处理器并行计算方法，其特征在于，步骤4)中当调用CUSPARSE库函数时，是调用CUSPARSE库函数中的csrmv格式函数进行dx₁←Ax_n计算。

4.权利要求1所述的一种基于矩阵指数的电磁暂态仿真图形处理器并行计算方法，其特征在于，步骤5)具体包括如下步骤：

(1)根据硬件配置参数，设置与硬件配置参数以及数据规模相对应的内核执行参数，线程数选为32的倍数；

(2)在每一个线程内计算当前步的f(x_n,t_n)内每一元素，为了实现指令优化，使用在中央处理器预先计算的f(x,t)所需的常量；

(3)将显存空间中存放当前步的中间计算向量dx₁与f(x_n,t_n)叠加，实现dx₂←dx₁+f(t_n,x_n)，将计算结果数据存放在图形处理器显存空间中对应的区域。

5.权利要求1所述的一种基于矩阵指数的电磁暂态仿真图形处理器并行计算方法，其特征在于，步骤6)中当调用CUBLAS库函数时，是调用CUBLAS库函数中的gemv格式函数进行计算。

6.权利要求1所述的一种基于矩阵指数的电磁暂态仿真图形处理器并行计算方法，其特征在于，步骤7)中当调用CUBLAS库函数时，是调用CUBLAS库函数中的axpy格式函数进行x_n+1←x_n+dx₃计算。