资讯详情

资讯详情

Matlab四旋翼无人机PID控制仿真:从零构建代码与调参实战

搞控制或者说搞机器人的朋友应该都经历过这个阶段看了无数篇四旋翼无人机 PID 控制的文章手里的代码却要么跑不起来要么一跑就发散对着曲线完全不知道问题出在哪。我当年入门的时候也在这种坑里趴了很久后来用 Matlab 从零写了一个完整的四旋翼无人机 PID 控制仿真把模型、控制律、参数整定、结果可视化全部串在一起才算真正把“PID 是怎么把飞机稳住”这件事吃透了。这篇教程就是按我当时的学习路径整理的。整体用 Matlab 脚本实现不依赖 Simulink也不需要额外工具箱完整代码直接贴在后面复制保存为.m文件就能运行。仿真对象是一个简化版四旋翼模型包含高度和俯仰两个通道控制律采用位置式 PID并加入了重力补偿。跑通之后你能清楚看到高度、姿态在 PID 控制下的动态响应也能直观理解 kp、ki、kd 这三个参数到底在干什么为后续做级联 PID 或者真机调参打个底。这个仿真适合三类读者刚接触控制理论、想找个具体对象练手的学生想快速验证 PID 调参思路、又不想一上来就碰硬件的入门工程师以及准备做真机实验但想先在低成本环境里把逻辑理清的小伙伴。接下来我按“从思路到代码再到排错”的顺序把整个搭建过程完整讲一遍。1. 仿真的整体思路新手为什么要从模型入手1.1 仿真闭环帮你建立控制直觉很多新手拿到一段控制代码第一件事就是到处调 Kp、Kd调了半天也不知道为什么发散、为什么稳态误差消不掉。归根结底是不清楚被控对象长什么样。仿真最大的价值不是真去复现一架无人机而是用数学方程把你对“物理系统如何运动”的理解变成一个可以随时观测的闭环给一个期望值看输出怎么响应控制量怎么变化误差怎么收敛。这个闭环一旦在脑子里成形你再去碰真机、看飞控日志、调参就不会全靠猜。仿真里的每一次发散都是在帮你建立“参数和响应之间的直觉”某个环节加大会有什么趋势某个环节减小会不会稳一点。这种直觉是任何文档都教不会的只能自己跑出来。我在这个项目里刻意没有用 Simulink原因有两点。一是脚本代码的每一步都显式可见误差怎么算、PID 怎么输出、模型怎么更新一目了然对新手来说没有“黑盒子”二是代码方便改造想换控制律、加扰动、改成级联 PID都是在文本层面直接操作比拖拽 Simulink 模块直观得多。1.2 简化模型是从哪来的完整四旋翼无人机模型是六自由度的三个平动位置x、y、z加三个姿态角roll、pitch、yaw方程之间耦合严重非线性项也不少。对新手来说一上来就把完整的牛顿-欧拉方程糊在脸上很容易被劝退就算代码写出来了一个小错误也会导致完全无法收敛排错难度很高。所以这个仿真里我做了一个经典简化只考虑悬停附近的运动把高度通道和俯仰通道解耦出来。换句话说我们先控制 z 方向高度和俯仰角 phi忽略横滚、偏航以及它们之间的耦合把四旋翼当成两个独立的单输入单输出系统来处理。高度通道很简单就是牛顿第二定律。记无人机质量为 m四个电机产生的总升力对应的等效加速度为 u1竖直方向运动方程为m * z u1 * m - m * g两边同时除以 m得到z u1 - g这里 u1 的单位已经是加速度单位不是力。好处是控制量跟质量解耦后面调参数的时候不用反复改质量。俯仰通道对应刚体转动定律。记俯仰转动惯量为 I_yy控制俯仰的力矩为 u2则phi u2 / I_yy这两个方程已经足够演示 PID 的核心控制逻辑给高度一个目标PID 算出需要的加速度再加上重力补偿项 g就是油门指令给俯仰角一个目标PID 算出需要的角加速度乘上转动惯量就是力矩指令。有读者可能会问简化掉横滚和偏航仿真还有意义吗我的看法是对于入门阶段先把两个通道吃透远比一次性面对六个通道的耦合要有效。而且等你理解了内环外环的分层思想后把横滚、偏航通道往代码里补只是同样的模式复制几遍而已。1.3 模型、控制器、可视化缺一不可一套完整的仿真至少包含三个部分被控对象模型、控制器算法、数据可视化。模型是你对物理世界的近似描述控制器是你要验证的算法而可视化是帮助你判断算法好坏的“眼睛”。三者缺一不可。这个仿真里还有第四个隐含部分就是主循环本身。它模拟了控制器的运行节拍每隔一个固定周期 dt读取一次状态计算一次控制量更新一次模型。这个“周期性采样”的概念对理解数字控制很关键真机飞控就是在固定的控制频率下不断重复这个过程和仿真主循环没有本质区别。2. PID控制的原理与参数设计逻辑2.1 比例、积分、微分到底在干什么PID 三个字母对应的三个环节可以用一句不太严谨但很好记的话概括比例管现在积分管过去微分管未来。比例项 Kp * e(t) 直接放大当前误差。误差越大控制量越大系统回正的趋势就越强。只有比例项时系统经常会出现稳态误差或者来回振荡。打个比方你开车看到偏离车道打方向盘的幅度只跟当前偏差成正比那速度一快就容易冲过头速度慢了又回不来。积分项 Ki * ∫e(t)dt 把历史上累计的所有误差加回来。它专门处理稳态误差。比如这个仿真里飞机悬停时重力始终向下如果控制量里只有比例项油门指令里没有多出来的重力补偿部分飞机就会持续往下掉最终停在某个高度再也上不去。积分项就是为了把这种“持续存在的小偏差”累积成足够的控制量去抵消它。缺点是积分容易过头积分项太大会导致超调明显。微分项 Kd * de(t)/dt 看的是误差变化趋势。误差在快速减小微分项就产生反向力提前“踩刹车”抑制超调。它会放大高频噪声所以真机上微分项通常要配低通滤波。把这三项合在一起就形成了完整的位置式 PID 控制律u Kp * e Ki * integral(e * dt) Kd * de/dt2.2 离散化与代码里的实现细节因为计算机和微控制器都是离散系统需要把连续 PID 离散化。我用的位置式离散形式是积分累加err_int err * dt; 微分近似err_dot (err - err_prev) / dt; 控制输出u Kp*err Ki*err_int Kd*err_dot;积分项直接累加误差与步长的乘积这是最简单的矩形积分微分项用一阶差分近似虽然对噪声敏感但在仿真环境下完全够用而且逻辑最直白。代码里我实际写的微分项有一点变化值得单独说明。因为期望高度和期望俯仰角是常数误差的导数就等于状态量的负导数e z_ref - z; e_dot 0 - zdot -zdot;所以微分项直接写成Kd_z * (-zdot)。这个写法和Kd * (err - err_prev) / dt在数学上等价但直接调用速度状态不会因为差分放大数值噪声数值上更干净。2.3 参数整定为什么是那个顺序网上一搜 pid 参数整定方法Ziegler-Nichols、临界比例度、衰减曲线法一大堆。但对新手来说最有用的还是先建立“参数和响应之间的直觉”。我的建议顺序永远是先调 P再调 D最后调 I。第一步把 Ki、Kd 全部设成 0只加 Kp。从小往大慢慢加观察曲线。Kp 太小响应慢稳态误差大Kp 太大系统开始等幅振荡。这一步能让你直观感受系统的开环特性和临界振荡点。第二步加 Kd。当曲线出现明显超调时开始加Kd 的作用是“刹车”让响应更快稳定下来。这里要注意 Kd 方向千万不能反方向反了会加速振荡甚至直接发散。第三步加 Ki。在 P、D 已经把系统基本稳定住的条件下如果发现最终的稳定值和期望值有偏差再加入 Ki 去消除稳态误差。我代码里给的高度环参数是 Kp8、Ki0.5、Kd5俯仰环是 Kp60、Ki10、Kd8。这些参数不是拍脑袋写的可以做个粗略推导。高度环忽略积分项后闭环特征方程近似为s^2 Kd_z * s Kp_z 0代入 Kp8、Kd5自然频率 wn sqrt(8) ≈ 2.83 rad/s阻尼比 zeta 5 / (2sqrt(8)) ≈ 0.88。这个阻尼比接近临界阻尼意味着响应快且几乎不超调收敛时间约 4 / (zetawn) ≈ 1.6 秒。俯仰环 Kp60、Kd8wn ≈ 7.75 rad/szeta ≈ 0.52姿态响应更快但会有少量超调实际调试时把 Kd 提高到 12 左右可以让姿态更稳。像这样先根据二阶系统公式估一个范围再微调参数比盲目试要高效得多。3. 从零搭建完整代码实现与调试3.1 环境准备与代码结构这个项目对 Matlab 版本没有特殊要求R2018a 以上都可以不需要安装额外工具箱。如果你机器上还没装 Matlab安装过程注意一下许可证激活装好以后新建一个脚本文件命名为quad_pid_sim.m把完整代码粘贴进去点击运行就行。整个代码分成四段第一段是参数初始化物理参数、仿真步长、期望值、PID 参数、初始状态全部集中在这里方便维护。第二段是主循环每过一个控制周期计算误差、PID 输出、更新模型、记录数据。第三段是绘图把高度、姿态、控制量画出来直观判断控制效果。第四段是终端输出直接打印最终的收敛结果方便快速验证。这样的组织方式对后期扩展非常友好。比如你想改期望高度为方波只需要定义z_ref_k放进循环想加扰动只需要在模型更新那行叠加一个外力项。3.2 主循环的四个关键步骤主循环是整个仿真器的核心先看骨架再看细节for k 1:N 计算当前误差 PID 计算控制量 模型积分更新状态 记录数据 end每一步都对应真实控制器的工作流程。dt0.01 就是控制周期对应 100Hz 的控制频率。真实飞控一般在 500Hz 到 1kHz100Hz 对入门仿真足够了。第一步计算误差。直接用期望值减当前状态err_z z_ref - z、err_phi phi_ref - phi。第二步PID 计算。这里高度通道和俯仰通道的写法略有不同高度通道因为要克服重力在 PID 输出之外必须叠加一个重力补偿项a_z_des Kp_z * err_z Ki_z * err_z_int Kd_z * (-zdot); u1 a_z_des g;新手最容易忽略的就是这个重力补偿。只把 PID 输出当成油门控制量里没有重力对应部分飞机永远飞不到目标高度最终会出现很大的稳态误差。加了 g 之后PID 输出的加速度全部用来纠正高度偏差物理含义更清晰。俯仰通道的控制量是力矩所以把 PID 输出的期望角加速度乘上转动惯量phi_ddot_des Kp_phi * err_phi Ki_phi * err_phi_int Kd_phi * (-phidot); u2 phi_ddot_des * I_yy;第三步模型更新。这一步本质是数值积分我用了最简单的欧拉法zdot zdot (u1 - g) * dt; z z zdot * dt; phidot phidot (u2 / I_yy) * dt; phi phi phidot * dt;欧拉法精度不算高但对这个仿真场景完全够用。如果后面你把模型改复杂了建议至少换成四阶 Runge-Kutta否则大步长下容易出现数值发散让你误以为是控制器的问题。第四步记录数据。把每个周期的高度、姿态、控制量存到数组里仿真结束后一次性绘图。3.3 完整代码复制保存直接跑下面是完整可运行的 Matlab 代码。新建脚本保存为quad_pid_sim.m直接运行即可%% 四旋翼无人机PID控制仿真简化模型高度俯仰 % 使用说明保存为 quad_pid_sim.m直接运行。 % 模型高度通道 z u1 - g俯仰通道 phi u2 / I_yy % 控制律位置式PID 重力补偿 clear; clc; close all; %% 1. 参数初始化 m 1.2; % 无人机质量 kg g 9.8; % 重力加速度 m/s^2 I_yy 0.015; % 俯仰转动惯量 kg*m^2 dt 0.01; % 仿真步长 s控制周期 t_end 20; % 仿真时长 s t 0:dt:t_end; N length(t); % 期望值 z_ref 10; % 期望高度 m phi_ref 0; % 期望俯仰角 rad悬停 % PID参数高度通道 Kp_z 8; Ki_z 0.5; Kd_z 5; % PID参数俯仰通道 Kp_phi 60; Ki_phi 10; Kd_phi 8; % 状态初始化 z 0; zdot 0; phi 0; phidot 0; % 误差积分 err_z_int 0; err_phi_int 0; % 数据存储 z_hist zeros(1, N); phi_hist zeros(1, N); u1_hist zeros(1, N); u2_hist zeros(1, N); %% 2. 主循环 for k 1:N % ---------- 高度通道 ---------- err_z z_ref - z; err_z_int err_z_int err_z * dt; a_z_des Kp_z * err_z Ki_z * err_z_int Kd_z * (-zdot); u1 a_z_des g; % 重力补偿 % ---------- 俯仰通道 ---------- err_phi phi_ref - phi; err_phi_int err_phi_int err_phi * dt; phi_ddot_des Kp_phi * err_phi Ki_phi * err_phi_int Kd_phi * (-phidot); u2 phi_ddot_des * I_yy; % ---------- 模型更新欧拉积分 ---------- zdot zdot (u1 - g) * dt; z z zdot * dt; phidot phidot (u2 / I_yy) * dt; phi phi phidot * dt; % ---------- 记录数据 ---------- z_hist(k) z; phi_hist(k) phi; u1_hist(k) u1; u2_hist(k) u2; end %% 3. 绘图 figure(Name, 四旋翼无人机PID控制仿真结果, Position, [100 100 1000 700]); subplot(2,2,1); plot(t, z_hist, b-, LineWidth, 1.5); hold on; plot(t, z_ref*ones(1,N), r--, LineWidth, 1.2); xlabel(时间 (s)); ylabel(高度 (m)); title(高度响应); legend(实际高度, 期望高度, Location, southeast); grid on; subplot(2,2,2); plot(t, phi_hist*180/pi, b-, LineWidth, 1.5); hold on; plot(t, phi_ref*ones(1,N), r--, LineWidth, 1.2); xlabel(时间 (s)); ylabel(俯仰角 (°)); title(姿态响应); legend(实际俯仰角, 期望俯仰角, Location, southeast); grid on; subplot(2,2,3); plot(t, u1_hist, b-, LineWidth, 1.5); xlabel(时间 (s)); ylabel(等效加速度控制量 u1 (m/s^2)); title(高度通道控制量); grid on; subplot(2,2,4); plot(t, u2_hist, b-, LineWidth, 1.5); xlabel(时间 (s)); ylabel(俯仰力矩 u2 (N·m)); title(俯仰通道控制量); grid on; %% 4. 终端输出 fprintf(仿真完成\n); fprintf(最终高度%.2f m期望 %.2f m\n, z_hist(end), z_ref); fprintf(最终俯仰角%.2f°期望 %.2f°\n, phi_hist(end)*180/pi, phi_ref*180/pi);3.4 运行结果怎么读跑完代码你会看到四张图和一行终端输出。高度响应那幅图蓝色的实际高度应该在 1 到 2 秒内平滑地爬到 10 米几乎没有超调这是因为高度环阻尼比接近 0.88接近临界阻尼。如果曲线爬得很慢多半是 Kp 太小如果冲过头再回来就是 Kd 不够。俯仰角响应那幅图初始可能有微小的角度偏差随后会快速收敛到 0 附近。这里动态要比高度快因为俯仰环的自然频率更高。高度通道控制量 u1起始值约等于 9.8也就是重力加速度。因为无人机要保持 10 米悬停控制量里必须先有重力补偿这部分然后 PID 才在它附近做小幅度调整修正高度误差。你如果看到 u1 远大于 g说明加速度指令里有很大一部分在爬升会对应高度快速上升的阶段。俯仰控制量 u2基本是小幅振荡后回到 0因为稳态时不需要额外的俯仰力矩修正角度。终端会输出最终高度和最终俯仰角。如果跟期望值差得比较多首先看仿真时间是不是不够长再把 Kp、Ki 适当加大。注意改参数之前先想清楚改的是哪一环、预期会发生什么变化再动手。改完一次只动一个参数这是调试的基本纪律。4. 常见问题与排查技巧实录4.1 仿真发散怎么办仿真发散的三个最常见原因按出现频率排序第一个是比例增益过大或微分方向错误。Kp 太大会让系统进入正反馈式的增幅振荡Kd 方向接反会让原本该“刹车”的环节变成“踩油门”。遇到发散先把 Ki、Kd 全设为 0只留一个小 Kp确认系统能稳定下来再逐项恢复。第二个是仿真步长太大。欧拉法本身是近似积分步长越大误差越大。判定方法很简单把 dt 从 0.01 缩小到 0.001再跑一遍如果两条曲线差异明显说明步长不够小需要减小。第三个是模型本身写错了。比如高度更新里忘了减重力或者俯仰力矩符号反了。这种情况需要回头逐行检查模型方程和 PID 计算的符号。判断发散原因有个经验看发散速度。一般 0.5 秒内直接冲到无穷大多数是符号问题缓慢增长然后振荡发散多半是参数问题高频抖动发散优先怀疑步长太小或微分项被噪声放大。4.2 稳态误差、振荡与参数微调稳态误差很典型高度停在 9.7 米或者 10.3 米再也不动。这通常说明积分项不够或者积分被限幅了。我在代码里没有给积分加限幅但真实系统里积分饱和是个大坑建议你加上一个简单的饱和函数比如err_z_int max(min(err_z_int, 10), -10);防止积分量过大导致大超调。振荡问题分两种低频大幅振荡周期比较长多半是 Kp 偏大或者 Kd 偏小。先把 Kd 按 1.5 到 2 倍往上调看超调是否明显减小。高频小幅振荡曲线毛刺多通常是微分项被噪声放大。仿真里没有噪声如果你自己加了传感器噪声就需要对微分项做低通滤波。高度环如果出现等幅振荡可以用我前面提到的公式反推参数。先固定 Kd调小 Kp 让阻尼比回到 0.7 到 1 之间振荡自然消失。下面整理了一个我调参时常用的速查表新手可以直接对照使用现象优先调整调整方向预期效果响应慢爬不到目标Kp增大更快接近目标超调大来回振荡Kd增大抑制超调更快稳定有小幅高频抖动Kd减小或加滤波降低噪声放大存在稳态误差Ki增大消除残余偏差积分导致过冲Ki减小或加限幅降低超调系统直接发散Kp / 符号调小 / 改正恢复稳定收敛但速度太慢Kp 和 Kd 同步同时增大提高响应速度4.3 从仿真到实机还有多远仿真跑通只是第一步。仿真和实机最大的差异在于仿真里没有传感器噪声、没有执行器饱和、没有气动扰动、也没有通信延迟。代码里 u1 没有做饱和限制但真机上电机的油门必然是 0% 到 100%不会出现负油门也不会无限大。所以到实机阶段至少要做三件事第一给 PID 输出加饱和限幅防止控制量超出执行器范围。第二微分项必须加低通滤波不然飞控的陀螺仪噪声会被放大到没法用。第三控制频率要提高真实飞控的内环角速度控制通常在 1kHz 左右外环姿态控制在 500Hz 左右。这也是为什么完整飞控普遍采用级联 PID而不是单级 PID。单级 PID 在简化仿真里够用但真机上的姿态控制通常分成角度环和角速度环内环角速度环响应更快外环角度环的输出是内环的期望角速度。你先把这个简化版跑熟后面再升级到级联 PID思路是连贯的。提示仿真里随便怎么改参数都不会烧东西但真机上调参前一定要保证桨叶周边空旷、飞机固定牢靠安全意识比任何算法都重要。4.4 我的个人调试习惯最后分享几个我自己实践下来很管用的习惯。改参数之前先把期望高度从阶跃改成方波比如 10 米和 5 米来回跳看系统能否来回跟踪。这个测试能同时暴露出响应速度、超调、稳态误差三方面问题比只跑一次阶跃响应信息量大多了。再一个是控制量曲线比状态曲线更能说明问题。有时候高度曲线看起来挺好但控制量在剧烈抖振说明参数其实偏临界了真机上根本飞不出来。我习惯每次仿真后先看 u1 和 u2 的曲线是否平滑再判断参数是否真的合适。还可以试着给自己制造麻烦比如在模型更新那行加入随机扰动或者阵风项看 PID 能不能扛住。这算是最简单的鲁棒性测试做一遍之后你对控制器的理解会深一层。这个代码我后来扩展过好几个版本加了横滚通道、加了积分限幅、加了方波跟踪和扰动测试都是在今天这个基础版上一点点补出来的。你跑通之后建议也按这个路子继续玩下去。
觉得有用,分享给同行:

为您的企业打造数字门面

稳重轻奢商务风格,端正雅致视觉,长效耐看不易过时。

立即咨询 →