1.被控模型和参考模型

顾名思义，首先我们要有参考模型，以二阶为例，参考模型如下：

被控二阶系统：

其中，d(t)为外部扰动，且|d(t)| ≤ D。

2.控制律推导

①定义跟踪误差e = y-ym，设计滑模面/滑模函数为：s = e'+c*e，则s' = e''+c*e';

②设滑模趋近律为s';

③将系统参数代入滑模导函数中，并联立滑模函数导数和趋近律，即

s' = e''+c*e' = (y''-ym'')+c*e'=（a*y'+b*u+d-(am*ym'+bm*r)）+c*e'；

所以：u = [s'-c*e'-d+(am*ym'+bm*r)-a*y']/b，若趋近律取为等速趋近律，即s' = -η*sign(s)，则

u = [-η*sign(s)-c*e'-d+(am*ym'+bm*r)-a*y']/b = [-η-c*|e'|-d+|am*ym'|+|bm*r|-|a*y'|]*sign(s)/b，因为

d ≤ D，不妨取d = D，则u = [-η-c*|e'|-D+|am*ym'|+|bm*r|-|a*y'|]*sign(s)/b

怎么推导和书上给的控制律不一样？加号那里符号有差别？

④使用Lyapunov函数证明闭环系统稳定性：

取Lyapunov函数V = 1/2*s^2

3.仿真分析

3.1s函数编写被控对象

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

switch flag,
case 0,
    [sys,x0,str,ts]=mdlInitializeSizes;
case 1,
    sys=mdlDerivatives(t,x,u);
case 3,
    sys=mdlOutputs(t,x,u);
case {2,4,9}
    sys=[];
otherwise
    error(['Unhandled flag = ',num2str(flag)]);
end

function [sys,x0,str,ts]=mdlInitializeSizes
sizes = simsizes;
sizes.NumContStates  = 2;
sizes.NumDiscStates  = 0;
sizes.NumOutputs     = 2;
sizes.NumInputs      = 1;
sizes.DirFeedthrough = 0;
sizes.NumSampleTimes = 1;
sys = simsizes(sizes);
x0 = [1.5;0];
str = [];
ts = [0 0];

function sys=mdlDerivatives(t,x,u)
a=25;
b=133;

sys(1)=x(2);
sys(2)=-a*x(2)+b*u+10*sin(t);
function sys=mdlOutputs(t,x,u)
sys(1)=x(1);
sys(2)=x(2);

3.2s函数编写控制器

function [sys,x0,str,ts] = spacemodel(t,x,u,flag)
switch flag,
case 0,
    [sys,x0,str,ts]=mdlInitializeSizes;
case 3,
    sys=mdlOutputs(t,x,u);
case {2,4,9}
    sys=[];
otherwise
    error(['Unhandled flag = ',num2str(flag)]);
end
function [sys,x0,str,ts]=mdlInitializeSizes
sizes = simsizes;
sizes.NumContStates  = 0;
sizes.NumDiscStates  = 0;
sizes.NumOutputs     = 1;
sizes.NumInputs      = 4;
sizes.DirFeedthrough = 1;
sizes.NumSampleTimes = 1;
sys = simsizes(sizes);
x0  = [];
str = [];
ts  = [0 0];
function sys=mdlOutputs(t,x,u)
a=25;b=133;
am=20;bm=100;
D=10;
c=10;
ym=u(1);y=u(3);
dym=u(2);dy=u(4);

e=y-ym;
de=dy-dym;
s=c*e+de;

r=sin(pi*t);
xite=0.02;

wt=1/b*(-c*abs(de)-abs(bm*r)-D-xite-abs(am*dym)-abs(a*dy));

M=2;
if M==1
    ut=wt*sign(s);
elseif M==2
    delta=0.02;
    if s>delta
        sats=1;
    elseif abs(s)<=delta
        sats=s/delta;
    elseif s<-delta
        sats=-1;
    end
    ut=wt*sats;
end
% ut = 1/b*(am*dym-bm*r-c*de-a*dy-(D+xite)*sign(s));%为什么要加绝对值？
sys(1)=ut;