
简介面向控制系统、信号处理与机器学习方向的学习者与工程人员这份资源提供一套最小二乘类参数辨识的MATLAB程序用于对系统参数进行在线估计与辨识。程序以递推最小二乘算法为主线构造极点位于单位圆内的稳定系统生成500点M序列输入并叠加白噪声扰动定义参数估计值Theta、估计方差P与增益矩阵K完成递推更新并借助MATLAB绘图功能呈现参数估计过程同时给出渐消记忆最小二乘的改写方式通过调整遗忘因子相关语句获得更准确的辨识结果。资源包共1个docx文件约78KB以文档形式集中呈现算法说明与完整代码便于对照理解与复现。已有2162人学习适合希望掌握递推最小二乘与渐消记忆法、并借助可视化结果排查参数收敛问题的读者参考。1. 从一组“收敛了但不对”的曲线说起跑过系统辨识的人大概都遇到过这种场面递推最小二乘RLS的曲线明明已经平了四条参数轨迹稳稳地趴在某个值上可你拿真值一对——1.8、-0.9、0.75、0.5——全都差那么一截而且这个偏差不随数据量增加而消失。这不是代码写错了是算法本身的“记忆”在作祟。这份最小二乘类的参数辨识 MATLAB 程序就是围绕这个现象展开的先用递推最小二乘把参数估出来看到偏差再换成渐消记忆最小二乘把老数据的影响压下去曲线才真正贴到真值上。它适合正在学系统辨识、需要一份能直接跑、能改参数、能对比两种算法差异的开发者也适合做控制或信号处理时想快速验证辨识思路的人。整套代码不长但把 M 序列激励、白噪声扰动、递推更新、遗忘因子这几个关键环节都串起来了拿来当二次开发的骨架很合适。2. 系统模型与激励设计为什么先造 M 序列再谈辨识辨识能不能做准一半取决于你喂给系统的输入够不够“有劲”。这份程序的第一步不是写辨识算法而是造激励信号和观测数据这个顺序是对的——没有足够丰富的输入参数根本激发不出来。2.1 被辨识系统的差分方程程序里设定的系统是二阶的Z(k2) 1.8*Z(k1) - 0.9*Z(k) 0.75*u(k1) 0.5*u(k) v(k)这里Z是系统输出u是输入v是白噪声。待辨识的参数就是四个系数a1 -1.8、a2 0.9、b1 0.75、b2 0.5注意程序里 Theta 的排列和符号约定后面会讲。特征方程z² - 1.8z 0.9 0的极点模长约为 0.949落在单位圆内系统稳定。这一点很关键如果系统本身发散递推过程会被数值淹没你看到的曲线爆炸不是算法问题而是模型问题。所以拿到任何一份辨识代码先看它的系统极点这是血泪经验。2.2 M 序列的生成9 级移位寄存器输入信号用的是 M 序列最大长度伪随机二进制序列程序里用 9 个寄存器、按xor(x(4), x(9))做反馈x [0 1 0 1 1 0 1 1 1]; % 9 级寄存器初值不能全 0 n 503; % 多跑 3 拍让序列进入稳定周期 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; % 反馈值填入首位 end逻辑说明这是一个标准的线性反馈移位寄存器。x(4)和x(9)是抽头位置异或结果回填到x(1)其余位依次右移末位x(9)输出。9 级寄存器对应周期2^9 - 1 511程序跑 503 拍、取前 500 点避开了初始暂态。参数上抽头位置决定了序列的周期和自相关特性改抽头就等于换了一个激励信号N决定数据长度500 点对二阶系统够用想看得更细可以加到 2000但递推次数也跟着涨。提示寄存器初值不能全 0否则异或永远是 0序列卡死。这是新手最容易翻车的地方之一。2.3 白噪声与观测序列的合成激励有了再叠上扰动才能模拟真实采集环境v randn(1, N); % 方差为 1 的高斯白噪声 z zeros(N2, 1); z(1) -1; z(2) 0; % 两个初始条件避免从零起步的死区 for i 3:N2 z(i) 1.8*z(i-1) - 0.9*z(i-2) ... 0.75*M(i-1) 0.5*M(i-2) v(i-2); end这里z的长度是N2因为差分方程用到k2前两个点作为初始条件。v(i-2)的下标对齐是为了让噪声和当前输出同步。randn默认均值 0、方差 1如果你想让信噪比更低、辨识更难把v乘个系数就行比如0.5*randn(...)提高信噪比2*randn(...)压低信噪比正好用来测试算法的鲁棒性。观测序列造好之后z和M就是后面递推的全部输入辨识算法不再碰真实参数这一点保证了仿真的“盲辨识”性质。3. 递推最小二乘核心一次参数更新的完整链路数据齐了进入正题。递推最小二乘的思想是每来一组新数据就修正一次参数估计而不是攒够所有数据再一次性求解。这份程序把更新过程写得非常紧凑四行代码就是全部核心。3.1 数据向量与参数向量的构造h [-z(i-1); -z(i-2); M(i-1); M(i-2)]; % 回归向量 Theta(:, i-1) Theta(:, i-2) K*(z(i) - h*Theta(:, i-2));h是回归向量把当前时刻能观测到的量拼成一列前两个是历史输出的负值后两个是历史输入。为什么输出要取负因为差分方程里a1*z(k1)移到等号左边变成-a1*z(k1)所以回归向量里带负号这样Theta估计出来的就是[a1; a2; b1; b2]的正值形式。z(i) - h*Theta(:, i-2)是预测误差新息它衡量的是“用旧参数预测当前输出差了多少”乘上增益K就是这次要修正的量。参数说明Theta是 4×1 列向量初值给[1;1;1;1]给多少不影响最终收敛只影响前几十拍的暂态。3.2 增益矩阵与协方差矩阵的更新K P*h*inv(h*P*h 1); % 增益矩阵 P (eye(4) - K*h)*P; % 协方差矩阵递推P是参数估计的协方差矩阵初值100*eye(4)表示“一开始我对参数一无所知方差给大”。h*P*h 1里的1对应噪声方差这里假设为 1如果噪声方差实际不是 1这个常数应该换成对应的值否则增益尺度会偏。K的维度是 4×1P是 4×4eye(4)保证维度对齐。这两行是整个算法的心脏K决定“信新数据多少”P决定“当前估计有多不确定”。随着递推进行P的元素逐渐衰减增益变小参数趋于稳定——这就是 RLS 收敛的直观解释。3.3 结果存储与可视化Pstore(:, i-1) [P(1,1), P(2,2), P(3,3), P(4,4)]; % 只存对角元 ... 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);Pstore只存对角元因为对角元代表各参数各自的方差够用来观察收敛趋势。作图时把四条估计曲线和四条真值直线画在一起一眼就能看出偏差。这里有个细节真值画的是-1.8和0.9而Theta估计的是a1、a2的正值形式所以图上Theta(1,:)应该收敛到1.8附近、Theta(2,:)收敛到-0.9附近——符号约定要和回归向量对齐否则你会以为算法错了。我第一次跑的时候就因为没对齐符号盯着图看了半天。跑完这段你会看到四条曲线都收敛了但和真值线之间留着一条稳定的缝隙。这不是数值误差是 RLS 的“无限记忆”特性老数据权重不衰减噪声累积导致估计有偏。要消掉这条缝就得请出渐消记忆。4. 渐消记忆最小二乘两行改动带来的精度跃升原文里最精彩的部分就在这里——只改两个语句辨识精度明显提升。改动虽小背后的逻辑值得掰开讲。4.1 遗忘因子的引入位置原始 RLSK P*h*inv(h*P*h 1); P (eye(4) - K*h)*P;渐消记忆版lambda 0.998; % 遗忘因子越接近 1 记忆越长 K P*h*inv(h*P*h lambda); P (eye(4) - K*h)*P / lambda;两处改动分母里的1换成lambdaP更新末尾除以lambda。lambda取 0.998 意味着每步旧信息权重乘以 0.998大约 500 步后权重衰减到0.998^500 ≈ 0.37也就是“半衰期”约 350 步。这个值不是随便定的太小比如 0.9会让算法过度健忘参数曲线抖动大、甚至发散太大比如 0.9999又退化成普通 RLS偏差依旧。0.998 是在 500 点数据量下的折中如果你把N加到 2000可以适当调大到 0.999。4.2 为什么除以 lambda 能纠偏普通 RLS 的P单调递减增益越来越小新数据的话语权越来越弱噪声的累积偏差就固化下来了。渐消记忆版在每步把P放大1/lambda倍相当于人为“注入”一点不确定性让增益不会衰减到零新数据始终有修正能力。这就是它能跟踪时变参数、也能消掉稳态偏差的原因。参数说明lambda是唯一需要调的旋钮建议从 0.995 到 0.999 之间扫一遍看哪条曲线最贴真值且不抖。4.3 两种算法的对比验证把两版结果画在同一张图上最直观。我一般会这样组织对比figure; subplot(2,1,1); plot(1:N1, Theta_rls(1,:), b, 1:N1, 1.8*ones(1,N1), r--); title(普通 RLSa1 估计); grid on; subplot(2,1,2); plot(1:N1, Theta_ff(1,:), b, 1:N1, 1.8*ones(1,N1), r--); title(渐消记忆 RLSa1 估计); grid on;对比时重点看三件事收敛速度多少拍进入真值 ±5% 带、稳态偏差最后 100 拍均值与真值之差、抖动幅度最后 100 拍的标准差。普通 RLS 通常收敛快但稳态有偏渐消记忆版稳态更准但前期可能略抖。这个对比表可以直接作为你二次开发时的回归测试基准。指标普通 RLS渐消记忆 RLSλ0.998进入 ±5% 带所需拍数约 5080约 80120稳态偏差a1明显约 0.050.1接近 0稳态抖动小略大对时变参数跟踪差好5. 避坑与排查跑不通时先看这几处这份代码不长但几个下标和符号的坑足够让新手卡半天。下面是我实际调试时踩过的几条按“现象 → 原因 → 解决”整理。现象一参数曲线一路发散到无穷大。原因多半是P初值给得过大或lambda太小增益爆炸。解决把P 100*eye(4)降到10*eye(4)试试lambda不要低于 0.99同时确认系统极点确实在单位圆内模型本身发散的话算法救不回来。现象二曲线收敛了但四条都收敛到 0 附近。这是回归向量和参数符号没对齐的典型表现。检查h里输出项是否带了负号以及作图时真值线的符号是否和Theta的约定一致。原文里真值画-1.8、0.9而Theta(1)收敛到1.8这个符号差必须心里有数。现象三报错“矩阵维度不一致”。常见于Theta(:, i-1)的下标越界——循环从i3开始但Theta只有N1列i-1最大到N1刚好卡边界。如果你改了N或循环范围记得同步检查Theta、Pstore的预分配列数。现象四渐消记忆版反而比普通版更差。先确认lambda是不是写成了 0.98 甚至 0.9遗忘太快会让估计追着噪声跑。其次检查P更新那行是不是漏了/lambda只改K不改P等于半吊子渐消效果可能更糟。现象五M 序列看起来像常数。寄存器初值全 0或者抽头位置选错导致反馈恒为 0。换一组非全零初值确认xor的两个抽头不同即可。注意调试时先把噪声v置零跑一遍看无噪声下能否精确收敛到真值。无噪声都收不准说明是代码逻辑问题无噪声准、加噪声偏才是算法特性问题。这个二分法能省掉大量瞎猜。6. 二次开发把这份脚本改成能复用的辨识工具这份代码是脚本形态变量全在基础工作区跑一次就完事。真要拿去做项目我建议第一步就是把它函数化把“造数据”和“辨识”拆开这样换系统、换算法、换噪声水平都不用改主逻辑。function [Theta, Pstore] rls_identify(z, u, lambda, P0) % z: 输出序列, u: 输入序列, lambda: 遗忘因子(1 表示普通RLS) % P0: 协方差初值, 返回参数估计轨迹和方差轨迹 N length(z) - 2; Theta zeros(4, N1); Theta(:,1) [1;1;1;1]; P P0 * eye(4); Pstore zeros(4, N1); for i 3:N2 h [-z(i-1); -z(i-2); u(i-1); u(i-2)]; K P*h / (h*P*h lambda); Theta(:,i-1) Theta(:,i-2) K*(z(i) - h*Theta(:,i-2)); P (eye(4) - K*h)*P / lambda; Pstore(:,i-1) diag(P); end end这样封装之后lambda1就是普通 RLSlambda0.998就是渐消记忆版一个函数覆盖两种算法对比实验只需改一个入参。注意我把inv换成了左除/数值上更稳也更快。参数说明P0建议做成可调噪声大的场景给 10100噪声小给 110。改造思路上还有几个方向值得试。一是把遗忘因子做成自适应监测新息z(i) - h*Theta的滑动方差方差大时减小lambda加快跟踪方差小时增大lambda稳住估计。二是把二阶模型推广到任意阶h的长度和Theta的维度由阶数na、nb决定循环里用拼接构造回归向量。三是加一个验证环节用辨识出的参数重建输出和真实输出算拟合度fit 1 - norm(z - z_hat)/norm(z - mean(z))这个指标比肉眼看曲线靠谱得多。验证方法上我习惯做三组对照无噪声、低噪声、高噪声每组都跑普通 RLS 和渐消记忆版记录稳态偏差和拟合度。如果渐消记忆版在高噪声下拟合度反而下降说明lambda调得太小往回加。这套流程跑顺了你手里就有了一份可复现的辨识基线后面换任何新算法都能拿它当标尺。说个自己的教训早些年我调这类递推算法总喜欢一上来就盯着最终曲线看结果被符号约定和下标对齐坑了好几回浪费的时间比写代码还多。从那以后我每次改辨识脚本都强制先跑无噪声、单参数、短数据的冒烟测试确认Theta能精确回到真值再往上叠噪声和渐消记忆。这个习惯帮我省掉了大量“曲线不对但不知道错在哪”的排查。希望这份拆解能帮你少走几步弯路把这份最小二乘参数辨识程序真正跑起来、改下去。本文还有配套的精品资源点击获取