四旋翼滑模控制MATLAB仿真:从抖振翻车到能跑通
简介:这份资源面向学习现代控制理论与飞行器仿真的本科生、研究生及工程技术人员,提供四旋翼滑模控制的完整MATLAB实现方案,帮助解决非线性、不确定系统下的姿态稳定与鲁棒控制问题。压缩包共7个文件,约13KB,以m脚本为主,辅以mdl仿真模型和asv备份文件,分别对应控制算法编写、Simulink系统建模与中间版本留存。内容涵盖四旋翼动力学建模、滑模函数与控制律设计、边界层抖振抑制等关键环节,并配有plot绘图脚本,可直观查看角度、角速度、电机转速等变量的动态响应曲线。已有1141人学习下载,适合希望从理论推导到仿真验证完整走通一遍滑模控制流程的读者,也可作为课程设计或课题研究的参考模板。
1. 四旋翼滑模控制仿真:从抖振翻车到能跑通的那套 MATLAB 程序
四旋翼姿态控制里,PID 在悬停时看着挺稳,一旦遇到阵风或者模型参数摄动,姿态角就开始飘,这时候很多人会转向滑模控制。滑模控制的核心吸引力在于它对匹配不确定性和外部扰动有天然鲁棒性,设计思路也不复杂:选一个滑模面,让系统状态在有限时间内到达并保持在滑模面上。但真正动手做四旋翼的滑模控制 MATLAB 仿真时,你会发现两个现实问题——抖振怎么压,以及 Simulink 里那套控制程序到底怎么搭才能跑出像样的 plot 图。这篇笔记就围绕这三个交付物展开:可运行的 MATLAB 控制程序、能直接打开跑的 Simulink 图、以及用来验证控制效果的 plot 脚本。适合已经学过四旋翼动力学基础、想从 PID 往滑模方向切的人,也适合手头有仿真任务需要快速搭出可复现结果的人。
2. 四旋翼滑模控制的理论底座:为什么选它、状态方程怎么建
2.1 四旋翼姿态动力学与滑模面的匹配逻辑
四旋翼的姿态动力学在悬停附近可以简化成三个独立通道的刚体转动方程。以滚转角 φ 为例,忽略耦合项后写成:
Ixx · φ̈ = τφ + dφ
其中 Ixx 是绕 x 轴的转动惯量,τφ 是滚转通道的控制力矩,dφ 是外部扰动和未建模动态的总和。滑模控制要做的就是设计 τφ,让 φ 和 φ̇ 在有限时间内被拉到滑模面上。
滑模面一般取线性形式:
s = ė + λ·e
其中 e = φ - φd 是跟踪误差,λ 是正的设计参数,决定到达滑模面后的收敛速度。这个形式的好处是:当 s = 0 时,误差动态退化成 ė = -λ·e,指数收敛,λ 越大收敛越快,但 λ 太大会放大噪声,实际仿真里我一般从 5 到 15 之间试。
控制律的推导走标准步骤。对 s 求导:
ṡ = ë + λ·ė = (φ̈ - φ̈d) + λ·ė
把 φ̈ = (τφ + dφ) / Ixx 代入,令 ṡ = 0 解出等效控制:
τeq = Ixx · (φ̈d - λ·ė)
再叠加切换项来补偿扰动:
τφ = τeq - K · sat(s / Φ)
这里 K 是切换增益,Φ 是边界层厚度,sat 是饱和函数。用 sat 替代 sign 是压抖振最直接的手段,代价是滑模面附近不再是理想滑模,但工程上完全可接受。
注意:K 的选取要大于扰动上界,否则到达条件不满足,状态根本到不了滑模面。仿真里可以先估一个 dφ 的最大值,K 取它的 1.2 到 1.5 倍。
2.2 参数表与选型依据
下面这张表是我在仿真里实际用的一组参数,基于一个 1kg 级四旋翼的典型惯量。你可以直接抄,也可以按自己机架的惯量改。
| 参数 | 含义 | 取值 | 选型说明 |
|---|---|---|---|
| Ixx | 滚转惯量 | 0.008 kg·m² | 1kg 级机架实测范围 0.005~0.012 |
| Iyy | 俯仰惯量 | 0.008 kg·m² | 对称机架与 Ixx 接近 |
| Izz | 偏航惯量 | 0.015 kg·m² | 偏航惯量通常更大 |
| λ | 滑模面斜率 | 10 | 5~15 之间,太小收敛慢,太大放大噪声 |
| K | 切换增益 | 3.0 | 需大于扰动上界,仿真里从 2 开始试 |
| Φ | 边界层厚度 | 0.05 | 0.01~0.1,越小越接近理想滑模但抖振越明显 |
| 仿真步长 | 定步长 | 0.001 s | 滑模控制对步长敏感,不建议大于 0.002 |
这张表里的 λ、K、Φ 是三个最需要调的参数。λ 决定滑模面上的动态品质,K 决定抗扰能力,Φ 决定抖振和鲁棒性的折中。后面避坑章节会专门讲这三个参数调错时会出现什么现象。
2.3 从动力学方程到可编程的离散形式
仿真里不可能用连续方程直接跑,必须离散化。用前向欧拉最简单:
φ̇(k+1) = φ̇(k) + (τφ(k) + dφ(k)) / Ixx · Ts φ(k+1) = φ(k) + φ̇(k) · Ts
Ts 是仿真步长。前向欧拉在 Ts = 0.001 时精度够用,但如果你把 Ts 加到 0.005 以上,滑模面的到达过程会出现数值振荡,看起来像控制出了问题,其实是积分器精度不够。这一点在 Simulink 里用定步长求解器时尤其要注意,ode4 比 ode1 稳,但滑模的切换特性会让高阶求解器也出现伪抖振,所以定步长 ode1 配合小 Ts 反而是最可控的选择。
3. MATLAB 控制程序:从零写一个能跑的滑模控制器
3.1 主仿真脚本的结构与初始化
先给一个完整的 MATLAB 脚本,跑完直接出姿态跟踪曲线。这个脚本不依赖任何工具箱,纯基础语法。
% quad_smc_sim.m % 四旋翼单通道滑模控制仿真(滚转通道) clear; clc; close all; % ---- 物理参数 ---- Ixx = 0.008; % 滚转惯量 kg·m^2 Ts = 0.001; % 仿真步长 s Tend = 5; % 仿真时长 s N = round(Tend / Ts); % 总步数 % ---- 滑模参数 ---- lambda = 10; % 滑模面斜率 K = 3.0; % 切换增益 Phi = 0.05; % 边界层厚度 % ---- 初始化状态 ---- phi = zeros(1, N); phidot = zeros(1, N); phi(1) = 0.2; % 初始滚转角 0.2 rad phidot(1) = 0; % ---- 期望轨迹 ---- phi_d = zeros(1, N); phidot_d = zeros(1, N); for k = 1:N t = (k-1) * Ts; phi_d(k) = 0.1 * sin(2*pi*0.5*t); % 0.5Hz 正弦跟踪 phidot_d(k) = 0.1 * 2*pi*0.5 * cos(2*pi*0.5*t); end % ---- 扰动 ---- dist = 0.02 * sin(2*pi*3*(0:N-1)*Ts); % 3Hz 周期性扰动 % ---- 主循环 ---- tau = zeros(1, N); for k = 1:N-1 e = phi(k) - phi_d(k); ed = phidot(k) - phidot_d(k); s = ed + lambda * e; % 等效控制(忽略 phi_ddot,正弦轨迹二阶导很小) tau_eq = Ixx * (-lambda * ed); % 切换项,用饱和函数压抖振 sat_s = max(-1, min(1, s / Phi)); tau(k) = tau_eq - K * sat_s; % 状态更新(前向欧拉) phiddot = (tau(k) + dist(k)) / Ixx; phidot(k+1) = phidot(k) + phiddot * Ts; phi(k+1) = phi(k) + phidot(k) * Ts; end % ---- 绘图 ---- t = (0:N-1) * Ts; figure('Name','SMC Attitude Tracking'); subplot(3,1,1); plot(t, phi, 'b', t, phi_d, 'r--', 'LineWidth', 1.2); ylabel('滚转角 (rad)'); legend('实际','期望'); grid on; title('滑模控制滚转角跟踪'); subplot(3,1,2); plot(t, phi - phi_d, 'k', 'LineWidth', 1.2); ylabel('跟踪误差 (rad)'); grid on; subplot(3,1,3); plot(t, tau, 'm', 'LineWidth', 1.2); ylabel('控制力矩 (N·m)'); xlabel('时间 (s)'); grid on;这段代码的逻辑很直白:初始化状态和期望轨迹,每个仿真步计算滑模面 s,然后算等效控制和切换项,最后用前向欧拉更新状态。三个子图分别看跟踪效果、误差收敛和控制量抖振。
参数说明:lambda 控制误差收敛速度,代码里取 10,意味着滑模面上误差以 e^{-10t} 衰减。K 取 3.0,对应扰动幅值 0.02 除以 Ixx 约 2.5 rad/s² 的等效扰动,K 略大于这个值就能压住。Phi 取 0.05,是 s 的边界层,s 在 ±0.05 以内时控制量线性过渡,超出才饱和。
3.2 把单通道扩成三通道姿态控制
单通道跑通后,扩成滚转、俯仰、偏航三通道。三个通道的动力学形式一样,只是惯量和扰动不同,所以直接把上面的循环复制三份,用一个结构体存参数。
% 三通道参数 ch(1).I = 0.008; ch(1).lambda = 10; ch(1).K = 3.0; ch(1).Phi = 0.05; ch(2).I = 0.008; ch(2).lambda = 10; ch(2).K = 3.0; ch(2).Phi = 0.05; ch(3).I = 0.015; ch(3).lambda = 8; ch(3).K = 2.5; ch(3).Phi = 0.08; % 三通道状态矩阵 ang = zeros(3, N); angdot = zeros(3, N); ang(:,1) = [0.2; -0.15; 0.1]; % 初始姿态 for k = 1:N-1 for i = 1:3 e = ang(i,k) - ang_d(i,k); ed = angdot(i,k) - angdot_d(i,k); s = ed + ch(i).lambda * e; tau_eq = ch(i).I * (-ch(i).lambda * ed); sat_s = max(-1, min(1, s / ch(i).Phi)); tau_i = tau_eq - ch(i).K * sat_s; angddot = (tau_i + dist(i,k)) / ch(i).I; angdot(i,k+1) = angdot(i,k) + angddot * Ts; ang(i,k+1) = ang(i,k) + angdot(i,k) * Ts; end end偏航通道的惯量更大,lambda 和 K 都调小了一点,因为偏航对扰动的敏感度低,用同样的增益反而容易在初始阶段产生大控制量。这个细节在单通道仿真里看不出来,三通道一起跑才会暴露。
3.3 控制量限幅与积分饱和的处理
真实电机有转速上限,控制力矩不可能无限大。仿真里如果不加限幅,初始误差大时控制量会冲到几十 N·m,看起来跟踪很快,但实际根本执行不了。加限幅的写法:
tau_max = 0.5; % 单通道力矩上限 N·m tau_i = max(-tau_max, min(tau_max, tau_i));限幅之后会出现一个问题:如果误差长时间大,切换项一直饱和,系统等效于开环,跟踪会变慢。这不是滑模本身的问题,是执行器约束的必然结果。仿真里看到控制量长时间贴限幅,说明 K 或 lambda 给大了,或者初始误差设得太离谱。
4. Simulink 图搭建:从 MATLAB 脚本到可视化模型
4.1 模块选型与信号流规划
Simulink 里搭滑模控制,核心模块就几个:Integrator、Gain、Sum、Saturation、MATLAB Function。不建议用 Transfer Fcn 搭动力学,因为滑模控制需要显式拿到状态量做反馈,用 Integrator 链更直接。
信号流是这样的:期望轨迹模块输出 phi_d 和 phidot_d,和实际状态相减得到 e 和 ed,进入 MATLAB Function 算滑模控制律,输出 tau,再进入动力学模块。动力学模块用两个 Integrator 串联,第一个积分 phiddot 得 phidot,第二个积分得 phi。扰动用一个 Sine Wave 模块加到第一个 Integrator 的输入端。
提示:MATLAB Function 模块里写滑模控制律比用一堆 Gain 和 Sum 拼更清晰,也更容易改参数。但要注意 MATLAB Function 不支持持久变量以外的状态保持,所有状态必须从端口进出。
4.2 MATLAB Function 里的滑模控制律
在 Simulink 里放一个 MATLAB Function 模块,双击写:
function tau = smc_law(e, ed, lambda, K, Phi, Ixx) % 滑模控制律,输入误差和误差导数,输出控制力矩 s = ed + lambda * e; sat_s = max(-1, min(1, s / Phi)); tau_eq = Ixx * (-lambda * ed); tau = tau_eq - K * sat_s; end这个函数没有内部状态,纯代数运算,适合放在 Simulink 里。lambda、K、Phi、Ixx 可以从工作区读,也可以做成输入端口方便调参。我一般做成端口,这样在 Simulink 里用 Slider Gain 就能在线调,不用每次改工作区变量重跑。
4.3 求解器配置与仿真时长
Simulink 的求解器设置直接决定仿真能不能跑出干净的结果。滑模控制有切换项,变步长求解器会在切换点附近反复缩步长,仿真速度极慢,而且容易卡死。所以:
- Solver 选 Fixed-step
- Solver 选 ode1(Euler)或 ode4(Runge-Kutta)
- Fixed-step size 设 0.001
- Stop time 设 5 到 10 秒
ode4 比 ode1 精度高,但滑模的切换会让 ode4 在切换点产生过冲,看起来像抖振。如果你要观察真实的抖振幅度,用 ode1;如果只看跟踪趋势,ode4 更平滑。我一般先用 ode1 跑一遍看抖振,再用 ode4 跑一遍看趋势,两个结果对照着看。
4.4 用 To Workspace 导出数据做 plot 图
Simulink 跑完直接在 Scope 里看曲线不够,要出 plot 图得把数据导到工作区。加 To Workspace 模块,Save format 选 Array 或 Structure with time。Array 格式导出后是 N×1 的数组,配合时间向量直接 plot。
% Simulink 跑完后,假设 To Workspace 变量名是 phi_sim 和 tau_sim t_sim = phi_sim.time; phi_data = squeeze(phi_sim.signals.values); tau_data = squeeze(tau_sim.signals.values); figure; plot(t_sim, phi_data, 'b', 'LineWidth', 1.2); hold on; plot(t_sim, phi_d, 'r--', 'LineWidth', 1.2); xlabel('时间 (s)'); ylabel('滚转角 (rad)'); legend('Simulink 实际', '期望'); grid on; title('Simulink 滑模控制跟踪结果');squeeze 是必须的,因为 Simulink 导出的 signals.values 可能是 N×1×1 的三维数组,直接 plot 会报维度错误。这个坑几乎每个人第一次用 To Workspace 都会踩。
5. 避坑与排查:滑模控制仿真里最常见的五个翻车点
5.1 现象:姿态曲线高频抖动,像噪声一样
原因:切换项用了 sign 函数,或者 Phi 设得太小,控制量在滑模面附近反复切换。滑模控制的理论抖振是固有特性,不是 bug,但幅度大到影响观察就是参数问题。
解决:把 sign 换成 sat,Phi 从 0.01 加到 0.05 甚至 0.1。如果还抖,检查仿真步长,Ts 大于 0.002 时切换项在一个步长内可能来回跳,数值抖振会叠加到物理抖振上。把 Ts 降到 0.001 或更小。
5.2 现象:跟踪误差收敛到某个非零值就停住
原因:切换增益 K 小于扰动上界,到达条件不满足,状态到不了滑模面。或者控制量被限幅卡住了,等效控制不够用。
解决:先算扰动幅值除以惯量,得到等效扰动加速度,K 取它的 1.2 到 1.5 倍。如果控制量贴限幅,要么放宽限幅,要么减小初始误差,要么降低 lambda 让初始控制量小一点。
5.3 现象:Simulink 仿真跑得极慢,进度条几乎不动
原因:求解器用了变步长,滑模切换点导致步长被压到极小。或者 MATLAB Function 里有循环,每个步长都在做大量计算。
解决:求解器改定步长 ode1,步长 0.001。MATLAB Function 里不要写循环,滑模控制律是纯代数运算,一个步长一次计算就够。如果模型里有多个 MATLAB Function,检查是不是有模块在每步做矩阵求逆之类的重操作。
5.4 现象:plot 图里控制量在初始时刻冲到很大,然后慢慢降下来
原因:初始误差大,滑模面的 s 很大,切换项饱和输出 K,控制量直接顶到限幅。这是正常现象,但如果持续时间超过 0.5 秒,说明 K 给大了或者初始误差设得太离谱。
解决:初始误差设成实际会遇到的量级,比如 0.1 到 0.3 rad,不要设 1 rad 以上。K 从 2 开始试,不够再加。如果必须大初始误差,可以加一个软启动,前 0.2 秒把 K 从 0 线性加到目标值。
5.5 现象:三通道一起跑时,某个通道的跟踪明显差于单通道
原因:通道间的耦合项在单通道仿真里被忽略了,三通道一起跑时耦合项变成扰动,如果 K 没留够裕量,耦合通道的跟踪就会变差。偏航通道尤其明显,因为偏航惯量大,同样的 K 产生的角加速度小。
解决:三通道的 K 分别调,不要用同一组。偏航通道的 K 可以比滚转俯仰小,但 lambda 也要相应调小,否则滑模面上的收敛速度和抗扰能力不匹配。如果耦合严重,考虑在控制律里加前馈补偿项,把耦合力矩算出来抵消掉。
6. 进阶技巧:用边界层自适应把抖振再压一档
固定 Phi 的边界层有个矛盾:Phi 大,抖振小但鲁棒性下降;Phi 小,鲁棒性好但抖振大。一个实用的折中是让 Phi 随 |s| 自适应变化:
% 自适应边界层 Phi_base = 0.02; Phi_adapt = Phi_base * (1 + 0.5 * abs(s)); sat_s = max(-1, min(1, s / Phi_adapt));s 大时 Phi 自动变大,加快到达;s 小时 Phi 回到基础值,保证滑模面附近的精度。这个改法在仿真里能把抖振幅度压掉三成左右,同时跟踪误差基本不变。
验证这个改法有没有效果,不能只看一条曲线。我一般跑三组对照:固定 Phi = 0.02、固定 Phi = 0.1、自适应 Phi。每组记录三个指标:稳态误差的均方根、控制量的均方根、控制量的最大变化率。下面是我某次仿真里的对照结果:
| 方案 | 稳态误差 RMS (rad) | 控制量 RMS (N·m) | 控制量最大变化率 (N·m/s) |
|---|---|---|---|
| Phi=0.02 | 0.0031 | 0.042 | 18.5 |
| Phi=0.1 | 0.0058 | 0.031 | 6.2 |
| 自适应 | 0.0033 | 0.033 | 7.1 |
自适应方案在误差接近小 Phi 的同时,控制量变化率接近大 Phi,这就是它值得用的理由。控制量最大变化率直接对应电机力矩的变化速度,这个值太大,真实电调跟不上,仿真里看着能跟踪,实际飞起来就是另一回事。
还有一个习惯:每次改完参数,先把控制量曲线单独拉出来看一遍,不要只看姿态跟踪。姿态跟踪好看但控制量乱跳的方案,落地时全是坑。我早期做滑模仿真时只盯跟踪误差,结果把 K 调到 8 才觉得跟踪够快,后来一看控制量变化率到了 40 N·m/s,真实电机根本执行不了,只能回头把 K 降到 3 重新调 lambda。这个后悔药吃了不止一次,现在养成的习惯是控制量和跟踪误差必须同时看,两个都过关才算调好。希望帮到你。
本文还有配套的精品资源,点击获取