简介一套参考模型自适应控制方法的无人机应用MATLAB代码目标读者是自动化、电子信息、数学及计算机等专业的大学生和研究生可用于课程设计、期末大作业或毕业设计也适合刚接触自适应控制的新手。压缩包内共5个文件包含1个主程序文件和4张流程图或结果示意图整体包体仅66KB结构简洁、便于快速查阅。代码采用参数化编程关键控制参数可灵活更改并附带可直接运行的案例数据极大降低了上手门槛同时注释详细逻辑清晰便于读者理解参考模型、自适应律和无人机控制律之间的关系。目前已有129人学习下载。借助这套代码读者可以完成算法建模、仿真运行、参数调整与结果可视化为后续深入无人机控制研究提供可扩展的基础。1. 参考模型自适应方法在无人机控制里到底解决了什么问题无人机控制中标称模型设计并不难难在飞行中不断变化的不确定性载荷变化、电机老化、阵风扰动都会让LQR或PID控制器的固定参数不再可靠。参考模型自适应方法直接把“希望无人机表现出什么动态”定义成一个参考模型再让控制器参数在线调整使真实飞行状态与参考模型的输出误差不断收敛。标题里的Matlab代码包就是把这一套方法做成可以仿真和修改的工程实现。这里从理论公式开始一路写到Matlab代码的模块拆分、参数调整和验证技巧适合读得懂理论但想快速上手复现的工程师。2. 参考模型自适应方法的理论基础与无人机模型搭建想跑通代码先要把几个矩阵的含义弄清楚否则后续调参只能乱试。2.1 从LQR控制到MRAC为什么需要在线修正很多无人机控制代码用LQR作为基线。LQR一般先建立线性状态方程[ \dot{x} A x B u ]然后用lqr函数算出固定的状态反馈增益K控制律是 (u-Kx)。这套方法在模型准确时性能很好但它的前提是A、B不能漂移。真实无人机系统却经常变为 (\dot{x} (A\Delta A)x (B\Delta B)u)其中 (\Delta A) 可能来自转动惯量变化(\Delta B) 可能来自执行机构增益衰减。固定K面对这种变化要么超调变大要么稳态误差逐渐明显。MRAC的做法完全不同。它不去精确辨识 (\Delta A) 和 (\Delta B)而是在控制器中布置可调参数用跟踪误差驱动这些参数收敛到理想值。典型结构包含参考模型、实际对象和自适应律模块作用在无人机中的对应参考模型定义期望动态期望的俯仰角/高度响应实际对象包含不确定性的无人机带风扰和负载变化的模型自适应律在线更新控制增益根据俯仰角误差更新反馈和前馈系数所以从LQR到MRAC改动并不是把一个K换成两个K而是把K从常值变成受误差驱动的时变参数。lqr不再出现在控制律里但可以参考它的最优极点来设计参考模型。打个比方你给无人机试飞员一份固定的操作手册飞机重量变了还是按旧手册控制自然容易出问题。LQR那份手册对应的就是名义A、BMRAC则是让操作手册随飞行状态自动修订。对于只有俯仰轴的测试平台这个差别可能只需要多写十几行代码。2.2 无人机姿态控制的状态空间表达式以四旋翼俯仰通道为例在临近悬停的小角度假设下俯仰角 (\theta) 和俯仰角速度 (\omega) 满足[ \dot{\theta} \omega,\qquad J \dot{\omega} M M_d ]其中 (J) 是俯仰转动惯量(M) 是控制力矩(M_d) 是扰动力矩。如果选用状态 (x[\theta,\omega]^T)、控制输入 (uM)那么标称模型为[ \dot{x} \begin{bmatrix} 0 1 \ 0 0 \end{bmatrix}x \begin{bmatrix}0\1/J\end{bmatrix} u ]这个模型并没有气动阻尼项说明你拿到的基础模型可以非常干净。自适应方法正好可以把未建模的阻尼、陀螺力矩和外部扰动打包成控制参数的变化去补偿。在Matlab里这个模型用很少的代码就能定义% 定义无人机俯仰通道标称模型 J 0.025; % 俯仰转动惯量单位kg·m^2 A [0 1; 0 0]; B [0; 1/J]; C eye(2); D 0; % 参考模型选用自然频率4rad/s、阻尼比0.7的二阶系统 Am [0 1; -16 -5.6]; Bm [0; 16];代码中的Am是参考模型的系统矩阵Bm是参考输入矩阵。这里16恰好等于自然频率的平方5.6等于两倍阻尼比乘以自然频率。建议你先把参考模型写好再用step(ss(Am,Bm,eye(2),0))看一眼阶跃响应确保上升时间和超调量符合预期再继续往下搭自适应律。2.3 参考模型选取规则与自适应律推导参考模型不是随便取的。它的极点决定了无人机在理想情况下的上升时间和超调但也必须匹配电机带宽。若把参考模型调得比物理系统能响应的快自适应误差会长期不收敛。我通常把参考模型自然频率设为目标闭环带宽的1.2到1.5倍阻尼比取0.8左右然后观察阶跃响应曲线确认无持续振荡。如果你的代码包里有lqr设计的闭环极点也可以直接把这些极点作为参考模型的极点这样MRAC至少不会比LQR更慢。设参考模型为[ \dot{x}_m A_m x_m B_m r ]实际系统为 (\dot{x}A x B u)控制律采用[ u \theta_1 x \theta_2 r ]这里的理想参数满足[ A B\theta_1^* A_m,\quad B\theta_2^* B_m ]若实际参数与理想参数不一致跟踪误差 (ex-x_m) 的动态中就会出现参数误差项。为了求更新律取Lyapunov函数[ V e^T P e \gamma_1^{-1}\tilde\theta_1^2 \gamma_2^{-1}\tilde\theta_2^2 ]其中P是Lyapunov方程 (A_m^T P P A_m -Q) 的解(\tilde\theta) 是参数误差。让 (\dot V) 负定可以得到标准更新律[ \dot\theta_1 -\gamma_1 B_m^T P e x^T,\qquad \dot\theta_2 -\gamma_2 B_m^T P e r ]这段推导在Matlab代码里最终只表现为几行矩阵运算% 计算P矩阵 Q eye(2); P lyap(Am, Q); % 自适应律更新e是当前时刻跟踪误差 theta1_dot -gamma1 * Bm * P * e * x(1); theta2_dot -gamma2 * Bm * P * e * r;注意第一个式子里的x(1)是我对俯仰角回路的特殊处理。对于高度控制或全状态控制这里应把x(1)换成完整状态向量x。使用完整状态向量时自适应参数的数量会成倍增加需要在仿真里仔细观察每个参数的收敛趋势。代码包中如果采用完整状态向量逻辑上更通用但调参时也更难判断是哪个状态在主导误差。如果你更习惯用Simulink也可以把参考模型放在连续模块里自适应律用S-Function实现。但纯m脚本的好处是每个中间变量都在工作区探索参数规律比Simulink方便尤其是做批量扫描时。3. 参考模型自适应方法的Matlab代码从RAR包到可运行的仿真闭环拿到RAR包后先不要直接双击运行。解压后应先分辨哪个是主脚本哪个是函数文件再理解数据流的走向。3.1 解压RAR后常见的文件组织方式我见过不少类似的代码包结构往往包含下面几类文件。为了让你快速对照我把常见文件名和它们的作用列成了一张表文件名作用运行时机init_params.m集中设置模型与自适应参数主脚本最前面调用drone_model.m无人机动力学模型含扰动接口每个仿真步调用mrac_controller.m参考模型与自适应律实现每个仿真步调用run_simulation.m主仿真脚本控制循环启动入口plot_results.m画图和计算指标仿真结束后调用实际解压后看到的目录可能长这样mrac_quadrotor/ ├── init_params.m ├── drone_model.m ├── mrac_controller.m ├── run_simulation.m └── plot_results.m你可以用Matlab当前文件夹视图查看。如果只有一个.m文件那一般就是把子函数写在同一个文件里。无论哪种布局重点都是找到循环或ode45所在的主脚本。RAR包里的代码通常需要Matlab搜索路径包含当前文件夹所以解压后最好在Matlab里右键设置当前文件夹或者用addpath把目录加进路径。如果你下载下来的文件被系统存成.txt记得开启“显示文件扩展名”把后缀改成.m后再交给Matlab打开否则编辑器里只看到一串纯文本而无法运行。如果你拿到的代码包里有.mat文件那多半是保存好的仿真结果或系统矩阵。可以在命令行里load xxx.mat看一眼变量列表避免被脚本里的名字误导。有些包还会放一个README说明运行顺序优先读它。3.2 无人机模型函数与控制器函数的代码实现工程上我会把模型和控制器拆开成两个函数。模型函数是真实系统的抽象里面可以故意加入J突变和扰动力矩用来考验控制器。function dotX drone_model(t, x, u, J_actual, tau_dist) % 四旋翼俯仰角动力学 % x(1)是俯仰角x(2)是俯仰角速度 % tau_dist是外部扰动力矩 dotX [x(2); u / J_actual tau_dist / J_actual]; end这里的J_actual对控制器不可见模拟实际飞行中转动惯量变化。控制器函数则负责计算控制力矩并同时更新两个自适应参数。function [u, theta1, theta2] mrac_controller(x, xm, r, theta1, theta2, params) e x - xm; P params.P; Bm params.Bm; % 控制律包含固定阻尼项 u theta1 * x(1) theta2 * r params.Kf * x(2); % 自适应律 theta1_dot -params.gamma1 * Bm * P * e * x(1); theta2_dot -params.gamma2 * Bm * P * e * r; theta1 theta1 theta1_dot * params.dt; theta2 theta2 theta2_dot * params.dt; end这个函数里的params.Kf * x(2)是我主动加的一个阻尼项它不属于教科书上最简MRAC。为什么要加因为自适应参数初始值不可能一次猜准如果初始误差太大纯靠自适应拉回来可能出现首拍控制量过大。固定阻尼项可以先托住系统让自适应算法在较小的误差范围内工作收敛也更平和。等你看到theta1和theta2曲线稳定后可以尝试逐步减小Kf观察系统能否靠纯自适应维持稳定。3.3 主仿真循环固定步长离散化与数据记录虽然ode45精度更高但MRAC的参数更新是时变的用固定步长循环反而更容易把所有中间量保存下来。仿真步长dt取1ms足够覆盖俯仰通道10Hz以内的动态。主循环里每个时刻的工作顺序是计算控制量、更新参考模型、更新真实模型、保存数据。参考模型更新必须使用当前时刻的参考状态否则误差序列会带入一阶滞后。init_params; dt 0.001; t 0:dt:10; N length(t); x zeros(2, N); xm zeros(2, N); u_rec zeros(1, N); theta1_rec zeros(1, N); theta2_rec zeros(1, N); x(:,1) [0.15; 0]; xm(:,1) [0; 0]; theta1 0.2; theta2 0.2; r deg2rad(8); % 参考输入8度 for k 1:N-1 % 扰动切换4秒后增加转动惯量 if t(k) 4 J_actual params.J; tau_dist 0.01 * sin(2*pi*0.8*t(k)); else J_actual params.J * 1.3; tau_dist 0.03 * sin(2*pi*0.8*t(k)); end [u, theta1, theta2] mrac_controller(x(:,k), xm(:,k), r, theta1, theta2, params); u_rec(k) u; theta1_rec(k) theta1; theta2_rec(k) theta2; % 参考模型更新 xm_dot params.Am * xm(:,k) params.Bm * r; xm(:,k1) xm(:,k) dt * xm_dot; % 真实模型更新 dotX drone_model(t(k), x(:,k), u, J_actual, tau_dist); x(:,k1) x(:,k) dt * dotX; end循环内有几个容易被忽略的细节。第一theta1_rec记录的是参数更新后的值也就是当前时刻的控制器参数画图时会看到明显的阶梯上升。第二步长改变后自适应律里的params.dt也应当同步改变不要只在两个地方分别写死。第三主循环中的r可以改成时变信号比如三角波或扫频指令用于测试参考输入变化时的跟踪能力。对于固定步长欧拉法步长不宜超过参考模型最小时间常数的十分之一否则离散误差会明显改变闭环极点。4. 参考模型自适应方法的增益参数怎么调才能让无人机不抖振MRAC的调参比LQR多一层因为没有现成的lqr函数直接给你算K。但反过来看可调参数通常只有两个自适应增益和一个参考模型调整规律还是有迹可循的。4.1 自适应增益gamma的调节策略与常见误区(\gamma_1)和(\gamma_2)分别控制参数(\theta_1)、(\theta_2)的更新速度。如果(\gamma)设得太小参数收敛慢跟踪误差会拖很久如果设得太大参数更新会高频振荡控制力矩也会跟着抖。常见误区是一味调大(\gamma)想加快响应结果却激发了未建模动态。下面的表格列出了我在仿真中常见的(\gamma)档位表现可以当成起点Gamma取值跟踪误差表现参数收敛表现推荐度0.05收敛慢2秒后误差才到0.1°参数缓慢爬坡保守0.50.8秒内误差进入±0.05°参数1秒内稳定首选2.00.3秒内误差减小但伴随轻微抖振参数出现小幅振荡谨慎10.0误差快速收敛但控制量饱和参数极限环振荡不可取我一般会先在init_params.m里把(\gamma_1)设为0.5(\gamma_2)按0.7倍设置跑一遍记录IAE。然后写一个参数扫描脚本用MATLAB优化工具箱里的fmincon或patternsearch自动搜索更优组合。不过要注意优化目标函数时要固定扰动场景否则优化出的参数在另一种扰动下可能表现很差。% 参数扫描脚本 gamma1_set [0.1 0.3 0.5 0.8 1.2]; for i 1:length(gamma1_set) params.gamma1 gamma1_set(i); params.gamma2 params.gamma1 * 0.7; run_simulation; e x(1,:) - xm(1,:); iae(i) trapz(t, abs(e)); end [best_iae, idx] min(iae); fprintf(最优gamma1%.2f, IAE%.4f\n, gamma1_set(idx), best_iae);这里用trapz求误差绝对值积分比简单看终值更能反映整个过程的跟踪质量。注意每次调用run_simulation后工作区变量会更新建议在脚本开头加clear避免上一次结果混进来。参数扫描脚本里run_simulation应当保持干净不要用plot阻塞不要在函数内部输入pause否则扫描到几十组参数时根本停不下来。4.2 参考模型带宽与执行机构饱和的匹配参考模型选择的中心问题是要匹配真实的执行机构带宽。以四旋翼为例电调和电机组合的闭环带宽通常在5到10Hz如果你把参考模型自然频率设为8rad/s基本接近这个范围上限。更高的频率会让控制器频繁请求大幅力矩变化而电机实际输出跟不上误差e长期存在自适应律陷入持续调整。一个比较稳妥的做法是先做参考模型的阶跃响应测试确认其上升时间不小于期望的物理极限。然后在主循环里给控制输入u加饱和模块u_max 0.5; % 最大控制力矩单位N·m u max(-u_max, min(u_max, u));加饱和之后自适应律仍然会对误差进行积分可能因此在饱和区间内参数增长。为了限制这种“饱和积分”给后期带来超调可在控制器里加入抗饱和项例如当u饱和时暂停更新参数if abs(u) u_max theta1_dot 0; theta2_dot 0; end这个小的保护逻辑在多数教科书代码里是没有的但在无人机真实飞行中非常关键。暂停更新不是唯一方法也可以把参数更新律减半但暂停最简单后续调参时也不会给系统引入额外相位滞后。4.3 扰动注入与不确定性建模让结果更可信如果仿真全程没有任何扰动MRAC和固定PID的差别很小。为了体现参考模型自适应方法的优势我习惯在仿真中加入三类扰动阶跃、周期和噪声。下面是组合示例% 阶跃扰动4秒后出现常值力矩 dist_step 0.02 * (cur_t 4); % 周期扰动模拟螺旋桨引起的振动和阵风 dist_sine 0.03 * sin(2*pi*0.8*cur_t) 0.02 * sin(2*pi*1.3*cur_t); % 量测噪声模拟传感器输出 dist_noise 0.005 * randn; tau_total dist_step dist_sine dist_noise;加入噪声后跟踪误差曲线不再平滑这时你才能测试自适应律在噪声环境里是否漂移。如果发现参数估计噪声很大可以在自适应律前面加一阶低通滤波或者减小(\gamma)。但不要完全去掉噪声否则代码在仿真里看起来稳妥真机一上的差距会很大。要留意的是噪声会让(\dot\theta)的跟踪产生有偏的梯度尤其是当噪声幅度与参考模型误差幅度相当时自适应增益增大会放大这个偏差所以在有噪声的代码包里优先降低(\gamma_2)因为(r)通常是常值(\theta_2)更新方程的噪声增益比(\theta_1)更高。5. 用Matlab的误差积分曲线验证参考模型自适应方法的收敛效果代码跑通后验证的重点不是看俯仰角是否跟踪参考输入而是看自适应参数是否有实质变化。如果参数曲线平直说明控制器的表现主要来自固定阻尼项MRAC部分没有真正发挥作用。我会用下面这段代码把轨迹、误差和参数放到同一个图窗里这也是MATLAB画图里最常用的布局figure; subplot(2,2,1); plot(t, rad2deg(x(1,:)), b, t, rad2deg(xm(1,:)), r--); xlabel(时间/s); ylabel(俯仰角/°); legend(真实角,参考模型角,Location,best); grid on; subplot(2,2,2); plot(t, rad2deg(x(1,:)-xm(1,:))); xlabel(时间/s); ylabel(跟踪误差/°); grid on; subplot(2,2,3); plot(t, theta1_rec, t, theta2_rec); xlabel(时间/s); ylabel(自适应参数); legend(theta1,theta2); grid on; subplot(2,2,4); plot(t, u_rec); xlabel(时间/s); ylabel(控制力矩/N·m); grid on;除了看图数值指标更客观。我常用IAE误差绝对值的积分来评估跟踪性能IAE trapz(t, abs(x(1,:) - xm(1,:)));在4秒扰动注入前后分别计算IAE如果差异不大说明自适应律确实在持续补偿不确定性。这个方法还能用来比较不同代码包之间的差异比单纯看轨迹更有说服力。实际项目里我会把每组仿真得到的IAE、(\theta_1)终值和(\theta_2)终值统一导出到表格文件边跑仿真边把新一行追加进去这样多组参数横向对比时不用重新翻命令历史窗口。本文还有配套的精品资源点击获取
简介一套参考模型自适应控制方法的无人机应用MATLAB代码目标读者是自动化、电子信息、数学及计算机等专业的大学生和研究生可用于课程设计、期末大作业或毕业设计也适合刚接触自适应控制的新手。压缩包内共5个文件包含1个主程序文件和4张流程图或结果示意图整体包体仅66KB结构简洁、便于快速查阅。代码采用参数化编程关键控制参数可灵活更改并附带可直接运行的案例数据极大降低了上手门槛同时注释详细逻辑清晰便于读者理解参考模型、自适应律和无人机控制律之间的关系。目前已有129人学习下载。借助这套代码读者可以完成算法建模、仿真运行、参数调整与结果可视化为后续深入无人机控制研究提供可扩展的基础。1. 参考模型自适应方法在无人机控制里到底解决了什么问题无人机控制中标称模型设计并不难难在飞行中不断变化的不确定性载荷变化、电机老化、阵风扰动都会让LQR或PID控制器的固定参数不再可靠。参考模型自适应方法直接把“希望无人机表现出什么动态”定义成一个参考模型再让控制器参数在线调整使真实飞行状态与参考模型的输出误差不断收敛。标题里的Matlab代码包就是把这一套方法做成可以仿真和修改的工程实现。这里从理论公式开始一路写到Matlab代码的模块拆分、参数调整和验证技巧适合读得懂理论但想快速上手复现的工程师。2. 参考模型自适应方法的理论基础与无人机模型搭建想跑通代码先要把几个矩阵的含义弄清楚否则后续调参只能乱试。2.1 从LQR控制到MRAC为什么需要在线修正很多无人机控制代码用LQR作为基线。LQR一般先建立线性状态方程[ \dot{x} A x B u ]然后用lqr函数算出固定的状态反馈增益K控制律是 (u-Kx)。这套方法在模型准确时性能很好但它的前提是A、B不能漂移。真实无人机系统却经常变为 (\dot{x} (A\Delta A)x (B\Delta B)u)其中 (\Delta A) 可能来自转动惯量变化(\Delta B) 可能来自执行机构增益衰减。固定K面对这种变化要么超调变大要么稳态误差逐渐明显。MRAC的做法完全不同。它不去精确辨识 (\Delta A) 和 (\Delta B)而是在控制器中布置可调参数用跟踪误差驱动这些参数收敛到理想值。典型结构包含参考模型、实际对象和自适应律模块作用在无人机中的对应参考模型定义期望动态期望的俯仰角/高度响应实际对象包含不确定性的无人机带风扰和负载变化的模型自适应律在线更新控制增益根据俯仰角误差更新反馈和前馈系数所以从LQR到MRAC改动并不是把一个K换成两个K而是把K从常值变成受误差驱动的时变参数。lqr不再出现在控制律里但可以参考它的最优极点来设计参考模型。打个比方你给无人机试飞员一份固定的操作手册飞机重量变了还是按旧手册控制自然容易出问题。LQR那份手册对应的就是名义A、BMRAC则是让操作手册随飞行状态自动修订。对于只有俯仰轴的测试平台这个差别可能只需要多写十几行代码。2.2 无人机姿态控制的状态空间表达式以四旋翼俯仰通道为例在临近悬停的小角度假设下俯仰角 (\theta) 和俯仰角速度 (\omega) 满足[ \dot{\theta} \omega,\qquad J \dot{\omega} M M_d ]其中 (J) 是俯仰转动惯量(M) 是控制力矩(M_d) 是扰动力矩。如果选用状态 (x[\theta,\omega]^T)、控制输入 (uM)那么标称模型为[ \dot{x} \begin{bmatrix} 0 1 \ 0 0 \end{bmatrix}x \begin{bmatrix}0\1/J\end{bmatrix} u ]这个模型并没有气动阻尼项说明你拿到的基础模型可以非常干净。自适应方法正好可以把未建模的阻尼、陀螺力矩和外部扰动打包成控制参数的变化去补偿。在Matlab里这个模型用很少的代码就能定义% 定义无人机俯仰通道标称模型 J 0.025; % 俯仰转动惯量单位kg·m^2 A [0 1; 0 0]; B [0; 1/J]; C eye(2); D 0; % 参考模型选用自然频率4rad/s、阻尼比0.7的二阶系统 Am [0 1; -16 -5.6]; Bm [0; 16];代码中的Am是参考模型的系统矩阵Bm是参考输入矩阵。这里16恰好等于自然频率的平方5.6等于两倍阻尼比乘以自然频率。建议你先把参考模型写好再用step(ss(Am,Bm,eye(2),0))看一眼阶跃响应确保上升时间和超调量符合预期再继续往下搭自适应律。2.3 参考模型选取规则与自适应律推导参考模型不是随便取的。它的极点决定了无人机在理想情况下的上升时间和超调但也必须匹配电机带宽。若把参考模型调得比物理系统能响应的快自适应误差会长期不收敛。我通常把参考模型自然频率设为目标闭环带宽的1.2到1.5倍阻尼比取0.8左右然后观察阶跃响应曲线确认无持续振荡。如果你的代码包里有lqr设计的闭环极点也可以直接把这些极点作为参考模型的极点这样MRAC至少不会比LQR更慢。设参考模型为[ \dot{x}_m A_m x_m B_m r ]实际系统为 (\dot{x}A x B u)控制律采用[ u \theta_1 x \theta_2 r ]这里的理想参数满足[ A B\theta_1^* A_m,\quad B\theta_2^* B_m ]若实际参数与理想参数不一致跟踪误差 (ex-x_m) 的动态中就会出现参数误差项。为了求更新律取Lyapunov函数[ V e^T P e \gamma_1^{-1}\tilde\theta_1^2 \gamma_2^{-1}\tilde\theta_2^2 ]其中P是Lyapunov方程 (A_m^T P P A_m -Q) 的解(\tilde\theta) 是参数误差。让 (\dot V) 负定可以得到标准更新律[ \dot\theta_1 -\gamma_1 B_m^T P e x^T,\qquad \dot\theta_2 -\gamma_2 B_m^T P e r ]这段推导在Matlab代码里最终只表现为几行矩阵运算% 计算P矩阵 Q eye(2); P lyap(Am, Q); % 自适应律更新e是当前时刻跟踪误差 theta1_dot -gamma1 * Bm * P * e * x(1); theta2_dot -gamma2 * Bm * P * e * r;注意第一个式子里的x(1)是我对俯仰角回路的特殊处理。对于高度控制或全状态控制这里应把x(1)换成完整状态向量x。使用完整状态向量时自适应参数的数量会成倍增加需要在仿真里仔细观察每个参数的收敛趋势。代码包中如果采用完整状态向量逻辑上更通用但调参时也更难判断是哪个状态在主导误差。如果你更习惯用Simulink也可以把参考模型放在连续模块里自适应律用S-Function实现。但纯m脚本的好处是每个中间变量都在工作区探索参数规律比Simulink方便尤其是做批量扫描时。3. 参考模型自适应方法的Matlab代码从RAR包到可运行的仿真闭环拿到RAR包后先不要直接双击运行。解压后应先分辨哪个是主脚本哪个是函数文件再理解数据流的走向。3.1 解压RAR后常见的文件组织方式我见过不少类似的代码包结构往往包含下面几类文件。为了让你快速对照我把常见文件名和它们的作用列成了一张表文件名作用运行时机init_params.m集中设置模型与自适应参数主脚本最前面调用drone_model.m无人机动力学模型含扰动接口每个仿真步调用mrac_controller.m参考模型与自适应律实现每个仿真步调用run_simulation.m主仿真脚本控制循环启动入口plot_results.m画图和计算指标仿真结束后调用实际解压后看到的目录可能长这样mrac_quadrotor/ ├── init_params.m ├── drone_model.m ├── mrac_controller.m ├── run_simulation.m └── plot_results.m你可以用Matlab当前文件夹视图查看。如果只有一个.m文件那一般就是把子函数写在同一个文件里。无论哪种布局重点都是找到循环或ode45所在的主脚本。RAR包里的代码通常需要Matlab搜索路径包含当前文件夹所以解压后最好在Matlab里右键设置当前文件夹或者用addpath把目录加进路径。如果你下载下来的文件被系统存成.txt记得开启“显示文件扩展名”把后缀改成.m后再交给Matlab打开否则编辑器里只看到一串纯文本而无法运行。如果你拿到的代码包里有.mat文件那多半是保存好的仿真结果或系统矩阵。可以在命令行里load xxx.mat看一眼变量列表避免被脚本里的名字误导。有些包还会放一个README说明运行顺序优先读它。3.2 无人机模型函数与控制器函数的代码实现工程上我会把模型和控制器拆开成两个函数。模型函数是真实系统的抽象里面可以故意加入J突变和扰动力矩用来考验控制器。function dotX drone_model(t, x, u, J_actual, tau_dist) % 四旋翼俯仰角动力学 % x(1)是俯仰角x(2)是俯仰角速度 % tau_dist是外部扰动力矩 dotX [x(2); u / J_actual tau_dist / J_actual]; end这里的J_actual对控制器不可见模拟实际飞行中转动惯量变化。控制器函数则负责计算控制力矩并同时更新两个自适应参数。function [u, theta1, theta2] mrac_controller(x, xm, r, theta1, theta2, params) e x - xm; P params.P; Bm params.Bm; % 控制律包含固定阻尼项 u theta1 * x(1) theta2 * r params.Kf * x(2); % 自适应律 theta1_dot -params.gamma1 * Bm * P * e * x(1); theta2_dot -params.gamma2 * Bm * P * e * r; theta1 theta1 theta1_dot * params.dt; theta2 theta2 theta2_dot * params.dt; end这个函数里的params.Kf * x(2)是我主动加的一个阻尼项它不属于教科书上最简MRAC。为什么要加因为自适应参数初始值不可能一次猜准如果初始误差太大纯靠自适应拉回来可能出现首拍控制量过大。固定阻尼项可以先托住系统让自适应算法在较小的误差范围内工作收敛也更平和。等你看到theta1和theta2曲线稳定后可以尝试逐步减小Kf观察系统能否靠纯自适应维持稳定。3.3 主仿真循环固定步长离散化与数据记录虽然ode45精度更高但MRAC的参数更新是时变的用固定步长循环反而更容易把所有中间量保存下来。仿真步长dt取1ms足够覆盖俯仰通道10Hz以内的动态。主循环里每个时刻的工作顺序是计算控制量、更新参考模型、更新真实模型、保存数据。参考模型更新必须使用当前时刻的参考状态否则误差序列会带入一阶滞后。init_params; dt 0.001; t 0:dt:10; N length(t); x zeros(2, N); xm zeros(2, N); u_rec zeros(1, N); theta1_rec zeros(1, N); theta2_rec zeros(1, N); x(:,1) [0.15; 0]; xm(:,1) [0; 0]; theta1 0.2; theta2 0.2; r deg2rad(8); % 参考输入8度 for k 1:N-1 % 扰动切换4秒后增加转动惯量 if t(k) 4 J_actual params.J; tau_dist 0.01 * sin(2*pi*0.8*t(k)); else J_actual params.J * 1.3; tau_dist 0.03 * sin(2*pi*0.8*t(k)); end [u, theta1, theta2] mrac_controller(x(:,k), xm(:,k), r, theta1, theta2, params); u_rec(k) u; theta1_rec(k) theta1; theta2_rec(k) theta2; % 参考模型更新 xm_dot params.Am * xm(:,k) params.Bm * r; xm(:,k1) xm(:,k) dt * xm_dot; % 真实模型更新 dotX drone_model(t(k), x(:,k), u, J_actual, tau_dist); x(:,k1) x(:,k) dt * dotX; end循环内有几个容易被忽略的细节。第一theta1_rec记录的是参数更新后的值也就是当前时刻的控制器参数画图时会看到明显的阶梯上升。第二步长改变后自适应律里的params.dt也应当同步改变不要只在两个地方分别写死。第三主循环中的r可以改成时变信号比如三角波或扫频指令用于测试参考输入变化时的跟踪能力。对于固定步长欧拉法步长不宜超过参考模型最小时间常数的十分之一否则离散误差会明显改变闭环极点。4. 参考模型自适应方法的增益参数怎么调才能让无人机不抖振MRAC的调参比LQR多一层因为没有现成的lqr函数直接给你算K。但反过来看可调参数通常只有两个自适应增益和一个参考模型调整规律还是有迹可循的。4.1 自适应增益gamma的调节策略与常见误区(\gamma_1)和(\gamma_2)分别控制参数(\theta_1)、(\theta_2)的更新速度。如果(\gamma)设得太小参数收敛慢跟踪误差会拖很久如果设得太大参数更新会高频振荡控制力矩也会跟着抖。常见误区是一味调大(\gamma)想加快响应结果却激发了未建模动态。下面的表格列出了我在仿真中常见的(\gamma)档位表现可以当成起点Gamma取值跟踪误差表现参数收敛表现推荐度0.05收敛慢2秒后误差才到0.1°参数缓慢爬坡保守0.50.8秒内误差进入±0.05°参数1秒内稳定首选2.00.3秒内误差减小但伴随轻微抖振参数出现小幅振荡谨慎10.0误差快速收敛但控制量饱和参数极限环振荡不可取我一般会先在init_params.m里把(\gamma_1)设为0.5(\gamma_2)按0.7倍设置跑一遍记录IAE。然后写一个参数扫描脚本用MATLAB优化工具箱里的fmincon或patternsearch自动搜索更优组合。不过要注意优化目标函数时要固定扰动场景否则优化出的参数在另一种扰动下可能表现很差。% 参数扫描脚本 gamma1_set [0.1 0.3 0.5 0.8 1.2]; for i 1:length(gamma1_set) params.gamma1 gamma1_set(i); params.gamma2 params.gamma1 * 0.7; run_simulation; e x(1,:) - xm(1,:); iae(i) trapz(t, abs(e)); end [best_iae, idx] min(iae); fprintf(最优gamma1%.2f, IAE%.4f\n, gamma1_set(idx), best_iae);这里用trapz求误差绝对值积分比简单看终值更能反映整个过程的跟踪质量。注意每次调用run_simulation后工作区变量会更新建议在脚本开头加clear避免上一次结果混进来。参数扫描脚本里run_simulation应当保持干净不要用plot阻塞不要在函数内部输入pause否则扫描到几十组参数时根本停不下来。4.2 参考模型带宽与执行机构饱和的匹配参考模型选择的中心问题是要匹配真实的执行机构带宽。以四旋翼为例电调和电机组合的闭环带宽通常在5到10Hz如果你把参考模型自然频率设为8rad/s基本接近这个范围上限。更高的频率会让控制器频繁请求大幅力矩变化而电机实际输出跟不上误差e长期存在自适应律陷入持续调整。一个比较稳妥的做法是先做参考模型的阶跃响应测试确认其上升时间不小于期望的物理极限。然后在主循环里给控制输入u加饱和模块u_max 0.5; % 最大控制力矩单位N·m u max(-u_max, min(u_max, u));加饱和之后自适应律仍然会对误差进行积分可能因此在饱和区间内参数增长。为了限制这种“饱和积分”给后期带来超调可在控制器里加入抗饱和项例如当u饱和时暂停更新参数if abs(u) u_max theta1_dot 0; theta2_dot 0; end这个小的保护逻辑在多数教科书代码里是没有的但在无人机真实飞行中非常关键。暂停更新不是唯一方法也可以把参数更新律减半但暂停最简单后续调参时也不会给系统引入额外相位滞后。4.3 扰动注入与不确定性建模让结果更可信如果仿真全程没有任何扰动MRAC和固定PID的差别很小。为了体现参考模型自适应方法的优势我习惯在仿真中加入三类扰动阶跃、周期和噪声。下面是组合示例% 阶跃扰动4秒后出现常值力矩 dist_step 0.02 * (cur_t 4); % 周期扰动模拟螺旋桨引起的振动和阵风 dist_sine 0.03 * sin(2*pi*0.8*cur_t) 0.02 * sin(2*pi*1.3*cur_t); % 量测噪声模拟传感器输出 dist_noise 0.005 * randn; tau_total dist_step dist_sine dist_noise;加入噪声后跟踪误差曲线不再平滑这时你才能测试自适应律在噪声环境里是否漂移。如果发现参数估计噪声很大可以在自适应律前面加一阶低通滤波或者减小(\gamma)。但不要完全去掉噪声否则代码在仿真里看起来稳妥真机一上的差距会很大。要留意的是噪声会让(\dot\theta)的跟踪产生有偏的梯度尤其是当噪声幅度与参考模型误差幅度相当时自适应增益增大会放大这个偏差所以在有噪声的代码包里优先降低(\gamma_2)因为(r)通常是常值(\theta_2)更新方程的噪声增益比(\theta_1)更高。5. 用Matlab的误差积分曲线验证参考模型自适应方法的收敛效果代码跑通后验证的重点不是看俯仰角是否跟踪参考输入而是看自适应参数是否有实质变化。如果参数曲线平直说明控制器的表现主要来自固定阻尼项MRAC部分没有真正发挥作用。我会用下面这段代码把轨迹、误差和参数放到同一个图窗里这也是MATLAB画图里最常用的布局figure; subplot(2,2,1); plot(t, rad2deg(x(1,:)), b, t, rad2deg(xm(1,:)), r--); xlabel(时间/s); ylabel(俯仰角/°); legend(真实角,参考模型角,Location,best); grid on; subplot(2,2,2); plot(t, rad2deg(x(1,:)-xm(1,:))); xlabel(时间/s); ylabel(跟踪误差/°); grid on; subplot(2,2,3); plot(t, theta1_rec, t, theta2_rec); xlabel(时间/s); ylabel(自适应参数); legend(theta1,theta2); grid on; subplot(2,2,4); plot(t, u_rec); xlabel(时间/s); ylabel(控制力矩/N·m); grid on;除了看图数值指标更客观。我常用IAE误差绝对值的积分来评估跟踪性能IAE trapz(t, abs(x(1,:) - xm(1,:)));在4秒扰动注入前后分别计算IAE如果差异不大说明自适应律确实在持续补偿不确定性。这个方法还能用来比较不同代码包之间的差异比单纯看轨迹更有说服力。实际项目里我会把每组仿真得到的IAE、(\theta_1)终值和(\theta_2)终值统一导出到表格文件边跑仿真边把新一行追加进去这样多组参数横向对比时不用重新翻命令历史窗口。本文还有配套的精品资源点击获取
