CN101382969A - 一种多步变步长电磁暂态仿真方法 - Google Patents

Abstract

本发明提供了一套电磁暂态计算中能考虑任意多重事件的多步变步长仿真方法——每计算一步都进行事件搜索，精确得到事件发生时刻，并重新计算至断点处，处理事件，解决多重事件同时或在一个仿真步长内发生的次序问题；然后采用多步可变步长的后退欧拉法彻底消除电磁暂态仿真中的数值振荡问题。该方法对事件发生时刻计算精确、可以完全消除数值振荡、能处理任意组合的事件序列，适合于具有大量快速开断的电力电子器件的系统的电磁暂态技术。

Description

一种多步变步长电磁暂态仿真方法

技术领域

本发明涉及电力系统电磁暂态过程的模拟与计算领域，具体涉及一种考虑任意重事件发生的多步变步长电磁暂态仿真方法。

背景技术

电磁暂态仿真中，由于开关器械动作、网络结构变化，如含有电感电路直接断开、电力电子器件的导通与截止等，引起的非状态变量在事件发生后不正常地摆动，这就是电磁暂态仿真中的数值振荡现象。

电磁暂态仿真计算一般采用隐式梯形积分将微分和偏微分方程转化为简单代数方程特性差分方程，具有相当高的准确性和数值稳定性(A稳定)。但从70年代开始，发现电磁暂态计算中使用隐式梯形积分可能导致数值振荡问题：例如，在开关断开电感支路电流时，电感两端电压会呈现围绕正确值数值振荡；同样，当某一电压源通过开关向电容器突然充电时，电容器电流也会呈现类似的数值振荡现象。此外，在非线性电感的工作状态发生变化(如从饱和至不饱和区，或者相反)也会出现数值振荡。

另外在电力电子仿真中，有许多电力电子器件的关断和闭合互为因果关系，例如，GTO的关断造成其他电力电子器件如二极管的闭合。它们虽为因果关系，但很可能是在同一瞬间完成，因此应把它看成同一瞬间的行为，完全有必要发展一种事件搜索方法，以检测不同电力电子器件之间的相互动作关系。

1、数值振荡现象、原因及现有的处理方法

以一个电感支路举例(参见图1)：已知在t-△t时刻的电压、电流，计算t时刻的电压和电流，其离散差分方程可以表示为：

$i_{L_{\mathrm{jk}}}\left(t\right)=\frac{\mathrm{\Δ}\left(t\right)(1+\mathrm{\α})}{2L_{\mathrm{jk}}}[u_j\left(t\right)-u_k\left(t\right)]+I_{\mathrm{histjk}}$

$=\frac{\mathrm{\Δt}(1+\mathrm{\α})}{2L_{\mathrm{jk}}}{u}_{\mathrm{jk}}\left(t\right)+I_{\mathrm{histjk}}$

$I_{{\mathrm{hist}}_{\mathrm{jk}}}=i_{L_{\mathrm{jk}}}(t-\mathrm{\Δt})\frac{\mathrm{\Δ}t(1-\mathrm{\α})}{2L_{\mathrm{jk}}}[u_j(t-\mathrm{\Δt})-u_k(t-\mathrm{\Δt})]$

$=i_{L_{\mathrm{jk}}}(t-\mathrm{\Δt})+\frac{\mathrm{\Δt}(1-\mathrm{\α})}{2L_{\mathrm{jk}}}{u}_{\mathrm{jk}}(t-\mathrm{\Δt})$

其中：

I_histjk——历史项，与上一步的状态量和非状态量都有关

α——系数，取值为0，计算采用隐式梯形积分；

取值为1，计算采用后退欧拉法；

取值0~1，算法介于后退欧拉法和隐式梯形积分之间

α＝0——隐式梯形积分公式，历史项和上一步电流和电压有关：

$i_{L_{\mathrm{jk}}}\left(t\right)=\frac{\mathrm{\Δt}}{2L_{\mathrm{jk}}}{u}_{\mathrm{jk}}\left(t\right)+i_{L_{\mathrm{jk}}}(t-\mathrm{\Δt})+\frac{\mathrm{\Δt}}{2L_{\mathrm{jk}}}{u}_{\mathrm{jk}}(t-\mathrm{\Δt})$

α＝1——后退欧拉法积分公式，历史项只和上一步电流有关：

$i_{L_{\mathrm{jk}}}\left(t\right)=\frac{\mathrm{\Δt}}{L_{\mathrm{jk}}}{u}_{\mathrm{jk}}\left(t\right)+i_{L_{\mathrm{jk}}}(t-\mathrm{\Δt})$

电磁暂态仿真计算中，在遇到网络发生突变时，非状态变量一般无法通过计算准确得到。考察上面电感仿真差分方程式，如果采用隐式梯形积分，由于历史项与上一步的非状态量(两端电压)有关，并且两端电压在0+瞬间计算不正确，这就造成了网络突变时隐式梯形算法产生数值振荡；而如果采用后退欧拉法，历史项与上一步的非状态量(两端电压)无关，就不会有数值振荡的问题，后退欧拉法的这个独特的优势，成为大量的电磁暂态程序在消除数值振荡问题时的首选算法。

从以上分析可以看出，电磁暂态程序中，数值振荡问题的产生是与数值算法的选择直接相关的。

一个简单的算例，参见图2。

直流电源1V，电感0.1H，开关合闸电阻1.0欧姆，断开时电阻无穷大，在0.5秒时刻，电感电流接近稳态值1A。开关SW1在t0＝0.5秒时刻强制开断，电流强制为0，电感两端电压UL将会发生突变，产生一个冲击响应(幅值和电感历史电流以及仿真计算的步长有关)。根据差分方程式，如果采用纯隐式梯形算法，在0.5秒后的计算中，电压将出现数值振荡现象：

$\frac{\mathrm{\Δt}}{2L}[u_L(t_0+\mathrm{\Δt})+u_L\left(t_0\right)]=0$

……

u_L(t₀+kΔt)＝(-1)^ku_L(t₀)

仿真结果如图3。

在时间段(0.49998秒~0.50010秒)，放大后的电压波形如图4。

为了解决数值振荡问题，中国电力科学研究院和加拿大大不列颠哥伦比亚大学(UBC)合作，共同研究开发了数值临界阻尼法(Critical Damp Adjustment，CDA)，并将该方法用于EMTP程序中。CDA法的主要原理就是利用后退欧拉法能避开非状态变量在突变时刻的值，从而不会产生数值振荡。

其主要的计算步骤为：

(1)在一般情况下仍然使用隐式梯形积分法，其时间步长为△t；

(2)若在t＝tz时刻网络发生突变，采用后退欧拉法，其步长改为△t/2，共进行两次步长为△t/2的后退欧拉法积分计算；

(3)然后在tz+△t后，继续采用隐式梯形积分法计算，其步长恢复到△t。

研究表明，二个半步长的后退欧拉法可以消除数值振荡，而且，也必须要经过二次后退欧拉法才能消除网络突变引起的数值振荡，主要原因是：

第一个半步长后退欧拉法，计算出了正确的状态变量，而非状态变量很有可能是一个冲击响应，并且在某些特殊算例中，需要根据这个非状态变量的冲击响应进行网络突变的判断和操作；

第二个半步长后退欧拉法才是真正消除数值振荡的有效步骤。

最初在进行CDA算法研究时，研究人员认为第一个半步后退欧拉法的作用只是为第二个半步后退欧拉法提供初值，其本身没有物理意义，因此在第一个半步长后退欧拉法结束时不进行网络突变判断，而将它留到整数步长时刻进行判断。后来才发现这个观点是错误的，第一个半步长后退欧拉法的结果反映了瞬变过程的冲激信号特性，其物理意义十分明显，结果不可忽略，必须进行网络突变的判断和操作。中国电科院林集明教授等在其完善开发的电磁暂态程序版本EMTPE中，应用该思路改进了电磁暂态程序，并将这种改进后的能反映冲激信号特性的CDA法称为ICDA(Improved CDA)。

CDA法(包括ICDA法)既保留了隐式梯形积分精度高、稳定性好、编程简单的优点，又基本能消除数值振荡，已经作为EMTP-RV和EMTPE缺省使用的算法；免费使用的程序ATP没有任何处理数值振荡的方法，在后面的测试算例可以看到，其数值振荡极其严重，ATP不适合具有快速开断的电力电子器件的FACTS和直流系统的电磁暂态仿真。

但是，CDA法还是有一定的问题，两次后退欧拉法并不能彻底地消除数值振荡，这个问题是本专利需要解决的重点之一。

2、事件搜索与处理

在电力电子仿真中，有许多电力电子器件的关断和闭合互为因果关系，例如，GTO的关断造成其他电力电子器件如二极管的闭合。它们虽为因果关系，但很可能是在同一瞬间完成，因此应把它看成同一瞬间的行为，完全有必要发展一种事件搜索方法，以检测不同电力电子器件之间的相互动作关系。

EMTPE程序中，提出了“同步响应法”(Simultaneous Response Procedure，简称SRP)，其主要步骤为：

(1)在t＝tz时刻，求解网络方程，得到tz时刻各元件的电压电流等物理量；

(2)根据tz时刻的电压、电流等物理量，找出是否有开关和电力电子器件需要改变其状态，如果没有就退出搜索过程，转向步骤(5)；

(3)如果需要改变状态，就根据其状态修改节点导纳阵、历史项以及外部电流源，保持t＝t1不变，再次求解方程，得到t＝tz时刻其他各个支路的电压、电流等物理量；

(4)重新搜索，转向步骤(2)，直至没有新的元件和新的状态发生改变为止；

(5)退出t＝tz计算，按正常方式进入下一个时间步的计算。

EMTPE中的“同步响应法”(SRP)，附加的计算工作量并不多，计算速度快，仿真结果准确，尤其处理电力电子仿真中复杂动作逻辑非常有效。

EMTP-RV和免费程序ATP都没有和“同步响应法”(SRP)类似的动作逻辑搜索的功能，因此在一些特殊算例中仿真不正确，例如图5的算例。

开关SW1常闭，1.1秒GTO导通，电感电流达到稳定运行值2A；1.2秒时刻，GTO断开，开断瞬间二极管D1导通，电感电流有通路，不应该会产生一个很大的脉冲。没有事件搜索功能的电磁暂态仿真程序EMTP-RV的结果如图6。

明显这个结果是错误的，原因是二极管没有在GTO断开的瞬间同时导通。

“同步响应法”的缺陷是不能精确计算出事件发生时间，在一个仿真步长上有多个事件先后发生时，不能处理这个先后次序，只能认为多个事件在一个时刻同时发生，因此“同步响应法”的计算步长不能过大。

本专利提出的仿真方法中，能精确计算事件发生时间，严格地按照事件发生的顺序正确动作，可以具有较大的仿真步长。

发明内容

本专利提出的技术方法主要包括两部分：(一)多步变步长后退欧拉法；(二)任意重事件搜索与处理功能。

1、多步变步长后退欧拉法

CDA法(或ICDA法)采用连续两次半步后退欧拉法来消除数值振荡，至于为什么要采用半步，主要原因是采用步长△t/2后，后退欧拉法的导纳阵和步长为△t的隐式梯形积分公式一致，因此在后退欧拉法转向隐式梯形积分公式时，不需要再对导纳阵修改分解，只需要修改历史项，计算速度非常快；但是，如果单纯从消除数值振荡角度来说，后退欧拉法采用多大步长都能消除数值振荡。

在图2的算例中，正常步长10us，如果事件发生后，第一个后退欧拉法采用0.35倍正常步长，第二个后退欧拉法采用0.48倍正常步长，计算结果如下：

虽然两次后退欧拉法的步长(分别为0.35△t和0.48△t)不相等，当然也不等于半个正常的步长，但是同样可以消除数值振荡。只要严格遵循“发生事件后一定要使用两次及以上后退欧拉法”这一原则，就不会影响消除数值振荡的效果。这个结论是本专利的重要基础。

如图7的例子，如果只采用两次后退欧拉法，GTO两端电压在断点处的计算结果为：

如图8所示，仍然具有稍微明显的数值振荡问题(电压振幅达到将近9.2V左右)

但是如果采用三步后退欧拉法，GTO两端电压在断点处的计算结果为：

从表中可以看出，数值振荡已经平息得很好，只有0.0026振幅

因此，本专利提出变步长的多步变步长后退欧拉法，其技术要点为：

(1)采用多步后退欧拉法的次数为3次及以上，次数可以由用户自行调节；

(2)步长可以根据计算的要求进行调整，不限制为正常仿真步长的一半；

(3)改变算法和仿真步长时，根据以下公式判断是否需要重新修改导纳阵：

Δt(1+α)

其中：Δt——仿真步长，

      α——系数，取值为0，计算采用隐式梯形积分；取值为1，计算采用后退欧拉法。

每一步计算前，根据算法和步长得到Δt(1+α)的值，如果该值发生改变，就需要重新修改导纳阵的值。

计算结果表明，多步变步长后退欧拉法，编程简单，计算灵活，能彻底消除数值振荡问题。

2、任意重事件搜索与处理功能

在本专利中，事件主要指元件的开断操作、网络故障以及变压器发电机等饱和段的改变等。根据发生的事件引起导纳阵改变的情况，事件可以分为两种类型：

(1)事件引起导纳矩阵的结构发生改变，例如元件的切除等；

(2)事件只影响导纳矩阵的数值发生改变，如开关动作，开关的电阻发生改变。

对于第一种类型的事件，电磁暂态程序PSModel处理包括重新节点编号、重新生产导纳阵、重新分解导纳阵等处理过程；而对于第二种类型的事件，PSModel只需要改变导纳阵的值，重新分解导纳阵，开始新的计算。

每一步计算完毕，PSModel都要检测是否有事件发生，如果在这一步上，没有事件发生，就表明计算成功，可以进入下一步计算；如果有事件发生，就需要对事件进行处理。

假设在时间段(t~t+△t)的计算中，有事件发生在中间时刻t+△t’，如图9所示。具体的处理步骤如下：

步骤①：时间段t～t+△t的计算后，判断在时刻t+△t’发生事件，这一步计算抛弃；

步骤②：回退到时刻t，以时刻t的网络状态，计算t～t+△t’时间段，由于这一步计算的步长由△t变为△t’，导纳阵的值也会发生改变，需要重新形成和分解方程式；

步骤③：在t+△t’时刻，处理事件，改变网络的状态，可能会改变导纳阵的值或结构；

步骤④：然后用新的网络状态和导纳阵的值，采用后退欧拉法和△t/2，计算时间段t+△t’～t+△t’+△t/2；

步骤⑤：再次采用后退欧拉法和△t/2，计算时间段t+△t’+△t/2~t+△t’+△t；

步骤⑥：计算成功后，采用后退欧拉法或隐式梯形公式，继续进行下一步计算。

当然在复杂系统的计算中，很可能会出现在应用后退欧拉法处理断点时，又发生了事件，PSModel只需简单地调整事件发生的时间，继续使用后退欧拉法，就可以应付任意复杂的连续事件发生情况，详细的事件搜索及处理流程参见图10。

3、仿真实例

下面是高压直流仿真计算中，采用本专利方法后，整流侧和逆变侧换流阀的阀电压，采用2us和2000us的计算结果对比，如图11所示。

从上面的测试结果可以看出，本方法计算结果都是正确的，彻底消除数值振荡问题；对于仿真步长有非常好的适应性，从2微秒到2000微秒，不同的步长，计算结果几乎没有区别。这就是采用本专利后的电磁暂态程序体现的巨大优势。

因此本发明提出了一种多步变步长电磁暂态仿真方法，采用三步及以上可变步长的后退欧拉法进行电磁暂态仿真计算，其特征在于包括以下步骤：

(1)在每一步电磁暂态计算完毕后进行事件搜索，判断是否发生了引起网络拓扑结构或发生的事件；

(2)如果有事件发生，根据搜索出的所有事件，并得到发生事件的最小精确时刻：断点时刻；

(3)抛弃上一次的计算结果，改变计算步长，重新计算至上述断点时刻；

(4)在断点时刻，处理相应的事件并改变网络的状态；

(5)在本步所有的事件处理完毕后，在后续三步及以上的计算过程中，都采用多步可变步长后退欧拉法处理数值振荡。

本发明的方法还包括：

(1)在每一步电磁暂态计算完毕后，判断是否发生了引起网络拓扑结构或发生的事件，如果发生了所述事件，在后面的连续三步及以上的计算中，严格采用后退欧拉法进行积分仿真；

(2)为了提高计算速度，后退欧拉法的步长采用1/2正常仿真步长或任意步长。

所述事件包括开关器械动作、网络故障、变压器饱和段的改变、发电机饱和段的改变和电力电子器件的导通与截止。

本发明具有以下优点：

1.在网络发生改变后，一般地，可以采用后退欧拉法，用来消除数值振荡；如果后退欧拉法的步长是正常步长的一半，在转向正常步长隐式梯形算法后，不需要修改导纳阵，计算速度快；

2.发生事件后，使用后退欧拉法计算次数可以由用户来调节，经过我们的研究，本专利认为使用三次及以上的后退欧拉法具有更好的消除数值振荡的效果；

3.在每一步计算后，不管是采用后退欧拉法还是隐式梯形算法，都要进行事件搜索和处理，该功能使得电磁暂态程序在处理同时发生的事件或间隔时间非常短的事件方面的能力非常强；

4.如果在使用后退欧拉法的过程中，继续有事件发生，后退欧拉法就不可能再严格保持半个步长，只要遵循“发生事件后一定要使用三次及以上后退欧拉法”这一原则，就不会影响电磁暂态程序消除数值振荡的效果，因此本专利的重要特点是采用可变步长的后退欧拉法，计算非常灵活；

发生事件后，PSModel不是通过插值的方法来求取断点处的状态量和非状态量的值，而是采用抛弃上一次计算结果，然后重新计算至断点处的方法，这种方法计算量稍大一些，但是对于事件发生的时刻搜索准确，计算结果精确；不需要对于情况更为复杂的控制部分进行插值处理，编程简单；并且对步长具有很大的适应性，不同步长的计算结果相差不大。

附图说明

图1是一个电感支路图；

图2是数值振荡简单例子电路图；

图3是图2例子仿真结果；

图4是图3中的曲线在时间段(0.49998秒~0.50010秒)放大后的电压波形图；

图5是事件搜索与处理的算例图；

图6是没有事件搜索功能的电磁暂态仿真程序EMTP-RV的结果图；

图7是两步后退欧拉法仍然不能完全消除数值振荡的算例图；

图8是只采用两步后退欧拉法，仍然有数值振荡曲线图；

图9是处理的具体步骤；

图10是事件搜索及处理流程图；

图11是高压直流仿真步长采用2us和2000us的计算结果对比图；

图12是图5算例的正确计算结果图。

具体实施方式

采用图五所示的例子，下面将通过实际的计算过程来描述该仿真方案具体的实施方式：

步骤1：GTO和二极管初始状态是不导通的，GTO两端初始电压为电源电压2000v，电感L1和二极管的电流都为0值；正常情况下的计算，积分算法采用隐式梯形算法，正积分步长设置为10us；

步骤2：从时刻0仿真到1.1秒，从时步(1.100000秒~1.100010秒)的仿真中，先试算，然后搜索事件，发现GTO检测在1.100000秒时刻有导通脉冲，且满足导通条件，应该导通，也就是有事件发生，添加该事件到事件序列中，这一步试算失败，抛弃；

步骤3：在1.100000秒时刻处理事件——改变GTO导通状态，改变积分算法为后退欧拉法(这是第一次使用后退欧拉法)，改变积分步长为正常积分步长的一半5us，在时步(1.100000秒~1.100005秒)上，重新形成仿真用的差分方程式并进行求解。同样先试算，然后搜索事件，无事件发生，试算成功，本次计算采纳；

步骤4：在时步(1.100005秒~1.100010秒)上，仍然采用后退欧拉法仿真(这是第二次使用后退欧拉法)，先试算，然后搜索事件，无事件发生，试算成功，本次计算采纳；

步骤5：在时步(1.100010秒~1.100015秒)上，第三次采用后退欧拉法仿真，同样的过程：先试算，然后搜索事件，无事件发生，试算成功，本次计算采纳；

步骤6：从下一步计算开始，三次后退欧拉法使用完毕，改为正常的隐式梯形算法，仿真步长改为正常步长10us，在时步(1.100015秒~1.100025秒)上采用相同的过程，先试算，然后搜索事件，无事件发生，试算成功，本次计算采纳；

步骤7：从1.100025秒开始，采用隐式梯形算法和10us的正常步长仿真计算直至1.199995秒；

步骤8：在时步(1.199995秒~1.200005秒)上，先试算(采用隐式梯形算法和10us的正常步长)，然后搜索事件，发现GTO检测在1.199996秒时刻有截止脉冲进行强制截止，也就是有事件发生，添加该事件到事件序列中，这一步试算失败，抛弃；

步骤9：改变仿真步长，在时步(1.199995秒~1.199996秒)上仿真，采用隐式梯形算法和1us步长，搜索事件，无事件发生，计算成功；

步骤10：在1.199996秒时刻处理事件——改变GTO导通状态为截止，改变积分算法为后退欧拉法(这是第一次使用后退欧拉法)，改变积分步长为正常积分步长的一半5us，在时步(1.199996秒~1.200001秒)上，重新形成仿真用的差分方程式并进行求解。同样先试算，然后搜索事件，发现由于GTO截止，二极管上电压出现正值，应该导通，插值出二极管的精确导通时刻应该为1.199996秒，添加该事件到事件序列中，这步计算失败，抛弃；

步骤11：在1.199996秒时刻再次处理事件——改变二极管的导通状态为导通，积分算法为后退欧拉法(这是第一次使用后退欧拉法)，改变积分步长为正常积分步长的一半5us，在时步(1.199996秒~1.200001秒)上，重新形成仿真用的差分方程式并进行求解。同样先试算，然后搜索事件，无新的事件发生，计算成功；

步骤12：在时步(1.200001秒~1.200006秒)上，仍然采用后退欧拉法仿真(这是第二次使用后退欧拉法)，先试算，然后搜索事件，无事件发生，试算成功，本次计算采纳；

步骤13：在时步(1.200006秒~1.200011秒)上，第三次采用后退欧拉法仿真，先试算，然后搜索事件，无事件发生，试算成功，本次计算采纳；

步骤14：从下一步计算开始，三次后退欧拉法使用完毕，改为正常的隐式梯形算法，仿真步长改为正常步长10us，在时步(1.200011秒~1.200021秒)上采用相同的过程，先试算，然后搜索事件，无事件发生，试算成功，本次计算采纳；

步骤15：从1.200021秒开始，采用隐式梯形算法和10us的正常步长仿真计算直至仿真结束时刻2.0秒；

图5中的算例，电感电流、GTO两端电压以及二极管电流的正确的仿真结果参见图12。

Claims (3)

1、一种多步变步长电磁暂态仿真方法，采用三步及以上可变步长的后退欧拉法进行电磁暂态仿真计算，其特征在于包括以下步骤：

(1)在每一步电磁暂态计算完毕后进行事件搜索，判断是否发生了引起网络拓扑结构或发生的事件；

(2)如果有事件发生，根据搜索出的所有事件，并得到发生事件的最小精确时刻：断点时刻；

(3)抛弃上一次的计算结果，改变计算步长，重新计算至上述断点时刻；

(4)在断点时刻，处理相应的事件并改变网络的状态；

(5)在本步所有的事件处理完毕后，在后续三步及以上的计算过程中，都采用多步可变步长后退欧拉法处理数值振荡。

2、如权利要求1所述的仿真方法，其特征在于：

(1)在每一步电磁暂态计算完毕后，判断是否发生了引起网络拓扑结构或发生的事件，如果发生了所述事件，在后面的连续三步及以上的计算中，严格采用后退欧拉法进行积分仿真；

(2)为了提高计算速度，后退欧拉法的步长采用1/2正常仿真步长或任意步长。

3、如权利要求1-2所述的仿真方法，所述事件包括开关器械动作、网络故障、变压器饱和段的改变、发电机饱和段的改变和电力电子器件的导通与截止。

Images