最近看了一下无刷直流电机的相关概念及仿真，看到大多数的文献仿真中都使用到了S函数，因此下了点资料看了一番，在本博文中简单地说一下S函数的概念及使用。

    S函数即系统函数System Function的意思，为什么要使用S函数呢？是因为在研究中，有时需要用到复杂的算法设计等，而这些算法因为其复杂性不适合用普通的Simulink模块来搭建，即matlab所提供的Simulink模块不能满足用户的需求，需要用编程的形式设计出S函数模块，将其嵌入到系统中。如果恰当地使用S函数，理论上，可以在Simulink下对任意复杂的系统进行仿真。

    S函数具有固定的程序格式，用matlab语言可以编写S函数，此外还允许用户使用C、C++、Fortran和Ada等语言进行编写，用非matlab语言进行编写时，需要采用编译器生成动态链接库DLL文件。

在主窗口中输入sfundemos，或者点击Simulink->User-Defined Functions->S-Function Examples，即可出现如图1所示的界面，可以选择对应的编程语言查看演示文件。

http://s7/middle/4b013fb1n984f8b1da8b6&690

图1 S函数范例库

Matlab为了用户使用方便，有一个S函数的模板sfuntmpl.m，一般来说，我们仅需要在sfuntmpl.m的基础上进行修改即可。在主窗口输入edit sfuntmpl即可出现模板函数的内容，可以详细地观察其帮助说明以便更好地了解S函数的工作原理。模板函数的定义形式为function[sys,x0,str,ts]=sfuntmpl(t,x,u,flag)，一般来说，S函数的定义形式为[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函数的使用过程中有2个概念值得注意：1、direct feedthrough，系统的输出是否直接和输入相关联，即输入是否出现在输出端的标志，若是为1，否则为0，一般可以根据在flag=3的时候，mdlOutputs函数是否调用输入u来判断是否直接馈通。2、dynamically sized inputs，主要给出连续状态的个数、离散状态的个数、输入数目、输出数目和直接馈通否。

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函数。

http://s4/middle/4b013fb1n985c6576a053&690

图2 不同flag情况下S函数执行情况

在实际仿真过程中，Simulink会自动将flag设置为0，进行初始化过程，然后将flag的数值设置为3，计算模块的输出，一个仿真周期后，Simulink将flag的数值先后设置为1和2，更新系统的连续和离散状态，再将其设置为3，计算模块的输出，如此循环直至仿真结束条件满足。 

在S函数的编写过程中，首先需要搞清楚模块中有多少个连续和离散状态，离散模块的采样周期是如何的，同时需要了解模块的连续和离散的状态方程分别是什么，输出如何表示。下面以实例说明S函数的具体应用。

实例一：在进行BLDC的速度电流双闭环控制时候，如果电流采用滞缓控制，那么速度环的输出作为参考电流的输入Is，A、B、C三相的参考电流分别为I_ar、I_br、I_cr，转子的位置记为Pos，那么转子位置和三相参考电流之间的关系表如图3所示。

http://s8/middle/4b013fb1n9850079a2f67&690

图3 转子位置和三相参考电流的关系表

分析系统情况：输入为两路，一是角度Angle，通过mod(angle,2*pi)转化为Pos，二是参考电流输入Is，输出为三路，分别为A、B、C的三相参考电流，离散和连续的状态均为0，输入端出现在输出端，将S函数命名为given_current，其程序代码如下所示：

%参考电流模块

%Author:dingqian12345@126.com

%很明显：2个输入分别为角度和速度控制器的输出、3个输出为三相电流的参考电流、为直通模型

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

clc;

switch flag

    case 0

        [sys,x0,str,ts]=mdlInitializeSizes;     %初始化

    case 1

        sys=[];                                 %连续状态的更新

    case 2

        sys=[];              %离散状态的更新

    case 3

        sys=mdlOutputs(u);                      %求取系统的输出信号

    case 4

        sys=[];                                 %计算下一时刻的仿真时间

    case 9

        sys=[];                                 %终止仿真

    otherwise

        error(['Unhandled flag=',num2str(flag)]);

end

%%%%%在flag=0的时候进行整个系统的初始化

function [sys,x0,str,ts]=mdlInitializeSizes

sizes = simsizes;            %读入初始化参数模板

sizes.NumContStates  = 0;    %连续状态个数

sizes.NumDiscStates  = 0;    %离散状态个数

sizes.NumOutputs     = 3;    %输出变量个数

sizes.NumInputs      = 2;    %输入信号个数

sizes.DirFeedthrough = 1;    %输入直接传入输出信号否

sizes.NumSampleTimes = 1;    % at least one sample time is neededWO  一般来说为1个

sys = simsizes(sizes);

x0=[];    %状态初始化

str=[];

ts=[-1 0];    %采样周期若写成-1表示继承其输入信号

%%%%%%在flag=1的时候进行连续系统状态的更新

function sys=mdlDerivatives(t,x,u)

sys = [];

%%%%%%在flag=2的时候进行离散系统状态的更新

function sys=mdlUpdates

%sys(1,1)=x(1)+T*x(2);                %为什么会写成这样呢

%sys(2,1)=x(2)+T*fst2(x,u,r,h);

%%%%%%在flag=3的时候进行系统输出信号的求取

function sys=mdlOutputs(u)

m=current_dq(u(1),u(2));

sys(1,1)=m(1);

sys(2,1)=m(2);

sys(3,1)=m(3);

%%%%%%在flag=4的时候进行下一时刻仿真时间的计算

function sys=mdlGetTimeOfNextVarHit(t,x,u)

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

%sys = t + sampleTime;

%%%%%%在flag=9的时候终止仿真过程

function sys=mdlTerminate(t,x,u)

%sys = [];

    如果某个flag对应的函数不起作用，我们将其所调用的函数设为空或者sys=[]，为了保持模块化以便思路清晰和方便扩充，一般不要使用case {1,4,9 }的形式，而是一个一个的列出。在以上代码中，在flag=3的时候调用了自编的子程序m=current_dq(u(1),u(2))，current_dq的代码如下：

%%%%%%用户自定义子程序

function x=current_dq(angle,current)

pos=mod(angle,2*pi);

x(1)=current;

x(2)=-current;

x(3)=0;

if 0<=pos && pos<pi/3

     x(1)=current;

     x(2)=-current;

     x(3)=0;

end

if pi/3<=pos && pos<2*pi/3

     x(1)=current;

     x(2)=0;

     x(3)=-current;

end

if 2*pi/3<=pos && pos<pi

    x(1)=0;

    x(2)=current;

    x(3)=-current;

end

if pi<=pos && pos<4*pi/3

    x(1)=-current;

    x(2)=current;

    x(3)=0;

end

if 4*pi/3<=pos && pos<5*pi/3

    x(1)=-current;

    x(2)=0;

    x(3)=current;

end

if 5*pi/3<=pos && pos<2*pi

   x(1)=0;

   x(2)=-current;

   x(3)=current;

end

    S函数的代码完成后，下面进行Simulink模块的搭建，如图4所示，其中的S-Function模块提取路径为Simulink->User-Defined Functions->S-Function，双击S-Function模块，在name行中输入函数名称given_current即可，仿真结果如图5所示。

http://s16/middle/4b013fb1498508de2ddcf&690

图4 参考电流模块Simulink测试图

http://s6/middle/4b013fb1n98505ca16cc5&690

图5 参考电流模块测试结果

    实例二：我觉得对于不做理论研究的我来说，如果用到S函数的话，一般的情况就如实例一的情况那样，没有离散和连续的状态，只需要对输入做某些运算或者判断就可。上学期学的一门交流调速的课程，问题稍微复杂一点，在此也说一下。

    对于含有齿隙环节的闭环系统，系统控制量反映了系统运动方向的变化，可以对系统控制量U加一个补偿量Ucomp以缩短系统在齿隙内的运行时间，补偿控制系统框图如图6所示。具体的补偿规则为：

if U(k)>0 and U(k-1)<=0 then flag=1
if U(k)<0 and U(k-1)>=0 then flag=1 Ucomp=flag*const*exp(-t/T0)

http://s7/middle/4b013fb1n98506a652686&690

图6 补偿控制系统框图

    分析系统情况：需离散状态变量，系统状态方程X1(k+1)=U(k)、X2(k+1)=f(X1(k)，U(k))，其中f代表一个如上等式所示的关系，输出Y= flag*const*exp(-t/T0)，因此得到如下的信息，输入为一路，即为图6中的u，输出为一路，连续状态个数为0，离散状态个数为2，输入端不出现在输出端，将S函数命名为comp，根据以上分析，其程序代码如下所示：

%Author:dingqian12345@126.com

%Note:实现基于规则的补偿算法

%2个离散状态、0个连续状态、1个输入信号、1个输出信号、输出不涉及输入信号

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

clc;

global Cflag;

global Cont;

global T0;

Cont=20;

T0=1.5;

switch flag

    case 0

        [sys,x0,str,ts]=mdlInitializeSizes; %初始化

    case 1

        sys=mdlDerivatives(t,x,u);          %连续状态的更新

    case 2

        sys=mdlUpdates(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

        error(['Unhandled flag=',num2str(flag)]);

end

%%%%%在flag=0的时候进行整个系统的初始化

function [sys,x0,str,ts]=mdlInitializeSizes

sizes = simsizes;            %读入初始化参数模板

sizes.NumContStates  = 0;    %连续状态个数

sizes.NumDiscStates  = 2;    %离散状态个数

sizes.NumOutputs     = 1;    %输出变量个数

sizes.NumInputs      = 1;    %输入信号个数

sizes.DirFeedthrough = 0;    %输入直接传入输出信号否

sizes.NumSampleTimes = 1;    % at least one sample time is neededWO  一般来说为1个

sys = simsizes(sizes);

x0=[0;0];        %状态初始化

str=[];

ts=[-1 0];    %采样周期若写成-1表示继承其输入信号

%%%%%%在flag=1的时候进行连续系统状态的更新

function sys=mdlDerivatives(t,x,u)

sys = [];

%%%%%%在flag=2的时候进行离散系统状态的更新

function sys=mdlUpdates(t,x,u)

global Cflag;

global Cont;

global T0;

Cflag=comp_decide(t,x,u);

sys(1,1)=u;

sys(2,1)=Cflag*Cont*exp(-t/T0);

%%%%%%在flag=3的时候进行系统输出信号的求取

function sys=mdlOutputs(t,x,u)

sys=x(2);

%%%%%%在flag=4的时候进行下一时刻仿真时间的计算

function sys=mdlGetTimeOfNextVarHit(t,x,u)

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

%sys = t + sampleTime;

%%%%%%在flag=9的时候终止仿真过程

function sys=mdlTerminate(t,x,u)

sys = [];

    可以看出，这里的S函数比实例1的S函数多用到了离散状态更新函数mdlUpdates，注意状态的更新方程的写法，sys=x(n+1)，可以通过help sfuntmpl查看。在以上代码中，在flag=2的时候调用了自编的子程序m=comp_decide(t,x,u)，comp_decide的代码如下：

%%%%%%用户自定义子程序

function cc=comp_decide(t,x,u)

cc=0;

if u>0 && x(1)<=0

    cc=1;

end;

if u<0 && x(1)>=0

    cc=-1;

end;

    最后，将S函数嵌入到Simulink中，所搭建的补偿仿真框图如图7所示，我们通过Scope4观察控制量U和补偿Ucomp，其中cont为20，To=1.5。图8为补偿系统中U和Ucomp仿真结果局部放大图，可见，系统可以捕捉到控制量过零点的时刻并给出相应的补偿量。

http://s16/middle/4b013fb1075a1aa8e211f&690

图7 补偿系统Simulink仿真框图

http://s14/middle/4b013fb149850ac79319d&690

图8 补偿系统中U和Ucomp仿真结果局部放大图

 

参考文献

[1] 薛定宇 著.控制系统仿真与计算机辅助设计[M].北京：机械工业出版社出版社，2009.

[2] 很好的S函数学习资料-Simulink的S-Function编写指导[EB/OL].

http://www.ilovematlab.cn/thread-71226-1-1.html.[2010-12-26].

[3] Simulink-s_function使用及应用举例[EB/OL].

http://www.ilovematlab.cn/thread-62299-1-1.html.[2010-12-26].

[4] 祁麟.交流伺服系统齿隙非线性控制及其网络化研究[D].硕士论文.南京理工大学，2006.

 

2010年12月26日晚上完成于njust 电工楼 310房间

 

CopyRight：版权所有若需转载或使用请联系作者
Email:dingqian12345@126.com