系统仿真

S函数控制器程序：controller.m
function [sys,x0,str,ts]=s_function(t,x,u,flag)
Байду номын сангаас
switch flag, %Initialization case 0, [sys,x0,str,ts]=mdlInitializeSizes; %Outputs case 3, sys=mdlOutputs(t,x,u); %Unhandled flags case {1, 2, 4, 9 } sys = []; %Unexpected flags otherwise error(['Unhandled flag = ',num2str(flag)]); end
在控制理论研究中，经常需要用复杂的算法设计 控制器，而这些算法经常因其复杂度又难以用模块搭 建。这样的系统如果需要在Simulink下进行仿真研究， 则需要用编程的形式设计出S-函数模块，将其嵌入到 系统中。成功使用M-函数和S-函数，则可以在 Simulink下对任意复杂的系统进行仿真。 介绍用MATLAB语言设计M-函数与S-函数的方法，
用MATLAB编写S-函数举例
以二阶线性传递函数为被控对象，进行模拟PID控制。
133 G( s) 2 s 25 s
在信号发生器中选择正弦信号，仿真时取kp=60，ki=1，kd=3， 输入信号为rin(t)=Asin(2πFt)，其中A=1.0，F=0.20Hz。
仿真方法一：
仿真方法二：
%mdlInitializeSizes function [sys,x0,str,ts]=mdlInitializeSizes sizes = simsizes; sizes.NumContStates = 0; sizes.NumDiscStates = 0; sizes.NumOutputs = 1; sizes.NumInputs = 3; sizes.DirFeedthrough = 1; sizes.NumSampleTimes = 0; sys=simsizes(sizes); x0=[]; str=[]; ts=[]; function sys=mdlOutputs(t,x,u) error=u(1); derror=u(2); errori=u(3); kp=60; ki=1; kd=3; ut=kp*error+kd*derror+ki*errori; sys(1)=ut;
Y=0; % Y为N×1阵，记录输出y，初始时为1×1阵，N为总步数 H=t; % H为N×1阵，记录时间t，初始时为1×1阵 while(t<tmax) xs=x+(A*x+B)*T; % 计算离散状态方程 y=xs(1); % 计算离散输出方程 t=t+T; Y=[Y;y];H=[H;t]; % 记录y和t的值，这时Y阵和T阵均增加1行 x=xs; end plot(H,Y); % 绘制输出曲线 grid; % 在“坐标纸”上画小方格
其中，fun2_3是保存在名为fun2_3的M文件中的M函数：
function xdot=fun2_3(t,x) xdot1(1)=-0.001*x(1)*x(2); xdot1(2)=0.001*x(1)*x(2)-0.072*x(2); xdot1(3)=0.072*x(2); xdot=xdot1'; % 第一个微分方程 % 第二个微分方程 % 第三个微分方程
•
S-函数的基本结构
用M语言编写的S-函数的引导语句为：
function [sys,x0,str,ts]=fun(t,x,u,flag,p1,p2,…)
其中： fun： S-函数的函数名； t,x,u：分别为时间、状态和输入信号； flag： 标志位。标志位取值不同，S-函数执行的任
务与返回数据是不同的。
并将通过例子介绍M-函数和S-函数的应用。
•
M-函数模块的基本结构
1 1
x y 1
1 x 1 x 1或x 1
用M-函数的形式描述该模块：(satur_non.m)
function y=satur_non(x) if abs(x)>=1,y=1*sign(x); else, y=x; end
或
function y = satur_ non (x) if x>=1 y=1; 选择User-Defined-Functions子 elseif x<= -1 库中的MATLAB Fcn模块, 并将 y= -1; 参数MATLAB Function设置为： else y=x; satur_non。 end
式中， x1表示可能受到传染的人数，x2 表示已经 被传染的病人数，x3表示已治愈的人数。试用ode45编 程，对其进行仿真研究，并绘制出对应的时间响应曲 线。
M文件：
clear x0=[620,10,70]; % 置状态变量初值 tspan=[0,30]; % 置仿真时间区间 [t,x]=ode45('fun2_3',tspan,x0); % 调用ode45求仿真解 plot(t,x(:,1),t,x(:,2),‘--’,t,x(:,3),‘:’); % 用不同的线型绘制仿真结果曲线 xlabel('time(天) t0=0, tf=30'); % 对t-x轴进行标注 ylabel('x(人):x1(0)=620,x2(0)=10;x3(0)=70'); legend('x1','x2','x3'); grid;
y[(n 1)T ] x1[(n 1)T ]
ⅳ 编程和调试
clear m=1;k=4; f=input('请输入阻尼系数f:'); t=0;T=0.01; A=[0 1;-k/m -f/m]; B=[0 1/m]'; tmax=10; x=[0,0] ';
% 质量系数m值，弹簧刚性系数m值 % 从键盘输入阻尼系数f值 % 置时间变量t和仿真步长T的初值 % 计算状态方程矩阵 % 置仿真总时间tmax的初值 % 置状态变量初值，其中x(i)代表xi(0)
133 G( s) 2 s 25 s
x Ax Bu y Cx Du
1 0 0 A , B , C 1 0, D 0 0 25 133
仿真方法三：
在仿真一的基础上，利用S函数实现对象的表达、控制器 的设计及仿真结果的输出。
2.控制系统仿真方法
编写M文件对系统进行仿真
【例】 考虑如右图所示的质量-弹 簧-阻尼器系统。当质量系数m=1， 弹簧弹性系数k=4时，为了使系统 的单位阶跃响应不发生振荡，阻尼 系数f（ 0 f 10 ）应在什么范围内 取值。
ⅰ问题的描述 要求研究的问题是：当m=1，k=4时，分析系统在外 力r(t)=1(t)的作用下，要使响应不发生振荡， f 应在什 么范围内取值，其中f的取值约束为
flag值为2时，将作离散状态变量的更新，将调用 mdlUpdate()函数。 flag值为3时，将求取系统的输出信号，将调用 mdlOutputs()函数。 flag值为4时，将调用mdlGetTimeOfNextVarHit()函 数，将计算下一步的仿真时刻。 flag值为9时，将终止仿真过程，将调用 mdlTerminate()函数，不返回任何变量。
采用M函数的形式，利用ODE45编程求解
[t , x]=ode45('方程函数名', tspan , x0, tol)
【例】 某地区某病菌传染的系统动力学模型为
x 0.001 x x 1 1 2 x2 0.001 x1 x2 0.072 x2 x3 0.072 x2 x1 (0) 620 x2 (0) 10 x3 (0) 70
0 f 10
ⅱ 建立系统的数学模型 描述该系统输入—输出关系的数学模型为
m fy ky r y
1 f m
0 1 x k x 2 m
x1 0 1 r x2 m
利用Simulink对系统进行仿真 首先将描述该系统输入—输出关系的数学模型改写成 传递函数形式 Y ( s) 1
G( s) R( s ) ms 2 fs k
从 Simulink模块库的信号源子库（Sources）中，选择 Step模块，拖入模型窗口，双击该模块，设置参数Step time为：0；从Simulink模块库的连续模块子库 （Continuous）中，选择Tranfer Fcn模块，拖入模型窗口, 双击该模块，设置参数Numerator为：[1]， Denominator 为：[1 f 4] ；从 Simulink模块库的接收器子库（Sinks） 中，选择Scope模块，拖入模型窗口。
flag值为0时，将启动S函数所描述系统的初始化 过程，这时将调用一个名为mdlInitializeSizes() 函数。 sizes=simsizes;获得默认的系统参数变量sizes sizes是结构体变量包括：
NumContStates；NumDiscStates； NumInputs； NumOutputs； DirFeedthrough；NumSampleTimes；
ⅲ 数学模型转换成仿真模型
x1 y [1 0] x2
1 f m 0 x1 (nT ) 1 r (nT )T x (nT ) 2 m
0 x1 (( n 1)T ) x1 (nT ) x (n 1)T ) x (nT ) k 2 2 m
S函数被控对象程序：plant.m
%mdlInitializeSizes function [sys,x0,str,ts]=mdlInitializeSizes %S-function for continuous state equation sizes = simsizes; function [sys,x0,str,ts]=s_function(t,x,u,flag) sizes.NumContStates = 2; sizes.NumDiscStates = 0; switch flag, sizes.NumOutputs = 1; %Initialization sizes.NumInputs = 1; case 0, sizes.DirFeedthrough = 0; [sys,x0,str,ts]=mdlInitializeSizes; sizes.NumSampleTimes = 0; case 1, sys=simsizes(sizes); sys=mdlDerivatives(t,x,u); x0=[0,0]; %Outputs str=[]; case 3, ts=[]; sys=mdlOutputs(t,x,u); %Unhandled flags function sys=mdlDerivatives(t,x,u) case {2, 4, 9 } sys = []; sys(1)=x(2); %Unexpected flags %sys(2)=-(25+5*sin(t))*x(2)+(133+10*sin(t))*u; otherwise sys(2)=-(25+10*rands(1))*x(2)+(133+30*rands(1))*u; error(['Unhandled flag = ',num2str(flag)]); function sys=mdlOutputs(t,x,u) end
1.2 Simulink仿真环境
Simulink是MATLAB下的数字仿真工具，其文件类型 为.mdl，是一个用来对动态系统进行建模、仿真和分 析的集成环境。 Simulink提供了用鼠标“画出”系统框图的方式，可 以进行图形建模，具有直观、方便、灵活的优点。 在MATLAB的命令窗运行指令simulink或单击MATLAB 操作桌面中的Simulink图标，便打开Simulink模块库浏 览器（Simulink Library Browser）窗口。
1.MATLAB/Simulink仿真环境
1.1 MATLAB语言简介
MATLAB是一种解释性程序设计语言，用户既可以在 MATLAB环境下键入一个命令，也可以用MATLAB语 言编写应用程序。
在Windows环境下，启动MATLAB后，就打开了一个 MATLAB环境的窗口（如图），在该Command Window窗内，可以键入各种送给MATLAB运行的指令、 函数、表达式；显示除窗口外的所有运行结果。
sys=simsizes(sizes); sys除外，还应设置系统的初始状态变量x0，说明变量 str和采样周期变量ts。 ts为双列的矩阵，每一行对应一个采样周期。对于连续
系统和有单个采样周期的系统来说，[t1,t2]。
flag值为1时，将作连续状态变量的更新，将调用 mdlDerivatives()函数。
利用M-函数、S-函数对系统仿真 在实际仿真中，如果模型中某个部分数学运算特 别复杂，则不适合用普通Simulink模块来搭建这样的 部分，而应该采用程序来实现。Simulink中支持两种 用语言编程的形式来描述这样的模块，即M-函数和S函数，它们的用途是不同的。
• M-函数适合于描述输出和输入之间为代数运算的模块。 • S-函数适合于动态关系的描述，所谓动态关系即是由 状态方程描述的关系。（系统函数）