系统辨识与自适应控制MATLAB仿真：程序拆解与工程迁移

简介北航《系统辨识与自适应控制MATLAB仿真》课程全套程序包面向正在学习控制理论、需要完成课程实验或论文仿真的本科生和研究生可作为自适应控制方向的上手参考。压缩包共114个文件仅约1007KB以46个m脚本为主另有pas、dpr、dcu、dfm等Delphi工程文件、exe演示程序以及vbw/vbp工程与jpg图片涵盖课程各章节的程序实现、界面示例和辅助文档。目前已有272人学习、下载常用于课程作业和控制系统仿真研究。程序内容覆盖最小二乘法、正规化最小二乘法等系统辨识方法也包含LQG自适应控制、自校正控制、鲁棒自适应控制等典型算法并配有仿真案例与数据方便读者按章节运行、修改参数并对比结果深入理解参数辨识、控制器设计与不确定性处理等关键环节。通过对这套代码的拆解与复现能明显提升MATLAB/Simulink编程能力也能为机械臂、飞行器、电力系统等实际控制对象建模和控制器调试提供直接借鉴。 总有人问我“北航那套《系统辨识与自适应控制MATLAB仿真》程序到底应该从哪儿开始看”问这个问题的多半是已经把脚本下载下来、双击运行、看到波形弹出来、然后陷入沉默的人。程序能跑和你真看懂了是两码事看懂了和你未来能拿这套方法做自己的仿真研究又是另一码事。这篇文章我想换个视角不讲“程序清单”也不贴整段整段的代码而是把这套程序包当作一个完整的学习系统来拆。我会先讲清楚系统辨识和自适应控制这两块知识是怎么通过MATLAB仿真串成一条线的然后分别拆解辨识程序、控制程序里的核心逻辑和那些教材上不会写的调试细节最后聊聊怎么把课程程序改造成能迁移到你自己研究对象上的工具。适用的人很明确正在学这门课、手头有程序包但觉得“信息过载”的本科生或研究生以及工作中需要做系统建模或自适应控制器设计、想快速上手MATLAB仿真验证的工程师。不保证你看完能背下代码但我能保证你看完再打开那个程序包时知道每一步在干什么、为什么要这么干。1. 拿到程序包先别急着跑先看懂这门课的知识闭环1.1 从系统辨识到自适应控制到底串的是哪条逻辑线很多人把“系统辨识”和“自适应控制”当成两门独立的课在学这是最大的误区。放在一起开课恰恰因为它们是一个闭环的两端。系统辨识解决的是“给我一段输入输出数据我怎么把这个对象的数学模型给反推出来”。这里的模型可以是传递函数、差分方程也可以是状态空间表达式。你可以把它理解成黑箱里有台未知的机器你往里面送信号记录它吐出来的信号然后逆推出这台机器内部大概长什么样。自适应控制解决的是“对象模型不够准、参数还一直在变我怎么让控制器自动跟着调整”。传统PID是出厂前调好参数就固定了但很多实际对象——比如飞行器的气动参数随高度速度变化、电机负载随工况变化——模型根本定不下来。自适应控制的思路就是既然模型不准那我就边运行边辨识用实时辨识出来的模型不断刷新控制器参数。这两者一拼闭环就出来了辨识给控制提供“当前时刻的模型认知”控制用这个认知算控制量系统运行后产生新的数据再喂给辨识算法更新模型。北航这套MATLAB仿真程序本质上就是把这条闭环里的每一个环节都用代码实现了一遍。1.2 程序包的目录结构里藏着课程主线打开程序包第一件事不是点开某个.m文件而是先看目录结构。我拿到的这套程序文件组织方式大致是按知识模块来的常见的是这样以chap或lec开头的文件夹对应教材章节比如最小二乘辨识、递推最小二乘、极大似然、自校正控制、模型参考自适应等每个模块里至少有三种类型的文件脚本文件用来跑通整个流程、函数文件封装算法、Simulink模型用于闭环控制系统仿真根目录一般有一个README或者主说明文档或者教材配套的“程序使用说明”这里有个非常实用的判断技巧先找“带数据生成的那段脚本”。系统辨识仿真一定要先有输入输出数据所以每个辨识模块的脚本开头几乎必然有一段用于生成仿真数据或读取实测数据的代码。你把这个入口找到了整个文件的阅读顺序基本就是从上往下走一遍流程。1.3 用一条仿真主线跑通所有知识点我自己总结了一个“五步主线”几乎所有实验都能套进去生成或导入输入输出数据包括设计激励信号、设置噪声调用辨识算法估计模型参数离线或递推验证辨识出的模型对比拟合效果、检查残差根据模型设计控制器自校正或模型参考自适应闭环仿真并分析控制效果这个顺序建议你先原封不动跑一遍哪怕每个细节都没看懂也要先建立“数据流”的感觉。这和调试程序一样——先跑通再理解最后改造。不要一开始就扎进某个算法的公式里出不来。2. 系统辨识仿真程序拆解怎么让MATLAB自己学出模型参数2.1 激励信号设计M序列和PRBS是辨识的起点辨识算法再漂亮喂给它的数据质量不行全都是白搭。这是最容易忽略但最要命的一点。系统辨识里有个概念叫“持续激励”。简单说输入信号要足够“丰富”各种频率的分量都得有一点模型参数才能被充分激发出来。如果输入是恒定值输出也是恒定值你再厉害的最小二乘也辨识不出动态模型。我见过很多同学仿真结果很差第一反应是改算法其实源头是激励信号设计不合理。课程程序里最常用的激励源是M序列也就是伪随机二进制信号PRBS。在MATLAB里可以直接用系统辨识工具箱的idinput命令生成% 生成一个幅值为1、周期为63、切换钟周期为2的PRBS信号 u idinput(1000, prbs, [0 0.2], [-1 1]);这里有个参数容易把人绕晕带宽那个[0 0.2]表示信号的频率覆盖范围大约是0到0.2倍采样频率。理论上信号切换越快高频成分越多但切换太快会让对象来不及响应实际激励效果反而差。课程程序里一般会给一个折中的切换速率你调试时不要盲目调小钟周期。实操上还有个小细节数据前一小段信号不要参与辨识。因为刚开始系统还没进入稳态初始条件的影响还在直接拿去辨识会给参数估计引入偏差。程序里如果有修剪数据的前处理步骤别嫌多余那是必要的。2.2 离线最小二乘与递推最小二乘代码脉络对比辨识部分最核心的算法就是最小二乘。课程程序里一般会给两套实现一套是离线批量计算一套是递推在线计算。离线最小二乘的代码逻辑非常清楚。假设你要辨识一个二阶差分方程模型y(k) a1·y(k-1) a2·y(k-2) b1·u(k-1) b2·u(k-2) e(k)那么程序的核心就是构造数据矩阵Phi和输出向量Y然后解法方程% 构造回归矩阵第k行对应k时刻的输入输出历史 Phi zeros(N-2, 4); for k 3:N Phi(k-2, :) [-y(k-1), -y(k-2), u(k-1), u(k-2)]; end Y y(3:N); % 最小二乘解theta_hat (Phi*Phi)^(-1) * Phi * Y theta_hat (Phi * Phi) \ (Phi * Y);注意这里用的是反斜杠运算符而不是显式求逆。MATLAB里反斜杠求解线性方程组用的算法更稳定数值上比直接算inv(Phi*Phi)靠谱得多。这个习惯建议你在写任何代码时都保持。递推最小二乘RLS的核心则是一套逐时刻更新的机制。程序里通常会有这样一个循环体% 初始化 theta zeros(4,1); % 参数向量初值 P 1e4 * eye(4); % 协方差矩阵初始给大值 lambda 0.98; % 遗忘因子 for k 3:N phi [-y(k-1); -y(k-2); u(k-1); u(k-2)]; K P * phi / (lambda phi * P * phi); theta theta K * (y(k) - phi * theta); P (P - K * phi * P) / lambda; theta_record(:, k) theta; % 记录参数变化轨迹 end这段代码是整套辨识程序里最值得反复咀嚼的。K是增益向量它决定了当前时刻的预测误差有多大比例会折算进参数修正里。P矩阵是参数估计的不确定度度量初始值取1e4量级的大数表示“一开始我根本不知道参数是多少”。遗忘因子lambda小于1作用是让过去的数据权重逐渐衰减适合跟踪慢时变参数。我调试这套程序时踩过的坑是P矩阵发散。如果你发现参数估计值震荡得厉害甚至变成NaN先检查每一轮循环里P矩阵是不是还保持正定。遗忘因子太小比如0.9以下、激励信号又不够丰富时P会迅速膨胀整个递推就崩了。一般课程实验里用0.95到0.99之间越接近1算法对参数变化的跟踪越慢但对噪声越不敏感。2.3 模型阶次与模型验证光看拟合还不够辨识算法把参数估计出来了不等于工作结束。你还要回答两个问题阶次选得对不对模型到底能不能用阶次辨识在课程程序里一般有专门的小节常见做法是不断升高模型阶次比较AIC准则或者F检验。AIC的核心思想很朴素模型阶次增加会让拟合残差下降但如果下降得不多那多出来的参数就不值当。程序里通常画一条“AIC值随阶次变化”的曲线取拐点处的阶次即可。模型验证的部分MATLAB系统辨识工具箱给了现成命令compare用于把辨识模型的预测输出和实际输出画在一起对比resid用于检查残差是否接近白噪声。一个合格的辨识结果残差的自相关函数应该在置信区间内如果残差还明显成结构性地波动说明模型结构没选对还有动态没抓住。课程程序里通常也会有这两步的绘图代码。我的建议是不要只盯着“拟合率97%”这种数字高兴重点看残差。有时候拟合率很高但残差在某个频率段有明显峰值说明系统里存在未建模的动态这在后续做控制器设计时会成为隐患。3. 自适应控制仿真程序拆解自校正与模型参考两条路线3.1 自校正控制在线辨识加控制器设计的双重循环自适应控制仿真程序里最常见的是自校正控制STR结构。它的核心思想是每个采样周期先递推辨识对象参数再用最新参数重新设计控制器。程序通常是一个大循环循环体里按顺序做四件事采集系统输出、递推更新模型参数、按最新模型计算控制器参数、施加控制量。最小方差控制是自校正里最经典的一个分支。它的目标是让系统输出的方差最小。课程程序里如果你看到代码里算了一个K控制器增益向量然后和参考输入、历史输出做线性组合算出控制量大概率就是最小方差控制器。它的推导起点是把系统输出用模型预测拆成“可控部分”和“不可预测的噪声部分”控制策略直接让可控部分抵消掉让输出只留下噪声。但最小方差控制有个毛病对控制量没有约束算出来的控制幅度可能特别大。课程程序一般在后面会引入广义最小方差控制就是在目标函数里加上控制量的加权项比如J E[y(kd)^2 rho·u(k)^2]。程序里增加的代码不会多就是在控制器增益计算里多了一个rho参数。调试时你会发现rho越大控制量越温和但输出方差也会变大。这个权衡是这类仿真实验最直观的训练点。3.2 模型参考自适应控制MIT律和参数调整的动态过程另一条路线是模型参考自适应控制MRAC。它和自校正的思路完全不同不是在线辨识对象模型而是先定义一个当对象的“参考模型”这个模型的动态特性是你期望的闭环特性。然后通过自适应律实时调整控制器参数让实际系统的输出尽量去跟踪参考模型的输出。课程程序里最典型的自适应律是MIT律。它的写法非常直接——用一个灵敏度的概念把输出误差对控制器参数做梯度下降% 参考模型ym am*ym bm*r ym am * ym bm * r; % 输出误差 e y - ym; % MIT律调整前馈增益kc kc kc - gamma * e * ym * Ts;这行代码里gamma是自适应增益也是整个仿真里最值得调的参数。gamma太小参数调整太慢输出要很久才能跟上参考模型gamma太大系统会震荡甚至发散。北航这套程序里通常会给一个基准值但实际跑的时候你会发现换一个参考模型或者换一组对象参数这个gamma很可能要跟着调。这就是这门课要你亲手体会的东西自适应控制不是“装上就完事”它的收敛性和稳定性强依赖自适应增益与信号幅值的关系。程序里如果涉及李雅普诺夫律你会看到调整律的形式更复杂一些里面会带系统的状态量而不是单纯误差信号。把MIT律和李雅普诺夫律放进同一个程序里对比调参你就能直观感受到MIT律程序简单但稳定性没有保证李雅普诺夫律推导繁琐但给出了稳定的保证。这也是很多学校作业喜欢让同学做对比的原因。3.3 Simulink与m脚本联合仿真时的关键配置自适应控制实验做到闭环仿真程序包里的文件往往就不仅仅是.m脚本了还会出现.slx文件。原因很自然闭环控制系统里包含对象模型、控制器、信号观测多个模块用Simulink搭出来结构更直观。但Simulink仿真有个老坑和课程程序本身没关系却会让程序跑出来结果完全不对——求解器的步长选择。如果你在模型里同时有连续对象比如传递函数和离散控制器比如用Z变换写的数字控制器仿真步长必须远小于离散控制器的采样周期。举个例子控制器采样周期是0.01秒仿真步长如果你用可变步长默认值它可能在某个区间跳得很大导致离散模块的更新时点错乱。我建议把求解器切换成固定步长步长设为采样周期的十分之一到二十分之一离散控制器的采样周期独立设置。另外Simulink里查看数组数据和查看波形不是一回事。很多课程会提到用Scope看波形但如果你想把某个数组信号比如参数估计轨迹、控制量序列导出来做进一步分析直接拖一个To Workspace模块更省事在m脚本里用sim函数跑完仿真后这些数据会存成结构体或者数组方便你后续绘图比较。4. 从“跑通程序”到“做成自己的仿真”调试与迁移经验4.1 程序跑不通时九成是这三类问题先把最常见的报错和异常归纳成一张表你再遇到问题可以直接对照现象最可能的原因排查方向矩阵维度不匹配回归矩阵Phi的列数和你选的模型阶次不一致检查模型结构确认滞后阶数d是否处理正确参数估计值振荡发散遗忘因子太小、激励信号不够丰富、数据含趋势项把lambda调回0.95以上检查输入信号是否有持续激励Simulink仿真结果跳变步长过大离散模块的更新时点错乱改成固定步长步长设为采样周期的1/10以下递推中途出现NaNP矩阵失正定数值溢出检查是否输入输出数据里有异常值考虑加正则化这里我要特别展开一下“数据含趋势项”这个坑。系统辨识里的最小二乘假设噪声是零均值的如果测量数据里叠加了一个直流偏移或缓慢漂移辨识结果里就会多出一个伪静态增益导致零点点位全错。课程程序的原始数据未必有这个处理步骤但你如果要跑真机采集的数据去均值这一步几乎必不可少。4.2 参数越改越乱试试“只动一个”的调试纪律自适应控制程序里参数多得吓人遗忘因子、自适应增益、参考模型极点、控制量加权、白噪声方差……很多人拿到手就忍不住同时调好几个参数结果仿真曲线变好了却根本不知道是谁起的作用。这个习惯一定要改。我的调试纪律是每次只动一个参数固定其他所有条件跑完一组对比曲线后记录结论再动下一个。比如想研究遗忘因子对辨识的影响就把lambda分别设为1、0.98、0.95、0.9其他全不动画四组参数收敛曲线叠加在一起。你会发现lambda接近1时参数轨迹平滑但跟踪慢lambda太小时参数轨迹毛刺多但能跟上突变。这种“参数扫描”的方法能让你很快建立起对算法的直觉比漫无目的地乱改有效一百倍。MATLAB里做参数扫描很方便写个外层for循环每次运行后把关键变量存到cell数组或者多维数组里最后统一画图。课程程序里如果已经有参数扫描代码建议你直接改参数范围跑一遍这是最快理解程序手感的方法。4.3 把课程案例迁移到自己的被控对象上最后说说是怎么把这套程序变成自己的工具的。最常见的问题是“我要辨识的对象不是差分方程是个传递函数该怎么办”其实很简单。任何一个线性传递函数只要给定采样时间MATLAB里可以用c2dm或者c2d命令离散化成差分方程形式。你先把自己的连续模型写成传递函数格式sys_c tf([1 2], [1 3 4]); % 假设你的对象是 (s2)/(s^23s4) sys_d c2d(sys_c, Ts, zoh); % 零阶保持器离散化然后从sys_d里把分子分母系数抠出来作为差分方程的已知真值再用程序里的辨识函数去反推。这样你就有了一个“真值在手”的验证环境能检验你的辨识代码在迁移过程中有没有写错。迁移时最容易忽略的是采样时间的选取。采样太快信号接近连续数据里相邻样本高度相关辨识矩阵会病态采样太慢连续系统的动态细节丢失辨识精度上不去。经验规则是取系统最快时间常数的十分之一到二十分之一。比如你的对象主导时间常数是0.5秒采样时间取0.025秒到0.05秒之间比较合适。另外一个实用的思路如果你需要反复做参数扫描研究可以把课程程序里的主循环封装成一个函数输入为对象参数、噪声方差、算法参数输出为性能指标比如跟踪误差的均方根、参数收敛时间。这样你就能跑几百组蒙特卡洛仿真统计算法的平均表现。这也是从“做作业”跨越到“做研究”最关键的一步。最后分享一个我自己用着很顺手的技巧无论跑什么实验在每个脚本开头都用rng固定随机种子否则每次运行的白噪声序列都不一样你很难判断一条曲线的变化是算法调整带来的还是纯粹随机波动。固定种子之后所有对比实验都在同一组噪声下进行公平性才有保证。这个细节做仿真的人都知道但新手往往要踩过几次“为什么两次运行结果不一样”的困惑后才能真正记住。本文还有配套的精品资源点击获取