【沧海拾昧】MATLAB\Simulink用S-函数仿真Lorenz系统

#C0403

---

沧海茫茫千钟粟，且拾吾昧一微尘

——《沧海拾昧集》@CuPhoenix

---

【阅前敬告】 沧海拾昧集仅做个人学习笔记之用，所述内容不专业不严谨不成体系 【如有问题必是本集记录有谬，切勿深究】

---

---

前言

    S-函数 是系统函数（System Function）的简称，是在 MATLAB 中仿真控制工程相关问题的常用方法。现如今，使用 S-函数 对复杂动态系统进行仿真已经是较通用、简便的做法，在 MATLAB 中也有成型的 S-函数 模板可供使用。

    本文将简要介绍在 MATLAB\Simulink 中使用 S-函数 进行仿真的基础内容和部分常见问题，并以 Lorenz 混沌系统为例，对仿真步骤进行基本介绍。

本文用到的软件环境是：

• MATLAB R2015b

---

一、S-函数

 
     S-函数 是采用非图形化的方式描述的一个功能块（非图形化方式，即使用计算机语言而非像 Simulink 一样拖动系统模块进行组装），用以描述并实现连续系统、离散系统、复合系统等动态系统。可以用 MATLAB、C\C++、Fortran等语言进行编写。

     S-函数 相比于 Simulink 的图形化仿真其优点在于：

• 便捷简单，对于复杂的系统可以用文本的形式更方便的写出来（并且有成型的模板可以方便套用）。
• 移植性好，能够如同一个模块一样在 Simulink 中使用。
• 功能丰富，对 Simulink 的功能进行扩展（如M文件S-函数能扩展图形能力，C语言S-函数能提供与操作系统的接口）

    在 MATLAB 命令行窗口中输入命令打开 MATLAB中 S-函数 的通用模板。

edit sfuntmpl

    这里给出一个简化注释后的纯净版模板便于理解与复制。

% 主函数，至少 4 个输出。其中：sys数组包含子函数的返回值，x0为状态初始化向量，str是空矩阵(保留参数)，ts为系统采样时间。
% 注意：编写 S-函数 时必须修改 sfuntmpl 为自定义函数名，且与M文件的文件名保持一致。
function [sys,x0,str,ts,simStateCompliance] = sfuntmpl(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

% 初始化例程子函数(FLAG=0)
function [sys,x0,str,ts,simStateCompliance]=mdlInitializeSizes
sizes = simsizes;
sizes.NumContStates  = 0;	% 连续状态的个数
sizes.NumDiscStates  = 0;	% 离散状态的个数
sizes.NumOutputs     = 0;	% 输出的个数
sizes.NumInputs      = 0;	% 输入的个数
sizes.DirFeedthrough = 1;	% 是否直接馈通（bool型变量，缺省则默认为1-存在）
sizes.NumSampleTimes = 1;   % 采样时间的个数(至少为1)
sys = simsizes(sizes);
x0  = [];	% 初始化状态
str = [];	% 保留变量，不得赋值
ts  = [0 0];	% 采样时间，[采样周期值 偏移量]，是 m*2 的矩阵
simStateCompliance = 'UnknownSimState';

% 计算导数子函数（FLAG=1）
function sys=mdlDerivatives(t,x,u)
sys = [];

% 状态更新子函数（FLAG=2）
function sys=mdlUpdate(t,x,u)
sys = [];

% 计算输出子函数（FLAG=3）
function sys=mdlOutputs(t,x,u)
sys = [];

% 计算下次一个采样时间子函数（FLAG=4）
% 注：仅在系统是变采样时间系统时调用
function sys=mdlGetTimeOfNextVarHit(t,x,u)
sampleTime = 1;
sys = t + sampleTime;

% 结束仿真子函数（FLAG=9）
function sys=mdlTerminate(t,x,u)
sys = [];

模板说明：

• 该模板是一个标准的 M文件S-函数 ，由 1 个主函数和 6 个子函数组成。
• 主函数使用一个 Switch-Case 结构，根据 Flag 标记的值将流程转移到对应的子函数。
• 状态向量的一阶导数是状态 $x$ 、输入 $u$ 、时间 $t$ 的函数。

M文件S-函数 流程

 

二、实例：Lorenz 混沌系统仿真

     本例以在 MATLAB/Simulink 中仿真某 Lorenz 混沌系统为例，对流程及步骤进行简要介绍。该系统的微分方程为：
$$\begin{gathered} {\dot{x}}_1(t)=-a{x}_1(t)+x_2(t)x_3(t) \\ {\dot{x}}_2(t)=-b{x}_2(t)+b{x}_3(t) \\ {\dot{x}}_3(t)=-x_1(t)x_2(t)+c{x}_2(t)-x_3(t) \end{gathered}$$其中，$a=8/3、b=10,c=28$ 是参数取值。

1、自定义 S-函数

     新建一个 M文件S-函数 文件（如：MySFunction1.m，注意修改模板时主函数名必须与文件名相同），并将模板粘贴到文件。然后对以下内容进行修改：

• 初始化子函数：修改输出输出个数，设置状态初始条件。

  % 初始化例程子函数(FLAG=0)
  function [sys,x0,str,ts,simStateCompliance]=mdlInitializeSizes
  % (未改动部分略)
  sizes.NumContStates  = 3;	% Lorenz 函数输出3个连续状态
  sizes.NumDiscStates  = 0;	
  sizes.NumOutputs     = 3;	% 3个状态均输出
  sizes.NumInputs      = 0;	% Lorenz 函数自激振荡，不需要输入
  % (未改动部分略)
  x0  = [1 1 1];	% 设置初始化状态
  % (未改动部分略)
  

• 计算导数子函数：定义系统微分方程（多个方程用分号隔开）。

  % 计算导数子函数（FLAG=1）
  function sys=mdlDerivatives(t,x,u)
  sys = [-8/3*x(1)+x(2)*x(3);
      -10*x(2)+10*x(3);
      -x(1)*x(2)+28*x(2)-x(3)];
  

• 计算输出子函数：定义系统输出，本例中输出状态 $x$ 。

  % 计算输出子函数（FLAG=3）
  function sys=mdlOutputs(t,x,u)
  sys = x;
  

将修改后的 MySFunction1.m 文件保存待用。

2、在 Simulink 中搭建仿真

注：关于 Simulink 的基础使用可参阅 MATLAB/Simulink仿真的基本操作 。

• 路径设置
       新建 Simulink Model 文件，保存在和 MySFunction1.m 相同的路径下，并在左侧当前文件夹中打开该路径。

路径设置

• 仿真结构搭建
       新建 Simulink Model 文件，保存在和自定义 S-函数 相同的路径下，并在左侧当前文件夹中打开该路径。在 Simulink 中使用 S-function 模块（位于 Library-Simulink-UserDefined Functions），按下图搭建仿真结构。
  

Lorenz系统的S-函数仿真

• S-函数设置
       双击 S-function 模块，将函数名调整为自定义S-函数的函数名。
  

修改函数名为 MySFunction1

• 设置求解时间与步长：求解时间设为 $100s$ ，步长设为 $auto$ 。

3、仿真结果

• 仿真结果查看

     在示波器中可以查看仿真结果，图像为 $x(1)\sim x(3)$ 在 $0\sim 100s$ 范围内的值。

Lorenz系统在时域的响应图像

• 工作空间绘图

     利用 ToWorkSpace 模块向工作空间记录数据后。可以利用如下代码绘图。（简单的绘图格式可参阅 用MATLAB画一张简单的图 ）

% 绘制Lorenz系统的三维相空间曲线
plot3(simout.Data(:,1),simout.Data(:,2),simout.Data(:,3))

% 调整其他绘图格式
grid on
title('\fontname{Times}Lorenz\fontname{宋体}系统的三维相空间曲线')
xlabel('x')
ylabel('y')
zlabel('z')
set(gca,'FontName','Times New Roman','FontSize',12,'FontWeight','bold') 

Lorenz系统的三维相空间曲线

 

三、其他资料

可参考的一些资料：

• 赵海滨;于清文;刘冲;陆志国;颜世玉.基于Matlab/Simulink的混沌同步控制实验[J].实验室研究与探索,2019,v.38;No.275,22-25.
       本文提供了一些其他的基于 MATLAB\Simulink 的仿真方法，如使用 ode 45() 函数、运用 Simulink 基本模块进行建模、M-函数模块等。

• 用一个简单的例子说明如何使用S函数进行Simulink建模 @鹿寨1314
       本文举例了一个更简单的 S-函数 仿真实例，并用推导计算的方法进行了一些说明。

敬谢诸君。

---

金陵钟山之阳。