后向离散状态事件驱动电力电子仿真方法、设备和介质

摘要

本发明公开了一种后向离散状态事件驱动电力电子仿真方法，包括：仿真初始化和第k步的仿真计算，k≥0，其中第k步的仿真计算包括以下步骤：S1：生成第k步的系统状态方程；S2：确定各状态变量在第k+1步的候选量化函数值；S3：建立有限状态机，逐个确定各状态变量在第k+1步的量化函数值；S4：基于所述第k+1步的量化函数值和所述系统状态方程计算导数向量；S5：基于事件发生的时刻，确定第k+1步仿真计算的时刻。本发明还公开了一种仿真计算设备和存储了仿真程序的计算机可读存储介质。本发明提出的基于有限状态机实现的仿真方法能高效地选取Q函数向量，在解算刚性系统时其仿真效率和数值稳定性明显优于量化状态系统QSS方法和传统的时间离散刚性解法。

说明书

技术领域

本发明涉及电力电子仿真技术领域，更具体地，涉及电力电子变换器系统的仿真方法、仿真计算设备、和存储仿真程序的计算机可读存储介质。

背景技术

计算机数值仿真是分析电力电子系统非理想特性重要的工具。然而，在考虑功率开关器件非理想模型和线路杂散参数的电力电子系统中，其仿真模型具有大时标跨度、不连续点多、非线性和刚性强等特点。传统的基于时间离散的数值计算方法(如梯形法、龙格库塔法等)在解算这样的系统时存在以下问(1)步长选择两难——定步长可能会导致计算误差的累积，变步长又会导致大量的计算时间浪费在寻找适当步长上，还有可能导致过小的步长使仿真难以推进；(2)仿真效率低——无论是变步长还是定步长算法，都不能在一个步长内跨过不连续状态，否则将使仿真结果出现错误，必须通过迭代等方法来定位不连续点，使得计算量大大增加，拖慢仿真速度；(3)算法收敛性差——在如此大的时标跨度下易出现数值不稳定现象。

针对传统时间离散数值算法的缺陷，Kofman E等人从数值分析的角度提出了一种采用事件离散的量化状态系统算法(Quantized-state systems，QSS算法)。参见Kofman E，Junco S.Quantized-state systems：a DEVS Approach for continuous systemsimulation[J]，Simulation Transactions of the Society for Modeling&SimulationInternational，2001，18(3)：123-132，此处引入其全文作为参考。QSS算法以各状态变量的Q函数对状态变量或其导数进行离散化切分，并以Q函数构成的向量代替状态变量计算导数向量，求解到达下一个切分点的时刻。该算法在收敛性、稳定性等方面具备优异的数学性质，其全局截断误差可被所选切分长度严格约束，截断误差不会随仿真步数叠加。同时，该算法计算过程简单，计算过程无迭代，判定不连续点不需要迭代和插值。另外，该算法本身天然具有变步长性质，变化平缓处天然以大步长解算，变化剧烈处天然以精细尺度解算，无需浪费时间反复调整步长，因而具备高效解算考虑非理想因素的电力电子变换器模型的潜力。

然而QSS算法在面对刚性系统时的效率和稳定性尚不尽如人意。电力电子系统是典型的刚性系统，其刚性一方面源于大量存在的开关器件，另一方面则源于其包含的多时间尺度过程。面对刚性系统时，如果不采用刚性解法，将导致最大时间步长与雅克比矩阵中变换速度最快的特征值而非最快解梯度相关，为了满足算法稳定域要求步长必须取得极小，否则将导致数值不稳定，极大地拖慢仿真速度。在时间离散算法中，所有的刚性解法都是后向隐式算法。QSS算法则是前向显式算法，基于QSS的方法在解算强刚性系统时会面临效率低、数值不稳定等问题。

Kofman E等人提出了一种后向QSS(BQSS)算法，即选取一组合适的量化函数值作为下一步的Q函数向量，使得据此计算出的导数正负符号能让各状态变量在下一步向各自的Q函数趋近。参见Migoni G，Kofman E，Cellier F.Quantization-based newintegration methods for stiff ordinary differential equations[J]，Simulation，2012，88(4)：387-407，此处引入其全文作为参考。BQSS作为隐式的事件驱动算法，最大的问题是如何高效实现，其困难之处在于：(1)状态变量的导数不仅与自身Q函数取值有关，还会由于别的状态变量的Q函数取值不同而发生变化，所以各状态变量的Q函数的确定是互相影响、互相制约的；(2)每一步计算中各状态变量的Q函数都有两种取值，对于高维复杂系统若采用枚举的方式将导致组合量巨大，难以实施。目前BQSS已有的实现方案未考虑Q函数确定时的相互制约，对于算法进程中已确定Q函数取值的部分，可能会由于待定部分的取值而使导数正负号发生变化，使其状态变量向Q函数趋近的性质遭到破坏，在多维复杂系统中，将无法找出正确的Q函数向量，甚至根本无法进行。有限的算例也都是简单的低维系统，对高维复杂系统尚缺乏有效的实现方案。

发明内容

本发明的目的是提供一种有效地解算强刚性电力电子变换器系统的仿真方法，提升在解算强刚性电力电子系统时的仿真速度和数值稳定性。

根据本发明第一方面的实施例，提供一种后向离散状态事件驱动(BackwardDiscrete State Event Driven，BDSED)电力电子仿真方法，包括：仿真初始化和第k步的仿真计算，k≥0，

其中第k步的仿真计算包括以下步骤：

S1：生成第k步的系统状态方程；

S2：确定各状态变量在第k+1步的候选量化函数值；

S3：建立有限状态机，逐个确定各状态变量在第k+1步的量化函数值；

S4：基于所述第k+1步的量化函数值和所述系统状态方程计算导数向量；

S5：基于事件发生的时刻，确定第k+1步仿真计算的时刻。

可选地，所述步骤S1由以下步骤实现：

(a)以各状态变量x_i的量化函数对x_i进行离散化，所述量化函数为：

其中：i＝1，2，...，n，n为状态变量的个数，节点间距Δq_i为离散区间的长度，ε_i是用于防止Q_i陷入跳变死循环的滞环宽度；以及

(b)基于各状态变量x_i的量化函数值所构成的向量Q^(k)、电路拓扑信息和外部控制信号判断系统状态，求取第k步时的系统状态方程。

可选地，所述步骤S2由以下步骤实现：

找出距离各状态变量x_i^(k)最近的相邻上下切分节点q_i^(k)+和q_i^(k)-，作为在第k+1步的候选量化函数值，即Q_i^(k+1)的候选取值：

其中，q_is为切分节点，m为离散区间的个数。

可选地，所述步骤S3由以下步骤实现：

建立包括四个状态的有限状态机，在每一个状态依据导数向量的极值信息、系统状态方程的雅克比矩阵、以及各自由变量对极值的影响程度，动态地对系统进行降维解耦，由此逐个确定各状态变量x_i在k+1步时的量化函数值Q_i^(k+1)。

根据本发明的基于有限状态机实现的仿真方法，充分利用了电力电子系统状态方程的稀疏特性及其数学性质，在各个状态都对尚未确定取值的状态变量进行检测，动态地对系统进行降维解耦，将系统状态变量间的耦合解开，逐个确定各状态变量在k+1步时的Q函数取值，并保证已确定部分的导数符号不会受到待确定部分的影响，从而在多维复杂系统中高效地实现BDSED。

可选地，所述步骤S4由以下步骤实现：

(a)用所述第k+1步的量化函数值代替状态变量，代入到所述系统状态方程计算导数向量；以及

(b)对步骤(a)中获得的导数向量进行后处理，其中根据电力电子系统的电应力对导数模值进行限幅修正。

可选地，在步骤(b)中，以上界对导数模值进行限幅：

其中|U|_max和|I|_max为系统电容电压和电感电流的稳态峰值，Lmin和Cmin为等效电感和等效电容的最小值，β为限幅系数。

通过根据电力电子变换器系统中的电应力确定每一个状态变量导数的最大幅值，对导数模值进行限幅修正，能够在有效抑制高频振荡的同时保证不破坏计算的准确性。

可选地，所述步骤S5由以下步骤实现：

确定各事件发生的时刻，并选取其中最早的时刻作为第k+1步仿真计算的时间戳，

其中，所述事件包括状态变量事件和非状态变量事件，

所述状态变量事件包括状态变量x_i到达第k+1步的量化函数值Q_i^(k+1)，其中基于所述导数向量计算所述状态变量事件发生的时刻，

所述非状态变量事件包括以下事件中的一个或多个：主动源和受控源激励改变、开关器件动作、器件内部暂态过程状态切换、和控制电路采样。

可选地，所述步骤S5还包括在发生事件的情况下更新状态变量的值。

通过扩展事件的内涵，不仅包括状态变量事件，即某状态变量的值到达切分节点，还包括非状态变量事件，例如主动源和受控源激励改变、开关器件动作、器件内部暂态过程状态切换、控制电路采样等，使仿真方法具备通用的解算电力电子系统的能力。

根据本发明第二方面的实施例，提供一种电力电子仿真计算设备，包括：

存储器；和

处理器，

其中，所述存储器中存储了程序代码，当该程序代码被执行时，使得所述处理器执行如上面所述的后向离散状态事件驱动电力电子仿真方法。

根据本发明第三方面的实施例，提供一种计算机可读的存储介质，其上存储了程序代码，所述程序代码能够被一个或多个处理器执行来实现如上面所述的后向离散状态事件驱动电力电子仿真方法。

总体而言，本发明的实施例具有以下优点中的一部分或全部：

1、可在高维系统中高效地进行Q函数向量选取：本方法采用有限状态机实现，充分利用状态方程的数学性质和电力电子变换系统雅克比矩阵的稀疏特性，动态地对系统进行降维解耦，避免枚举的同时确保已确定部分的导数符号不会受到待确定部分的影响，无需反复地对已确定部分的取值进行调整，解决了后向算法在高维系统中难以实现的问题。

2、抑制数值震荡能力强：本发明提出的BDSED仿真方法为后向隐式方法，数值稳定性较前向算法更强，在等长度的离散区间下可有效抑制数值震荡，仿真精度更高。

3、解算刚性系统效率更高：本发明提出的BDSED方法具备离散事件驱动算法固有的无需定位不连续点、无需反复调整步长、收敛性好的性质；另一方面由于数值稳定性强，可将切分区间放宽，同时又能抑制数值震荡避免无谓事件的产生，因而可大大减少计算步数。和传统的时间离散解法和QSS相比，BDSED在仿真速度上具备明显优势，解算刚性系统时的效率更高，总耗时更少。

4、对电力电子系统有更好的适应性：本发明的一些实施例提供了完备的仿真方法，包含了仿真分析的各个流程，在事件拓展、实现方案、后处理等方面都充分地与电力电子模型相结合，从而能通用地解算各类电力电子系统。

附图说明

通过结合附图对实施例的描述，可以更好地理解本发明的上述和其他特征及优点。

图1示出根据本发明实施例的后向离散状态事件驱动电力电子仿真方法的流程图；

图2示出一个仿真算例的电路拓扑；

图3是图2所示的仿真算例的简化拓扑；

图4是图2所示的仿真算例所用到的器件模型；

图5示出各种仿真方法所得到的TB+管波形及耗时；

图6示出根据本发明实施例的仿真计算设备的示意性框图。

具体实施方式

下面结合附图详细描述本发明的实施例。这些实施例旨在用于解释本发明，而不应该理解为对本发明的限制。在附图中，相同或相似的附图标记表示相同或类似的部件。

图1示出根据本发明一个实施例的后向离散状态事件驱动电力电子仿真方法的流程图。

电力电子变换器系统既不是连续系统，也不是完全离散系统，而是典型的连续与离散混杂系统，可视为由不连续事件分割开来的连续系统，在不连续点之间系统是连续的，系统可由如下状态方程描述：

其中雅克比矩阵A和与输入向量对应的矩阵B可以均为常系数矩阵，x为系统模型和开关器件模型中的状态变量，u为系统的输入向量。在电力电子变换器系统中，输入向量u例如可以包括主动源和受控源，状态变量x例如可以包括电容电压和电感电流。

参照图1，在仿真初始化阶段，k＝0，输入初值x⁽⁰⁾，Q⁽⁰⁾，和Δq。x⁽⁰⁾是状态变量的初始值，Q⁽⁰⁾是状态变量的离散化初始值，Δq是离散区间的长度。

第k步的仿真计算包括以下步骤。

步骤S1，生成第k步的系统状态方程。

在运算开始前，以各状态变量x_i的量化函数即Q函数对x_i进行离散化，其中离散区间的切分例如可以通过一系列等间距的切分节点完成：

其中：i＝1，2，...，n，n为状态变量的个数，ε_i是为防止Q_i陷入跳变死循环而引入的滞环宽度，节点间距Δq_i为离散区间的长度，则Q_i是分段常函数。

在所讨论的实施例中，节点间距Δq_i是固定的。但本领域技术人员可以理解，本发明的原理也适用于节点间距Δq_i可变的情况。

由第k步各状态变量x_i的Q函数值所构成的向量Q^(k)、电路拓扑信息和外部控制信号判断系统状态，求取第k步时的系数矩阵A、B和u^(k)里的受控源部分，u^(k)的主动源部分可直接由电路信息提取，从而生成第k步的系统状态方程。

步骤S2，确定各状态变量x_i在第k+1步的候选量化函数值。

具体地，找出距离状态变量x_i^(k)最近的相邻上下切分节点q_i^(k)+和q_i^(k)-，作为Q_i^(k+1)在k+1步的候选取值。

若状态变量x_i在第k步时的值为x_i^(k)，则在第k+1步其Q函数的取值Q_i^(k+1)有两种可能，分别为上下两个距离x_i^(k)最接近的切分节点，记为q_i^(k)+和q_i^(k)-，可由x_i^(k)和Δq_i根据式(3)计算而得。

其中，q_is为切分节点，m为离散区间的个数。

步骤S3，建立有限状态机，逐个确定各状态变量x_i在第k+1步的Q函数值Q_i^(k+1)。

该步骤应使得据此计算出的导数尽量能够让x_i^(k)向Q_i^(k+1)趋近。

本发明的后向离散状态事件驱动(BDSED)算法通过选取k+1步时的Q函数向量Q^(k+1)＝(Q₁^(k+1)，Q₂^(k+1)，......Q_n^(k+1))^T来代替状态变量计算导数向量，其中Q_i^(k+1)的取值为q_i^(k)+和q_i^(k)-之间的一个。Q^(k+1)函数向量应使得据此计算出的导数尽量能够让x_i^(k)向Q_i^(k+1)趋近，即满足：

定义此性质为方向一致性。其中导数向量各分量可表示为：

f_i是根据状态方程得到的x_i导数计算表达式。对于部分状态变量来说必定不满足方向一致性，即无论怎样选取Q函数组合总是有：

对于这部分状态变量可强制令其导数为0。注意到雅克比矩阵A是常系数矩阵，在每一步计算前输入向量已经确定，因此f_i(Q^(k+1)，U^(k))是一个关于Q^(k+1)的多元线性函数，表达式为：

f_i(Q^(k+1)，U^(k))＝a_ilQ_l^(k+1)+a_i2Q₂^(k+1)......+a_inQ_n^(k+1)+C_i>

其中，a_ij为雅克比矩阵A第i行j列的元素，C_i为与输入向量有关的常数，在输入确定后已知。每一个自变量Q_i^(k+1)只有两种取值q_i^(k)+或q_i^(k)-，多元线性函数极大值和极小值易于求得，且每固定一个变量的取值，都会使极大值和极小值发生变化。本发明利用这一性质设计一个有限状态机，建立四个状态，在每一个状态依据极值信息、雅克比矩阵以及各自由变量对极值的影响程度，动态地对系统进行降维解耦，逐个确定各状态变量x_i的Q函数在k+1步时的取值Q_i^(k+1)，保证已确定部分的导数符号不会受到待确定部分的影响，并找出确定不满足方向一致性的状态变量。

步骤S4，基于所述第k+1步的量化函数值和所述系统状态方程计算导数向量。

具体地，利用Q^(k+1)代替x^(k)代入到状态方程计算导数向量。并且，根据一个优选实施例，对导数向量进行后处理。

首先将Q^(k+1)代入到状态方程中计算导数向量：

在步骤S3已找出确定不满足方向一致性的状态变量，相应地在(8)的结果中将对应的分量调整为0。同时由于在状态变量量化长度固定的条件下，某状态变量导数绝对值过大将直接导致该步计算时间间隔过小，如果这种状态不断重复，会出现高频振荡。因此，根据一个优选实施例，在计算中限制导数幅值来抑制异常振荡，由此得到最终导数向量。电力电子系统中的状态变量为电容电压和电感电流，其导数的实际物理意义与电应力互相对应，因此依据系统中的最大电应力进行模值限幅处理。模值限幅处理的原则是在有效抑制高频振荡的同时保证不破坏计算的准确性。

步骤S5，基于事件发生的时刻，确定第k+1步仿真计算的时刻t^(k+1)。

利用步骤S4所得到的最终导数向量，计算各状态变量到达下一个切分节点，也就是到达Q_i^(k+1)所需时间。

另外，根据本发明的优选实施例，对“事件”内涵进行拓展。除了上述状态变量事件(即，状态变量到达Q_i^(k+1))以外，将主动源和受控源激励改变、开关器件动作、器件内部暂态过程状态切换、控制电路采样等非状态变量事件，一并纳入事件驱动体系中来。通过扩展事件的内涵，使仿真方法具备通用的解算电力电子系统的能力。

选取其中最先发生的事件时刻作为下一步计算的时间戳，更新系统时间以及状态变量的值，完成第k步计算。

如果所确定的第k+1步仿真计算的时刻未达到设定的仿真终止时间，则进入第k+1步计算，否则仿真终止。

从上述描述可以理解，每一步运算只在事件发生时进行，仿真计算完全由事件驱动。

为了进一步说明本发明实施例的采用有限状态机实现的BDSED电力电子仿真方法并验证其效果，使用考虑开关器件非理想模型和线路杂散参数的带星形阻感负载三相两电平变换器电路进行算例分析，仿真算例电路拓扑如图2所示。

使用的IGBT三端电容行为模型如图4A所示，根据IGBT的开关特性将其开通、关断瞬态过程各分为5个阶段，在不同阶段改变参数R_g、C_gc、R_PN的值以及I_T的表达式。二极管模型则等效为可变电阻r_d并联可变RC以模拟其正反向恢复特性，如图4B所示。仿真电路如图3所示，其中L_S1+、L_S2+、L_S3+、L_S1-、L_S2-、L_S3-为母排杂散电感，设其杂散电感值均相等，负载为三相对称阻感负载。

仿真过程为单个开关周期，其中T_A+管处于关断稳态，T_C+管处于开通稳态，T_B+管经历关断、开通两个过程，开关周期为0.2ms，占空比为50％。由于负载电感值相对杂散参数较大，可以认为在瞬态过程中负载电流不会改变方向，故而可只保留负载电流通过的换流电路，仿真模型可简化为图3。其换流通路为A、B相桥臂上管IGBT及其下管反并联二极管以及C相桥臂上管的反并联二极管及其下管IGBT。该仿真模型中，输入变量为各IGBT的门极驱动电压U_{g_x}以及各电流源I_{T_x}的电流；模型可变参数包括各IGBT的门极电阻R_{g_x}、基区电阻R_{PN_x}、gc极间杂散电容C_{gc_x}、等效电流源I_{T_x}，各IGBT二极管的动态电容C_{d_x}、动态电阻R_{d_x}、静态电阻r_{d_x}。上述各量的下标“_x”表示其对应的IGBT或二极管编号。系统状态变量为独立的电容端电压和电感电流，该仿真模型中：每一个IGBT含有3个独立的电容电压，每一个二极管含有1个独立的电容电压；三相杂散电感中L_S1+与L_S1-、L_S2+与L_S2-、L_S3+与L_S3-的电流分别相等因此只有3个独立的电感电流，而三相负载电感中由于其电流之和为零故而只有2个独立的电感电流。整个电路合计包含12个独立电容电压和5个独立电感电流，因此是一个17维系统。

仍然参照图1，第k步(k≥0)的仿真计算包括以下步骤。

步骤S1，求取第k步的系统状态方程。

运算开始前，按照上述(2)式对系统状态变量离散化。根据第k步的Q函数向量Q^(k)以及外部控制信号求解上述可变参数，分析电路拓扑是否发生改变，随后更新电路节点信息，由计算机程序自动列写系统状态方程。

步骤S2，找出距离x_i^(k)最近的相邻上下切分节点q_i^(k)+和q_i^(k)-，作为Q_i^(k+1)的候选取值。

对17个状态变量，均按照式(3)计算其上下相邻的切分节点q_i^(k)+和q_i^(k)-，易知q_i^(k)+、q_i^(k)-和x_i^(k)之间有如下性质：

(q_i^(k)--x_i^(k))·(q_i^(k)+-x_i^(k))＜0>

将本步的计算结果储存在数组中用于后续步骤的计算。

步骤S3，建立有限状态机，逐个确定各状态变量x_i在第k+1步的Q函数值Q^(k+1)。

构造17维的初始极大值方阵maxQ和初始极小值方阵minQ。maxQ第i个行向量maxq_i是使得第i个状态变量导数取到最大值的Q函数组合；minQ第i个行向量minq_i是使得第i个状态变量的导数取为最小值的Q函数组合：

初始化一个状态变量索引集合Xset＝{1，2，3......n}和两个空集X1和X2，这里n等于17。每确定一个Q_i^(k+1)的值，则maxQ和minQ矩阵第i列也要动态调整，并将相应的下标索引i从Xset中移除，f_i(Q^(k+1)，U^(k))的自由变量的维数降低1，直到Xset为空集为止，17个状态变量的Q函数取值全部确定。

状态1：对每一个i∈Xset，分别计算并检查其对应的f_i(maxq_i，U^(k))和f_i(minq_i，U^(k))的正负号情况：若f_i(maxq_i，U^(k))≤0，则可确定Q_i^(k+1)＝q_i^(k)-，将maxQ和minQ两个矩阵第i列的值全部替换为q_i^(k)-并从Xset中移除元素i；若f_i(minq_i，U^(k))≥0，则可确定Q_i^(k+1)＝q_i^(k)+，将maxQ和minQ两个矩阵第i列的值全部替换为q_i^(k)+，同样从Xset中移除元素i；。检查结束后，如果有i被移除，则下一状态仍进入状态1，对Xset中剩余的索引代入的maxQ和minQ重新检查；否则进入状态2。

状态2：对每一个i∈Xset，若矩阵A中对应的对角线元素a_ii≥0，按如下方式构造两个向量和其第j个分量为：

分别计算并检查其对应的和的正负号情况：若则可确定Q_i^(k+1)＝q_i^(k)-，将maxQ和minQ两个矩阵第i列的值全部替换为q_i^(k)-并从Xset中移除元素i；若f_i(，U^(k))≥0，则可确定Q_i^(k+1)＝q_i^(k)+，将maxQ和minQ两个矩阵第i列的值全部替换为q_i^(k)+，同样从Xset中移除元素i。检查结束后，如果有i被移除，则下一状态仍进入状态1，对Xset中剩余的索引代入的maxQ和minQ重新检查；否则进入状态3。

状态3：对每一个i∈Xset，若矩阵A中对应的对角线元素a_ii＜0，则如同状态2一样构造两个向量并检查和向正负号情况，若同时有：

f_i(Q^(k+1)，U^(k))·(Q_i^(k+1)-x_i^(k))＜0>

则说明对于状态变量x_i，不论如何选取Q^(k+1)，都无法满足方向一致性，对于这样的状态变量强制其导数为0即可。令Q_i^(k+1)＝Q_i^(k)，更新矩阵并移除相应地元素。检查结束后，如果有i被移除，则下一状态跳回到状态1；否则进入状态4。

状态4：从当前Xset剩余的元素中，选取下标s为：

s＝arg min(f_i(maxq_i，U^(k))-f_i(minq_i，U^(k)))>

其中i∈Xset。构造向量该向量第j个分量为：

若则令Q_s^(k+1)＝q_s^(k)+，否则令Q_s^(k+1)＝q_s^(k)-。将maxQ和minQ两个矩阵第i列的值全部替换为Q_s^(k+1)，从Xset中移除s，并将s添加到集合X2中。下一状态跳回状态1。

算法停机条件为Xset变为空集，此时已全部确定Q_i^(k+1)的取值，Q^(k+1)已经得到。

步骤S4，计算导数向量

首先由状态机得到Q^(k+1)后据此按照(8)式计算向量得到后对X1与X2集合里的下标所对应的分量进行调整。优选地，为了进一步抑制异常震荡，以上界对导数模值进行限幅，的大小由变换器系统电应力决定，具体和电压、电流峰值和等效电容、电感值的大小有关：

其中|U|_max和|I|_max为整个系统电容电压、电感电流的稳态峰值，β为限幅系数，为保证导数限幅不破坏仿真精度，β必须大于1。β取值越大，导数限幅越不会破坏仿真精度，但进一步抑制振荡的能力越弱，需要选取合适的值。最终的导数向量为

步骤S5，基于各事件发生的时刻，确定下一次计算发生的时刻，完成第k步计算，推动仿真进行。

由步骤S4所得的导数向量，计算各状态变量到达其Q函数的下一个切分节点，也就是到达Q_i^(k+1)所需时间：

当有x_i值达到下一个切分节点，则系统产生一次事件，需要更新状态变量的值。除状态变量事件外，对事件的内涵进行拓展，将输入改变、开关状态改变、控制采样等事件一并纳入进来，任何事件发生都要更新电路方程，仿真进程完全由事件驱动推进，只在每次有事件发生的时刻进行计算。设第k步计算时的时间戳为t^(k)，则两次计算的时间间隔Δt^(k)、第k+1步的时间戳t^(k+1)，以及k+1步的状态变量x^(k+1)分别为：

其中Δt_e表示最先发生的非状态变量事件距当前仿真时刻的时间间隔。至此BDSED算法第k步运算结束，进入第k+1步运算，直到到达设定的仿真终止时间。

仿真结果分析：

在MATLAB平台上，分别用MATLAB自带的4种刚性ODE(常微分方程)算法(经试验相对误差设为0.4％时即可获得较准确的仿真波形)、QSS算法以及本发明提出的后向BDSED方法对上述仿真模型算例进行解算。为在同一环境下比较，所有数值计算统一用MATLAB的m语言实现。为使得运算精确而稳定地进行，QSS和BDSED中各状态变量的离散区间大小与其稳态下的幅值须成等比例关系，用k表示此比例系数：

离散区间相当于时间离散中的步长。k的值越小，则离散区间越小，对状态变量的切分越精细，算法精度和数值稳定性越高，但也会产生更多的事件，仿真耗时更长。图5分别为传统ODE算法中速度最快的ode23tb以及QSS和BDSED算法所得到的T_B+管压降、电流波形。由图5可知本文提出的有限状态机的实现方案可在多维复杂系统中高效地选取Q函数向量，得到准确的仿真波形。在k＝0.8％时QSS由于离散区间过大而产生严重的数值震荡，导致仿真结果失真，直到将离散区间缩小四倍后才能得到基本准确的波形，但稳态部分仍存在小幅度的数值震荡；而BDSED算法数值稳定性好，在k＝0.8％下仍能有效地抑制数值震荡，得到非常准确的波形。对比图5(e)和图5(f)，BDSED算法虽然单步计算较QSS算法更加复杂，但由于其数值稳定性强，可将算法的离散区间放宽，同时又抑制震荡避免无谓的事件，可减少计算步数，总耗时更少。各数值算法的仿真步数和耗时如表1所示。

表1BDSED与QSS及传统时间离散解法仿真速度对比

由表1可知，有限状态机的实现方案可在多维复杂系统中正确而高效地选取Q函数向量，BDSED的平均单步运算耗时与TR-BDF2等传统刚性ODE算法接近；而由于能够大大减少仿真步数，同传统的刚性ODE算法和QSS相比，BDSED在总耗时上具备明显优势，解算刚性系统时的效率和数值稳定性更高。

综上所述，本发明借鉴BQSS数值算法的思想，面向电力电子强刚性系统，提出了一种后向离散状态事件驱动(BDSED)仿真方法，针对BQSS算法以及基于BQSS算法的BDSED方法在高维系统中难以实现的问题，给出一种采用有限状态机的实现方案。另外，本发明的优选实施例将数值算法与电力电子变换器系统本身更多地结合起来，拓展事件驱动的内涵；同时根据电力电子变换器的参数对导数进行限幅。在一个17维的算例模型上进行验证的结果表明，本发明提出的基于有限状态机实现的BDSED仿真方法能高效地选取Q函数向量，在解算刚性系统时其仿真效率和数值稳定性明显优于QSS方法和传统的时间离散刚性解法。

图6示出根据本发明一种实施例的仿真计算设备600的示意性框图。该仿真计算设备600包括存储器601和处理器602。存储器601中存储了仿真程序，当该仿真程序被执行时，使得处理器602执行如上文中所描述的后向离散状态事件驱动(BDSED)电力电子仿真方法。

本发明还可以实现为一种计算机可读存储介质。该计算机可读存储介质可以存储程序指令，后者可以被一个或多个处理器执行来实现本发明实施例的BDSED电力电子仿真方法。该计算机可读存储介质例如可以包括磁或光介质，例如固定或可移动盘，磁带，CD-ROM，DVD-ROM，CD-R，CD-RW，DVD-R，DVD-RW或者蓝光盘。该存储介质还可以包括易失性或非易失性存储介质，例如RAM，ROM，闪存，USB存储器，等等。该计算机可读存储介质还可以包括微机电系统(MEMS)，以及可通过通信介质例如网络或无线链接访问的存储介质。

以上参考附图描述了本发明的实施例。这些实施例仅仅是示例性的，而并非旨在对本发明的范围进行限制。本领域技术人员可以在不偏离本发明的精神和原则的情况下做出各种修改、等同替换和改进。本发明的保护范围仅由权利要求书限定。