什么是S函数？
S-函数是系统函数（System Function）的简称，在 Simulink 中用非图形化的方式来描述一个模块。一个完整的S-函数结构体系包含了描述一个动态系统所需要的全部能力。使用S-函数用户可以向 Simulink 模型中添加自己的模块，可以自由选择使用 MATLAB、C、C++ 等语言来创建自己的模块。Matlab 语言编写的S-函数可以充分利用 MATLAB 所提供的丰富资源，方便地调用各种工具箱函数和图形函数；使用 C 语言编写的S-函数可以实现对操作系统的访问，如实现与其他进程的通信和同步等。
Simulink 已经提供了大量的内置的系统模块，并且允许用户自定义模块，那么为何还要使用S-函数呢？诚然，对于大多数动态系统仿真分析语言，使用 Simulink 提供的模块即可实现，而无需使用S-函数。
但是，当需要开发一个新的通用的模块作为一个独立的功能单元时，使用S-函数实现则是一种相当简便的方法。另外，由于S-函数可以使用多种语言编写，因此可以将已有的代码结合进来，而不需要在 Simulink 中重新实现算法，从而在某种程度上实现了代码移植。此外，在S-函数中使用文本方式输入公式、方程，非常适合复杂动态系统的数学描述，并且在仿真过程中可以对仿真进行更精确的控制。如果恰当地使用S函数，理论上，可以在 Simulink 下对任意复杂的系统进行仿真。
基础知识
Simulink仿真过程
Simulnk 仿真分为两步：初始化、仿真循环。仿真是由求解器控制的，求解器主要作用是：计算模块输出、更新模块离散状态、计算连续状态。求解器传递给系统的信息包括：时间、输入和当前状态。系统的作用：计算模块的输出、更新状态、计算状态导数，然后将这些信息传递给求解器。求解器和系统之间的信息传递是通过不同标志来控制的。
Simulink 仿真过程：
S函数工作原理
Matlab 为了用户使用方便，有一个S函数的模板 sfuntmpl.m，一般来说，我们仅需要在 sfuntmpl.m 的基础上进行修改即可。在主窗口输入 edit sfuntmpl 即可出现模板函数的内容，可以详细地观察其帮助说明以便更好地了解S函数的工作原理。
S函数的定义形式为function[sys,x0,str,ts]=sfunc(t,x,u,flag,p1,…Pn)，其中的 sfunc 为自己定义的函数名称，以上参数中，t、x、u 分别对应时间、状态、输入信号，flag 为标志位，其取值不同，S函数执行的任务和返回的数据也是不同的，Pn 为额外的参数，sys 为一个通用的返回参数值，其数值根据 flag 的不同而不同，x0 为状态初始数值，str 在目前为止的 matlab 版本中并没有什么作用，一般 str=[] 即可，ts 为一个两列的矩阵，包含采样时间（整个模型的基础采样时间，各个子系统或模块的采样时间，必须以这个步长为整数倍）和偏移量两个参数，如果设置为 [0 0]，那么每个连续的采样时间步都运行，[-1 0] 则表示按照所连接的模块的采样速率进行，[0.25 0.1] 表示仿真开始的 0.1s 后每 0.25s 运行一次，采样时间点为 TimeHit = n*period + offset。S函数的使用过程中有1个概念值得注意：
Direct Feedthrough，系统的输出是否直接和输入相关联，即输入是否出现在输出端的标志，若是为 1，否则为 0，一般可以根据在 flag=3 的时候，mdlOutputs 函数是否调用输入 u 来判断是否直接馈通。
S函数中目前支持的 flag 选择有 0、1、2、3、4、9 几个数值，下面说一下在不同的 flag 情况下S函数的执行情况。
1.flag=0。进行系统的初始化过程，调用 mdlInitializeSizes 函数，对参数进行初始化设置，比如离散状态个数、连续状态个数、模块输入和输出的路数、模块的采样周期个数、状态变量初始数值等。
2.flag=1。进行连续状态变量的更新，调用 mdlDerivatives 函数。
3.flag=2。进行离散状态变量的更新，调用 mdlUpdate 函数。
4.flag=3。求取系统的输出信号，调用 mdlOutputs 函数。
5.flag=4。调用 mdlGetTimeOfNextVarHit 函数，计算下一仿真时刻，由 sys 返回。
6.flag=9。终止仿真过程，调用 mdlTerminate 函数。
在实际仿真过程中，Simulink 会自动将 flag 设置为 0，进行初始化过程，然后将 flag 的数值设置为 3，计算模块的输出，一个仿真周期后，Simulink 将 flag 的数值先后设置为 1 和 2，更新系统的连续和离散状态，再将其设置为 3，计算模块的输出，如此循环直至仿真结束条件满足。
S函数执行顺序：
在S函数的编写过程中，首先需要搞清楚模块中有多少个连续和离散状态，离散模块的采样周期是如何的，同时需要了解模块的连续和离散的状态方程分别是什么，输出如何表示。下面以实例说明S函数的具体应用。
sfuntmpl.m模板代码内容详细解析
function [sys,x0,str,ts,simStateCompliance] = sfuntmpl(t,x,u,flag) % 主函数 % 主函数包含四个输出： % sys数组包含某个子函数返回的值 % x0为所有状态的初始化向量 % str是保留参数，总是一个空矩阵 % Ts返回系统采样时间 % 函数的四个输入分别为采样时间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 % 初始化回调子函数 % 提供状态、输入、输出、采样时间数目和初始状态的值 % 初始化阶段，标志变量flag首先被置为0，S-function首次被调用时 % 该子函数首先被调用，且为S-function模块提供下面信息 % 该子函数必须存在 sizes = simsizes; % 生成sizes数据结构，信息被包含在其中 sizes.NumContStates = 0; % 连续状态数，缺省为0 sizes.NumDiscStates = 0; % 离散状态数，缺省为0 sizes.NumOutputs = 0; % 输出个数，缺省为0 sizes.NumInputs = 0; % 输入个数，缺省为0 sizes.DirFeedthrough = 1; % 是否存在直馈通道，1存在，0不存在 sizes.NumSampleTimes = 1; % 采样时间个数，至少是一个 sys = simsizes(sizes); % 返回size数据结构所包含的信息 x0 = []; % 设置初始状态 str = []; % 保留变量置空 ts = [0 0]; % 设置采样时间 simStateCompliance = 'UnknownSimState'; function sys=mdlDerivatives(t,x,u) % 计算导数回调子函数 % 给定t,x,u计算连续状态的导数，可以在此给出系统的连续状态方程 % 该子函数可以不存在 sys = []; % sys表示状态导数，即dx function sys=mdlUpdate(t,x,u) % 状态更新回调子函数 % 给定t、x、u计算离散状态的更新 % 每个仿真步内必然调用该子函数，不论是否有意义 % 除了在此描述系统的离散状态方程外，还可以在此添加其他每个仿真步内都必须执行的代码 sys = []; % sys表示下一个离散状态，即x(k+1) function sys=mdlOutputs(t,x,u) % 计算输出回调函数 % 给定t,x,u计算输出，可以在此描述系统的输出方程 % 该子函数必须存在 sys = []; % sys表示输出，即y function sys=mdlGetTimeOfNextVarHit(t,x,u) % 计算下一个采样时间 % 仅在系统是变采样时间系统时调用 sampleTime = 1; % 设置下一次采样时间是在1s以后 sys = t + sampleTime; % sys表示下一个采样时间点 function sys=mdlTerminate(t,x,u) % 仿真结束时要调用的回调函数 % 在仿真结束时，可以在此完成仿真结束所需的必要工作 sys = [];
S函数实例
使用 M 文件编写S函数形式有两种，Level 1 和 Level 2，二者的包装模块是不同的。
连续有限输出的积分器
% 实现一个连续有限积分器，其中输出受下限和上限限制，并包括初始条件。 function [sys,x0,str,ts,simStateCompliance]=limintm(t,x,u,flag,lb,ub,xi) switch flag case 0 [sys,x0,str,ts,simStateCompliance] = mdlInitializeSizes(lb,ub,xi); case 1 sys = mdlDerivatives(t,x,u,lb,ub); case {2,9} sys = []; % do nothing case 3 sys = mdlOutputs(t,x,u); otherwise DAStudio.error('Simulink:blocks:unhandledFlag', num2str(flag)); end function [sys,x0,str,ts,simStateCompliance] = mdlInitializeSizes(lb,ub,xi) sizes = simsizes; sizes.NumContStates = 1; sizes.NumDiscStates = 0; sizes.NumOutputs = 1; sizes.NumInputs = 1; sizes.DirFeedthrough = 0; sizes.NumSampleTimes = 1; sys = simsizes(sizes); str = []; x0 = xi; ts = [0 0]; simStateCompliance = 'DefaultSimState'; function sys = mdlDerivatives(t,x,u,lb,ub) if (x <= lb & u < 0) | (x>= ub & u>0 ) sys = 0; else sys = u; end function sys = mdlOutputs(t,x,u) sys = x;
一连续系统
蹦极跳是一种挑战身体极限的运动，蹦极者系着一根弹力绳从高处的桥梁或山崖向下跳。如果蹦极者系在一个弹性系数为 k kk 的弹力绳索上。定义绳索下端的初始位置为 0，则蹦极者受到的弹性力是
其中 m mm 为物体的质量，g gg 为重力加速度取 10，x xx 为物体的位置，第二项为物体受到的弹性力，第三项和第四项表示空气的阻力。桥梁距离地面的高度为 50m,绳子的原长为 30m，弹性系数 k kk 为 50，a 1 = a 2 = 1 a_1 = a_2 = 1a 1 =a 2=1。试判断质量为 70kg 的人是否能够安全地享受此游戏带来的乐趣。下图为整个系统的示意图。
function [sys,x0,str,ts,simStateCompliance] = jumping(t,x,u,flag,len,m,height,k) switch flag, case 0, [sys,x0,str,ts,simStateCompliance]=mdlInitializeSizes(len,m,height,k); case 1, sys=mdlDerivatives(t,x,u,len,m,height,k); case 2, sys=mdlUpdate(t,x,u); case 3, sys=mdlOutputs(t,x,u,len,m,height,k); 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(len,m,height,k) sizes = simsizes; sizes.NumContStates = 2; sizes.NumDiscStates = 0; sizes.NumOutputs = 1; sizes.NumInputs = 0; sizes.DirFeedthrough = 1; sizes.NumSampleTimes = 1; sys = simsizes(sizes); x0 = [-len;0]; str = []; ts = [0 0]; simStateCompliance = 'UnknownSimState'; function sys=mdlDerivatives(t,x,u,len,m,height,k) if x(1)>0 b=-k*x(1); else b=0; end sys = [x(2);10+b/m-1/m*x(2)-1/m*abs(x(2))*x(2)]; function sys=mdlUpdate(t,x,u) sys = []; function sys=mdlOutputs(t,x,u,len,m,height,k) sys = [height-len-x(1)]; function sys=mdlGetTimeOfNextVarHit(t,x,u) sampleTime = 1; sys = t + sampleTime; function sys=mdlTerminate(t,x,u) sys = [];
s-function 参数的配置：
搭建 simulink 框图观察上述系统输出
综上所述，此蹦极系统并不能让 70kg 的人嗨。
一离散系统s-function 参数的配置：
搭建 simulink 框图观察上述系统输出
综上所述，此蹦极系统并不能让 70kg 的人嗨。
一离散系统
function [sys,x0,str,ts] = dsfunc(t,x,u,flag) A=[-1.3839 0.5097 1.0000 0]; B=[-2.5559 0 0 4.2382]; C=[ 0 2.0761 0 7.7891]; D=[ -0.8141 -2.9334 1.2426 0]; switch flag, case 0, [sys,x0,str,ts] = mdlInitializeSizes(A,B,C,D); case 2, sys = mdlUpdate(t,x,u,A,B,C,D); case 3, sys = mdlOutputs(t,x,u,A,C,D); case 9, sys = []; % do nothing otherwise DAStudio.error('Simulink:blocks:unhandledFlag', num2str(flag)); end function [sys,x0,str,ts] = mdlInitializeSizes(A,B,C,D) sizes = simsizes; sizes.NumContStates = 0; sizes.NumDiscStates = size(A,1); sizes.NumOutputs = size(D,1); sizes.NumInputs = size(D,2); sizes.DirFeedthrough = 1; sizes.NumSampleTimes = 1; sys = simsizes(sizes); x0 = ones(sizes.NumDiscStates,1); str = []; ts = [1 0]; function sys = mdlUpdate(t,x,u,A,B,C,D) sys = A*x+B*u; function sys = mdlOutputs(t,x,u,A,C,D) sys = C*x+D*u;
