简介面向控制系统、信号处理与机器学习方向的学习者与工程人员这份资源提供一套最小二乘类参数辨识的MATLAB程序用于对系统参数进行在线估计与辨识。程序以递推最小二乘算法为主线构造极点位于单位圆内的稳定系统生成500点M序列输入并叠加白噪声扰动定义参数估计值Theta、估计方差P与增益矩阵K完成递推更新并借助MATLAB绘图功能可视化参数收敛过程随后改用渐消记忆最小二乘算法通过调整增益与方差更新语句获得更准确的辨识结果。资源包共1个docx文件约78KB以文档形式整理程序代码、仿真结果与算法说明便于对照理解递推最小二乘与渐消记忆法的差异。已有2161人学习适合希望掌握参数辨识基本流程、复现仿真实验并理解算法实现细节的读者参考。1. 从一组“收敛了但始终差一点”的曲线说起做系统辨识的人大多遇到过这种场面递推最小二乘跑完五百个点四条参数曲线稳稳地收敛了可把估计值和真值摆在一起看每条都差那么一截不随数据量增加而消失。这不是代码写错了而是标准递推最小二乘在有色噪声或模型阶次不匹配时的固有偏差。这份最小二乘类的参数辨识 MATLAB 程序正是围绕这个痛点展开先用递推最小二乘把参数辨识流程跑通再改用渐消记忆遗忘因子最小二乘把偏差压下去两套结果对照着看偏差从哪来、怎么补一目了然。它适合正在做控制类课程设计、系统辨识毕设或者要给某个被控对象做在线参数估计的读者拿来就能改对象、改阶次、改遗忘因子。2. 需求与模型先把被辨识对象钉死2.1 辨识任务到底要估什么参数辨识的任务不是“拟合一条曲线”而是从输入输出数据里反推出差分方程的系数。这份程序设定的对象是一个二阶差分方程Z(k2) 1.8*Z(k1) - 0.9*Z(k) 0.75*u(k1) 0.5*u(k) v(k)待估参数共四个a1 -1.8、a2 0.9、b1 0.75、b2 0.5。写成回归形式就是 z(k) h·θ v(k)其中回归向量 h [-z(k-1); -z(k-2); u(k-1); u(k-2)]参数向量 θ [a1; a2; b1; b2]。这里有个容易看花眼的地方程序里 Theta 的第一分量对应的是 -1.8 而不是 1.8因为回归向量里 z 项带了负号画图时参考线也画在 -1.8 和 0.9 上对照时别把符号搞反。特征多项式 z² - 1.8z 0.9 的极点模长为 √0.9 ≈ 0.949落在单位圆内系统稳定。这一点很关键只有稳定系统输出才不会发散递推过程才有意义。如果换成极点跑到圆外的对象输出指数增长协方差矩阵 P 很快数值爆炸辨识曲线会直接飞出去那不是算法问题是对象选错了。2.2 输入信号为什么选 M 序列辨识的激励信号必须“足够丰富”能把系统的动态模态都激发出来。程序用 9 级移位寄存器生成 500 点 M 序列伪随机二进制序列初值x[0 1 0 1 1 0 1 1 1]反馈抽头取第 4 位和第 9 位异或x [0 1 0 1 1 0 1 1 1]; % 9 级寄存器初值不能全零 n 503; % 多跑 3 拍保证取到 500 个有效点 N 500; for i 1:n temp xor(x(4), x(9)); % 反馈抽头4 与 9 异或 M(i) x(9); % 从末级输出 for j 9:-1:2 x(j) x(j-1); % 寄存器右移 end x(1) temp; endM 序列的周期是 2⁹ - 1 511取 500 点基本覆盖一个完整周期自相关接近冲激频谱平坦对二阶对象来说激励充分。这里有个血泪经验寄存器初值不能全零全零状态下异或永远输出零序列直接死掉另外抽头位置决定了序列的周期和统计特性换抽头等于换了一套激励辨识结果的可比性就没了做对照实验时别随手改。2.3 噪声与观测序列的生成白噪声用v randn(1,N)生成方差为 1。观测序列递推生成时要注意索引对齐z 预分配 N2 个元素前两个给初值z(1)-1, z(2)0从 i3 开始按差分方程算输入取 M(i-1) 和 M(i-2)噪声取 v(i-2)。这种“多留两拍”的写法是为了让回归向量 h 在 i3 时就能凑齐 z(i-1)、z(i-2)、u(i-1)、u(i-2) 四项避免开头几拍索引越界。很多同学在这里翻车把 z 的长度设成 N循环里一取 z(i-1) 就报错或者错位一拍导致辨识出的 b1、b2 互换。3. 递推最小二乘的实现从公式到可跑代码3.1 三个核心矩阵的初始化递推最小二乘每拍做三件事算增益、更新参数、更新协方差。初始化直接决定前几十拍的收敛形态P 100*eye(4); % 协方差初值大数表示“对初值极不信任” Pstore zeros(4, N1); % 记录 P 对角线观察收敛 Pstore(:,1) [P(1,1), P(2,2), P(3,3), P(4,4)]; Theta zeros(4, N1); % 参数估计历史 Theta(:,1) [1; 1; 1; 1]; % 参数初值故意给偏离真值的值 K [1; 1; 1; 1]; % 增益初值第一拍会被覆盖P 取 100·I 是常规做法数值越大表示初始不确定性越大算法前期步子迈得大、收敛快取太小则前期几乎不动看起来像“算法没生效”。Theta 初值故意设成 [1;1;1;1]和真值 [-1.8;0.9;0.75;0.5] 差得很远正好用来观察收敛过程。如果初值直接给真值曲线一条直线什么信息都看不出来答辩时也讲不出东西。3.2 递推主循环与增益更新主循环从 i3 跑到 N2每拍构造回归向量、算增益、更新参数和协方差for i 3:N2 h [-z(i-1); -z(i-2); M(i-1); M(i-2)]; % 回归向量z 项带负号 K P*h*inv(h*P*h 1); % 增益分母的 1 是噪声方差 Theta(:,i-1) Theta(:,i-2) K*(z(i) - h*Theta(:,i-2)); % 参数更新 P (eye(4) - K*h)*P; % 协方差更新 Pstore(:,i-1) [P(1,1), P(2,2), P(3,3), P(4,4)]; end分母里的1对应噪声方差 σ²1如果噪声方差改成 0.1这里要同步改成 0.1否则增益尺度不对收敛速度会明显偏慢。inv(h*P*h1)在四维情况下是标量求逆开销可以忽略但如果把阶次提到十阶以上建议改成左除K P*h/(h*P*h1)数值更稳。参数更新项z(i) - h*Theta(:,i-2)是新息也就是当前观测与预测之差新息越大说明模型越不准增益就把它按比例分配到四个参数上。3.3 结果可视化与偏差的读法画图部分把四条参数曲线和四条真值参考线叠在一起i 1:N1; figure(1) plot(i, Theta(1,:), i, Theta(2,:), i, Theta(3,:), i, Theta(4,:), ... i, -1.8, i, 0.9, i, 0.75, i, 0.5) grid on; title(参数估计过程);跑完会看到曲线在前 50 拍快速收敛之后基本走平但四条都停在离参考线一小段距离的位置而且这段距离不随 N 增大而缩小。这就是标准递推最小二乘的偏差问题——当噪声与回归向量相关有色噪声、或模型阶次偏低时最小二乘估计不再是无偏的。读图时重点看两件事收敛速度前段斜率和稳态偏差后段与参考线的间距前者反映 P 初值和增益设计后者反映算法本身的统计性质。4. 渐消记忆最小二乘两行改动换来的精度4.1 遗忘因子改的是哪两个式子原文给出的改法非常克制只动两行% 原式 K P*h*inv(h*P*h 1); P (eye(4) - K*h)*P; % 渐消记忆版 K P*h*inv(h*P*h 0.998); P (eye(4) - K*h)*P/0.998;第一处把分母的 1 换成 0.998第二处在协方差更新后除以 0.998。0.998 就是遗忘因子 λ取值在 0 到 1 之间越接近 1 记忆越长、越接近标准递推越小则越“健忘”对近期数据加权更重。这里 λ0.998 对应的有效记忆长度约为 1/(1-λ) 500 拍和 500 点数据量刚好匹配所以稳态偏差被明显压下去。4.2 遗忘因子为什么能消偏差标准递推最小二乘等价于对全部历史数据做等权最小二乘旧数据和新数据权重一样。当噪声有色时旧数据里的噪声分量会持续“拖住”参数估计形成系统性偏差。渐消记忆相当于给每拍数据乘上一个指数衰减权重越老的数据权重越低参数估计跟着系统当前的统计特性走偏差自然被削弱。代价是估计方差会略微变大曲线尾部抖动比标准版明显一点这是精度和稳定性之间的正常交换不是程序出错。4.3 参数怎么调、边界在哪λ 的选取没有万能值通常的做法是按数据长度和系统时变程度折中数据 500 点、系统定常λ 取 0.995 到 0.999 都合理系统参数缓慢时变λ 要往小取比如 0.98让算法跟得上变化λ 取到 0.99 以下时有效记忆只剩一百拍左右曲线抖动会明显加剧答辩时容易被问“你这估计方差是不是太大了”。另外 λ 不能取 1取 1 就退化成标准递推改了等于没改也不能取太小太小相当于只用最近几拍数据协方差矩阵 P 会周期性膨胀数值上容易出问题。我一般会先固定 λ0.998 跑一遍看稳态偏差再上下各调一档做对比把三条曲线放进同一张图里偏差和抖动的权衡关系就直观了。5. 避坑与排查这类辨识程序最容易翻车的五处5.1 曲线收敛但四条参数整体错位现象是曲线形状正常、收敛也快但估出来的 b1、b2 和真值对不上甚至互换。原因多半是观测序列生成时输入索引错了一拍M(i-1)和M(i-2)写反或者 z 的循环起点不是 i3。解决方法是拿前五拍数据手工代入差分方程验算一遍确认 z(3) 确实由 z(2)、z(1)、M(2)、M(1) 算出索引对齐后再跑。5.2 前几拍参数剧烈跳变甚至出现 NaN现象是曲线开头几十拍上下乱窜偶尔冒出 Inf 或 NaN。原因通常是 P 初值给得过大比如 1e6叠加回归向量元素量级差异大h*P*h数值条件变差求逆失稳。解决办法是把 P 初值降到 100 量级或者把inv换成左除必要时对输入输出做归一化让回归向量各分量量级接近。5.3 改了遗忘因子结果反而更差现象是把 λ 从 0.998 调到 0.95 后稳态偏差没降、抖动却大幅增加。原因是 λ 太小导致有效记忆长度远小于系统收敛所需的数据量参数还没稳住旧信息就被丢光了。解决方法是按 1/(1-λ) 估算有效记忆长度保证它至少是系统时间常数的若干倍先试 0.995 再逐步下调每次只动一档并保留对照图。5.4 换一个被辨识对象就发散现象是原对象跑得好好的换成三阶或极点靠近单位圆的系统后输出发散、P 矩阵数值爆炸。原因是对象本身不稳定或接近临界稳定开环输出指数增长任何辨识算法都救不回来。解决办法是先检查特征多项式根的模长确认全部落在单位圆内且留有余量必要时先做闭环采集数据再辨识。5.5 噪声方差改了但增益没跟着改现象是把randn乘上 0.1 把噪声方差改成 0.01 后收敛速度肉眼可见地变慢。原因是增益公式分母里的1仍按方差为 1 写尺度和实际噪声不匹配。解决办法是把分母常数同步改成噪声方差或者干脆把它写成一个变量sigma2在程序开头统一定义改噪声时只改一处。6. 验收与答辩用三张图把方案讲圆验收这类辨识程序不看代码行数看三样东西能不能自圆其说。第一张图是标准递推最小二乘的收敛曲线要能指出前段收敛速度和后段稳态偏差的具体数值比如“第 50 拍基本收敛稳态偏差约 0.05”第二张图是渐消记忆版的对照曲线把两版画在同一坐标系里偏差缩小了多少、尾部抖动增加了多少用数字说话第三张图是协方差对角线 Pstore 的衰减曲线它证明算法确实在收敛而不是靠初值撑着。三张图摆出来从“有问题”到“怎么改”再到“改的代价”逻辑链就完整了。答辩时最容易被追问的是参数选择依据提前把下面这张表准备好被问到直接报数参数取值选取依据改动影响P 初值100·I表示初始不确定性大加快前期收敛过大易数值失稳过小前期不动θ 初值[1;1;1;1]故意偏离真值便于观察收敛给真值则曲线平直失去演示意义遗忘因子 λ0.998有效记忆约 500 拍匹配数据长度减小则跟踪快但抖动大数据长度 N500覆盖 M 序列一个完整周期过短则统计特性不足噪声方差1与增益分母常数一致改动需同步改分母还有一个常被忽略的细节答辩演示时别只放最终收敛图把前 20 拍的局部放大图也备一张。评审最爱问“你怎么证明它是在辨识而不是在拟合初值”局部放大图里参数从初值一路爬向真值的过程比任何解释都有说服力。我早年第一次讲这个项目只放了稳态图被追问“曲线平了是不是因为算法停了”当场卡壳从那以后我每次做辨识类演示都强制把前段收敛过程、稳态偏差数值、协方差衰减三样一起放再留一张参数对照表兜底。希望这套结构能帮到你把手里那份最小二乘辨识程序讲成一份站得住的完整方案。本文还有配套的精品资源点击获取