简介面向控制工程与自动化领域学习者这份压缩包围绕滑模变结构控制的核心理论与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 thetax2 theta_dot控制量u tau则状态方程为x1_dot x2x2_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_dot1s_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行为没有本质区别但原因在于——当状态反复穿越s0时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(:,k1) xk dt*k2; end关键参数是dt取0.001秒是底线离散控制系统仿真经验法则是dt要小于控制系统带宽对应时间常数的五分之一。滑模控制的切换频率非常高dt太大会让抖振幅度被数值放大仿真结果失真。3. 切换面设计c参数是滑模面的灵魂3.1 c与系统特征根的关系切换函数s c·x1 x2 0定义了相平面中的一条直线——滑模面。当s0时有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(c2, c5, c10);注意直接调节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.5310大于机械带宽3倍以上小于驱动饱和对应值永磁同步电机速度环820电流环带宽约100Hz速度滑模面要在电流环响应时间内飞行器姿态角回路26执行机构带宽低c过大会与舵机延迟耦合惯性较大的双连杆臂14结构频率低c太大会激起柔性模态注意c取值必须低于控制频率对应奈奎斯特频率的十分之一否则直接发散。3.3 调c时的典型误判把跟踪误差大归因于K不足一个常见的误操作系统状态始终不收敛初学者会习惯性增大切换增益K。但滑模控制的误差从两个源头来一是到达阶段误差而这由趋近律参数决定二是滑模面收敛误差由c决定。如果c太小即便K足够大系统到达滑模面之后滑动得很慢误差曲线会呈现“缓坡”而非快速归零。我判断c是否合适的方法是看仿真输出的相轨迹图figure; plot(x(1,:), x(2,:)); xlabel(角度(theta)); ylabel(角速度(dtheta));如果相轨迹先是迅速冲向相平面直线s0然后在s0附近小幅度振荡滑向原点说明c设置基本合适。如果相轨迹绕圈不贴近s0那就是c过小。4. 趋近律与抖振源头参数ε、K的仿真标定4.1 符号函数在数值仿真里的三种实现方式s_dot -K·sign(s)是最基本的等速趋近律sign(s)在连续域中是一个理想开关但在数值积分中会导致状态在s0附近来回穿越。MATLAB里实现切换项有几种选择每种都会影响仿真结果实现方式表达式数值特性适用场景signsign(s)切换最尖锐抖振最明显理论演示、验证滑模存在性sats/Δ (sΔ); sign(s) (其他)tanhtanh(s/ε)全程光滑但控制增益被软化需要连续控制量的场景仅用sign会导致两个数值问题第一状态在滑模面两侧以最高频率来回切换ode45会不断缩小步长仿真时间可能膨胀数倍第二真正的物理系统是不存在无限频率切换的仿真会高估实际系统的高频分量。4.1.1 工程上怎么在仿真里处理sign工程定位是通过给连续模型做一个执行机构带宽限制来反映现实比如加一个一阶惯性环节1/(tau_act*s1)。但如果只是验证控制器本身用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.010.05。4.2 指数趋近律的仿真脚本与参数标定等速趋近律的问题是靠近滑模面时s_dot仍很大会产生过冲。工程上更常用指数趋近律s_dot -ε·sign(s) - kq·s其中ε0保证到达条件kq0系数指数收敛项。把趋近律代入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.52视为合理的起始区间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.1msepsilon和K的合适参数区间会明显下移。原因是步长越细数值积分对高频切换的解析越真实系统越容易识别出真实的抖振频率步长太粗时抖振被平均化了epsilon可以取更大而不发散。由此建议的标定顺序是先固定dt为期望控制周期如1ms。用更大的epsilon观察抖振幅值是否在可接受范围。再调kq观察到达滑模面的时间是否满足指标。所有参数调完后将dt降到0.2ms重新仿真如果抖振明显变大适当减小epsilon。这也是为什么仿真参数不能直接搬到实际硬件上跑的深层原因数字仿真的步长和控制器的采样时间是两个独立变量混淆它们会导致设计失察。5. 用MATLAB S函数把滑模控制器做进Simulink5.1 S函数的结构与滑模控制器的S函数实现纯脚本仿真适合调参和验证但如果你最终要在Simulink里搭整个机电系统伺服驱动、负载模型、观测器把控制器写成一个S函数是最干净的方式。S函数本质上是一个MEX或MATLAB函数其回调函数被Simulink执行。一个完整的两输入三输出S函数框架输入是状态x输出是控制量ufunction [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的输入接MuxMux的输入接被控对象的状态S-Function输出接被控对象的输入端。用Scope观察角度theta和控制量u。S参数块设置的关键项表格参数设置值说明S-Function namesmc_sfun必须与.m文件名一致文件名不含中文S-Function parameters8, 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然后扫一组Delta0.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和趋近律参数一组不带观测器一组带一个线性扩张状态观测器观测器带宽设为控制带宽的510倍看抖振幅值是否有量级上的改善。这一步能直接回答“抖动到底是切换引起的还是模型不确定性引起的”这个核心问题。最后留下一个验证滑模面是否确实存在的指标画出s(t)曲线若s在趋近到达段结束后保持在零附近很小的邻域内而不是大幅振荡或缓慢漂移就说明控制律设计成立。这是滑模变结构控制MATLAB仿真程序的最后一道验收线——比看输出波形更严格因为它直接检验了理论前提是否被数值实现真实复现。本文还有配套的精品资源点击获取