简介:一套用于二自由度机械臂关节角度控制的MATLAB滑模控制仿真源码包,面向机器人控制、自动化相关专业学生与工程师。资源以滑模面设计与控制律构成为核心,滑模面选取目标姿态与当前状态之差并减去角速度加权项,控制扭矩则分成线性项与含滑模参数的切换项,便于理解滑模控制在机械臂轨迹跟踪中的实现方式。包内含8个文件,以m脚本、slx仿真模型、png图像及pdf文档为主,脚本覆盖主程序、绘图与动力学计算,slx适合可视化建模仿真,png为控制输入、相平面、位置响应等结果图,pdf可作原理对照,压缩包整体仅278KB,下载后便于快速部署与二次修改。目前已有91人学习下载,适合用来做课程设计、毕业设计或入门滑模控制的辅助参考,能帮助读者结合源码和仿真图直接看到控制效果,缩短调参时间。
1. 二自由度机械臂滑模控制,先从误差和参数失配说起
当一台二自由度机械臂在MATLAB仿真里做轨迹跟踪时,真正让人改参数改到烦的往往不是参考轨迹给得太快,而是模型内的惯量、连杆长度和实际对象总有偏差。用PID,增益低了跟不上高速轨迹,增益高了测量噪声又被放大;用计算力矩法,则要求M(q)和C(q,qd)完全已知。滑模控制换了一个思路:承认模型不准确,把不匹配的那部分当作有上界的扰动,用一个不连续切换项把它压下去。这个思路简单直接,代价是控制力矩会抖振。下面按“被控对象建模、控制器设计、参数调优、验证”的顺序,把滑模控制完整落地成MATLAB源码。
2. 二自由度机械臂动力学模型:先把被控对象写成MATLAB能算的形式
2.1 为什么控制器设计之前要重新写一遍动力学
很多入门者拿到滑模控制论文后,先把控制器代码写好,再回头补动力学,这样容易埋雷。滑模控制律里的等效控制项是从名义模型推出来的,如果被控对象方程里的惯量矩阵符号或者科氏力方向不对,控制器里每一项都会带着同样的错误,调参时只会觉得“怎么都不收敛”。
常见做法是把机械臂放在水平平面内,只保留两个旋转关节,这时重力项G(q)=0,摩擦和负载变化全部归入扰动项。这样做的好处是动力学方程足够简洁,能看出M(q)和C(q,qd)的结构;等仿真通了,再往G(q)里加重力项,思路并不会变。
2.2 拉格朗日方程给出的动力学方程
对水平面内刚性二连杆机械臂,用拉格朗日方程能得到如下标准形式:
M(q) * qdd + C(q, qd) = tau
状态取x = [q1; q2; q1d; q2d],其中q1和q2是两个关节角。惯量矩阵和科氏力向量为:
M(q) = [ (m1+m2)l1^2 + m2 l2^2 + 2 m2 l1 l2 cos(q2), m2 l2^2 + m2 l1 l2 cos(q2); m2 l2^2 + m2 l1 l2 cos(q2), m2 l2^2 ]
C(q,qd) = [ -m2 l1 l2 sin(q2) (2 qd1 qd2 + qd2^2); m2 l1 l2 sin(q2) qd1^2 ]
这里m1/m2是连杆质量,l1/l2是连杆长度。M矩阵是对称正定的,只要q2不等于使矩阵奇异的位置,M^{-1}存在;C中有sin(q2)项,说明关节2的转动会产生作用在关节1上的科氏力,这是二自由度机械臂和单关节模型最大的区别。
| 符号 | 含义 | 单位 | 典型初值 |
|---|---|---|---|
| m1, m2 | 连杆1、2质量 | kg | 1.0, 1.2 |
| l1, l2 | 连杆1、2长度 | m | 1.0, 0.9 |
| q1, q2 | 关节角 | rad | 0, 0 |
| tau1, tau2 | 关节力矩 | N·m | — |
参数具体数值不影响控制律结构,只影响仿真曲线。
2.3 参数不确定性从哪里进入滑模控制
实际对象用真实参数p_true描述,控制器只能用名义参数p_nom。真实方程与名义方程之差,最终会以M_true^{-1} delta的形式直接叠加在控制加速度上,和输入tau在同一个通道。这类影响称为匹配不确定性,滑模控制能完全抑制匹配不确定性,只要切换增益大于其幅值上界。这一点是滑模控制器抗扰能力的来源,也是后面调K参数时为什么总有“一定要大于某个值,否则误差不收敛”现象的根源。
2.4 MATLAB源码:动力学与惯量矩阵函数
先写一个公共函数arm_dynamics,返回M矩阵和C向量。控制器和仿真对象都调用它,避免两处公式不一致。
function [M, C] = arm_dynamics(q, qd, p) % 二自由度平面机械臂动力学:惯量矩阵M和科氏/离心向量C c2 = cos(q(2)); s2 = sin(q(2)); m1 = p.m1; m2 = p.m2; l1 = p.l1; l2 = p.l2; M = [(m1 + m2)*l1^2 + m2*l2^2 + 2*m2*l1*l2*c2, m2*l2^2 + m2*l1*l2*c2; m2*l2^2 + m2*l1*l2*c2, m2*l2^2]; C = [-m2*l1*l2*s2*(2*qd(1)*qd(2) + qd(2)^2); m2*l1*l2*s2*qd(1)^2]; end说明:这个函数把M和C放在一起,是因为C向量由M的偏导得到,分开写反而容易写错符号。注意M11里cos(q2)前的系数是2m2l1*l2,M12和M21相同,这样矩阵才是对称的。
被控对象函数把状态求导结果返回给ode求解器:
function dstate = two_dof_arm(t, state, tau, p) % t和state是ode45的标准接口,tau由控制律算好再传进来 q = state(1:2); qd = state(3:4); [M, C] = arm_dynamics(q, qd, p); ddq = M \ (tau - C); % 注意是矩阵左除,不是求逆 dstate = [qd; ddq]; end代码里用M \ (tau - C)而不是inv(M)*(tau - C),数值上更稳,MATLAB对2x2矩阵也会走高斯消元。state的布局决定了后续滑模面计算的下标,建议全篇统一用state(1:2)做角度、state(3:4)做角速度,不要一套代码里换布局。
提示:这段动力学要求MATLAB基础环境即可运行,不依赖Simulink或额外工具箱。如果用的是较早的MATLAB版本,注意把主脚本里的局部函数拆成独立m文件,避免脚本局部函数语法不支持。
3. 滑模控制律设计:等效控制加切换项就是MATLAB源码的全部核心
3.1 滑模面不是误差,而是误差的一阶组合
对轨迹跟踪问题,位置误差e = q - qd,速度误差de = qd - dqd。如果只用s = de,系统只在速度上锁定,位置误差可能漂泊;如果只用s = e,控制器退化成纯位置反馈,无法约束速度。二自由度机械臂的控制目标同时包含角度和角速度,因此最常见的选择是:
s = de + lambda .* e
lambda是正数,控制了滑模面的带宽。lambda越大,等效误差动力学de + lambda*e = 0收敛越快,但同时会把高频测量噪声放大,后面调参表里会再说明。对两个关节分别选lambda,可以不是相同值,但一般先设成一致。
3.2 等效控制:让系统轨迹保持在滑模面上
滑模控制的第一个步骤是求等效控制u_eq,它的含义是:假设名义模型完全准确,并且系统已经在滑模面上(s=0),那么要让sdot也等于0,控制力矩应该取多少。
令sdot = dde + lambda*de = 0,代入动力学qdd = M0^{-1}(u - C0),可以解出:
u_eq = M0(q) * (ddqd - lambda*de) + C0(q, qd)
注意这里用的是名义参数下的M0和C0,不是真实参数。所以就算模型失配,控制律也能写出确定数值。等效控制在理想情况下把非线性机械臂解耦成两个独立的二阶误差系统,剩下的模型误差和外部扰动交给切换项处理。
推导中我用的是水平面模型,重力项为零;如果目标机械臂在垂直平面内,需要在u_eq里额外加G0(q),公式变成u_eq = M0*(ddqd - lambda*de) + C0 + G0,其他不变。
3.3 趋近律:如何把状态拉回滑模面
等效控制只能在滑模面上维持滑模,没法把不在滑模面上的轨迹拉回来。所以在u_eq后面必须加上切换项u_sw = -K*sat(...)或-K*sign(s)。负号来自李雅普诺夫条件,要让sdot*s < 0,切换方向必须和s相反。
工程上常用的趋近律有几种,差别在于到达滑模面速度和到达后的抖振表现:
| 趋近律 | 数学形式 | 参数作用 | 特点 |
|---|---|---|---|
| 等速趋近 | sdot = -eps*sign(s) | eps决定切换项幅度 | 结构最简单,但恒幅切换,抖振明显 |
| 指数趋近 | sdot = -epssign(s) - ks | eps保证到达,k加速收敛 | 工程最常用,k项让误差大时趋近更快 |
| 幂次趋近 | sdot = -kabs(s)^alphasign(s), 0<alpha<1 | 越靠近滑模面切换越小 | 接近s=0时更平滑,但鲁棒性稍弱 |
二自由度机械臂仿真里最稳妥的是指数趋近。它比等速趋近多一个k*s项,误差大时趋近速度快,误差小时切换幅值不变,好调。下面的源码采用指数趋近的离散实现,但在符号处换成饱和函数以减少数值抖振。
3.4 MATLAB源码:命名为smc_controller的滑模控制函数
控制器的输入是机械臂状态state和参考轨迹ref(包含qd、dqd、ddqd),输出是两个关节力矩tau。该函数内部嵌套一个饱和函数sat,避免直接用sign导致ode45步长收缩。
function tau = smc_controller(state, ref, p) % 二自由度机械臂滑模控制器:等效控制 + 指数趋近 + 边界层饱和 q = state(1:2); qd = state(3:4); e = q - ref.qd; % 位置误差 de = qd - ref.dqd; % 速度误差 s = de + p.lambda .* e; % 滑模面 [M0, C0] = arm_dynamics(q, qd, p); % 名义模型 u_eq = M0 * (ref.ddqd - p.lambda .* de) + C0; u_sw = -p.K .* sat(s, p.phi); % 边界层内的等效线性反馈 tau = u_eq + u_sw; function y = sat(s, phi) y = s / phi; out_a = abs(y) > 1; y(out_a) = sign(y(out_a)); end end解释:第一段取状态并计算s,第二段用名义参数算M0和C0,第三段合成为力矩。这里u_eq里的ref.ddqd - lambda*de就是期望的“加速度修正项”,当系统在滑模面上时,靠近该项的轨迹会满足de + lambda*e的动态。u_sw的K如果设成对角矩阵,可以给两个关节不同的切换增益;如果设成标量,则两个关节同样处理。边界层饱和函数的含义是当|s|>phi时输出±1,当|s|<=phi时输出s/phi,这段线性区域等价于一个比例控制器,是抑制抖振的关键。
需要特别注意的是参数失配问题:smc_controller内部调用的arm_dynamics用的是控制器参数结构体p,而仿真对象two_dof_arm用的是真实参数结构体。两者可以相同,也可以人为失配,失配后只要K的幅值还能覆盖误配差,跟踪误差就会维持在边界层宽度量级。
4. 二自由度机械臂滑模控制的仿真主程序与参数调优
4.1 主程序:参考轨迹、闭环函数和ode45
我用m脚本而不是Simulink来跑仿真,原因是调参和断点排错更直接。参考轨迹选五次多项式,起止位置、速度和加速度都连续,避免在仿真启动时就给控制器一个跳变误差。
% main_smc.m % 真实对象参数 p_true.m1 = 1.0; p_true.m2 = 1.2; p_true.l1 = 1.0; p_true.l2 = 0.9; % 控制器名义参数,m2人为失配 0.2kg p_ctrl = p_true; p_ctrl.m2 = 1.0; p_ctrl.lambda = [6; 6]; p_ctrl.K = [18; 22]; p_ctrl.phi = 0.05; p_ctrl.t0 = 0; p_ctrl.tf = 2; p_ctrl.q0 = [0; 0]; p_ctrl.qf = [0.6; -0.35]; x0 = [p_ctrl.q0; 0; 0]; opts = odeset('RelTol', 1e-7, 'AbsTol', 1e-8); [t, X] = ode45(@(t, x) closed_loop(t, x, p_true, p_ctrl), [0 4], x0, opts); plot(t, X(:, 1:2)); xlabel('t/s'); ylabel('q/rad'); legend('q1', 'q2', 'Location', 'best');脚本末尾放两个局部函数:closed_loop把控制器和对象串起来,ref_traj生成五次多项式轨迹。如果读者用的MATLAB版本R2016a以下,就把这两个函数拆成单独m文件。
function dstate = closed_loop(t, x, p_true, p_ctrl) ref = ref_traj(t, p_ctrl.t0, p_ctrl.tf, p_ctrl.q0, p_ctrl.qf); tau = smc_controller(x, ref, p_ctrl); dstate = two_dof_arm(t, x, tau, p_true); end function ref = ref_traj(t, t0, tf, q0, qf) if t <= t0 tau_n = 0; elseif t >= tf tau_n = 1; else tau_n = (t - t0) / (tf - t0); end ref.qd = q0 + (qf - q0) * (10*tau_n^3 - 15*tau_n^4 + 6*tau_n^5); ref.dqd = (qf - q0) * (30*tau_n^2 - 60*tau_n^3 + 30*tau_n^4) / (tf - t0); ref.ddqd = (qf - q0) * (60*tau_n - 180*tau_n^2 + 120*tau_n^3) / (tf - t0)^2; end说明:closed_loop中的p_true和p_ctrl分开传递,控制器只知道名义参数,被控对象用真实参数,这样仿真结果才有说服力。ref返回一个结构体,控制器直接取qd、dqd、ddqd三个字段;因为仿真时间超过tf后tau_n保持为1,参考速度自动归零,机械臂会停在目标位置。
4.2 三个必调参数:lambda、K和边界层phi
直接在脚本里把lambda、K、phi改成不同值跑几轮,比看任何理论分析都直观。给出建议的调节方向和初始范围:
| 参数 | 初值范围 | 调大后 | 调小后 |
|---|---|---|---|
| lambda | 2到8 | 误差收敛更快,但对测量噪声更敏感,等效控制中的ddqd项会被放大 | 轨迹跟踪显得“软”,滞后明显 |
| K | 10到30 | 抵抗失配和扰动的能力更强,抖振幅值同步上升 | 模型失配稍大时就压不住,稳态误差变大 |
| phi | 0.01到0.1 | 抖振明显下降,稳态误差按比例增大 | 更接近纯sign,抖振回升,但精度高 |
K的选择有个实用下界:让模型失配为零时,误差能快速收敛,然后逐步加大失配或外加阶跃扰动,直到误差曲线仍保持有界。如果发现只要p_ctrl.m2改动5%误差就开始发散,说明K不够,不是控制器结构错了。lambda的上界由执行器带宽决定:机械臂控制频率50Hz时,lambda取到8已经会让u_sw在切换瞬间产生明显冲击,需要配合phi来吸收。
4.3 抖振怎么量化:比看示波器更可靠的做法
抖振不能只看力矩曲线“毛不毛”,更好的办法是统计控制力矩的高频符号变化次数。在仿真后处理里,用闭环输出的状态重新计算控制力矩,然后数diff发生符号跳变的次数:
% chattering_eval.m tau_rec = zeros(numel(t), 2); for i = 1:numel(t) x = X(i, :).'; ref = ref_traj(t(i), p_ctrl.t0, p_ctrl.tf, p_ctrl.q0, p_ctrl.qf); tau_rec(i, :) = smc_controller(x, ref, p_ctrl); end flip_cnt = sum(sum(diff(sign(tau_rec)) ~= 0));每次仿真结束记录一版flip_cnt,和最大跟踪误差一起对比。这样做的好处是把“抖振感”变成数字;调phi时能同时看两个指标,就不会只凭感觉往一个方向调。需要注意的是这个统计受ode45输出步长影响,步长越密统计越准,比较时保持同一套opts。
4.4 典型失败现象和排查顺序
仿真跑挂有几个高概率原因,按检查顺序列出:
第一,ode45报步长接近0或无法继续积分。这通常是因为sign(s)的切换让系统右侧不连续,solver在每个切换点不断缩步长。先把phi设置成0.05以上的饱和函数,再把RelTol加到1e-7,大多数情况立刻恢复正常。
第二,跟踪曲线在启动瞬时就飞掉。先检查ref_traj在t0时刻的ddqd是否为0,如果参考轨迹给的是阶跃或者正余弦起点,初始加速度跳变会直接灌进u_eq。停在t=t0附近打成断点,打印u_eq数值,看看有没有数量级异常大的成分。
第三,两个关节都收敛,但关节2稳态误差一直降不到边界层范围内。这多半不是控制律问题,而是m2作为控制输入时,模型失配超出了K的覆盖能力;把K加大或者把phi改小即可。
提示:如果不需要Simulink的可视化反馈,用ode45跑m文件比Simulink更容易定位源码问题。滑模控制的抖振在Simulink中会表现为仿真步长异常,排查成本更高。
5. 用相平面和自适应边界层验证二自由度机械臂滑模控制的边界
5.1 相平面看到达条件和滑模维持情况
滑模控制理论上保证s会趋向0,但实际仿真里s只在边界层内振荡。要验证控制器是否真的把系统拉到了滑模面,光看q的跟踪曲线不够,因为误差小可能只是lambda大,未必说明滑模条件满足。更直接的方法是画s和sdot的相轨迹:
s_all = zeros(numel(t), 2); ds_all = zeros(numel(t), 2); for i = 1:numel(t) x = X(i, :).'; ref = ref_traj(t(i), p_ctrl.t0, p_ctrl.tf, p_ctrl.q0, p_ctrl.qf); e = x(1:2) - ref.qd; de = x(3:4) - ref.dqd; s_all(i, :) = de + p_ctrl.lambda .* e; tau = smc_controller(x, ref, p_ctrl); [M_true, C_true] = arm_dynamics(x(1:2), x(3:4), p_true); ddq = M_true \ (tau - C_true); dde = ddq - ref.ddqd; ds_all(i, :) = dde + p_ctrl.lambda .* de; end plot(s_all(:, 1), ds_all(:, 1), '.'); xlabel('s1'); ylabel('s1dot');这段代码用后处理重新计算每个输出点的s和sdot,不干扰ode45的积分过程。如果相轨迹从初始点出发后能很快进入±phi的窄带区域并在其中打转,说明到达条件成立;如果s1dot和s1符号相反方向的时间很多,说明切换增益K还不足以弥补扰动,先增大K而不是急着改lambda。
5.2 让边界层厚度随时间自适应
最后一个常用技巧是把固定的phi改成时变函数:初期误差大时允许较宽的边界层,避免大切换幅值带来的剧烈抖振;后期系统接近稳态时再收缩边界层,把稳态误差压下去。在smc_controller里把sat的第二参数改为:
phi_t = p.phi_inf + (p.phi - p.phi_inf) * exp(-p.alpha * t); u_sw = -p.K .* sat(s, phi_t);其中phi是初始边界层宽度,phi_inf是稳态边界层下限,alpha控制收缩速度。主脚本里补上p.phi_inf = 0.02; p.alpha = 3;即可。例如phi=0.1、phi_inf=0.02、alpha=3,误差大时等效边界层为0.1量级,2秒后收缩到0.02。这样可以在不显著牺牲鲁棒性的前提下把抖振降低一半以上,也是工程上把滑模控制从论文搬到实物前最常做的一步。把alpha分别设成2、5、10跑三次,观察饱和函数在边界层内部的翻转频率,这个趋势比任何固定参数的选择都更适合作为二自由度机械臂滑模控制的最终调参依据。
本文还有配套的精品资源,点击获取