昨天和 TC 聊了聊求解器的事,就顺便记录一下。
这里不涉及数值求解算法的细节,我只从用户角度来谈谈 Simulink 求解器的基本概念。
我们先来看这个模型:
运行模型,随意打开一个 Scope,比如下面的这个。这时,你会看到一条熟悉的线,这就是我们常查看的“信号”,也就是求解器计算出来的输出。
如果我们改变线型。
你会看到,我们曾经查看的曲线,实际上是一个个时间轴上的点。
也就是说,模型求解,实际上就是在“各个时间点”上,计算模块,更新各个“模块的输出”。
那么它是按照什么顺序来进行更新的呢。
下图右上角红色数字:求解器按照红色数字标识的顺序更新计算模块输出。
这里有两个问题:
我们来看一个简单例子,现在目的是要更新右边那个输出1模块。
如果是你,你怎么确定当前时刻的模块更新顺序?
好,我们来看结果
这就是模型求解的基本回溯原则。
那我们再来看这个例子,如果是你,你如何确定它们的执行顺序呢?
想想之前的基本原则。
所以,类似于增益这种“当前时刻的输出一定要根据当前时刻输入才能计算”的特质,叫做 Direct feedthrough,我没翻译过它的专用名词,要么叫“直通”也行。
而像Unit Delay这类模块就可以打断这种“回溯”链条。
有没有发现,上面的这个序号跳过了一个 4 。这是因为:Simulink 模型在求解的时候会根据情况添加各种隐藏“模块”。
我们用“原始的Debug”模式查看一下这个简单模型里到底隐藏了什么。
下图列表,看到第四个了么 0:4。
发现求解器在 Gain 后面的位置添加了一个模块,用来输出。因为是隐藏模块,所以我们在模型框图界面上是看不到的。
我们把这个信号log输出取消。
重新刷新模型,那个隐藏编码消失了。
我们现在把中间两个模块做成子系统。
如果还记得我们一开始关于大括号的问题?
现在就明白了,大括号其实就是方便显示“虚拟子系统”内部的模块序号。之所以叫虚拟子系统,就是因为它是为了布局好看,在模型计算上不产生任何影响。
内部和之前一样:
如果我们把它变成原子子系统
再排序,你会看到排序只显示一个1。
也就是说求解器把它当做一个整体来看待。轮到计算它时,会跳进去把它内部的所有模块都计算完毕,再跳出来计算其它模块。
可能是我这个例子太简单,好像模块执行顺序一样的嘛,看不出来。
那我们再来看最开始那个f14。
看那个面积最大的 Aircraft Dynamics Model模块,它右上角{1,5,18...}。这就意味着它会跳进去算一部分,再跳出来算其它模块,再跳进去算另一部分。
如果把它设置成 Atomic Subsystem,它右上角显示一个8,这也就意味着把它当做一个整体来看,跳进去把里面的模块更新完毕,在跳出来算其它的地方。
所以如果你的子系统是个严格且独立具有某个特定功能的单元,原子子系统可以很好的保证你所测试的单元,在任何一个地方被使用时,它的计算过程都保持不变。
当然,它内部模块求解还是有顺序的。
最后再来看这个模型
现在似乎能很容易知道求解器是怎么确定它的计算顺序以及计算过程了。假如把它变化一下。
当我们求解这个模型,回溯到加法的时候,发现加法需要输出信号才能计算。
这些 Direct feedthrough 模块组成的闭环是代数环。
要么用一些具有内部状态量的模块来强行打断。
要么,改变模型的打发,去掉或许本来就可以不存在的代数环。
可以看到此刻求解器把代数环涉及到的回路,当做一个隐形的原子子系统,单独拿出来,通过迭代的方式求解。
实际上确定模块计算顺序还有很多其它的规则,但到目前为止,对于理解后面的部分也足够了。暂时也还没有需要研究茴香豆有几种写法的需求,于是到此为止。
定步长与变步长
好我们现在知道了,模型求解其实就是在各个时刻去更新每个模块的输出。
那是哪些时刻呢?
定步长求解器,从输入步长h的那一刻起,Simulink 模型求解器要工作的时刻点就已经被确定下来了,0,h,2*h,3*h,4*h......
求解器只需要在指定的时刻点根据各自的算法计算模型就行了。
而变步长求解器还多一个工作,就是要确定步长。
那它怎么确定呢?先看一下求解器的名字,ode就是Ordinary Differential Equations,后边的数字是阶数。
举个例子,比如 ode45,这个求解器怎么确定步长的呢?
其实就是先假设一个步长h',用 ode4 算一遍,再用 ode5 算一遍,如果他俩的差值很小很小了,那么就认为这步长h'合适了,否则就继续减少步长。
那,怎么评价这个“很小很小”呢,就是根据你在这里设置的误差:
多速率系统
到目前为止,是不是感觉对基本的求解规则已经了解了?
那么再来看个离散模型的例子,一个模块设置的采样时间 0.75,另一个设置为采样时间0.5. 看求解器怎么确定模型的步长的。
有时候我们用 Simulink Test等工具自动生成 Test Harness 的时候,会看到有这样奇怪的采样时间的设置:[1 0].
比如下图,手动设置为 [1,0.1]也就是[采样时间,偏移量],就表示这个模块的刷新时刻变成了 1.1,2.1,3.1
那为什么不直接设置为1.1呢,因为这样采样时间就变成1.1,2.2,3.3了,和原来不一样。
从 Simulink 模型仿真输出的tout,能看出来求解时刻。
从信号图上可以看到,(看对应的横坐标是时间)上面那个模块(蓝线)的更新时刻就是3.1,而下面的模块(红线)更新时刻1.4。看右侧的时刻表更清楚。
模型引用
求解时刻需要匹配。通常我们在做系统集成的时候,会用到 Model Reference 模块。而做代码生成的时候,通常又要配置成定步长求解器,因此经常会遇到各种关于采样时间不匹配的报错,会看不懂。
1如果 Top Model 是定步长求解器,则被引用模型必须也是定步长求解器。求解器步长要匹配,Top Model永远都是基础求解步长。
因此,但具体是哪一个定步长求解器,没关系。比如下图一个ode1,一个ode8.
但求解器步长要匹配,各个被引用模型的求解步长必须是基础步长的整数倍。
比如下图,反过来设置被引用模型步长更短为 0.05,而Top Model求解步长反而为0.1,就不行。
3.如果 Top Model 是变步长求解器,一般来说没什么特别的约束。
但是保不齐会有很混乱的子系统,比如左图一个采样时间为 0.7 的离散模块且同时还有一个连续模块。
这个子模型自己是可以跑没问题的,但是如果把它集成到变步长求解器的 Top Model,就会报错。
就是刚才表格里的第三种情况,这时候根据提示修改一下模型求解器,匹配就行了。
模型求解
我们最后再来看模型求解的基本含义,可能又有不同的体会了。
其实就是算出当前状态X(n),再计算步长h,导数 X(n),由此算出下一个状态X(n+1);
重复这个过程,算出每一个时间点上的信号。