SMILELAND

善始者实繁,克终者盖寡。


  • Home

  • About

  • Categories

  • Archives

  • Tags

  • Search

2020-02-27-关于克服拖延症的5个建议

Posted on 2020-02-27

关于克服拖延症的5个建议

  1. 改变想法

    一个事情我今天不去做,就不会再去做了,想到一个事情,若必须去做,那么今天就去做,可以做任务分解,然后分解到每一天,若任务太多,将任务分级,做更重要的。

  2. 每天进行任务划分

    每天做任务计划,然后按照计划去完成,每天都要完成,最开始可以划的粗糙一点,任务量轻一点,然后逐步加大。

  3. 仪式感

    每天做学习的的任务之前要给自己定一个仪式,就和睡眠仪式一样的,有助于进入状态

  4. 每日打卡

    对每天的完成情况打卡,有助于积累成就感。

  5. 对自己加油!

    加油吧!

Read more »

2020-02-25-微分方程(3)

Posted on 2020-02-25

The Calculus You Need

学习微分方程需要用到一些前置的知识,主要有以下几方面

1. 特殊的一些函数的微分

$x^n$ $sin(x)$ $cos(x)$ $e^x$ $ln(x)$

清楚以上特殊函数的微分。

2. 微分操作定律

$f+g$ $f(x)+g(x)$ $\frac { f ( x ) } { g ( x ) }$ $f(g(x))$

清楚以上方程式的微分。

3. 一个求解公式

用$y ( t ) = \int _ { 0 } ^ { t } e ^ { t - s } 2 ( s ) d s$用来求解$\frac { d } { d x } \int _ { 0 } ^ { x } y ( t ) d t = y ( x )$

我们可以对其证明,对左边的式子左右两端进行微分。

$\frac { d y ( t ) } { d t } = \frac { d } { d t } \left( e ^ { t } \cdot \int _ { 0 } ^ { t } e ^ { - s } q ( u ) d s \right)$

得到

$\frac { d y ( t ) } { d t } = e ^ { t } \int _ { 0 } ^ { t } e ^ { - s } g ( x ) d s + g(t)$

即

$\frac { d } { d x } \int _ { 0 } ^ { x } y ( t ) d t = y ( x )$

4. Taylor级数

我们用一阶导数来形容斜率,我们用二阶导数来形容bending,我们用三阶导数以及更高阶导数去修正曲线,在$t+\Delta t$这段,也就是$\frac{df}{dt}$段更加拟合原函数$f(x+\Delta t)$

$f(t+\Delta t) = f(t) + \Delta t \frac {df(t)}{dt} + \frac{1}{2}(\Delta t)^2 \frac{d^{2}t}{dt^{2}}+ \cdots + \frac { 1 } { n ! } ( \Delta t ) ^ { n } \frac { d ^ { n } f } { d t ^ { n } } + \cdots$

Read more »

2020-02-25-微分方程(2)

Posted on 2020-02-25

Overview of Differential Equations

https://www.youtube.com/watch?v=ghjOS7Q82s0&list=PLUl4u3cNGP63oTpyxCMLKt_JmB0WtSZfG&index=2

一个典型的一阶的微分方程

$\frac { d y } { d t } = a y + q ( t )$

该方程描述的是y的变化率由y自己的解来决定,同时会在等式右端引入一个激励。这是一个线性的y的方程,因为y是一次的。

$\frac { d y } { d t } = f ( y )$

这可能是一个非线性的,因为$f(y)$可以是y的任何函数,可以是$y^2$、$sin(y)$、$e^y$等等可能形式。

这里可以形象的理解,y是弹簧的长度,$\frac { d y } { d t }$是弹簧长度的变化率。

二阶微分方程

$\frac { d ^ { 2 } y } { d t ^ { 2 } } = - k y$

二阶微分方程可理解为加速度,也是区间的弯折率(beding of a curve),一阶微分方程给的是曲线的方向,是上升还是下降,二阶微分方程给的是曲线向上弯还是向下弯,也就是凸函数还是凹函数。也可以理解为力的变化加速度,而线性的与y有关。

更一般的形式

$m y ^ { \prime \prime } + b y ^ { \prime } + k y = f(t)$

$m y ^ { \prime \prime }$理解为加速度乘以质量,第二项$b y ^ { \prime }$表示阻尼,速度的变化量,第三项是一些强制条款,取决于y自己,$f(t)$是一些外部因素。

我们可解的是线性常微分方程,Linear Constant Coefficients,这样的微分方程我们认为是good equations。

对于不好的方程,比如非线性方程,比如不是常微分方程,即系数(m/k/b)是不变的方程,则需要用数值方法求解。

一个系统经常不止是一个方程,是Systems of n equations.

$\frac { d y } { d t }= A y$其中A是矩阵 $\frac { d ^ { 2 } y } { d t ^ { 2 } } = - S y$

特征值和特征向量。

数值解法很重要Numerical Solutions

目前很好用的方法是Ode45.

而Euler方法可看做是Ode1,Ode45具有更高的精度。“Ode 4 and 5”

偏微分方程,有2个变量,以下是热力学方程,时间以及x是空间向量。

$\frac { \partial u } { \partial t } = \frac { \partial ^ { 2 } u } { \partial x ^ { 2 } }$

以上是一个很重要常微分方程。(PDE,Partial differential equations)偏微分方程

$\frac { \partial ^ { 2 } u } { \partial t ^ { 2 } } = \frac { \partial ^ { 2 } u } { \partial x ^ { 2 } }$

上面是波动方程。

拉普拉斯方程。

$\frac { \partial ^ { 2 } u } { d x ^ { 2 } } + \frac { \partial ^ { 2 } u } { d y ^ { 2 } } = 0$

Read more »

2020-02-23-Cassie电弧模型离散化

Posted on 2020-02-23 | In 解问题

对于Cassie电弧方程的离散化。

Read more »

2020-02-22-如何在matlab sfunction 函数中调用自己写的函数?

Posted on 2020-02-22

引用

enter description here

自己编写了一个s函数,有几个参数引用了自己写的几个函数,在脚本中可以正确运行,但在写成s函数,进行 simulink 仿真的时候,已知提示“too many input auguments”,不知道怎么回事。

经过调试发现,把那几个参数换成常值,simulink 可以正常运行。说明:s函数中,不能直接调用自己写的函数。

那么问题来了:

在用 sfuntmpl 模板用m语言写的s函数中,不能直接调用自己写的函数,怎么办?

具体的解决办法是:

1、打开主页上的‘set path’选项,把自己写的函数所在的文件夹加到搜素路径中。

2、在matlab的安装文件路径 X:\Program Files\MATLAB\R2012\toolbox\ 下,新建一个文件夹 ‘myTool’。把自己写的函数复制到此文件夹下,然后把次文件夹加到路径中。

任选一种方法即可。个人比较推荐2。

Read more »

2020-02-21-常微分方程解法(1)

Posted on 2020-02-21

在进行电弧模型建模时候,传统的cassie电弧模型的写法都是连续的,但是连续状态的采样率设置上就出现了问题,无法按照自己期望的采样率进行仿真,所以需要对模型进行离散化处理。

转自csdn

enter description here

建立微分方程只是解决问题的第一步通常需要求出方程的解来说明实际现象,并加以检验。如果能得到即系形式的解固然是便于分析和应用的,但是我们知道,只有线性常系数微分方程,并且自由项是某些特殊类型的函数时,才可以肯定得到这样的解,而绝大多数变系数方程、非线性方程都是所谓“解不出来”的,即使看起来非常简单的方程如$\frac{d y}{d x}=y^{2}+x^{2}$,对于用微分方程解决实际问题来说,数值解法就是一个十分重要的手段。

1. 常微分方程离散化

下面主要讨论一阶常微分方程的初值问题,其一般形式是

\[\left\{\begin{array}{l} {\frac{d y}{d x}=f(x, y) \quad a \leq x \leq b} \\ {y(a)=y_{0}} \end{array}\right.\tag{1}\]

在下面讨论中,我们总假定函数$f(x,y)$连续,且关于y满足Lipschitz条件,即存在常数L,使得

\[|f(x, y)-f(x, \bar{y})| \leq L|y-\bar{y}|\]

这样,由常微分方程理论知,处置问题(1)的解必定存在唯一。

数值解法

所谓数值解法,就是求问题(1)的解$y(x)$在若干点

$a=x_{0}<x_{1}<x_{2}<\cdots<x_{N}=b$处的近似值$y_n(n=1,2,…,N)$,$y_n(n=1,2,…,N)$称为问题(1)的数值解,$h_n=x_{n+1}-x_n$称为由$x_n$到$x_{n+1}$的步长,今后如无特别说明,我们总取步长为常量h。

建立数值解法,首先要将微分方程离散化,一般采用以下几种方法:

(i)用差商近似导数——缠粉方程初值问题

若用向前差商$\frac{y\left(x_{n+1}\right)-y\left(x_{n}\right)}{h}$代替$y^{\prime}\left(x_{n}\right)$带入(1)中的微分方程,则得

\[\frac{y\left(x_{n+1}\right)-y\left(x_{n}\right)}{h} \approx f\left(x_{n}, y\left(x_{n}\right)\right) \quad(n=0,1, \cdots)\]

简化得

\[y\left(x_{n+1}\right) \approx y\left(x_{n}\right)+h f\left(x_{n}, y\left(x_{n}\right)\right)\]

如果用$y(x_n)$的近似值$y_n$带入上式右端,所得结果作为$y(x_{n+1})$的近似值,记为$y_{n+1}$,则有

\[y_{n+1}=y_{n}+h f\left(x_{n}, y_{n}\right) \quad(n=0,1, \cdots)\tag{2}\]

这样,问题(1)的近似解可通过求解下述问题

\[\left\{\begin{array}{l} {y_{n+1}=y_{n}+h f\left(x_{n}, y_{n}\right) \quad(n=0,1, \cdots)} \\ {y_{0}=y(a)} \end{array}\right.\tag{3}\]

按式(3)由初值$y_0$可逐次算出$y_1,y_2…$。式(3)是个离散化的问题,成为差分方程初值问题。

(ii)用数值积分方法

将问题(1)的解表成积分形式,用数值积分方法离散化。例如,对微分方程两端积分,得

\[y\left(x_{n+1}\right)-y\left(x_{n}\right)=\int_{x_{n}}^{x_{n+1}} f(x, y(x)) d x \quad(n=0,1, \cdots)\tag{4}\]

右边的积分用矩形公式或梯形公式计算。

(iii)Taylor多项式近似

将函数$y(x)$在$x_n$处展开,取一次Taylor多项式近似,则得

\[y\left(x_{n+1}\right) \approx y\left(x_{n}\right)+h y^{\prime}\left(x_{n}\right)=y\left(x_{n}\right)+h f\left(x_{n}, y\left(x_{n}\right)\right)\]

再将$y(x_n)$的近似值$y_n$代入上式右端,所得结果作为$y(x_{n+1})$的近似值$y_{n+1}$,得到离散化的计算公式

\[y_{n+1}=y_{n}+h f\left(x_{n}, y_{n}\right)\]

以上三种方法都是将微分方程离散化的常用方法,每一类方法又可导出不同形式的计算公式。其中Taylor展开法,不仅可以得到求数值解的公式,而且容易估计截断误差。

2 欧拉(Euler)方法

2.1 向前Euler公式、向后Euler公式

Euler方法就是用查分方程初值问题(3)的解来近似微分方程初值问题(1)的解,即由公式(3)依次算出$y(x_n)$的近似值$y_n$,($n=1,2,…$)。这组公式求问题(1)的数值解成为向前Euler公式。

如果在微分方程离散化时,用向后差商代替导数,即$y^{\prime}\left(x_{n+1}\right) \approx \frac{y\left(x_{n+1}\right)-y\left(x_{n}\right)}{h}$,则得计算公式

\[\left\{\begin{array}{l} {y_{n+1}=y_{n}+h f\left(x_{n+1}, y_{n+1}\right) \quad(n=0,1, \cdots)} \\ {y_{0}=y(a)} \end{array}\right.\tag{5}\]

用这组公式求问题(1)的数值解成为向后Euler公式。

向后Euler法与Euler法形式上相似,但实际计算时却复杂的多。向前Euler公式是显式的,可直接求解。向后Euler公式的右端含有$y_{n+1}$,因此是隐式公式,一般要用迭代法求解,迭代公式通常为

\[\left\{\begin{array}{l} {y_{n+1}^{(0)}=y_{n}+h f\left(x_{n}, y_{n}\right)} \\ {y_{n+1}^{(k+1)}=y_{n}+h f\left(x_{n+1}, y_{n+1}^{(k)}\right) \quad(k=0,1,2, \cdots)} \end{array}\right.\tag{6}\]

2.2 Euler方法的误差估计

对于向前Euler公式(3)我们看到,当n = 1,2,……时公式右端的$y_n$都是近似的,所以用它计算的$y_{n+1}$会有累积误差,分析累积误差比较复杂,这里先讨论比较简单的所谓局部截断误差。

假定用(3)式右端的$y_n$没有误差,即$y_n=y(x_n)$,那么由此算出

\[y_{n+1}=y\left(x_{n}\right)+h f\left(x_{n}, y\left(x_{n}\right)\right)\tag{7}\]

局部截断误差指的是,按(7)式计算由$x_n$到$x_{n+1}$这一步的计算值$y_{n+1}$与精确值$y(x_{n+1})$之差$y(x_{n+1})-y_{n+1}$。为了估计它,由Taylor展开得到的精确值$y(x_{n+1})$是

\[y\left(x_{n+1}\right)=y\left(x_{n}\right)+h y^{\prime}\left(x_{n}\right)+\frac{h^{2}}{2} y^{\prime \prime}\left(x_{n}\right)+O\left(h^{3}\right)\tag{8}\]

(7)、(8)两式相减(注意到$\left.y^{\prime}=f(x, y)\right)$)得

\[y\left(x_{n+1}\right)-y_{n+1}=\frac{h^{2}}{2} y^{\prime \prime}\left(x_{n}\right)+O\left(h^{3}\right) \approx O\left(h^{2}\right)\tag{9}\]

即局部截断误差是$h^2$阶的,而数值算法的精度定义为:

若一种算法的局部截断误差为$O(h^{p+1})$,则称该算法具有p阶精度。

显然p越大,方法的精度越高。式(9)说明,向前Euler方法是一阶方法,因此它的精度不高。

3 改进的Euler方法

3.1 梯形公式

利用数值积分方法将微分方程离散化时,若用梯形公式计算式(4)中之右端积分,即。

enter description here

这就是求解初值问题(1)的梯形公式。

直观上容易看出,用梯形公式计算数值积分要比矩形公式好。梯形公式为二阶方法。 梯形公式也是隐式格式,一般需用迭代法求解,迭代公式为

enter description here

如果实际计算时精度要求不太高,用公式(10)求解时,每步可以只迭代一次,由此导 出一种新的方法—改进 Euler 法。

3.2 改进Euler法

按式(5)计算问题(1)的数值解,如果每步只迭代一次,相当于将Euler公式与梯形公式结合使用:先用Euler公式求$y_n+1$的一个初步近似值$\bar{y}_{n+1}$,称为预测值,然后用梯形公式校正求得近似值$y_n+1$,即

\[\left\{\begin{array}{l} {\bar{y}_{n+1}=y_{n}+h f\left(x_{n}, y_{n}\right)} \\ {y_{n+1}=y_{n}+\frac{h}{2}\left[f\left(x_{n}, y_{n}\right)+f\left(x_{n+1}, \bar{y}_{n+1}\right)\right]} \end{array}\right.\tag{11}\]

式(11)称为由Euler公式和梯形公式得到的预测——校正系统,也叫改进Euler法。

为了便于编制程序上机,式(11)常改写成

\[\left\{\begin{array}{l} {y_{p}=y_{n}+h f\left(x_{n}, y_{n}\right)} \\ {y_{q}=y_{n}+h f\left(x_{n}+h, y_{p}\right)} \\ {y_{n+1}=\frac{1}{2}\left(y_{p}+y_{q}\right)} \end{array}\right.\tag{12}\]

改进Euler法是二阶方法。

4 龙格——库塔(Runge-Kutta)方法

这部分以后再补,暂时用不上。

https://blog.csdn.net/qq_29831163/article/details/89703598

5 Matlab解法

https://blog.csdn.net/qq_29831163/article/details/89703911

5.1 非刚性常微分方程的解法

Matlab 的工具箱提供了几个解非刚性常微分方程的功能函数,如 ode45,ode23, ode113,其中 ode45 采用四五阶 RK 方法,是解非刚性常微分方程的首选方法,ode23 采用二三阶 RK 方法,ode113 采用的是多步法,效率一般比 ode45 高。 Matlab 的工具箱中没有 Euler 方法的功能函数。

(1)对简单的一阶方程的初值问题

enter description here

自己编写改进的Euler方法函数eulerpro.m如下:

function [x,y]=eulerpro(fun,x0,xfinal,y0,n);
if nargin<5,n=50;end 
h=(xfinal-x0)/n;
x(1)=x0;y(1)=y0;
for i=1:n
    x(i+1)=x(i)+h;
    y1=y(i)+h*feval(fun,x(i),y(i));
    y2=y(i)+h*feval(fun,x(i+1),y1);
    y(i+1)=(y1+y2)/2;
end 
Read more »

2020-02-18-Level 2 M S函数常见问题

Posted on 2020-02-18
  1. 可能提示没有对应的function

看一看对应的function在setup模块中是否被引用。

  1. 可能提示Dwork不能有空的。。。

看一下Dwork设定参数的个数,与实际设定的Dwork参数个数是否相符。

  1. 不同的function中的参数无法互相通用,也是程序无法正确使用的原因。

Dwork有效性试错

如果Dwork可以看做是全局变量,那么在InitializeConditions中的Dwork的值是可以带入到Output子函数的,现在试试。

设置0.1s,step应该变值,但是很显然没有变化,所以这样是不行的。

enter description here

也就是说InitializeConditions并不是执行一次就初始化一次么?

调整将Dwork放到Derivatives模块中,这样就相当于原来的了,应该就好用了。

function Derivatives(block)

tau = block.DialogPrm(1).Data;
Uc = block.DialogPrm(2).Data; 

if block.CurrentTime <= block.DialogPrm(3).Data
    block.Dwork(3).Data = 0;
else
    block.Dwork(3).Data = 1;
end
    
block.Derivatives.Data = block.Dwork(3).Data * 1/tau * (block.InputPort(1).data^2 / Uc^2 - 1);
%endfunction

enter description here

这样是可以的。

Read more »

2020-02-16-Level 2 M S函数模板

Posted on 2020-02-16

官方模板写的比较啰嗦,找起来比较费力,重新做一个中文的模板。

function msfunCassie(block)
% smileland作品,免费使用,但请不要删除此句话
% 以下设定参数并不一定全部需要,不需要的注释掉就好
setup(block);

%endfunction

function setup(block)
  %block.CurrentTime 只读对象,调用当前时间

  % 初始化对话框参数,即对话框内给定的初始化参数.
  block.NumDialogPrms     = 2;%给定2个参数
  %block.DialogPrmsTunable = {'Tunable','Nontunable','SimOnlyTunable'};%参数是否可变(此句话可选是否需要,默认不可变)

  %%%%%%%%初始化输入输出端口的属性%%%%%%%
  
  % 初始化输入输出端口个数.
  block.NumInputPorts  = 1;%输入端口个数
  block.NumOutputPorts = 1;%输出端口个数
  
  % 设定输入输出端口时继承还是动态,以指示输入和输出端口从模型继承其已编译属性(尺寸、数据类型、复杂度和采样模式).
  block.SetPreCompInpPortInfoToDynamic;%设定输入端口为动态
  block.SetPreCompOutPortInfoToDynamic;%设定输出端口为动态

  % 设定输入端口属性.(也可以再Dwork中设定,在PostPropagationSetup子方法中)
  block.InputPort(1).DatatypeID  = 0;  % (-1 inherited;0 double;1 single;2 int8;3 uint8;4 int16;5 uint16;6 int32;7 uint12;8 boolean)
  block.InputPort(1).Complexity  = 'Real';%是否为复数
  
  % 设定输出端口属性.(也可以再Dwork中设定,在PostPropagationSetup子方法中)
  block.OutputPort(1).DatatypeID  = 0; % double
  block.OutputPort(1).Complexity  = 'Real';
  
  block.InputPort(1).Dimensions = 1;%输入口1的维数是1,最高支持到2维
  % block.InputPort(1).DimensionsMode = 'Fixed';%Fixed表示端口维数固定,Variable表示端口维数可变
  block.OutputPort(1).Dimensions = 1;%输出口1的维数是1,最高支持到2维
  % block.OutputPort(1).DimensionsMode = 'Fixed';%Fixed表示端口维数固定,Variable表示端口维数可变
  
  block.InputPort(1).DirectFeedthrough = true;%输入口有无直接馈通(true/false)
  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  
  % 设定有多少连续状态变量.
  block.NumContStates = 1;
  
  % 设定有多少离散状态变量.
  % block.NumDworkDiscStates = 1;

  % 设定采样时间.
  %  [0 0.5]            : 从0.5s开始连续采样
  %  [0.1 0.5]          : 从0.5s开始,每0.1s采样一次
  %  [-1, 0]            : 继承采样频率
  %  [-2, 0]            : 动态变化的采样频率
  block.SampleTimes = [0 0];%连续采样率
  
  %设定离散工作向量个数
  % block.NumDworks = 1;
  
  % -----------------------------------------------------------------
  % Options,若使用TLC则需要以下语句,否则注释掉即可
  % -----------------------------------------------------------------

  % block.SetAccelRunOnTLC(false);
  
  % block.OperatingPointCompliance = 'Default';
  
  % -----------------------------------------------------------------
  % 引用使用的子方法,主要有PostPropagationSetup;InitializeConditions;
  % Start;Outputs;Updates;Derivatives;Terminate.以下子方法仅保留着几个。
  % -----------------------------------------------------------------
  
  %{
  %包含PostPropagationSetup子方法,用于初始化DWork工作向量
  block.RegBlockMethod('PostPropagationSetup', @DoPostPropSetup);
  %}
  
  %包含InitializeConditions子方法,用于初始化状态变量或DWork工作向量的值
  block.RegBlockMethod('InitializeConditions', @InitializeConditions);
  
  %与InitializeConditions子方法功能一致,但只初始化一次,Ini方法每次调用子方法时均会初始化
  % block.RegBlockMethod('Start', @Start);

  %包含输出模块
  block.RegBlockMethod('Outputs', @Outputs);

  %包含离散状态更新模块
  %block.RegBlockMethod('Update', @Update);

  %包含连续状态微分模块
  block.RegBlockMethod('Derivatives', @Derivatives);
  
  %包含终止模块
  %block.RegBlockMethod('Terminate', @Terminate);

  %还有其他模块,仅将最主要的列出,若需要则使用S函数自带模板。
%endfunction

% -------------------------------------------------------------------
% 以下为子方法模块,同样主要有DoPostPropSetup;InitializeConditions;
  % Start;Outputs;Updates;Derivatives;Terminate.以下子方法仅保留这几个。
% -------------------------------------------------------------------
%{
function DoPostPropSetup(block)
  block.NumDworks = 4;
  
  block.Dwork(1).Name            = 'tau';
  block.Dwork(1).Dimensions      = 1;
  block.Dwork(1).DatatypeID      = 0;      % double
  block.Dwork(1).Complexity      = 'Real'; % real
  block.Dwork(1).UsedAsDiscState = false;
  
  block.Dwork(2).Name            = 'Uc';
  block.Dwork(2).Dimensions      = 1;
  block.Dwork(2).DatatypeID      = 0;      % double
  block.Dwork(2).Complexity      = 'Real'; % real
  block.Dwork(2).UsedAsDiscState = false;
  
  block.Dwork(3).Name            = 'step';
  block.Dwork(3).Dimensions      = 1;
  block.Dwork(3).DatatypeID      = 0;      % double
  block.Dwork(3).Complexity      = 'Real'; % real
  block.Dwork(3).UsedAsDiscState = false;
  
  block.Dwork(4).Name            = 'u';
  block.Dwork(4).Dimensions      = 1;
  block.Dwork(4).DatatypeID      = 0;      % double
  block.Dwork(4).Complexity      = 'Real'; % real
  block.Dwork(4).UsedAsDiscState = false;
  
  % Register all tunable parameters as runtime parameters.
  block.AutoRegRuntimePrms;

%endfunction
%}
%

function InitializeConditions(block)

block.ContStates.Data = log(1e-6);
    
%{
function Start(block)

  tau = block.DialogPrm(1).Data;
  Uc = block.DialogPrm(2).Data; 

%}

%{
  function Update(block)
  block.Dwork(1).Data = [];
%}

function Outputs(block)
  
  block.OutputPort(1).Data = block.InputPort(1).data * exp(block.ContStates.Data);
  
  
function Derivatives(block)

tau = block.DialogPrm(1).Data;
Uc = block.DialogPrm(2).Data; 
  
if block.CurrentTime <= 0.02
    step = 0;
else
    step = 1;
end
    
block.Derivatives.Data = step * 1/tau * (block.InputPort(1).data^2 / Uc^2 - 1);
%endfunction

%endfunction
  
%function Terminate(block)

%disp(['Terminating the block with handle ' num2str(block.BlockHandle) '.']);

%endfunction

Read more »

2020-02-13-S-Function学习笔记

Posted on 2020-02-13

10. S函数

10.1 概述

作者分析,整个Simulink引擎调用通用库的模块时候使用的方法应该和S函数差不多。

10.2 S函数的类型

类型实际上可以分成2种,一种是语言类型,最主要的有M语言和C语言。支持功能多少分为Level1和Level2。M语言能实现的功能少,但相对简单,C Mex S函数能实现非常多的功能。Level1能实现单输入和单输出,目前已不建议使用,应使用Level2函数。

C Mex S函数有以下几个优点

  1. 实现功能多
  2. 仿真速度快
  3. 可与外部设备进行交互
  4. 可以通过Matlab为设备生成C代码

10.3 S函数的要素

S函数的原理无论哪种语言,都是一样的,一个模块包括输入、输出以及模型内部变量,以及一个时间变量。

enter description here

状态量x:状态量为模块内部的计算量或缓存量。所谓状态量,根据系统性质分为连续系统中的微分量和离散系统中的差分量。

enter description here

10.4 S函数的组成及执行顺序

一个S函数有几个子方法构成,表征了S函数甚至整个Simulink引擎的工作过程。这些子方法包括模块初始化、采样时刻计算、模块输出的计算方法,模块离散状态量的更新方法、连续状态变量的微分方法和仿真结束前的终止方法等。

enter description here

采样时间:即采样间隔。对于Simulink模型来说,解算器中的一个步长决定了整个模型最小的采样时间间隔,模型中的模块的步长,可以设置-1集成系统或父层模块的采样时间,也可以设置明确的数字决定执行间隔,但是此时这个间隔应该是系统采样时间间隔的整数倍。采样时间也可以是inf,表示不进行采样,即使常数。对于连续模块,设置采样时间为0。

运行顺序优先度:对于有输入/输出的模块来说,当输入/输出是直接馈入时,输入/输出需要再同一个采样时刻计算完毕,那么要求输入值准备好之后,才能进行输出值的计算。这样的模块必然依赖于连线的前一个模块的输出。对于输入/输出之间不是直接馈入的模块,输出值与状态值相关,与当前输入值无关,可以先行根据状态值计算出输出值,再处理输入值以更新状态值。对于没有输入/输出的模块,它本身不依赖前序模块,也不是后续模块的依赖,其计算顺序则是另外一种情况。模块的执行顺序也收到模块输入/输出以及直接馈入的影响。

S函数的执行顺序:

  • 模型初始化之后,如果是采样时间变化(下称变采样时间)的模块,需要在每个采样时刻计算出下一个采样时刻;
  • 如果不是变采样时间的模块,则无需计算下一个采样时刻,直接计算模块的输出Outputs。
  • 紧接着再进行更新子方法(Update)执行,主要用于离散状态变量等的更新。
  • 接下来计算积分环节(integration),由多个子环节构成、主要用于处理连续状态的计算和更新及非采样过零检测的实施。
  • 对于存在连续状态的S函数,Simulink将会以小步长进行输出方法Outputs与微分计算方法derivative的调用。
  • 对于非采样过零检测,Simulink将会调用输出方法Outputs及过零检测子方法zero-crossings以精确定位过零的位置。
子方法 作用
初始化 在第一个采样时间的仿真之前运行的函数,用来初始化模块,包括设定输入/输出端口的个数和维数,输入是否直接馈入,参数个数设定采样时间,当使用工作向量时还需要为其分配存储空间
下一个采样点时间计算 根据模型解算器的算法求得下一个采样时刻点,通常用于变步长模块
输出函数计算 在每一个major step计算模型所有的输出口的输出值
离散状态更新 在每一个major step都进行一次离散状态的更新计算,在输出函数之后
积分计算 如果模型具有连续状态,才采取此方法,将major step分隔为数个minor step,在每一个minor step里进行一次输出函数与积分计算。积分计算主要用来更新连续状态。当模型存在非采样过零检测时,还会在minor step中进行过零检测
模型终止 当仿真终止时调用的子函数,用于清除不用的变量,释放内存空间等动作

minor step与major step是变步长解算器相关的知识点,后者表示连个采样点之间的步长;前者表示major step内为了精确计算积分,将次步长再划分为更小的多个步长。

10.5 编写S-Function

Level1 M : 支持简单的MATLAB接口及少数S函数API。

Level2 M : 支持扩展S函数API及代码生成功能,使用场合更加广泛。

C MEX 提供更灵活的变成方式,既可手写C代码也可以调用既存的C/C++代码或Fortran代码。要求掌握很多C MEX S函数 API 用法及TLC代码编写方法,才能够定制具有代码生成功能的C MEX S-Function。

10.5.1 Level1 M S-Function

函数原型为:[sys,x0,str,ts]=f(t,x,u,flag,p1,p2,…),其中f是S函数的函数名。Simulink会在仿真过程的每个步长内多次调用f。flag的值随着仿真过程自动变化,其值对应的S函数子方法如表所示。

flag值 Level1 M S函数子方法名 说明
0 mdlInitializeSizes 定义S函数的基本属性,如输入/输出维数、连续/离散状态变量个数、采样时间,以及输入是否为直接馈入
1 mdlDerivatives 连续状态变量的微分函数,这里通常通过给定的微分计算表达式,通过积分计算得到状态变量的值
2 mdlUpdate 更新离散状态变量
3 mdlOutputs 计算S函数的输出
4 mdlGetTimeOfNextVarHit 尽在变离散采样时间情况下(ts = [-2 0]时)使用,用于计算下一个采样时刻的绝对时间,若模块不是变步长此函数不会执行
9 mdlTerminate 在仿真结束时执行一些必要的动作,如清除临时变量,或显示提示信息等

注:模块是否直接馈入有一个简单的判断方法,就是查看mdlOutputs和mdlGetTimeOfNextVarHit两个子方法中有没有用到输入u。如果用到了,这个输入端口的直接馈入direct feedthrough就必须设置为1。

Level 1 M S函数的采样时间ts,适用于Level 2 M S函数和C MEX S 函数。有一个包含两个元素的数组[m,n]表示,m为模块的采样时间周期,或连续或离散或继承父层模型的采样时间,表示每隔多长时间采样一次;n为采样时间的偏移量,表示与采样时间周期所表示的时刻点的偏差时间。常见的设置见表。

采样时间表示 意义
[0 0] 连续采样时间
[-1 0] 继承S函数输入信号或父层模型的采样时间
[0.5 0.1] 离散采样时间,从0.1秒开始每0.5s采样一次

一个S函数可以同时存在多个采样速率,可以通过设置多维数组来实现多个速率采样时间的配置,如[0.25 0; 1 0.1],在这个S函数中,采样时刻为[0 0.1 0.25 0.5 0.75 1 1.1…]。

例子,通过Level1 M S函数来实现这个方程。

enter description here

可以很容易得知:

A = [-0.5572 -0.7814;0.7814 0]; B = [1 -1; 0 2]; C = [1.9691 6.4483]。

为S函数设置3个参数,用来传递矩阵A、B、C到S函数中,在Update和Output子方法中分别使用其中部分参数,Update的参数列表中仅需要使用A、B,Outputs子方法中仅需要使用C。Level1 M S函数构成的模块最多仅支持单个输入口,单个输出口。所以输入口的维数要设置为2.状态变量采用离散状态变量,在S函数中按列排序。然后需要判断直接馈入,由于输出y由状态变量x决定,不由输入u决定,所以不存在直接馈入。根据分析,可以编写代码如下。

function [sys,x0,str,ts,simStateCompliance] = sfun_state01(t,x,u,flag,A,B,C)

switch flag,
  case 0,
    [sys,x0,str,ts,simStateCompliance]=mdlInitializeSizes;
  case 1,
    sys=mdlDerivatives(t,x,u);
  case 2,
    sys=mdlUpdate(t,x,u,A,B);
  case 3,
    sys=mdlOutputs(t,x,u,C);
  case 4,
    sys=mdlGetTimeOfNextVarHit(t,x,u);
  case 9,
    sys=mdlTerminate(t,x,u);
  otherwise
    DAStudio.error('Simulink:blocks:unhandledFlag', num2str(flag));
end

function [sys,x0,str,ts,simStateCompliance]=mdlInitializeSizes
sizes = simsizes;

sizes.NumContStates  = 0;
sizes.NumDiscStates  = 2;
sizes.NumOutputs     = 1;
sizes.NumInputs      = 2;
sizes.DirFeedthrough = 0;
sizes.NumSampleTimes = 1;   % at least one sample time is needed

sys = simsizes(sizes);

x0  = [0 0]';
str = [];
ts  = [0 0];

simStateCompliance = 'UnknownSimState';

function sys=mdlDerivatives(t,x,u)

function sys=mdlUpdate(t,x,u,A,B)

sys = A * x + B * u;

function sys=mdlOutputs(t,x,u)

sys = C * x;

function sys=mdlGetTimeOfNextVarHit(t,x,u)

sampleTime = 1;    %  Example, set the next hit to be one second later.
sys = t + sampleTime;

function sys=mdlTerminate(t,x,u)

sys = [];

% end mdlTerminate

在S函数对话框中填入参数。

enter description here

不建议使用全局变量的原因

  1. 全局变量会破坏模块化设计,增加模块之间的耦合性。当系统复杂之后,依赖全局变量过多,会造成难以管理的混乱局面。建议尽量不适用全局变量,以养成好习惯。
  2. 不出现多次的global声明,代码看起来清爽简约并充分发挥了S函数的作用。

离散变量通过明确的表达式来进行值的更新。而连续状态下,mdlDerivatives子方法的计算方式有所不同。mdlDerivatives子方法中编写的表达式计算出来的值会自动进行一次积分,得到的值才更新到连续状态变量中去,作为下一个采样时刻的连续状态值进行下一次的计算。积分的S函数代码为:

function [sys,x0,str,ts,simStateCompliance] = sfun_hyo(t,x,u,flag)

switch flag,
  case 0,
    [sys,x0,str,ts,simStateCompliance]=mdlInitializeSizes;
  case 1,
    sys=mdlDerivatives(t,x,u);
  case 2,
    sys=mdlUpdate(t,x,u);
  case 3,
    sys=mdlOutputs(t,x,u);
  case 4,
    sys=mdlGetTimeOfNextVarHit(t,x,u);
  case 9,
    sys=mdlTerminate(t,x,u);
  otherwise
    DAStudio.error('Simulink:blocks:unhandledFlag', num2str(flag));
end

function [sys,x0,str,ts,simStateCompliance]=mdlInitializeSizes
sizes = simsizes;

sizes.NumContStates  = 0;
sizes.NumDiscStates  = 1;
sizes.NumOutputs     = 1;
sizes.NumInputs      = 1;
sizes.DirFeedthrough = 0;
sizes.NumSampleTimes = 1;   % at least one sample time is needed

sys = simsizes(sizes);

x0  = [0];
str = [];
ts  = [0 0];

simStateCompliance = 'UnknownSimState';

function sys=mdlDerivatives(t,x,u)

sys = u;

function sys=mdlUpdate(t,x,u,A,B)

sys = [];

function sys=mdlOutputs(t,x,u)

sys = x;

function sys=mdlGetTimeOfNextVarHit(t,x,u)

sampleTime = 1;    %  Example, set the next hit to be one second later.
sys = t + sampleTime;

function sys=mdlTerminate(t,x,u)

sys = [];

% end mdlTerminate

这个函数中,初始化输入/输出端口为一维,一个初始值为0的连续状态变量,连续采样时间,输入端口没有直接馈入。

多然由于对S函数运行原理、参数传递、直接馈入等方法或概念不太熟悉,导致编写S函数时遇到各种错误。

  1. flag = 3时仿真出错,往往由于输出维数配置不当或使用了输入u但是没有在初始化时将输入端口配置为直接馈入。
  2. 参数未定义,命名S函数参数列表中设置了参数却没有添加到每个子方法的参数列表中导致子方法中参数未定义错误。

10.5.2 Level2 M S-Function

Level 2 M S函数调用模块如图所示。

enter description here

Level 2 M S函数使得用户能够使用MATLAB语言来编写支持多个输入/输出端口的自定义模块,并且每个端口都能够支持包括矩阵在内的Simulink支持的所有数据类型。

Level 2 M S函数提供了一系列API设置模块属性和定义各个子方法,其中Setup和Outputs两个方法是必不可少的。

1.Setup子方法

Setup子方法是Level 2 M S函数整体中唯一调用的语句,对模块的属性和其他子方法进行初始化,Setup子方法类似Level 1 M S函数中mdlInitializeSizes子方法的功能,并且相比之下功能更强大。在Setup中可以设置多输入多输出,并且每个输出的端口信号的维数可以是标量或矩阵以及可变维数,S函数的其他子方法也是通过Setup子方法进行注册的,可以成为Level 2 M S函数的根本。Setup子方法实现以下功能:

  • 设定模块输入输出端口的个数;
  • 设定每一个端口的数据类型、数据维数、实数复数性和采样时间等;
  • 规定模块的采样时间;
  • 设定S函数参数的个数;
  • 注册S函数的子方法(将子方法函数的句柄传递到实时对象的RegBlockMethod函数的对应属性中)。

Setup子方法的参数为block,是一个Level 2 M S函数的实时对象,包括了一些属性和方法。其属性成员如表

Level 2 M S函数实时对象的属性列表

实时对象属性成员 说明 实时对象属性成员 说明
NumDialogPrms 模块GUI参数个数 NumContStates 连续状态变量数目
NumInputPorts 输入端口数 NumDworkDiscStates 离散状态变量个数
NumOutputPorts 输出端口数 NumRuntimePrms 运行时参数个数
BlockHandle 模块句柄,只读 SampleTimes 产生输出的模块的采样时间
CurrentTime 当前仿真时间,只读 NumDworks 离散工作向量个数

通过block变量的点操作符来访问上表中属性并且可以直接使用等号赋值:

block.NumDialogPrms = 0;

模块的输入/输出端口又包含自己的属性,其中常用的属性列表如表所示。

Level 2 M S函数端口的属性列表

端口属性名 说明
Dimensions 端口数据维数
DatatypeID/Datatype 端口数据类型,可通过ID号制定也可以直接制定数据类型名
Complexity 端口数据是否为复数
DirectFeedthrough 端口数据是否直接馈入
DimensioinsMode 端口维数是固定或可变(fixed/variable)

当模块中存在多个端口时,需要对每一个端口进行属性的设定,使用端口访问方法(InpuPort或OutputPort成员)及端口号索引(正整数)设定输入/输出端口,如:

block.InputPort(1).Dimensionis = 1;

Setup子方法中,出来能够通过block实时对象域操作符访问模块属性之余,还可以通过子方法获取模块对象的信息,可以调用的方法列表如表所列。

Level 2 M S 函数实时对象的方法列表

实时对象方法 说明 实时对象方法 说明
ContStates 获取模块的连续状态 Dwork 获取Dwork向量
DataTypelsFixedPoint 判断数据类型是否为固定点数 FixedPointNumericType 获取固定点数据类型属性
DatatypeName 获取数据类型的名称 InputPort 获取输入端口
DatatypeSize 获取数据类型大小 OutputPort 获取输出端口
Derivatives 获取连续状态的微分 RuntimePrm 获取运行时参数
DialogPrm 获取GUI中的参数    

Setup子方法中还需要使用RegBlockMethod函数来注册S函数中需要用到的其他各个子方法,如表所列

其他子方法列表

子方法名 说明
PostPropagationSetup 设置工作向量机状态变量的函数(可选)
InitializeConditions 在仿真开始时被调用的初始化函数(可选)
Start 在模型运行仿真时调用一次,用来初始化状态变量和工作向量(可选)
Outputs 在每个步长里计算模型输出
Update 在每个步长里更新离散状态变量的值(可选)
Derivatives 在每个步长里更新连续状态变量的微分量(可选)
Terminate 在仿真结束时调用,用来清除变量内存

表中“可选”表示不是必须存在于S函数中,根据需要进行选取。

下面例子,使用Setup子函数初始化一个带有1个参数、1个输入端口、1个输出端口的模块,输入/输出端口维数都为1,没有直接馈入;模块的采样时间为0.1s,注册4个子方法函数:DoPostPropSetup,InitConditions,Output和Update。

function msfuntmpl(block)
%MSFUNTMPL A Template for a MATLAB S-Function
%   The MATLAB S-function is written as a MATLAB function with the
%   same name as the S-function. Replace 'msfuntmpl' with the name
%   of your S-function.  

%   Copyright 2003-2018 The MathWorks, Inc.
  
%
% The setup method is used to setup the basic attributes of the
% S-function such as ports, parameters, etc. Do not add any other
% calls to the main body of the function.  
%设置方法用于设置S功能的基本属性,例如端口,参数等。请勿在功能主体上添加其他任何调用。
%   

setup(block);%设置模块
  
%endfunction

% Function: setup ===================================================
% Abstract:
%   Set up the S-function block's basic characteristics such as:%设定基础参数
%   - Input ports%输入端口
%   - Output ports%输出端口
%   - Dialog parameters%对话参数
%   - Options%选项
% 
%   Required         : Yes
%   C MEX counterpart: mdlInitializeSizes
%
function setup(block)

  % Register the parameters.初始化对话框参数
  block.NumDialogPrms     = 1;%初始化三个S函数对话框参数
  block.DialogPrmsTunable = {'Tunable','Nontunable','SimOnlyTunable'};%参数是否可变

  % Register the number of ports.设置输入输出端口数
  block.NumInputPorts  = 1;
  block.NumOutputPorts = 1;
  
  % Set up the port properties to be inherited or dynamic.将端口属性设置为继承或动态。以指示输入和输出端口从模型继承其已编译属性(尺寸,数据类型,复杂度和采样模式)
  block.SetPreCompInpPortInfoToDynamic;%设定输入端口属性为动态
  block.SetPreCompOutPortInfoToDynamic;%设定输出端口属性为动态

  % Override the input port properties.覆盖输入端口属性
  block.InputPort(1).DatatypeID  = 0;  % double
  block.InputPort(1).Complexity  = 'Real';
  
  % Override the output port properties.覆盖输出端口属性
  block.OutputPort(1).DatatypeID  = 0; % double
  block.OutputPort(1).Complexity  = 'Real';

  % Hard-code certain port properties
  block.InputPort(1).Dimensions = 1;
  block.InputPort(1).DirectFeedthrough = false;
  block.OutputPort(1).Dimensioins = 1;
  
  % Set up the continuous states.设定具有连续状态变量个数
  % block.NumContStates = 1;

  % Register the sample times.设定采样时间
  %  [0 offset]            : Continuous sample time连续的采样时间
  %  [positive_num offset] : Discrete sample time离散的采样时间
  %
  %  [-1, 0]               : Inherited sample time继承的采样时间
  %  [-2, 0]               : Variable sample time可变采样时间
  block.SampleTimes = [0.1 0];%设定上面的采样时间,设置为[-1 0]即制定S函数具有继承的采样时间
  
   % -----------------------------------------------------------------
   % Register the methods called during update diagram/compilation.
   % -----------------------------------------------------------------
   
  %
  % PostPropagationSetup:
  %   Functionality    : Set up the work areas and the state variables. You can
  %                      also register run-time methods here.
  %   C MEX counterpart: mdlSetWorkWidths
  %
  block.RegBlockMethod('PostPropagationSetup', @DoPostPropSetup);
  
  % 
  % InitializeConditions:
  %   Functionality    : Call to initialize the state and the work
  %                      area values.
  %   C MEX counterpart: mdlInitializeConditions
  % 
  block.RegBlockMethod('InitializeConditions', @InitializeConditions);
  
  % 
  % Outputs:
  %   Functionality    : Call to generate the block outputs during a
  %                      simulation step.
  %   C MEX counterpart: mdlOutputs
  %
  block.RegBlockMethod('Outputs', @Outputs);
  
  % 
  % Update:
  %   Functionality    : Call to update the discrete states
  %                      during a simulation step.
  %   C MEX counterpart: mdlUpdate
  %
  block.RegBlockMethod('Update', @Update);
  
  % 
  % Derivatives:
  %   Functionality    : Call to update the derivatives of the
  %                      continuous states during a simulation step.
  %   C MEX counterpart: mdlDerivatives
  %
  block.RegBlockMethod('Derivatives', @Derivatives);

2.PostPropagationSetup子方法

PostPropagationSetup子方法是用来初始化Dwork工作向量的方法,规定Dwork向量的个数及每个向量的维数、数据类型、离散状态变量的名字和虚实性,以及是否作为离散状态变量使用。

DWork向量时SImulink分配给模型中每个S函数实例的存储空间块。当不同S函数模块之间需要通过全局变量或者静态变量进行数据交互时,必须在S函数模块存在多个拷贝时,每一个模块中全局变量所使用的空间都需要十分小心的分配、修改和释放,特别是在使用C语言进行函数编写时更要进深。因为一旦管理不慎,可能会造成S函数的一个实例的数据被另一个实例的数据所覆盖,导致计算出错并最终使得仿真失败。使用Dwork向量可以有效防止这些情况发生,可重复性(reentrancy)有效地保证了S函数多个实例的可跟踪性,当S函数中使用Dwork向量时,这个S函数的实例就具备了可重入性,Simulink为这个S函数的每一个实例分别管理数据。

Dwork向量既可以作为全局变量在S函数模块之间传递数据,也可以作为某个S函数内部离散状态变量来使用。举一个例子说明在PostPropagationSetup子方法中使用Dwork向量作为离散状态变量的使用方法。在下面代码中,DoPostPropSetup函数作为PostPropagationSetup子方法的回调函数,程序中定义了一个Dwork向量,命名为x0,维数为1,数据类型是double型实数,此Dwork向量作为离散状态使用。

 function DoPostPropSetup(block)
  block.NumDworks = 1;
  
  block.Dwork(1).Name            = 'x0';
  block.Dwork(1).Dimensions      = 1;
  block.Dwork(1).DatatypeID      = 0;      % double
  block.Dwork(1).Complexity      = 'Real'; % real
  block.Dwork(1).UsedAsDiscState = true;

  % Register all tunable parameters as runtime parameters.
  block.AutoRegRuntimePrms;

%endfunction

Dwork向量的数据类型,通常使用DatatypeID来表示,每一个数据类型对应一个整数。

Level 2 M S函数实时对象支持的数据类型列表

数据类型 代表数据类型的整数 数据类型 代表数据类型的整数
inherited -1 int16 4
double 0 uint16 5
single 1 int32 6
int8 2 uint32 7
uint8 3 boolean或定点类型 8

上例中定义Dwork向量数据类型为double,故使用整数0代表。

3. InitializeConditions/Start子方法

InitializeConditions子方法可以用来初始化状态变量或DWork工作向量的值,下面例子中InitConditions函数作为InitializeConditions子方法的回调函数,其内容是通过S函数对话框获取DWork的初始值,连续状态变量的值初始化为1.0。

function InitializeConditions(block)

block.Dwork(1).Data = block.DialogPrm(1).Data;
block.ContStates.Data = 1;

%endfunction

Start子方法跟initializeConditions子方法功能一致,都可以用来初始化离散、连续状态变量和Dwork向量,但是二者也存在区别。Start子方法仅仅在仿真执行开始时初始化一次,而S函数模块放置在使能子系统中时,其InitializeConditions子方法在每次子系统被使能时都会被调用。

4.Output子方法

Output子方法跟Level 1 M S函数的mdlOutputs子方法作用一样,用于计算S函数的输出。下面的代码将Dwork向量值赋给输出端口作为当前采样时刻的输出值。

function Outputs(block)
  
  block.OutputPort(1).Data = block.Dwork(1).Data;
  
%endfunction

5.Update子方法

Update子方法跟Level 1 M S函数中的mdlUpdate子方法作用相同,用于计算离散状态变量的值。下面代码将当前采样时刻的输入端口值幅值给Dwork向量:

function Update(block)
  
  block.Dwork(1).Data = block.InputPort(1).Data;
  
%endfunction

6.Derivatives子方法

Derivatives子方法跟Level 1 M S函数中的mdlDerivatives子方法作用相同用于计算并更新连续状态变量的值。下面代码将当前采样时刻的输入端口值赋给连续状态变量:

function Derivatives(block)

block.Derivatives(1).Data = block.InputPort(1).Data;

%endfunction

7.Terminate子方法

S函数工作的首位工作放在Terminate子方法中进行,如存储空间的释放,变量的删除等。Level 2 M S函数不要求必须包含Terminate子方法。


例子

融合2张图像。

enter description here

两张图片的尺寸是[375×500×3]的uint8类型,即使是Level 2 M S函数也不能处理3维数据矩阵,对多支持二维数据。需要使用M脚本进行预处理,然后再S函数内部实现数据融合的计算。图片的第三维存储信息是每个像素点的RGB信息,若仅取RGB其中一组信息,那么图片就变成灰度图片,同时图片的第三维数据就被取出了从而变成二维矩阵,适合Level 2 M S函数进行处理。模块有两个输入,每个输入的维数是二维数据,375×500的矩阵,数据类型均为uint8.模块没有状态变量和输出端口,故不需要实现S函数的Update,Output子方法。模块在Terminate子方法中实现两个端口输入的矩阵进行融合的代码,并拖过figure展示融合后的图片。

首先在工作空间中将图片的信息读入并取出R维色彩信息以转化为二维矩阵。

m1 = imread('阿布罗狄.jpg');
m1 = m1(:,:,1);
m2 = imread('艾欧里亚.jpg');
m2 = m1(:,:,1);

在模型中可以使用Constant模块将图片转化之后的数据读入,在Constant Value处使用imread函数读取后得到的变量。

enter description here

编写名为sfun_image_merge的S函数:

setup(block);%设置模块
  
%endfunction

% Function: setup ===================================================
% Abstract:
%   Set up the S-function block's basic characteristics such as:%设定基础参数
%   - Input ports%输入端口
%   - Output ports%输出端口
%   - Dialog parameters%对话参数
%   - Options%选项
% 
%   Required         : Yes
%   C MEX counterpart: mdlInitializeSizes
%
function setup(block)

  % Register the number of ports.设置输入输出端口数
  block.NumInputPorts  = 2;
  block.NumOutputPorts = 0;
  
  % Set up the port properties to be inherited or dynamic.将端口属性设置为继承或动态。以指示输入和输出端口从模型继承其已编译属性(尺寸,数据类型,复杂度和采样模式)
  block.SetPreCompInpPortInfoToDynamic;%设定输入端口属性为动态

  % Override the input port properties.覆盖输入端口属性
  block.InputPort(1).DatatypeID  = 3;  % uint8
  block.InputPort(1).Complexity  = 'Real';
  block.InputPort(2).DatatypeID  = 3;  % uint8
  block.InputPort(2).Complexity  = 'Real';

  % Hard-code certain port properties
  block.InputPort(1).Dimensions = [375 500];
  block.InputPort(1).DirectFeedthrough = false;
  block.InputPort(2).Dimensions = [375 500];
  block.InputPort(2).DirectFeedthrough = false;
  
  %Register parameters
  block.NumDialogPrms = 0;
  
  % Set up the continuous states.设定具有连续状态变量个数
  % block.NumContStates = 1;

  % Register the sample times.设定采样时间
  %  [0 offset]            : Continuous sample time连续的采样时间
  %  [positive_num offset] : Discrete sample time离散的采样时间
  %
  %  [-1, 0]               : Inherited sample time继承的采样时间
  %  [-2, 0]               : Variable sample time可变采样时间
  block.SampleTimes = [0 0];%设定上面的采样时间,设置为[-1 0]即制定S函数具有继承的采样时间
  
  block.SimStateCompliance = 'DefaultSimState';
  
   % -----------------------------------------------------------------
   % Register the methods called during update diagram/compilation.
   % -----------------------------------------------------------------
   
  %
  % PostPropagationSetup:
  %   Functionality    : Set up the work areas and the state variables. You can
  %                      also register run-time methods here.
  %   C MEX counterpart: mdlSetWorkWidths
  %
 % block.RegBlockMethod('PostPropagationSetup', @DoPostPropSetup);
  
  % 
  % InitializeConditions:
  %   Functionality    : Call to initialize the state and the work
  %                      area values.
  %   C MEX counterpart: mdlInitializeConditions
  % 
%  block.RegBlockMethod('InitializeConditions', @InitializeConditions);
  
  % 
  % Outputs:
  %   Functionality    : Call to generate the block outputs during a
  %                      simulation step.
  %   C MEX counterpart: mdlOutputs
  %
%  block.RegBlockMethod('Outputs', @Outputs);
  
  % 
  % Update:
  %   Functionality    : Call to update the discrete states
  %                      during a simulation step.
  %   C MEX counterpart: mdlUpdate
  %
%  block.RegBlockMethod('Update', @Update);
  
  % 
  % Derivatives:
  %   Functionality    : Call to update the derivatives of the
  %                      continuous states during a simulation step.
  %   C MEX counterpart: mdlDerivatives
  %
%  block.RegBlockMethod('Derivatives', @Derivatives);

  % 
  % Terminate:
  %   Functionality    : Call at the end of a simulation for cleanup.
  %   C MEX counterpart: mdlTerminate
  %
  block.RegBlockMethod('Terminate', @Terminate);
  
  function Terminate(block)

  imshow((block.InputPort(1).data + block.InputPort(2).data) / 2);

  %endfunction

使用Level 2 M S函数模块,将S函数名sfun_image_merge填入S-Function Name栏,构成如图所示模型。

enter description here

仿真得到图像。

enter description here

编写自己的Level 2 M S函数时可以使用Simulink提供的模板作为基础,然后进行内容增加或删减。msfun_test.m就是其中一个常用的模板。

明天自己编写一个Level 2 M S函数的模板吧。

Read more »

2020-02-12-使用C语言写简单S-Function

Posted on 2020-02-12

1. 例子一

转自csdn博客 enter description here

功能描述:

实现一个双输入,双输出,使用两个参数,使得输出为:

y1 = para1 * u1 + 5 y2 = para2 * u2 + 3

代码实现

代码如下:

#define S_FUNCTION_NAME  test
#define S_FUNCTION_LEVEL 2
#define INPUT_NUM 2
#define OUTPUT_NUM 2
 
#include <stdio.h>
#include "simstruc.h"
#include <math.h>
 
static void mdlInitializeSizes(SimStruct *S)
{
    ssSetNumSFcnParams(S, 2);  /* Number of expected parameters */
    if (ssGetNumSFcnParams(S) != ssGetSFcnParamsCount(S)) {
        /* Return if number of expected != number of actual parameters */
        return;
    }
 
    ssSetNumContStates(S, 0);
    ssSetNumDiscStates(S, 0);
 
    if (!ssSetNumInputPorts(S, 2)) return;
    ssSetInputPortWidth(S, 0, 1);
    ssSetInputPortRequiredContiguous(S, 0, true); /*direct input signal access*/
    ssSetInputPortWidth(S, 1, 1);
    ssSetInputPortRequiredContiguous(S, 1, true);
 
    ssSetInputPortDirectFeedThrough(S, 0, 1);
    ssSetInputPortDirectFeedThrough(S, 1, 1);
 
    if (!ssSetNumOutputPorts(S, OUTPUT_NUM)) return;
    ssSetOutputPortWidth(S, 0, 1);
    ssSetOutputPortWidth(S, 1, 1);
 
    ssSetNumSampleTimes(S, 0.001);
    ssSetNumRWork(S, 0);
    ssSetNumIWork(S, 0);
    ssSetNumPWork(S, 0);
    ssSetNumModes(S, 0);
    ssSetNumNonsampledZCs(S, 0);
 
    /* Specify the sim state compliance to be same as a built-in block */
    ssSetSimStateCompliance(S, USE_DEFAULT_SIM_STATE);
 
    ssSetOptions(S, 0);
}
 
static void mdlInitializeSampleTimes(SimStruct *S)
{
    ssSetSampleTime(S, 0, CONTINUOUS_SAMPLE_TIME);
    ssSetOffsetTime(S, 0, 0.0);
 
}
 
 
static void mdlOutputs(SimStruct *S, int_T tid)
{
  real_T *para1 = mxGetPr(ssGetSFcnParam(S, 0));
  real_T *para2 = mxGetPr(ssGetSFcnParam(S, 1));
 
  const real_T *u1 = (const real_T*) ssGetInputPortSignal(S,0);
  const real_T *u2 = (const real_T*) ssGetInputPortSignal(S,1);
  real_T       *y1 = ssGetOutputPortSignal(S,0);
  real_T       *y2 = ssGetOutputPortSignal(S,1);
  y1[0] = para1[0]*u1[0] + 5;
  y2[0] = para2[0]*u2[0] + 3;
}
 
static void mdlTerminate(SimStruct *S)
{
}
 
#ifdef  MATLAB_MEX_FILE    /* Is this file being compiled as a MEX-file? */
#include "simulink.c"      /* MEX-file interface mechanism */
#else
#include "cg_sfun.h"       /* Code generation registration function */
#endif

代码解析

1. 初始化

static void mdlInitializeSizes(SimStruct *S)
{
    ssSetNumSFcnParams(S, 2);  /* Number of expected parameters */
    if (ssGetNumSFcnParams(S) != ssGetSFcnParamsCount(S)) {
        /* Return if number of expected != number of actual parameters */
        return;
    }
 
    ssSetNumContStates(S, 0);
    ssSetNumDiscStates(S, 0);  //连续和离散状态个数都为0
 
    if (!ssSetNumInputPorts(S, 2)) return;//2个输入
    ssSetInputPortWidth(S, 0, 1);  //1号输入的宽度为1,即维数为1
    ssSetInputPortRequiredContiguous(S, 0, true); /*direct input signal access*/
    ssSetInputPortWidth(S, 1, 1);   //2号输入的维数也为1
    ssSetInputPortRequiredContiguous(S, 1, true);
 
    ssSetInputPortDirectFeedThrough(S, 0, 1);
    ssSetInputPortDirectFeedThrough(S, 1, 1);
 
    if (!ssSetNumOutputPorts(S, OUTPUT_NUM)) return;
    ssSetOutputPortWidth(S, 0, 1);  //1号输出的维数为1
    ssSetOutputPortWidth(S, 1, 1);  //2号输出的维数为1
 
    ssSetNumSampleTimes(S, 0.001);   //采样时间
    ssSetNumRWork(S, 0);
    ssSetNumIWork(S, 0);
    ssSetNumPWork(S, 0);
    ssSetNumModes(S, 0);
    ssSetNumNonsampledZCs(S, 0);
 
    /* Specify the sim state compliance to be same as a built-in block */
    ssSetSimStateCompliance(S, USE_DEFAULT_SIM_STATE);
 
    ssSetOptions(S, 0);
}

2. 控制算法

static void mdlOutputs(SimStruct *S, int_T tid)
{
  real_T *para1 = mxGetPr(ssGetSFcnParam(S, 0));  //获得参数1
  real_T *para2 = mxGetPr(ssGetSFcnParam(S, 1));  //获得参数2
 
  const real_T *u1 = (const real_T*) ssGetInputPortSignal(S,0);  //获得输入u1
  const real_T *u2 = (const real_T*) ssGetInputPortSignal(S,1);  //获得输入u2
  real_T       *y1 = ssGetOutputPortSignal(S,0);  //输出y1
  real_T       *y2 = ssGetOutputPortSignal(S,1);  //输出y2
  y1[0] = para1[0]*u1[0] + 5;  //输出与输入关系
  y2[0] = para2[0]*u2[0] + 3;
}

最后,使用mex 文件名.c编译,使用simulink s-function模块,将保存在同一文件夹下的模型名填入,仿真即可。

enter description here

对于输入都为1,参数分别为1,2,其输出结果:

enter description here

2. 例子二

转自博客园

enter description here

2.1 MATLAB配置

1.先配置一下 mex -setup

mex -setup Please choose your compiler for building external interface (MEX) files:

Would you like mex to locate installed compilers [y]/n? mex -setup

Select a compiler: [1] Borland C++Builder version 6.0 in D:\soft\c++builer [2] Lcc C version 2.4 in D:\MATLAB\sys\lcc [3] Microsoft Visual C/C++ version 6.0 in D:\soft\viual c++

[0] None

Compiler: 2

Please verify your choices:

Compiler: Lcc C 2.4 Location: D:\MATLAB\sys\lcc

Are these correct?([y]/n): y

Try to update options file: C:\Documents and Settings\zjj\Application Data\MathWorks\MATLAB\R14\mexopts.bat From template: D:\MATLAB\BIN\WIN32\mexopts\lccopts.bat

Done . . .

2.打开viual c++ ,新建文本文件,编写s-function,保存Text1.c

3.在matlab中mex Text1.c

4.在simulink中拖出一个s-function模块,改名为Text.c

补充:封装和测试步骤

编好S - 函数后,对相应的模块进行封装和测试: 1)向模型窗口中加入S-function 模块,双击该模块,打开参数设置对话框(Block Parameters),输入 源文件名和用户定义的参数表; 2)选择该模块,按Ctrl + M,打开封装编辑器(Mask Editor)。在Mask type 编辑框中输入模块类型 名;在Icon 页中用MATLAB 的有关绘图命令绘制模块图标;在Initialization 页中添加用户定义的变量 参数,必要时对变量初始化;在Documentation 页中添加模块的说明和帮助文档。 3)给模块命名; 4)建立一个简单仿真模型框图,测试模块功能。

2.2 MATLAB中的S-Function的用法(C语言)

  1. S-Function简介
    S-Function是system-function的缩写。说得简单,S-Function就是用MATLAB所提供的模型不能完全满足用户,而提供给用户自己编写程序来满足自己要求模型的接口。

  2. MEX函数与M文件的区别 第一, MEX 函数能实现的回调函数比M-文件能实现的回调函数要多得多;

第二, MEX 函数直接访问内部数据结构SimStruct,SimStruct 是Simulink 用来保存关于S-function 信息的一个数据结构;

第三, MEX 函数也可使用MATLAB MEX 文件API 直接来访问MATLAB 的工作空间。

如果一个C MEX文件与一个M文件具有相同的名字,则C MEX文件被优先使用,即在S-Function块中使用的是C MEX文件。

  1. 基础知识 3.1 直接馈通(direct feedthrough) 直接馈通表示系统的输出或可变采样时间是否受到输入的控制。

a. 输出函数(mdlOutputs或flag==3)是输入u的函数。即,如果输入u在mdlOutputs中被访问,则存在直接馈通。

b. 对于一个变步长S-Function的“下一个采样时间”函数(mdlGetTimeOfNextVarHit或flag==4)中可以访问输入u。

例如,一个需要其输入的系统(也就是具有直接馈通)是运算y=kXu,其中,u是输入,k是增益,y是输出。

又如,一个不需要其输入的系统(也就是没有直馈通)是一种简单的积分运算:

输出:y=x;

导数:dx/dt=u

其中,x是状态,dx/dt是状态对时间的导数,u是输入,y是输出。

正确设置直接馈通标志是十分重要的,因为它影响模型中块的执行顺序,并可用检测代数环。

3.2 dynamically sized inputs 主要是给出:输入连续状态数目(size.NumContStates),离散状态数目(size.NumDiscStates) ,输出数目(size.NumOutputs),输入数目(size.NumInputs),Direct Feedthrough(size.Dir Feedthrough)。

3.3 setting sample times and offsets setting smaple times and offsets主要设置采样时间.

3.4 Level-1 和Level-2 Level 1 提供一个简单的接口,可与少部分的S函数API交互。Matlab对于这种方式的支持更多的是为了保持与以前版本的兼容,现在推荐采用的是Level 2 S函数。

  1. S-Function实例

enter description here

S-Function的仿真流程

例如要创建一个有1输入(2维),2输出(1维),3个参数,还有全局变量的S-Function。 过程如下:

a. 新建sfunction的C语言文件

打开simulink,点击User-Defined Functions里面的S-Function Examples。这个里面有多个语言版本的模板,有C,C++,Ada,Fortran和M语言的版本,其实都大同小异,只要了解几个函数就很容易使用了。 选择C语言的版本:从S-function模块中选择C-file S-functions里面的Basic C-MEX template。打开后,另存为自己的模块名字,如test.c 。下面我们来分析代码:

#define S_FUNCTION_NAME  test//这里把文件名sfuntmpl_basic修改为test
#define S_FUNCTION_LEVEL 2
#include "simstruc.h"
//程序里面要用到的头文件在这里引用,如“math.h”等。
float global_var; //定义全局变量
static void mdlInitializeSizes(SimStruct *S)
{
 //这个函数用来设置输入、输出和参数的。
    ssSetNumSFcnParams(S, 3);  /*设置参数个数,这里为3 */
    if (ssGetNumSFcnParams(S) != ssGetSFcnParamsCount(S)) {
        return;
    }
    ssSetNumContStates(S, 0);//设置连续状态的个数,缺省为0;
    ssSetNumDiscStates(S, 0);//设置离散状态的个数,缺省为0;
    if (!ssSetNumInputPorts(S, 1)) return;//设置输入变量的个数,这里为1
    ssSetInputPortWidth(S, 0, 2); //设置输入变量0的维数为2
ssSetInputPortRequiredContiguous(S, 0, true); //设置input0的访问方式,true就是临近访问,这样指针的增量后就可以直接访问下个input端口了。
ssSetInputPortDirectFeedThrough(S, 0, 1);// 设置输入端口的信号是否mdlOutputs函数中使用,这儿设置为true。
    if (!ssSetNumOutputPorts(S, 2)) return;//设置输出变量的个数
ssSetOutputPortWidth(S, 0, 1);//设置输出变量0的维数为1维
    ssSetOutputPortWidth(S, 1, 1);//设置输出变量1的维数为1维
ssSetNumSampleTimes(S, 1); //设置采样时间,此处为1s。
    ssSetNumRWork(S, 0);//不管
    ssSetNumIWork(S, 0);
    ssSetNumPWork(S, 0);
    ssSetNumModes(S, 0);
    ssSetNumNonsampledZCs(S, 0);
ssSetOptions(S, 0);
//下面可以写全局变量的初始化程序
global_var=1;
}
static void mdlInitializeSampleTimes(SimStruct *S)//暂时不管
{
    ssSetSampleTime(S, 0, CONTINUOUS_SAMPLE_TIME);
    ssSetOffsetTime(S, 0, 0.0);
 
}
#define MDL_INITIALIZE_CONDITIONS   /* Change to #undef to remove function */
#if defined(MDL_INITIALIZE_CONDITIONS)
 
  static void mdlInitializeConditions(SimStruct *S)//暂时不管
  {
  }
#endif /* MDL_INITIALIZE_CONDITIONS */
#define MDL_START  /* Change to #undef to remove function */
#if defined(MDL_START) 
  static void mdlStart(SimStruct *S)//暂时不管
  {
  }
#endif /*  MDL_START */
static void mdlOutputs(SimStruct *S, int_T tid)//这里填入相关的运算、算法等
{
real_T *para1 = mxGetPr(ssGetSFcnParam(S,0));
real_T *para2 = mxGetPr(ssGetSFcnParam(S,1));
real_T *para3 = mxGetPr(ssGetSFcnParam(S,2));
const real_T *u = (const real_T*) ssGetInputPortSignal(S,0);
real_T       *y1 = ssGetOutputPortSignal(S,0);
real_T       *y2 = ssGetOutputPortSignal(S,1);
y1[0]=u[0]*para1[0]+u[1]*para2[0];
y2[0]=u[1]*para3[0]+u[0]*para1[0];
}
#define MDL_UPDATE  /* Change to #undef to remove function */
#if defined(MDL_UPDATE)
 
  static void mdlUpdate(SimStruct *S, int_T tid)
  {
  }
#endif /* MDL_UPDATE */
#define MDL_DERIVATIVES  /* Change to #undef to remove function */
#if defined(MDL_DERIVATIVES)
  static void mdlDerivatives(SimStruct *S)
  {
  }
#endif /* MDL_DERIVATIVES */
static void mdlTerminate(SimStruct *S)//这里需要把global变量全部初始化,否则下次运行程序时,全局变量还是之前的值。
{
}
 
#ifdef  MATLAB_MEX_FILE    /* Is this file being compiled as a MEX-file? */
#include "simulink.c"      /* MEX-file interface mechanism */
#else
#include "cg_sfun.h"       /* Code generation registration function */
#endif

b. 编译

在matlab的command window 里面输入“mex test.c”,即可将test.c编译为mex文件。

c.调用sfunction

在simulink空间里面拉入sfunction,在s-function name里面填入test,参数里面填入要设定的参数,然后仿真即可。

总结

整体结构并不复杂,但是,如何使用C语言进行算法编程是难点。

比如

real_T *para1 = mxGetPr(ssGetSFcnParam(S,0));

mxGetPr这样的命令如何获取,指针的操作如何进行,数据类型如何定义,如何求微分方程,如何调用.mat文件。

Read more »
1 … 7 8 9 10
月白

月白

Thinking can do nothing, but action will !!!

97 posts
5 categories
157 tags
RSS
GitHub Email
Friends
  • 111qqz
  • VectorLu
  • Yalv
  • YJK
  • 大黄菌
  • Liadbiz
  • CharsTech
  • 五斤
  • YangHan
  • 一之笔
  • DaiChao1997
  • Oubindo
  • Danielhu
  • hhuysqt
  • 瓠樽
  • Leon-wtf
  • Lighting Studio
0%
© 2016 - 2023 月白
Powered by Jekyll
Theme - NexT.Muse
知识共享许可协议