针对电力电子混杂系统仿真的离散状态事件驱动仿真方法与流程

本发明涉及电力电子系统仿真技术领域，特别是涉及针对电力电子混杂系统仿真的离散状态事件驱动仿真方法。

背景技术：

电力电子系统本质上是一种包含连续状态与离散事件的混杂动态系统。基于电力电子系统的混杂特性，一个准确与高效的电力电子系统仿真方法应该具有以下特点：

1)连续状态的灵活积分：阶数与步长均可以根据系统状态的变化自适应地调整，且该自适应调整的机制应能够完全嵌入于数值积分的过程中，从而避免产生额外的计算量；

2)离散事件的高效定位：高效的事件定位机制应能够有效区分电力电子系统中的主动事件与被动事件，并根据其特点分别采取不同的定位方法，以尽可能降低计算量。

目前虽然已有不少电力电子系统的商业仿真软件，例如simulink与plecs，但它们仍然在以下方面表现出了其低效性。在连续状态的积分方面，这些商业软件仅提供了步长可变但阶数固定的数值积分算法，例如：经典的dormand-prince算法(即simulink中所采用的ode45算法与plecs中所采用的dopri算法)。由于这些算法的阶数固定，因此其数值积分的灵活性与高效性有限。在离散事件的定位方面，由于仿真计算仍然基于连续状态的时间离散数值积分，所有的离散事件均被视为被动事件进行处理，忽略了提前计算主动事件发生时刻从而避免迭代计算的可能性。此外，由于系统状态不能被表示为时间的显式函数，因此被动事件的定位过程会涉及整个数值积分的过程，导致该过程的计算量较大。

因此希望有一种对电力电子混杂系统仿真的离散状态事件驱动仿真方法以解决现有技术中的问题。

技术实现要素：

本发明公开了一种针对电力电子混杂系统仿真的离散状态事件驱动仿真方法，其基于电力电子系统的混杂特性，针对电力电子系统中共存的连续状态与离散事件分别设计了仿真方法，所述仿真方法包括：阶数与步长均可变的灵活自适应算法和由离散事件自动触发仿真计算的事件驱动机制；其中连续状态以状态变量法进行描述，离散事件被定义为仿真计算的触发源。

优选地,所述连续状态的建模包括：在第k个仿真间隔内，电力电子系统被建模为一个线性时不变系统，该系统通过公式(1)的状态方程，公式(2)的输出方程以及公式(3)的阈值方程建模，

y(t)＝ckx(t)+dku(t)tk＜t≤tk+1(2)

y2(t)＝cth(3)

其中，x为系统的独立状态变量，y为系统的输出变量，u为系统输入变量，a,b,c,d是由系统拓扑与参数决定的系统方程的系数矩阵，输出变量y包含用于仿真波形显示的显示向量y1和能触发被动事件的触发向量y2，当y2的任意一个分量越过其在阈值向量cth中对应的阈值分量时，会检测到一次被动事件，系统方程的系数矩阵被写成开关器件开关状态向量swk和系统参数向量pk的函数，如公式(4)所示：

mk＝m(swk,pk)m＝a,b,c,d(4)。

优选地,所述离散事件包括四类：

(i)状态事件：当某一步的数值积分误差随着步长的增大而达到用户所设定的误差上限时，发生一次状态事件；

(ii)控制事件：在电力电子系统中，由数字控制器触发的离散事件为控制事件；

(iii)外部事件：用户仿真开始前预定义的离散事件为外部事件；

(iv)被动事件：发生时刻不可提前获知的离散事件为被动事件；

其中(i)、(ii)和(iii)类事件为主动事件，其发生时刻可以在其发生之前提前获知。

优选地,所述仿真方法包括以下步骤：

步骤(1)仿真准备工作与系统状态初始化：完成系统参数输入、方程列写、外部事件设置与系统状态初始化，该步骤(1)结束后进入步骤(2)-(5)的仿真架构主循环；

步骤(2)系统状态更新与系统事件处理：首先更新系统状态，然后判断是否发生系统事件，若发生控制事件与/或外部事件，则根据事件内容更新器件开关状态与/或电路元件参数；

步骤(3)系数矩阵更新与仿真波形输出：若器件的开关状态与电路参数均未变化，则保持系数矩阵不变，计算输出方程并进行波形输出；若器件的开关状态或电路参数变化，则在更新系数矩阵之前计算输出方程并完成第一次波形输出，然后根据最新的器件开关状态更新系数矩阵，随后再次计算输出方程并完成第二次波形输出；

步骤(4)系统状态的灵活自适应数值积分：首先通过事件排序给出积分步长的初始估计，随后计算状态向量的各阶导数与增量，其次进行精度检测，然后进行自适应的阶数选取与步长调整，最后进行数值积分；

步骤(5)被动事件定位：若有触发变量越过阈值条件，则发生被动事件，此时采用割线法定位被动事件发生的时刻，同时缩短积分步长以满足误差要求，并重新进行数值积分；

步骤(6)仿真进程判断与推进：判断下一仿真时刻是否已经超过仿真终止时刻，若未超过，则跳转至步骤(2)重复步骤(2)-(5)的仿真架构主循环；若已超过，则以仿真终止时刻为下一仿真时刻，完成该时刻的系统状态计算与波形输出后，结束仿真进程。

优选地，所述步骤(1)具体包括以下步骤：

步骤(1.1)所述系统参数输入是输入待仿真电路的主电路元件参数与控制器参数；

步骤(1.2)所述系统方程列写是根据电路拓扑、元件参数与器件开关状态组合列写电路的状态方程、输出方程与阈值方程，如公式(1)-(4)所示；

步骤(1.3)所述外部事件设置是预先设置外部事件发生的时间轴及其相应的事件内容；

步骤(1.4)所述系统状态初始化对状态变量、器件开关状态与控制器存储器进行初始化，步骤(1.4)结束后进入所述步骤(2)-(5)的仿真架构主循环。

优选地,所述步骤(2)具体包括以下步骤：

步骤(2.1)所述系统状态更新的方式如公式(5)所示：

步骤(2.2)系统事件检测与处理，若发生所述控制事件，则根据控制器存储的信息进行采样或开关操作，更新主动器件的开关状态向量swact,k，并计算下一次控制事件的内容与其发生的时间间隔δtctr；若发生所述外部事件，则根据其内容更新所述系统参数向量pk，并计算下一次外部事件发生的时间间隔δtext。

优选地，所述步骤(3)具体包括以下步骤：判断所述器件开关状态向量swk与所述系统参数向量pk是否发生变化，若swk与pk均未变化，则保持系数矩阵不变；若swk或pk变化，则首先利用第k-1步的系数矩阵c1,k-1与d1,k-1计算输出方程，并完成第一次波形输出，然后根据最新的器件开关状态向量swk与系统参数向量pk更新系数矩阵，如公式(4)所示，随后再次计算输出方程，并完成所述第二次波形输出。

优选地，所述步骤(4)具体包括以下步骤：

步骤(4.1)所述积分步长的初始估计以当前时刻距下一次控制事件发生的时间间隔δtk,ctr与距下一次外部事件发生的时间间隔δtk,ext中的较小者为积分步长的初始估计δtk，同时，δtk不应超过用户所设置的最大步长δtmax；

步骤(4.2)所述导数与增量计算，在第k步数值积分过程中，可将被仿真系统视为一个线性时不变系统，其系统方程的系数矩阵ak,bk,ck,dk为与仿真时间t和状态向量x均无关的常矩阵，状态向量的各阶导数按照公式(6)所示的递归格式依次求得：

x⁽ⁱ⁾(t)＝akx^(i-1)(t)+bku^(i-1)(t)i≥1(6)

其中，x⁽ⁱ⁾(t)和u⁽ⁱ⁾(t)分别为状态向量与输入向量的i阶导数；一般情况下，在电力电子系统中，由于输入向量一般由显式的独立电源组成(例如：含谐波的交流电源与直流电源)，因此输入向量u⁽ⁱ⁾(t)(i≥0)的各阶导数容易显式获取。

利用状态向量的各阶导数，基于taylor展开的数值积分方法，t＝tk+1时刻的系统状态数值解如公式(7)所示：

定义数值积分过程中，状态向量的第i阶增量为公式(8)

步骤(4.3)所示精度检测，判断公式(9)是否成立：

δixk≤abstolk(9)

其中，abstolk为所设定的绝对误差上限，若公式(9)成立，则为情况1，跳至步骤(4.4)；若公式(9)不成立则将所分析的阶数提高一阶，跳回至步骤(4.2)；一般情况下，为防止所分析的阶数升至过高，会要求用户在仿真开始前设定所分析阶数的上限值qmax。若所分析的阶数已达到qmax，则为情况2，跳至步骤(4.4)；

步骤(4.4)所述阶数选取与步长调整，阶数选择方法如公式(10)所示，步长调整方法如式(11)所示：

步骤(4.5)所述数值积分，数值积分方法如式(12)所示：

对于一个包含n个分量的状态向量xk＝[xk,1,xk,2,…,xk,n]^t，步骤(4.3)-(4.4)中的操作将被作用于状态向量的每一个分量，只有当每一个分量均满足精度要求时，才能通过步骤(4.3)中的精度检测。

优选地，所述步骤(5)具体包括以下步骤：

步骤(5.1)阈值条件判断，被动事件是否发生，可以根据是否存在触发变量越过阈值条件进行判断，定义如公式(13)所示的阈值函数gth(δt)：

gth(t)＝y2(t)-cth(13)

则某一触发变量在该积分步长内越过阈值条件，等价于其对应的阈值函数在该积分区间两端处的值异号，若存在触发变量越过阈值条件，则可知有被动事件发生，跳至步骤(5.2)，反之则可知被动事件未发生，跳至步骤(5.5)；

步骤(5.2)采用割线法定位被动事件的发生时刻，由步骤(5.1)可知触发向量中的某一个或某几个分量达到了阈值条件，假设积分步长足够短，每个触发变量在一个积分区间内仅能达到一次阈值条件，由于状态向量和输入向量在t＝tk时刻的各阶导数已在步骤(4)中求出，被检测到的触发变量的各阶导数由公式(14)求出：

基于公式(14)，阈值函数gth(δt)可被表示为时间显示多项式函数公式(15)：

采用割线法定位阈值函数gth(δt)的零点若触发向量中有多个分量j1,j2,…在同一积分区间内达到阈值条件，则取其中的最小者为当前时刻距下一次被动事件发生时刻的时间间隔δtk,pas；

步骤(5.3)确定开关状态。根据事件内容改变发生被动事件的电路元件的状态；

步骤(5.4)数值积分，若有被动事件发生，则意味着此时的数值积分步长δtk相较步骤(4)中所确定的积分步长被缩短，因此步骤(4)中的数值积分结果已不能继续使用，需要重新进行数值积分；

步骤(5.5)计算下一仿真时刻。按照公式(16)计算下一仿真时刻：

tn,k＝tp,k+δtk,pas(16)。

本发明公开的针对电力电子混杂系统仿真的离散状态事件驱动仿真方法在保证精度的前提下，解决了复杂电力电子系统仿真计算量大、耗时长的问题，本发明的有益效果为：

1、所提dsed仿真方法基于电力电子系统的混杂特性，针对系统中共存的连续状态与离散事件分别设计了不同的仿真方法，既能够实现连续状态的灵活积分，又能够实现离散事件的高效定位，可在保证仿真精度的前提下显著降低仿真计算量，极大地提高了电力电子系统的仿真效率。

2、所提dsed仿真方法中的fa算法能够同时自适应地调整数值积分的阶数与步长，以最小的计算代价满足用户所设定的精度要求，且该自适应调整方法完全嵌入于数值积分的过程中，不产生额外的计算量。

3、所提dsed仿真方法中的ed机制能够区分电力电子系统中的主动事件与被动事件的定位问题。对于主动事件采用提前计算其发生时刻的方式进行定位，从而避免了迭代计算，显著降低了主动事件定位的计算量。

4、所提dsed仿真方法中的ed机制采用所提的割线法定位被动事件。由于系统状态的各阶导数已经在fa算法中计算出，用于事件定位的迭代计算仅涉及一个关于时间的显式多项式函数(阈值函数)，而不涉及整个数值积分过程，显著降低了被动事件定位的计算量。

附图说明

图1是本发明针对电力电子混杂系统仿真的离散状态事件驱动仿真方法流程示意图。

图2是本发明中被动事件定位原理图((a)二极管的关断事件定位原理图，(b)二极管的开通事件定位原理图)示意图。

图3是用于仿真算例研究的功率放大器的电路拓扑图与测试工况图((a)被测功率放大器的电路拓扑图，(b)测试过程中发生突变的参考波)。

图4是图3所示仿真算例的仿真结果与实验波形对比图。

图5是图3所示仿真算例的仿真性能测试曲线图。

具体实施方式

为使本发明实施的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行更加详细的描述。在附图中，自始至终相同或类似的标号表示相同或类似的元件或具有相同或类似功能的元件。所描述的实施例是本发明一部分实施例，而不是全部的实施例。下面通过参考附图描述的实施例是示例性的，旨在用于解释本发明，而不能理解为对本发明的限制。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

针对电力电子系统仿真的离散状态事件驱动(dsed)仿真方法，基于电力电子系统的混杂特性，针对系统中共存的连续状态与离散事件分别设计了不同的仿真方法：1)为实现系统连续状态的灵活积分，设计了一种阶数与步长均可变的灵活自适应算法(flexibleadaptivealgorithm,faalgorithm)；2)为实现系统离散事件的高效定位，设计了一套由离散事件自动触发仿真计算的事件驱动机制(event-drivenmechanism,edmechanism)。

所述仿真方法如图1所示包括以下步骤：

步骤(1)仿真准备工作与系统状态初始化：完成系统参数输入、方程列写、外部事件设置与系统状态初始化，该步骤(1)结束后进入步骤(2)-(5)的仿真架构主循环；

步骤(2)系统状态更新与系统事件处理：首先更新系统状态，然后判断是否发生系统事件，若发生控制事件与/或外部事件，则根据事件内容更新器件开关状态与/或电路元件参数；

步骤(3)系数矩阵更新与仿真波形输出：若器件的开关状态与电路参数均未变化，则保持系数矩阵不变，计算输出方程并进行波形输出；若器件的开关状态或电路参数变化，则在更新系数矩阵之前计算输出方程并完成第一次波形输出，然后根据最新的器件开关状态更新系数矩阵，随后再次计算输出方程并完成第二次波形输出；

步骤(4)系统状态的灵活自适应数值积分：首先通过事件排序给出积分步长的初始估计，随后计算状态向量的各阶导数与增量，其次进行精度检测，然后进行自适应的阶数选取与步长调整，最后进行数值积分；

步骤(5)被动事件定位：若有触发变量越过阈值条件，则发生被动事件，此时采用割线法定位被动事件发生的时刻，同时缩短积分步长以满足误差要求，并重新进行数值积分；

步骤(6)仿真进程判断与推进：判断下一仿真时刻是否已经超过仿真终止时刻，若未超过，则跳转至步骤(2)重复步骤(2)-(5)的仿真架构主循环；若已超过，则以仿真终止时刻为下一仿真时刻，完成该时刻的系统状态计算与波形输出后，结束仿真进程。

所述步骤(1)具体包括以下步骤：

步骤(1.1)所述系统参数输入是输入待仿真电路的主电路元件参数与控制器参数；

步骤(1.2)所述系统方程列写是根据电路拓扑、元件参数与器件开关状态组合列写电路的状态方程、输出方程与阈值方程，如公式(1)-(4)所示；

步骤(1.3)所述外部事件设置是预先设置外部事件发生的时间轴及其相应的事件内容；

步骤(1.4)所述系统状态初始化对状态变量、器件开关状态与控制器存储器进行初始化，步骤(1.4)结束后进入所述步骤(2)-(5)的仿真架构主循环。

优选地,所述步骤(2)具体包括以下步骤：

步骤(2.1)所述系统状态更新的方式如公式(5)所示：

步骤(2.2)系统事件检测与处理，若发生所述控制事件，则根据控制器存储的信息进行采样或开关操作，更新主动器件的开关状态向量swact,k，并计算下一次控制事件的内容与其发生的时间间隔δtctr；若发生所述外部事件，则根据其内容更新所述系统参数向量pk，并计算下一次外部事件发生的时间间隔δtext。

所述步骤(3)具体包括以下步骤：判断所述器件开关状态向量swk与所述系统参数向量pk是否发生变化，若swk与pk均未变化，则保持系数矩阵不变；若swk或pk变化，则首先利用第k-1步的系数矩阵c1,k-1与d1,k-1计算输出方程，并完成第一次波形输出，然后根据最新的器件开关状态向量swk与系统参数向量pk更新系数矩阵，如公式(4)所示，随后再次计算输出方程，并完成所述第二次波形输出。

所述步骤(4)具体包括以下步骤：

步骤(4.1)所述积分步长的初始估计以当前时刻距下一次控制事件发生的时间间隔δtk,ctr与距下一次外部事件发生的时间间隔δtk,ext中的较小者为积分步长的初始估计δtk，同时，δtk不应超过用户所设置的最大步长δtmax；

步骤(4.2)所述导数与增量计算，在第k步数值积分过程中，可将被仿真系统视为一个线性时不变系统，其系统方程的系数矩阵ak,bk,ck,dk为与仿真时间t和状态向量x均无关的常矩阵，状态向量的各阶导数按照公式(6)所示的递归格式依次求得：

x⁽ⁱ⁾(t)＝akx^(i-1)(t)+bku^(i-1)(t)i≥1(18)

其中，x⁽ⁱ⁾(t)和u⁽ⁱ⁾(t)分别为状态向量与输入向量的i阶导数；一般情况下，在电力电子系统中，由于输入向量一般由显式的独立电源组成(例如：含谐波的交流电源与直流电源)，因此输入向量u⁽ⁱ⁾(t)(i≥0)的各阶导数容易显式获取。

利用状态向量的各阶导数，基于taylor展开的数值积分方法，t＝tk+1时刻的系统状态数值解如公式(7)所示：

定义数值积分过程中，状态向量的第i阶增量为公式(8)

步骤(4.3)所示精度检测，判断公式(9)是否成立：

δixk≤abstolk(21)

其中，abstolk为所设定的绝对误差上限，若公式(9)成立，则为情况1，跳至步骤(4.4)；若公式(9)不成立则将所分析的阶数提高一阶，跳回至步骤(4.2)；一般情况下，为防止所分析的阶数升至过高，会要求用户在仿真开始前设定所分析阶数的上限值qmax。若所分析的阶数已达到qmax，则为情况2，跳至步骤(4.4)；

步骤(4.4)所述阶数选取与步长调整，阶数选择方法如公式(10)所示，步长调整方法如式(11)所示：

步骤(4.5)所述数值积分，数值积分方法如式(12)所示：

对于一个包含n个分量的状态向量xk＝[xk,1,xk,2,…,xk,n]^t，步骤(4.3)-(4.4)中的操作将被作用于状态向量的每一个分量，只有当每一个分量均满足精度要求时，才能通过步骤(4.3)中的精度检测。

所述步骤(5)具体包括以下步骤：

步骤(5.1)阈值条件判断，被动事件是否发生，可以根据是否存在触发变量越过阈值条件进行判断，定义如公式(13)所示的阈值函数gth(δt)：

gth(t)＝y2(t)-cth(25)

则某一触发变量在该积分步长内越过阈值条件，等价于其对应的阈值函数在该积分区间两端处的值异号，若存在触发变量越过阈值条件，则可知有被动事件发生，跳至步骤(5.2)，反之则可知被动事件未发生，跳至步骤(5.5)；

步骤(5.2)采用割线法定位被动事件的发生时刻，由步骤(5.1)可知触发向量中的某一个或某几个分量达到了阈值条件，假设积分步长足够短，每个触发变量在一个积分区间内仅能达到一次阈值条件，由于状态向量和输入向量在t＝tk时刻的各阶导数已在步骤(4)中求出，被检测到的触发变量的各阶导数由公式(14)求出：

基于公式(14)，阈值函数gth(δt)可被表示为时间显示多项式函数公式(15)：

采用割线法定位阈值函数gth(δt)的零点(即被动事件发生的时刻)，以一个二极管为例，其被动事件定义及其检测过程如图2所示，其中图2(a)所示为二极管关断事件的定位过程，图2(b)所示为二极管开通事件的定位过程，图2中所示定位过程依次从①-④，④号点对应的时刻为最终求得的被动事件发生时刻，其中，t4被选取为被动事件发生时刻准确解tsw的数值近似，若触发向量中有多个分量j1,j2,……在同一积分区间内达到阈值条件，则取其中的最小者为当前时刻距下一次被动事件发生时刻的时间间隔δtk,pas；

步骤(5.3)确定开关状态。根据事件内容改变发生被动事件的电路元件的状态；以二极管为例，若其当前处于开通状态，则发生被动事件后应将其关断，反之则应将其开通，由此确定第k+1步不可控器件的开关状态向量swpas,k+1。

步骤(5.4)数值积分，若有被动事件发生，则意味着此时的数值积分步长δtk相较步骤(4)中所确定的积分步长被缩短，因此步骤(4)中的数值积分结果已不能继续使用，需要重新进行数值积分；

步骤(5.5)计算下一仿真时刻。按照公式(16)计算下一仿真时刻：

tn,k＝tp,k+δtk,pas(28)。

下面给出一个具体的实施例，对如图3(a)所示的功率放大器进行仿真，通过将dsed仿真结果与实验结果进行比较验证dsed仿真的准确性，通过对比dsed与现有商业仿真软件的仿真的计算点数与仿真耗时验证dsed仿真的高效性。该功率放大器的具体参数为：网侧三相电源有效值ea,b,c(rms)＝87v，网侧三相电阻rgrid＝0.2ω，网侧三相电感lgrid＝250nh，直流母线电容cbus＝2695μf，桥臂电感l1,2,3,4＝190μh，桥臂电阻rl1,2,3,4＝0.2ω，输出电容co＝60nf，负载电阻rload＝30ω，负载电感lload＝105μh，开关频率fs＝100khz，参考波频率fr＝1khz，调制比m＝0.9。该算例的仿真时长为0.2s，在仿真过程中，参考波由三角波突变为正弦波，如图3(b)所示。

仿真过程具体包括以下步骤：

步骤(1)仿真准备工作与系统状态初始化，具体分为以下步骤：1)将上述主电路与控制器参数输入至仿真程序；2)根据电路拓扑与参数列写系统的状态方程、输出方程与阈值方程；3)设置参考波由三角波突变为正弦波这一外部事件发生的时刻。

本算例中，状态方程为公式(17)

其中，状态向量输入向量

为便于表示系数矩阵a,b，定义开关函数sa,b,c分别表示三相二极管桥式整流三个桥臂的开关状态，开关函数sl1,l2,l3,l4分别表示功率放大器4个并联桥臂的开关状态，开关函数的含义为：s＝1表示桥臂上管导通，s＝0表示桥臂下管导通，s＝-1表示桥臂阻断(即上下管均关断)。基于上述开关函数定义，系数矩阵的具体形式为公式(18)：

其中，当sa＝-1,sb＝-1,sc＝-1时，a11＝03×3,b11＝03×2；当sa＝-1,sb≠-1,sc≠-1时，当sa≠-1,sb＝-1,sc≠-1时，当sa≠-1,sb≠-1,sc＝-1时，当sa≠-1,sb≠-1,sc≠-1时，

输出方程为公式(19)：

y＝cx+du(31)

其中，

阈值方程为公式(20)：

其中，id与vd分别为二极管的管电流与管电压。

步骤(2)系统状态更新与系统事件处理：首先更新系统状态，然后判断是否发生系统事件。若发生控制事件与/或外部事件，则根据事件内容更新器件开关状态与/或电路元件参数。具体流程参见发明内容的步骤(2)。

步骤(3)系数矩阵更新与仿真波形输出：若器件的开关状态与电路参数均未变化，则保持系数矩阵不变，计算输出方程并进行波形输出；若器件的开关状态或电路参数变化，则在更新系数矩阵之前计算输出方程并完成第一次波形输出，然后根据最新的器件开关状态更新系数矩阵，随后再次计算输出方程并完成第二次波形输出。具体流程参见发明内容的步骤(3)。

步骤(4)系统状态的灵活自适应(fa)数值积分：首先通过事件排序给出积分步长的初始估计，随后计算状态向量的各阶导数与增量，其次进行精度检测，然后进行自适应的阶数选取与步长调整，最后进行数值积分。具体流程参见发明内容的步骤(4)。

步骤(5)被动事件定位：若有触发变量越过阈值条件，则发生被动事件，此时采用割线法定位被动事件发生的时刻，同时缩短积分步长以定位被动事件，并重新进行数值积分。具体流程参见发明内容的步骤(5)。

步骤(6)仿真进程判断与推进：判断下一仿真时刻是否已经超过仿真终止时刻，若未超过，则令k＝k+1并跳转至步骤(2)重复步骤(2)-(5)所述主循环；若已超过，则以仿真终止时刻为下一仿真时刻，完成该时刻的系统状态计算与波形输出后，结束仿真进程。具体流程参见发明内容的步骤(6)。

应用本发明的仿真方法得到仿真结果，并与实验结果进行对比。测试结果如图4所示。可以看到，本发明所提dsed仿真方法所给出的仿真波形与实验波形吻合良好，验证了dsed仿真的准确性。

在不同的仿真精度下，将所提dsed仿真方法的计算点数与仿真耗时与现有商业软件的计算点数与仿真耗时进行对比。参与仿真性能对比的仿真器与算法如表1所示。为了单独评估本发明所提fa数值积分算法对于仿真的提速效果，除fa数值积分算法外，ode45算法也在本发明所提dsed仿真方法中实现作为对比的对象。嵌入fa数值积分算法与ode45算法的dsed仿真方法均以c++编程实现，二者除数值积分算法不同外，其余部分(主要包括事件定位机制)完全相同。

表1参与仿真性能对比的仿真器与算法

为了衡量仿真波形的精度，定义某一输出变量仿真波形的相对误差为公式(21)：

其中，ysim为仿真波形向量，yref为参考波形向量，||·||2为向量2-范数。如果仿真的输出变量不止一个，则综合的相对误差为所有输出变量相对误差的几何平均值。本算例中，参考波形由simulink在以下条件下提供：仿真算法为ode45，最大仿真步长为1ns。

本算例的仿真性能测试结果如图5所示。可以看到，在所有被比较的仿真器与算法中，在达到同等仿真精度的条件下，本发明所提出的带有fa数值积分算法的dsed仿真方法表现出了最高的仿真效率(即最短的仿真耗时)，且占用的存储资源最少(即仿真点数最少)。与plecs和simulink相比，本发明所提出的带有fa数值积分算法的dsed仿真方法分别实现了约10倍与100倍的仿真速度提升，其中fa数值积分算法贡献了3-6倍的仿真速度提升。

最后需要指出的是：以上实施例仅用以说明本发明的技术方案，而非对其限制。尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的精神和范围。