简介：面向控制工程与自动化领域学习者，这份压缩包围绕滑模变结构控制的核心理论与MATLAB实现，汇集了从滑模面设计、切换函数选择到控制律搭建、Simulink仿真分析的完整示例程序，尤其适合正在研读《滑模变结构控制MATLAB仿真（第3版）》并希望动手验证的读者。资源共357个文件，大小约979KB，以259个.m脚本为主，配合68个.mdl模型、20个r2011a版本模型文件以及.mat数据文件、.fis模糊推理文件，覆盖多种被控对象与控制器设计场景，便于按章节对照学习与二次开发。已有1234人浏览学习。通过运行这些源码，读者可以直观观察系统响应曲线，理解抖振抑制与参数摄动下的鲁棒性表现，并借助现有框架调节滑模面参数或切换增益，提升自身控制器设计能力与MATLAB/Simulink仿真水平。
1. 滑模变结构控制MATLAB仿真在做什么——从抖振开始的工程叙事
一个连续系统在数值仿真里被符号函数来回切换，相轨迹会在滑模面附近来回穿越，形成锯齿状的高频抖动——这是滑模变结构控制最精妙也最折磨人的地方。MATLAB仿真程序要解决的不是“能不能跑”，而是“跑得对不对”：切换逻辑在离散步长下是否真实反映连续系统的行为，趋近律参数取多少才能让系统既快速收敛又不至于振荡发散。这篇内容面向两类人：一类是写课程设计或研一课题的学生，需要把滑模控制从公式变成可运行的脚本；另一类是做电机控制、机械臂或飞行器姿态的工程师，想在Simulink里搭一个能挂在被控对象上的控制器原型，并搞清楚抖振到底从哪来、用哪个参数去压。
2. 先立住滑模控制的数学骨架：状态方程、切换面与等效控制
2.1 二阶非线性系统怎么改写成MATLAB能处理的状态空间
滑模变结构控制面向的是仿射非线性系统，形如
x_dot = f(x) + g(x) · u
在MATLAB仿真里，第一步不是写控制器，而是把被控对象整理成“输入u到状态x_dot”的一阶微分方程组。以单连杆机械臂绕关节转动为例，动力学方程是
J · theta_ddot + b · theta_dot + m·g·l·sin(theta) = tau
取状态x1 = theta，x2 = theta_dot，控制量u = tau，则状态方程为：
x1_dot = x2
x2_dot = (u - b·x2 - m·g·l·sin(x1)) / J
这两行必须写在MATLAB函数句柄或单独的函数文件里，作为所有仿真脚本的共同底座。
% plant_dynamics.m % 单连杆机械臂的动力学模型，给仿真脚本和被控对象函数共用 function dx = plant_dynamics(x, u, p) q = x(1); dq = x(2); dx = zeros(2,1); dx(1) = dq; dx(2) = (u - p.b*dq - p.m*p.g*p.l*sin(q)) / p.J; end这段代码里，p是结构体参数，集中存放J、b、m、g、l，便于修改。x是当前的二维状态向量，u当前时刻的控制力矩。把动力学隔离成独立函数的好处是：后面做控制器仿真、S函数、参数扫描时不用反复改方程。
2.1.1 参数怎么给：别用魔法数字
直接在脚本里写dx(2) = (u - 0.1*dq - 1.2*9.8*0.3*sin(q))/0.5看起来很省事，但等你调趋近律参数时就会后悔。我一般把物理参数放在脚本顶部的结构体里：
p.J = 0.5; % 关节转动惯量 kg·m^2 p.b = 0.1; % 粘性摩擦系数 N·m·s/rad p.m = 1.2; % 连杆质量 kg p.g = 9.8; % 重力加速度 m/s^2 p.l = 0.3; % 质心到关节距离 m参数集中管理之后，替换模型对象（比如把机械臂换成永磁同步电机的d-q轴方程）只需要改p和plant_dynamics函数，控制器代码可以完全不碰。
2.2 切换函数与控制律推导：先求等效控制，再加切换项
滑模控制器的设计分两步。第一步，设计切换函数s(x) = c·x1 + x2，这里c > 0 是滑模面系数。第二步，令滑模面不变性条件成立：s_dot = 0，反解出的控制量称为等效控制u_eq。
对上述机械臂模型：
s = c·x1 + x_dot1
s_dot = c·x2 + x2_dot = c·x2 + (u - b·x2 - m·g·l·sin(x1))/J
令 s_dot = 0，得
u_eq = -J·c·x2 + b·x2 + m·g·l·sin(x1)
这体现了等效控制的物理含义：它刚好抵消系统的非线性项和阻尼项，使状态在滑模面上滑动。实际控制律还需要加一个切换项，保证系统从滑模面外被“拉”回滑模面，也就是满足可达条件 s·s_dot < 0。
完整控制律是 u = u_eq + u_sw，其中 u_sw 通常取 -K·sign(s) 或带趋近律的形式。sign(s)在MATLAB里就是符号函数，但仿真里直接用sign会带来数值振荡，后面第4章会专门展开。
2.2.1 为什么仿真程序里要区分“模型函数”和“控制器函数”
很多初学者把被控对象和控制器写进同一个for循环，这样跑几次是没问题，但一旦要对比不同趋近律、不同滑模面的效果，就得复制粘贴整段代码。更好的做法是让控制律也成为一个独立的MATLAB函数，输入状态和参数，输出控制量：
function u = smc_controller(x, p, c, K, s_type) q = x(1); dq = x(2); % 等效控制：抵消模型中的标称项 u_eq = -p.J*c*dq + p.b*dq + p.m*p.g*p.l*sin(q); % 切换项 s = c*q + dq; switch s_type case 'sign' u_sw = -K * sign(s); case 'sat' Delta = 0.02; % 边界层厚度，见第6章 u_sw = -K * sat(s, Delta); case 'tanh' eps_sw = 0.2; % 光滑化系数 u_sw = -K * tanh(s/eps_sw); end u = u_eq + u_sw; end function val = sat(s, Delta) % 边界层饱和函数 if abs(s) < Delta val = s/Delta; else val = sign(s); end end注意两点：第一，u_eq里的-c·x2项来自滑模面不变性条件，不能省；第二，切换项u_sw的增益K必须大于系统不确定性和扰动的上界，否则到达条件不满足。后面调参时，K数值偏小时状态到不了滑模面，直接表现为跟踪误差不收敛。
2.3 仿真的整体架构：ode45还是for循环，这是个问题
滑模控制仿真有两种主流的数值求解方式。常见的是用MATLAB的ode45配合事件函数，让步长自适应地跟随非光滑的切换行为；我自己在工程上更推荐用固定步长显式欧拉或ode4跑，原因后面解释。
先给出用ode45的标准脚本结构：
% run_smc_ode45.m % 外层用ode45积分，控制器在odefun内部被调用 p.J = 0.5; p.b = 0.1; p.m = 1.2; p.g = 9.8; p.l = 0.3; c = 8; K = 5; x0 = [0.5; 0]; % 初始角度0.5rad，初始角速度0 tspan = [0 5]; opts = odeset('RelTol',1e-6,'AbsTol',1e-8); [t, x] = ode45(@(t,x) closed_loop(t,x,p,c,K), tspan, x0, opts); function dx = closed_loop(t, x, p, c, K) u = smc_controller(x, p, c, K, 'sign'); dx = plant_dynamics(x, u, p); end这段代码的逻辑是：每一时刻，控制函数根据当前x计算u，再把u和x一起交回plant_dynamics求导数，ode45根据该导数推进积分步。自适应步长下，symbolic的switch行为没有本质区别，但原因在于——当状态反复穿越s=0时，ode45会被forces去细分步长，导致仿真时间成倍变长。
2.3.1 为什么固定步长更适合滑模仿真
控制系统的仿真要尊重一个事实：实际的控制器是按固定周期运行的，比如1kHz或500Hz，中间的状态演化是零阶保持的。用变步长ode45反而违背了这个物理时序。所以更贴近工程的做法是用ode4（固定步长四阶Runge-Kutta）跑：
dt = 0.001; % 控制周期1ms t = 0:dt:5; x = zeros(2,length(t)); x(:,1) = [0.5; 0]; for k = 1:length(t)-1 xk = x(:,k); u_k = smc_controller(xk, p, c, K, 'sign'); k1 = plant_dynamics(xk, u_k, p); xh = xk + dt*k1; u_h = smc_controller(xh, p, c, K, 'sign'); k2 = plant_dynamics(xh, u_h, p); % 二阶中点法简化示意 x(:,k+1) = xk + dt*k2; end关键参数是dt：取0.001秒是底线，离散控制系统仿真经验法则是dt要小于控制系统带宽对应时间常数的五分之一。滑模控制的切换频率非常高，dt太大会让抖振幅度被数值放大，仿真结果失真。
3. 切换面设计：c参数是滑模面的灵魂
3.1 c与系统特征根的关系
切换函数s = c·x1 + x2 = 0定义了相平面中的一条直线——滑模面。当s=0时，有x2 = -c·x1，即
x1_dot = x2 = -c·x1
这是一个一阶线性方程，解为 x1(t) = x1(0)·e^{-c·t}。由此可见，c越大，滑模面上的状态收敛越快，但代价是滑模面的高频增益和切换项的震荡幅度都会增大。c本质上是滑模面上的一个极点，和线性系统极点配置里的期望闭环极点直接对应。
在MATLAB中验证这个关系很简单：
% 观察c对滑模面上收敛速度的影响 c = [2, 5, 10]; t = 0:0.01:3; x0 = 1; for i = 1:3 plot(t, x0*exp(-c(i)*t), 'LineWidth', 2); hold on; end grid on; legend('c=2', 'c=5', 'c=10');注意，直接调节c时系统在滑模面上的运动还受被控对象本身动态的影响，机械臂的电机延迟、结构谐振如果被c激励起来，仿真出来的高频段噪声会显著增大。
3.2 从极点配置角度确定c的初值
更规范的做法是把滑模面系数当作线性系统的极点来设计。对二阶系统，期望滑模面多项式是 λ + c = 0，即期望的一个闭环极点为 -c。对于更高阶系统（如机械臂的关节空间模型是二阶，但带电机时会是四阶），用place函数解：
% 以电机+关节的二阶线性化模型为例 A = [0 1; 0 -p.b/p.J]; B = [0; 1/p.J]; % 期望滑模面对应的极点 desired_pole = -5; c_design = place(A, B, desired_pole); % 对二阶系统，c_design是1x2向量 % 取滑模面系数为第一个元素（x2的系数归一化后）这里说一下place的坑：place对系统可控性极为敏感，如果A是奇异阵，place可能报错。在滑模控制里其实没必要用place——手动给定c-(-c是特征根)更直接。我通常在0.5到15之间扫几个值，观察跟踪误差的收敛时间。
3.2.1 c参数调参速查表
| 被控对象 | 推荐c区间 | 依据 |
|---|---|---|
| 单连杆机械臂（J≈0.5） | 3～10 | 大于机械带宽3倍以上，小于驱动饱和对应值 |
| 永磁同步电机速度环 | 8～20 | 电流环带宽约100Hz，速度滑模面要在电流环响应时间内 |
| 飞行器姿态角回路 | 2～6 | 执行机构带宽低，c过大会与舵机延迟耦合 |
| 惯性较大的双连杆臂 | 1～4 | 结构频率低，c太大会激起柔性模态 |
注意，c取值必须低于控制频率对应奈奎斯特频率的十分之一，否则直接发散。
3.3 调c时的典型误判：把跟踪误差大归因于K不足
一个常见的误操作：系统状态始终不收敛，初学者会习惯性增大切换增益K。但滑模控制的误差从两个源头来：一是到达阶段误差，而这由趋近律参数决定；二是滑模面收敛误差，由c决定。如果c太小，即便K足够大，系统到达滑模面之后滑动得很慢，误差曲线会呈现“缓坡”而非快速归零。
我判断c是否合适的方法是看仿真输出的相轨迹图：
figure; plot(x(1,:), x(2,:)); xlabel('角度(theta)'); ylabel('角速度(dtheta)');如果相轨迹先是迅速冲向相平面直线s=0，然后在s=0附近小幅度振荡滑向原点，说明c设置基本合适。如果相轨迹绕圈不贴近s=0，那就是c过小。
4. 趋近律与抖振源头：参数ε、K的仿真标定
4.1 符号函数在数值仿真里的三种实现方式
s_dot = -K·sign(s)是最基本的等速趋近律，sign(s)在连续域中是一个理想开关，但在数值积分中会导致状态在s=0附近来回穿越。MATLAB里实现切换项有几种选择，每种都会影响仿真结果：
| 实现方式 | 表达式 | 数值特性 | 适用场景 |
|---|---|---|---|
| sign | sign(s) | 切换最尖锐，抖振最明显 | 理论演示、验证滑模存在性 |
| sat | s/Δ ( | s | <Δ); sign(s) (其他) |
| tanh | tanh(s/ε) | 全程光滑，但控制增益被软化 | 需要连续控制量的场景 |
仅用sign会导致两个数值问题：第一，状态在滑模面两侧以最高频率来回切换，ode45会不断缩小步长，仿真时间可能膨胀数倍；第二，真正的物理系统是不存在无限频率切换的，仿真会高估实际系统的高频分量。
4.1.1 工程上怎么在仿真里处理sign
工程定位是通过给连续模型做一个执行机构带宽限制来反映现实，比如加一个一阶惯性环节1/(tau_act*s+1)。但如果只是验证控制器本身，用sat就足够了。
这里把sat函数定义成独立函数，方便在多个脚本中复用：
function val = sat(s, Delta) % 饱和函数：边界层内线性，边界层外饱和 if abs(s) <= Delta val = s / Delta; else val = sign(s); end end注意sat中Delta的量纲与s一致，s的量纲取决于状态变量的单位。对机械臂角度，s单位是rad/s，所以Delta取0.01～0.05。
4.2 指数趋近律的仿真脚本与参数标定
等速趋近律的问题是靠近滑模面时s_dot仍很大，会产生过冲。工程上更常用指数趋近律：
s_dot = -ε·sign(s) - kq·s
其中ε>0保证到达条件，kq>0系数指数收敛项。把趋近律代入s_dot的展开式，可解出控制律：
function u = smc_controller_exp(x, p, c, epsilon, kq, Delta) q = x(1); dq = x(2); s = c*q + dq; % 等效控制：来自 s_dot = 0 u_eq = -p.J*c*dq + p.b*dq + p.m*p.g*p.l*sin(q); % 指数趋近律的切换项 u_sw = -p.J*(epsilon*sat(s,Delta) + kq*s); u = u_eq + u_sw; end这里的epsilon取0.5～2视为合理的起始区间，kq取10到50。epsilon的作用是在状态到达滑模面附近时提供恒定推力保证穿越，kq的作用是在远离滑模面时加快趋近速度。显然，如果kq过大，系统在到达段会非常急促，对执行机构造成冲击。
跑一组不同参数得到时间响应对比图：
% 参数对比脚本 epsilon_list = [0.5, 1.0, 2.0]; kq_list = [10, 20, 50]; colors = {'r','b','k'}; for i = 1:3 x = simulate_closed_loop(p, c, epsilon_list(i), kq_list(i)); plot(t, x(1,:), colors{i}, 'LineWidth', 1.5); hold on; end仿真结果通常会说明一个规律：epsilon对抖振幅值的影响远大于对收敛时间的影响，而kq影响的是到达时间但几乎不改变滑模面上的动态。记住这一点，调参时就能少做很多盲试。
4.3 仿真步长与控制周期对参数的影响
这一步特别容易被忽略。实测经验：把仿真步长从1ms改成0.1ms，epsilon和K的合适参数区间会明显下移。原因是步长越细，数值积分对高频切换的解析越真实，系统越容易识别出真实的抖振频率；步长太粗时，抖振被平均化了，epsilon可以取更大而不发散。
由此建议的标定顺序是：
- 先固定dt为期望控制周期（如1ms）。
- 用更大的epsilon观察抖振幅值是否在可接受范围。
- 再调kq，观察到达滑模面的时间是否满足指标。
- 所有参数调完后，将dt降到0.2ms重新仿真，如果抖振明显变大，适当减小epsilon。
这也是为什么仿真参数不能直接搬到实际硬件上跑的深层原因：数字仿真的步长和控制器的采样时间是两个独立变量，混淆它们会导致设计失察。
5. 用MATLAB S函数把滑模控制器做进Simulink
5.1 S函数的结构与滑模控制器的S函数实现
纯脚本仿真适合调参和验证，但如果你最终要在Simulink里搭整个机电系统（伺服驱动、负载模型、观测器），把控制器写成一个S函数是最干净的方式。S函数本质上是一个MEX或MATLAB函数，其回调函数被Simulink执行。
一个完整的两输入三输出S函数框架，输入是状态x，输出是控制量u：
function [sys, x0, str, ts] = smc_sfun(t, x, u, flag, c, K, epsilon, Delta) switch flag case 0 sizes = simsizes; sizes.NumContStates = 0; % 无连续状态 sizes.NumDiscStates = 0; % 无离散状态 sizes.NumOutputs = 1; sizes.NumInputs = 2; % 角度和角速度 sizes.DirFeedthrough = 1; % 输出直接依赖输入，必须置1 sizes.NumSampleTimes = 1; sys = simsizes(sizes); x0 = []; str = []; ts = [-1 0]; % 继承被控对象采样时间 case 3 q = u(1); dq = u(2); s = c*q + dq; % 被控对象参数写死为便于演示，工程上应作为参数传入 J = 0.5; b = 0.1; m = 1.2; g = 9.8; l = 0.3; u_eq = -J*c*dq + b*dq + m*g*l*sin(q); u_sw = -J*(epsilon*sat(s,Delta) + K*s); sys = u_eq + u_sw; case 4 sys = []; case 9 sys = []; otherwise sys = []; end end这里的DirFeedthrough必定是1，因为控制律里直接含有当前输入u（这里是状态x）。如果把状态直接接到S函数输入，或者S函数里不需要当前输入，才可以把该位置为0。ts = [-1 0] 表示继承前一个块的采样时间，对S函数接入连续被控对象来说通常没问题。
5.2 Simulink模型搭建与参数设置
在Simulink里搭这个闭环需要的块有：被控对象模型（State-Space或自定义的连续传递函数）、S-Function块、Mux、Scope。
搭建步骤：
- 拖入一个S-Function块，双击后填写S-Function名称：smc_sfun。参数列表依次输入c、K、epsilon、Delta，用英文逗号分隔。
- 拖入Integrator两个，用一个Sum和Gain搭出机械臂的连续模型，或者直接用Continuous库里的State-Space，状态矩阵A和B从2.1节的plant_dynamics线性化得到。
- S-Function的输入接Mux，Mux的输入接被控对象的状态；S-Function输出接被控对象的输入端。
- 用Scope观察角度theta和控制量u。
S参数块设置的关键项表格：
| 参数 | 设置值 | 说明 |
|---|---|---|
| S-Function name | smc_sfun | 必须与.m文件名一致，文件名不含中文 |
| S-Function parameters | 8, 5, 1.0, 0.02 | 与S函数输入参数顺序一一对应 |
| Allow direct feedthrough | 由编码自动设为1，不勾选会报错 | |
| 采样时间 | -1（继承） | 若设为0表示连续，设其它则为离散控制周期 |
如果出现“Error evaluating parameter”类型错误，多半是S-function里的sizes结构少了字段，例如忘了NumSampleTimes，或ts设置出了问题。熟悉MATLAB Simulink仿真的可以查一下doc simsizes的字段说明。
5.3 用Simulink与纯脚本对比验证控制器实现一致性
S函数写完之后，务必做一次脚本/S函数的对比验证。方法是：在Simulink里跑出总时间为1s的仿真，把数据导到工作区，同时用固定步长跑一次纯脚本版本，然后画在一起看是否重合。
% 将simulink输出加载到工作区 (假设输出名为theta_sim) theta_sim = yout(:, 1); % 纯脚本运行 [t_script, x_script] = simulate_script(p, c, K, epsilon, Delta); plot(t_script, x_script(1,:)); hold on; plot(t_sim, theta_sim, 'r--'); legend('脚本', 'Simulink S函数');两条曲线如果不重合，优先检查三点：
- S函数里是否有隐藏的全局变量污染。
- 控制周期是否一致，Simulink里的固定步长和脚本dt是否相同。
- 被控对象的初始状态是否一致，脚本里初始状态是[0.5; 0]，Simulink里积分器初始值也是这个值。
6. 抖振抑制的最后一公里：边界层加验证对比
抖振是滑模控制最直观的工程障碍。前面的内容里已经用了饱和函数sat，但边界层Delta到底取多少才能把抖振幅值压下去，同时不牺牲稳态精度，这需要专门的对比验证。
我常用的方法是以等速趋近律为主，把sign换成sat，然后扫一组Delta：0.005、0.01、0.02、0.05。每个Delta下记录控制量的高频分量RMS值和稳态跟踪误差的均方根：
% 对比不同边界层厚度的抖振抑制效果 Delta_list = [0.005, 0.01, 0.02, 0.05]; for i = 1:4 [t, x] = simulate_with_sat(p, c, K, epsilon, Delta_list(i)); u = compute_control_sequence(x, p, c, K, epsilon, Delta_list(i)); u_rms(i) = rms(u(ceil(end/2):end)); % 后半段控制量RMS e_rms(i) = rms(x(1,ceil(end/2):end)); % 稳态角度误差RMS end % 画柱状图对比 bar([u_rms; e_rms]'); legend('控制量RMS', '角度误差RMS');仿真结果大多会显示：Delta从0.005加到0.02，抖振幅值显著下降，而跟踪误差只增加一两个数量级的微小量。但Delta超过0.05后，控制量RMS不再明显下降，误差却开始上升，说明边界层已经大到把系统“滑模性质”都磨掉了。所以工程上的折中值通常在系统状态量程的1%～5%之间。
除了边界层，另一个在仿真实战中有效的做法是用观测器或扩张状态观测器估计总扰动，在等效控制里前馈补偿掉扰动，这样切换增益K可以大幅降低。仿真的对比方式依然是对照组实验——同样的K和趋近律参数，一组不带观测器，一组带一个线性扩张状态观测器，观测器带宽设为控制带宽的5～10倍，看抖振幅值是否有量级上的改善。这一步能直接回答“抖动到底是切换引起的，还是模型不确定性引起的”这个核心问题。
最后留下一个验证滑模面是否确实存在的指标：画出s(t)曲线，若s在趋近到达段结束后保持在零附近很小的邻域内，而不是大幅振荡或缓慢漂移，就说明控制律设计成立。这是滑模变结构控制MATLAB仿真程序的最后一道验收线——比看输出波形更严格，因为它直接检验了理论前提是否被数值实现真实复现。
本文还有配套的精品资源，点击获取
