状态事件的处理

我们已经介绍了时间事件状态事件。现在让我们来看看与状态事件相关的一些重要的困难。出人意料的是,这些困难甚至会在简单的模型里出现。

基本衰减模型

请考虑以下简单得几近不值一提的模型:

model Decay1
  Real x;
initial equation
  x = 1;
equation
  der(x) = -sqrt(x);
end Decay1;

如果试图对此模型进行5秒的仿真,我们会发现仿真大概用2秒后就停止了。结果是如下的轨迹:

../../../_images/Decay1.png

再一次,数值问题混了进来。尽管在数学上x的值应该不可能的降到低于零,但在使用数值积分方法时,少量错误可能混入结果,使得x低于零。发生这种情况时,sqrt(x)的表达式产生浮点异常,结果仿真终止。

保护表达式

为了避免这种情况,我们可以引入一个if表达式防止计算负数的平方根,如下:

model Decay2
  Real x;
initial equation
  x = 1;
equation
  der(x) = if x>=0 then -sqrt(x) else 0;
end Decay2;

这个模型进行仿真,我们得到如下的轨迹[1]:

../../../_images/Decay2.png

再一次地,仿真失败。但是为什么呢?它出于相同的原因失败了:对负数的求平方根而导致的数值异常。

大多数人在看到像这样(或者例如被零除)的浮点异常错误信息时,他们都非常困惑,毕竟他们都像我们一样写下了一个保护表达式。他们自然认为是当x小于零时sqrt(x)根本不可能会被计算。但这个假设是错误的。

事件和条件表达式

事件在行为上的作用

若有如下if表达式:

der(x) = if x>=0 then sqrt(x) else 0;

sqrt被以负参数调用是完全有可能的。其原因与这是状态事件有关。请记住,时间事件的发生时间是可以预知的。但是,状态事件则不然。为了确定事件何时会发生,我们必须在解的轨迹寻找以确定条件(例如:x>=0)何时为假。

要了解的重要一点是,直到事件发生时,行为都不会改变。换句话说,这个if表达式的两个值代表了两种类型的行为,der(x)=sqrt(x)der(x)=0。由于x最初大于零,最初的行为便是der(x)=sqrt(x)求解器将继续使用该方程,直到它确定了x>=0所表示的事件发生的时间。而为了确定该事件发生的时间,求解器必须越过条件表达式的值发生变化的点。这意味着,在试图准确确定条件x>=0从真假的变化的时间这一过程中,虽然x已经为负,求解器仍会继续使用公式der(x)=sqrt(x)

大多数用户开始是认为每次der(x)被计算时,if表达式也会被计算(特别是在if表达式中的条件表达式)。但愿前面的章节已经很清楚地阐释了,这并非事实。

我们可以减少在尝试和重试积分步长中所花费的时间。这是得益于Modelica可以从if表达式提取出所谓的“过零”(zero-crossing)函数。这些函数被之所以被称为过零函数,因为它通常被构造成为事件发生的位置有根存在的形式。例如,如果我们有以下if表达式:

y = if a>b then 1 else 0;

这里的过零函数是\(a-b\)。之所以选择这个函数的原因是,它在a>b成立的时刻正好是由正转为负。

回想我们在前面的公式

der(x) = if x>=0 then sqrt(x) else 0;

在此过零函数显然是\(x\)。这是因为事件在\(x\)自己穿过零时发生。

Modelica编译器收集在模型中所有的过零函数以给予积分器使用。在积分时,积分器检查是否有任何的过零函数变了号。如果这些过零函数变了号,那积分器便使用在该步长中计算出的解去内插过零函数,以此计算交叉函数变号的时刻。 而这就是事件发生的时间点。这种方法非常有效率。因为求根算法可以使用更多的信息来帮助他们识别为根的位置(如过零函数的导数)。加上由于它不涉及采取额外的积分步骤, 而仅从触发事件的积分步长开始计算内插函数的值,这种求解算法本身非常花销非常低。

事件的抑制

但经过这一切,如何避免我们在Decay1Decay2模型中所看到的问题,仍然是不明确的。答案是一个特殊的操作符,noEventnoEvent操作符抑制了这种特殊的事件处理方式。相反,它提供大多数用户一开始预期中的行为。也就是,在每一个x的取值都会进行条件表达式的计算。我们可以从下面的模型中看到noEvent操作符的作用:

model Decay3
  Real x;
initial equation
  x = 1;
equation
  der(x) = if noEvent(x>=0) then -sqrt(x) else 0;
end Decay3;

结果如下:

../../../_images/Decay3.png

现在仿真毫无问题地完成了。这是因为noEvent确保了sqrt(x)不会在x值为负时被调用。

这似乎有点奇怪,我们必须显式地加入noEvent操作符才能得到我们认为最直观的行为。为什么不把默认行为设为最直观的一个呢?答案是性能原因。使用条件表达式生成事件提高了仿真的性能。因为这给予求解器何时要为行为突变做准备的相关线索。大多数时候,这种方法不会导致任何问题。我们在本章中提出了的例子都是为了强调这个问题,但这些例子并不能代表大多数情况。出于这个原因,noEvent不是默认的,而是必须显式地使用。但是应当指出的是,noEvent操作符应该只用于行为变化是平滑过渡的时候。否则它会带来性能问题。

抖振

用Modelica语言时,你早晚会遇到一个被称为“抖动”(chattering)常见现象。考虑下面的模型:

model WithChatter
  Real x;
initial equation
  x = 2;
equation
  der(x) = if x>=1 then -1 else 1;
end WithChatter;

实际上,这个模型的行为是,对x的任何初值,状态x都会线性渐近到1。在数学上讲,一旦x值到达1,它的值便不再会变更。因为,任何从1的偏移,不论正负,都将立即导致它返回到1。

但我们不会用一个严格的数学方法去求解这些方程。相反,我们将使用浮点数表示法表示变量,而用数值积分器计算结果。因此,我们只能满足于有限的精度和积分误差。这些误差的净效应是,x的轨迹不会完全保持1而将偏离略微上方和下方。而每次发生这种情况都会生成一个事件。

对此模型进行仿真给了我们如下结果:

../../../_images/WC.png

这种模型可能会引起一种被称为“抖动”的现象。抖动其实只是由于求解器在大量事件产生时,由于过度地缩短时间步长所引致的仿真速度下降。如果我们观察在仿真过程中所使用的CPU时间,抖动对仿真性能的影响其实是显而易见的。当x接近1时,CPU时间便开始急剧上升。原因是在幕后事件造成大量的非常小的时间步长而显着地增加所执行的计算的数量。这是因为在幕后,事件造成大量非常小的时间步长,从而让所执行的计算量也急剧增加。WithChatter这个例子之所以重要,在于本例有一个显而易见的解析解,但即使如此,本例仍然受到事件的高频率产生而减低了仿真性能。

这是另一种noEvent运算符有所助益的情况。因为我们知道,这if表达式并不引入任何行为上的突然变化,我们可以用noEvent操作符把条件表达式包裹起来。如下:

model WithoutChatter
  Real x;
initial equation
  x = 2;
equation
  der(x) = if noEvent(x>=1) then -1 else 1;
end WithoutChatter;

这样做以后,我们会得到几乎相同的解,但仿真性能更佳:

../../../_images/WOC.png

请注意x到达1的前后,CPU占用率的斜率没有明显的变化。这与在WithChatter的情况对比下有着显著差异。

在现实中,像这样的方程是罕见的。在本例中,我们使用一个极端情况试图清楚地显示抖动所带来的影响。这里所描述的行为并不是特别现实或物理的。在这种情况下,我们夸大了抖动产生的效果以清楚地表明仿真性能受到的影响。

抖动在现实世界中的更典型的例子包括,在某个稳定点周围生效的条件表达式(Decay2就是一个很好的例子)。在这种情况下就会发生抖动,因为该系统往往会自然稳定到其中在条件表达式的发生点或者其附近。而由于精确度和数值计算上的考虑,与条件式相关联的事件会被频繁触发。倘若系统有许多部件在这样的平衡点附近,这个效果便会被加剧了。

速度与准确性

希望讨论至今大家已经清楚,为什么在某些情况下有必要抑制事件的产生。不过,大家可能会问,为什么不直接跳过事件,每次都计算条件表达式?那么,让我们花一些时间来探讨这个问题,并解释为何就整体而言把条件表达式与事件进行关联是个非常好的主意[2]

如果没有事件检测,积分器只会直接跳过事件。而如果这种情况发生,积分器会错过重要的行为改变。这将会对仿真的准确性有显著的影响。这是因为大多数积分程序的准确性是基于对被积函数及其导数连续性的假设。倘若这些假设不再成立了,我们需要让积分程序知道这点。如此,积分器便可以考虑这些行为的改变。

这就到事件出场的时候了。事件强迫积分器在行为发生变化的时刻停止,然后在发生行为变化后的时点重新启动。其结果是在牺牲模拟速度的前提下获得更高的精度让我们看一个具体的例子。请考虑下面这个简单的Modelica模型:

model WithEvents "Integrate with events"
  parameter Real freq = 1.0;
  Real x(start=0);
  Real y = time;
equation
  der(x) = if sin(2*Modelica.Constants.pi*freq*time)>0 then 2.0 else 0.0;
end WithEvents;

我们可以从这个系统看到在一半的时间里x的微分为2。而在另外一半时间里x的微分则为0。所以,在每一个周期中,``x``平均的微分应该是1。这意味着在每一个周期结束时,x``y``应该相等。

倘若我们使用WithEvents对模型进行仿真,我们会得到如下的结果:

../../../_images/WE.png

请注意在每个周期结束时,xy的轨迹如何在一点汇合。这是潜在的积分精度的视觉指示。即使我们增加基本周期的频率,我们可以看到这个特点仍然存在:

../../../_images/WEf.png

不过,在使用完全相同的积分参数的情况下,现在让我们考虑通过noEvents操作符以抑制事件产生的情况:

model WithNoEvents "Integrate without events"
  parameter Real freq = 1.0;
  Real x(start=0);
  Real y = time;
equation
  der(x) = noEvent(if sin(2*Modelica.Constants.pi*freq*time)>0 then 2.0 else 0.0);
end WithNoEvents;

在这种情况下,积分器不能得知行为变化。尽管它会尽量准确地进行积分,但如果积分器没有明确知道这些行为改变发生,那么它会继续盲目地使用错误的微分值,并在行为变化后很长时间继续进行外推。如果使用相同的积分器设置对使用了WithNoEvents的模型进行仿真,我们可以看到两者的结果有如何显著的不同:

../../../_images/WNE.png

请注意积分器何其迅速地引入了些非常显著的错误。

../../../_images/WNEf.png

在这些例子中使用的积分器设定是为了说明noEvent操作符对精度可能造成的影响而设置的。然而,选择这些设置的目的确实是为了强调这些差异。如果使用更典型的设置,结果上的差异很可能不这么会引人注目。此外,使用noEvent的影响无法预测或量化。因为,这些影响在使用不用的求解器时会显著有所不同。但核心的问题名曲:使用noEvent操作符可能对仿真结果的精度有显著影响。

[1]

这个模型不会总会仿真失败。失败的出现取决于仿真引入多少积分误差。而积分误差则又取决于所使用的数值容差。

[2]

特别感谢Lionel Belmon对我原来讨论的质疑,以及找出了我的几个未经证实的假设。因为如此,现在的解释有了很大的改进,而且还加上了仿真结果以支持我的结论。