递推最小二乘参数辨识MATLAB编程：二阶系统与渐消记忆实现-北京尧图网络

简介这是一篇面向控制系统、信号处理与机器学习等领域学习者的最小二乘参数辨识程序说明文档。资源共包含一个docx格式文档压缩包大小七十八KB内容围绕递推最小二乘算法和渐消记忆最小二乘算法展开先构造一个极点位于单位圆内的稳定系统生成五百点M序列作为输入并叠加白噪声扰动随后通过增益矩阵与估计方差矩阵的递推更新实现系统参数的在线估计与辨识同时借助绘图功能将参数估计过程可视化。文档还针对渐消记忆法给出关键改进通过引入遗忘因子削弱旧数据影响使估计结果更接近真实值其中包含完整的递推公式、参数初值设定、循环迭代流程以及仿真图形输出便于读者直接对照学习或改造成其他辨识对象。已有两千一百五十八人学习下载适合需要掌握最小二乘类参数辨识原理、递推算法实现及仿真验证方法的初学者和相关技术开发人员。1. 最小二乘参数辨识一套能直接改着用的 MATLAB 递推程序做控制系统或信号处理的人大概率都遇到过这种事系统结构大致清楚但里面的系数拿不准。传递函数写不出来状态空间更是无从下手。这时候最直接的办法就是给系统一个已知的激励信号记录输出再用参数辨识把模型里的未知系数一点一点“逼”出来。这篇笔记要拆的是一份操作型很强的 MATLAB 程序——它用最小二乘递推算法RLS辨识一个二阶系统的四个参数然后给出渐消记忆最小二乘的改进版本顺带把 M 序列激励、白噪声扰动、增益矩阵更新这些细节全部摊开来讲。适合刚接触辨识理论但想尽快跑通一个完整例子的从业者也适合那些已经在用工具箱、但想搞明白内部迭代逻辑的人。这程序里辨识的对象是Z(k2)1.8Z(k1)-0.9Z(k)0.75u(k1)0.5u(k)v(k)。四个待辨识参数分别是 -1.8、0.9、0.75、0.5注意程序里 h 向量带负号所以实际迭代目标是 -1.8 和 0.9 这两个带符号的量。极点都在单位圆内系统稳定所以数据是收敛可用的。M 序列做输入、白噪声做扰动这是辨识里最常见的测试组合。下面从算法原理到每一段代码怎么改一步步来过。2. 递推最小二乘的数学内核为什么是 P、K、Theta 三个量在滚2.1 最小二乘的核心思想与递推化的必要性最小二乘参数辨识的基本思路很简单把系统写成 y(k) h(k) * θ v(k) 的形式其中 h(k) 是回归向量由过去的输入输出组成θ 是待辨识参数向量v(k) 是噪声。我们的目标是找一组 θ让残差平方和最小。批量最小二乘可以直接用矩阵求逆算出 θ(HH)^(-1)HZ但问题是如果数据是一批一批到的或者系统参数本身可能缓慢漂移每次都重新做一次矩阵求逆计算量会随数据量增长而爆炸而且老数据会把新数据的信息“稀释”掉。递推最小二乘就是为了解决这个矛盾。它把批处理公式拆成一个迭代式子每来一个新数据就用增益矩阵 K(k) 对当前的参数估计 Theta(k-1) 做一个修正。修正量等于 K(k) 乘以预测误差实际输出减去当前参数预测的输出。这个过程中有三个量在同步滚参数估计值 Theta、误差协方差矩阵 P、增益矩阵 K。P 矩阵本质上记录了参数估计的不确定度K 则是把 P 和当前回归向量 h 组合起来算出的修正步长。三者配合理论上可以从任意初值出发逐步逼近真实参数。这个程序用的系统方程是 Z(k2)1.8Z(k1)-0.9Z(k)0.75u(k1)0.5u(k)v(k)。把它改写成最小二乘的标准形式回归向量 h(k)[-z(k-1); -z(k-2); u(k-1); u(k-2)]程序中循环从 i3 开始所以 h 里取的是 z(i-1)、z(i-2)、M(i-1)、M(i-2)对应前面的 z(k1)、z(k)、u(k1)、u(k)待辨识参数 θ[a1; a2; b1; b2]其中 a1-1.8、a20.9、b10.75、b20.5。注意这里回归向量里的 z 和 u 都带了负号或正号符号问题在后面看代码时要特别留神——很多人自己写的时候就是在这里把符号搞反导致辨识出来的参数符号全错。2.2 RLS 迭代公式与程序中的变量一一对应RLS 的标准迭代公式有三个增益矩阵K(k) P(k-1) * h(k) / (h(k) * P(k-1) * h(k) λ)。λ 是遗忘因子标准 RLS 取 λ1渐消记忆法取 λ1如 0.998。协方差更新P(k) (I - K(k) * h(k)) * P(k-1) / λ。参数更新Theta(k) Theta(k-1) K(k) * (z(k) - h(k) * Theta(k-1))。程序里的P100*eye(4)对应 P 的初始值Theta(:,1)[1;1;1;1]对应参数初值K[1;1;1;1]只是初始化一个列向量在循环里会被真正计算的值覆盖。h[-z(i-1);-z(i-2);M(i-1);M(i-2)]就是当前时刻的回归向量。循环体里三行代码和上面三个公式一一对应。P 的初值取得比较大100 倍单位阵表示一开始对参数估计非常不确定这样 K 的初始修正步长比较大参数能快速从初值 [1;1;1;1]这是一个故意给偏的初值收敛到真值附近。这里有个值得注意的细节程序里 hPh 1 中的 1 是噪声方差假设。如果实际系统噪声方差不是 1这个 1 应该相应调整。但在这个例子里vrandn(1,N) 产生的白噪声方差恰好是 1所以不用改。很多人直接把工具箱里的 RLS 拿过来用根本没意识到这个数字和噪声方差挂钩等换了系统换了噪声辨识结果变差了却找不到原因——这就是其中一条隐蔽的坑。3. 程序逐段拆解M 序列、白噪声、观测数据和递推主循环3.1 M 序列产生移位寄存器与异或反馈的配合程序开头这段是产生 M 序列伪随机二进制序列作为输入激励x[0 1 0 1 1 0 1 1 1]; % 寄存器初值9级 n503; % 脉冲总数 N500; % 实际使用的点数 for i1:n tempxor(x(4),x(9)); % 第4级与第9级异或作为反馈 M(i)x(9); % 输出取第9级 for j9:-1:2 x(j)x(j-1); % 寄存器右移一位 end x(1)temp; % 反馈值写入第1级 end逻辑说明这是一个 9 级线性反馈移位寄存器LFSR反馈多项式取第 4 级和第 9 级的异或。每来一个时钟脉冲寄存器内容右移一位最高位丢弃最低位补上反馈值。输出的 M 序列取的是第 9 级的值。这个序列是周期性的周期为 2^9 - 1 511。程序产生 503 个点取前 500 个使用所以实际用的是一整个周期内的连续片段。参数说明寄存器初值 x 不能全零否则序列永远输出 0这一点要特别注意。n 设为 503 而 N 设为 500多出来的 3 个点是留给后面生成观测数据时用的——观测数据 z 需要从 i3 开始循环所以输入序列必须多准备几个点才能让下标对齐。反馈抽头选择第 4 和第 9 级不是随意的9 级寄存器用这对抽头能保证序列周期达到最大长度 511。如果你要改寄存器级数抽头位置也要跟着变否则周期会缩短激励信号的频谱覆盖就不够均匀。3.2 白噪声与观测数据生成验证模型时的关键路径vrandn(1,N); % 方差为1的白噪声 zzeros(N2,1); z(1)-1; z(2)0; for i3: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 是系统的观测输出序列。z(1)-1、z(2)0 是两个初始状态不是为了凑数而是为了让循环从 i3 开始时有前两个输出值可用。第 i 个观测值由前两个输出、前两个输入和当前噪声共同决定。公式里的系数 1.8、-0.9、0.75、0.5 就是我们要辨识的真值。注意噪声项是 v(i-2)和输出 z(i) 对齐——因为回归向量里的输入输出都是滞后项噪声也要对应滞后两个时刻才合理。参数说明randn(1,N) 产生标准正态分布随机数均值 0 方差 1正好匹配 RLS 公式里分母上的“1”。如果你把噪声方差调大比如改成 2*randn(1,N)那么 RLS 公式里的那个 1 也要相应改成 4方差值否则增益矩阵的步长和噪声水平不匹配参数估计的波动会明显变大。程序里 z 的长度是 N2比 M 序列多两个点也是为了让循环下标从 3 开始后能恰好用到 M 序列的全部 500 个点。3.3 递推主循环P、K、Theta 三个量的更新顺序是成败关键P100*eye(4); % 估计方差初值给大 Pstorezeros(4,N1); Pstore(:,1)[P(1,1),P(2,2),P(3,3),P(4,4)]; Thetazeros(4,N1); % 参数估计值存储 Theta(:,1)[1;1;1;1]; % 故意给偏的参数初值 K[1;1;1;1]; % 增益矩阵初值占位 for i3:N2 h[-z(i-1);-z(i-2);M(i-1);M(i-2)]; % 回归向量 KP*h*inv(h*P*h1); % 增益矩阵 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逻辑说明循环里三步依次是算 K、更新 Theta、更新 P顺序是固定的——必须先用旧的 P 算 K再用 K 修正 Theta最后才更新 P。换顺序会导致 P 和 K 不匹配参数收敛轨迹会乱掉。存储数组 Pstore 和 Theta 的列数都是 N1从 1 到 501但循环体里 Theta(:,i-1) 的下标范围是 2 到 501所以初值放在第 1 列最后一个估计值放在第 501 列。画图时 i1:N1 正好覆盖全部估计点。参数说明P 的初值取 100eye(4)表示初始协方差为 100即参数不确定度很高。这个值如果取太小比如 1参数收敛会明显变慢取太大比如 10^6则前几步的参数估计会剧烈震荡。Theta 初值取 [1;1;1;1] 离真值 [-1.8;0.9;0.75;0.5] 有不小距离这样能直观看到收敛过程。回归向量 h 里的符号很关键因为差分方程是 z(k2)1.8z(k1)-0.9*z(k)...移到等式左边后z 项的系数分别变成 -1.8 和 0.9所以 h 里对应位置是 -z(i-1) 和 -z(i-2)这样辨识出来的 θ(1) 和 θ(2) 会收敛到 -1.8 和 0.9。4. 渐消记忆最小二乘遗忘因子 0.998 背后的原理和改法4.1 为什么标准 RLS 会有稳态偏差标准 RLS 的推导里P 矩阵会随着迭代次数增加而不断缩小。这意味着后期增益矩阵 K 的模越来越小参数的修正步长也随之变小。理论上如果系统是线性定常的、噪声是零均值的参数会收敛到真值附近但实际仿真中你往往会发现估计值和真值之间始终存在一个小的偏差而且这个偏差不随迭代次数增加而消失。原因在于RLS 的代价函数对所有历史数据等权重看待越老的数据积累越多在总代价里占比越大新数据对参数修正的贡献相对被稀释。如果系统本身有轻微时变或者噪声不是纯白噪声比如有色噪声这个偏差会更明显。把系统输出和参数估计代入程序画出的图上你会发现四条参数曲线都收敛了但每条曲线和真值之间都有一条肉眼可见的“缝”。这个缝的大小取决于噪声水平、P 初值和数据长度。P 初值取 100、数据取 500 点时偏差通常较小但不为零——这正是程序注释里写的“辨识结果会与其精确值始终存在着偏差”。4.2 改成渐消记忆法的关键两个语句而已渐消记忆法的思路是给历史数据加一个指数衰减的权重越老的数据权重越小这样新数据能在代价函数中占据更大比例从而持续修正参数。实现上只需要改两处% 原语句 KP*h*inv(h*P*h1); P(eye(4)-K*h)*P; % 改为渐消记忆法遗忘因子 λ0.998 KP*h*inv(h*P*h0.998); P(eye(4)-K*h)*P/0.998;逻辑说明分母上的 1 改成 0.998对应遗忘因子 λ。P 更新时除以 0.998即乘以 1/0.998 ≈ 1.002效果是让 P 矩阵不会像标准 RLS 那样快速收缩到接近零。这样每步结束后 P 维持在一个“有弹性”的量级下一个数据点到来时 K 仍然有足够大的模参数修正步长不会衰减到微乎其微。综合效果就是估计值能持续跟随真值稳态偏差明显缩小。参数说明遗忘因子 λ 的取值范围通常是 0.951.0。λ 越接近 1记忆越长参数越稳定但跟踪越慢λ 越小比如 0.95记忆越短跟踪时变参数的能力越强但稳态波动也越大。0.998 是一个比较保守的值适合参数基本不变但希望消除稳态偏差的场景。如果你的系统参数是缓变的可以考虑取 0.980.99如果参数变化较快可能需要 0.950.97但这时要考虑噪声放大的副作用。注意在这里噪声方差假设从 1 变成了 λ0.998所以分母的含义不只是噪声方差而是噪声方差与遗忘因子的组合——这个细节很多人会忽略。4.3 两种算法的结果对比怎么看对比标准 RLS 和渐消记忆法的参数估计曲线最直观的差异在稳态段。标准 RLS 的曲线虽然收敛但每条曲线和真值之间始终有偏差渐消记忆法的曲线几乎贴近真值线。P 矩阵的对角元变化也解释了这一点标准 RLS 的 P 对角元持续衰减到接近零渐消记忆法的 P 对角元稳定在一个较小的非零值附近。这个差异在你需要高精度参数值时非常关键——比如你辨识出来的参数接下来要用于设计控制器偏差可能导致闭环极点偏出期望位置。5. 避坑指南从初值设定到噪声匹配四个最容易翻车的地方5.1 现象P 初值给太小导致收敛极慢现象把 P100eye(4) 改成 Peye(4) 或 P0.1eye(4)参数估计曲线需要很长时间才能逼近真值500 个数据点可能不够用。原因P 初值代表对参数初始不确定度的估计给小了等于默认“初始参数已经很准”增益矩阵 K 的初始步长就小参数只能慢慢爬过去。解决保持相对较大的初值如 100 倍单位阵或者根据参数可能的变化范围手动估计一个合适的数量级。我一般的做法是先跑一次标准 RLS 看 P 对角元收敛到多少再反推一个合理的初值一般取收敛值的 10100 倍。5.2 现象遗忘因子太小导致参数估计波动大现象把遗忘因子从 0.998 改成 0.95参数曲线虽然在真值附近波动但波动的幅度明显变大看起来“毛刺”很多。原因遗忘因子越小有效数据窗口越短等效于用更少的数据做估计方差自然增大。解决先试 0.990.998 这个区间确认稳态偏差确实被消除后再尝试适当减小。如果偏差还在先检查是不是噪声方差没匹配而不是一味调小遗忘因子。我自己的习惯是只在对参数有明确时变需求时才把遗忘因子调到 0.98 以下否则尽量保持在 0.99 附近。5.3 现象回归向量的符号写反辨识出的参数符号全错现象把 h[-z(i-1);-z(i-2);M(i-1);M(i-2)] 写成 h[z(i-1);z(i-2);M(i-1);M(i-2)]结果 θ(1) 收敛到 1.8 而不是 -1.8θ(2) 收敛到 -0.9 而不是 0.9。原因差分方程移项后系数符号变了但回归向量没有同步调整。解决先把差分方程写成标准回归形式 y(k)h(k)*θv(k)逐项确认符号再写代码。这一步最基础但最容易错而且错的隐蔽——最终参数的绝对值是对的正负号却是反的。5.4 现象M 序列初值全零导致激励信号为零辨识完全失败现象寄存器初值改成 x[0 0 0 0 0 0 0 0 0]程序运行后参数曲线完全不动或者发散到异常值。原因全零初值的 LFSR 永远输出 0系统完全没有被激励观测数据里不含任何参数信息。解决确保寄存器初值至少有一个 1最省事的就是保持程序原来的初值 [0 1 0 1 1 0 1 1 1]。如果你自己重新生成 M 序列记得在代码里加一个判断如果 sum(x)0就报错提示重新初始化。5.5 现象噪声方差与 RLS 公式中的常数不匹配估计值偏离真值现象把噪声改成 v2*randn(1,N)但 RLS 公式里的常数保持 1最后参数估计偏离真值且曲线波动变大。原因公式里的常数同时承担了噪声方差设定和遗忘因子两个角色噪声实际方差和假设值不一致时增益矩阵的步长不是最优的。解决如果你改了噪声方差记得同步修改inv(h*P*h常数)里的常数改成噪声方差的数值。这是最容易被忽略的一个参数因为它藏在矩阵求逆的表达式里不是单独一句话。6. 把辨识结果落地残差校验利用带遗忘因子的 RLS 做在线跟踪当你把渐消记忆法的程序跑通、参数曲线也收敛到真值附近之后还有一个值得做的验证步骤检查模型输出和实际观测之间的残差。做法很简单把辨识出的参数代回系统方程用同一段输入重新仿真一遍然后算残差% 辨识完成后用估计参数计算模型输出 z_est zeros(N2,1); z_est(1) -1; z_est(2) 0; for i3:N2 z_est(i) -Theta(1,end)*z_est(i-1) - Theta(2,end)*z_est(i-2) ... Theta(3,end)*M(i-1) Theta(4,end)*M(i-2); end residual z(3:N2) - z_est(3:N2); % 残差序列 figure(2) subplot(2,1,1) plot(1:N, residual) title(模型残差序列) subplot(2,1,2) histogram(residual, 30) title(残差直方图)逻辑说明残差等于实际输出减模型输出。如果辨识效果好残差应该接近白噪声——均值接近零、直方图呈钟形。如果残差里还有明显的周期性或趋势说明模型结构没对比如阶次不够高或者输入信号没有充分激励出系统的所有模态。参数说明注意这里模型输出用 Theta(1,end)、Theta(2,end) 取的是最后一列的估计值。Theta(1,end) 是辨识出的 -a1因为回归向量里 h(1)-z(i-1)所以计算模型输出时要写成-Theta(1,end)*z_est(i-1)。这个符号细节和前面回归向量里的符号是呼应的别在这里搞反。一个我自己常用的习惯是在线辨识场景里不只要画参数曲线还要同时画 P 矩阵的迹或对角元和因为它能直观反映“当前参数可信度”。如果迹在迭代后期还在持续下降说明系统还没进入稳态如果迹已经稳住但有波动说明遗忘因子在起作用。从那以后我每次跑辨识实验都会强制走一遍“观察 P 迹 → 检查残差 → 确认符号”这三步再把结果用于后续设计。这一步看着多花几分钟实际能挡掉大量后期的低级返工。希望这份程序和拆解思路能帮到你尤其是那些刚开始接触递推辨识、想尽快把理论落到代码上的朋友——按上面的步骤改参数、看曲线、检查残差你很快就能摸清这套算法的手感。本文还有配套的精品资源点击获取

递推最小二乘参数辨识MATLAB编程：二阶系统与渐消记忆实现

相关推荐

requests+BeautifulSoup：爬虫抓取h1标题的入门实战指南

全国计算机二级C语言真题PDF高效刷题指南

OpenGL面试核心考点：渲染管线、矩阵变换与性能排查全解析

更多资讯

IEEE Transactions论文模板全解析：格式规范与投稿避坑指南

H3CNE-MSA存储认证：DAS/NAS/SAN与RAID快照考点

Python词云从出图到答辩：7个案例覆盖完整交付链

CAT6网线真伪鉴别：200MHz PS-ACR≥0dB才是千兆稳定运行硬门槛

GPTQ量化原理：用Cholesky分解破解Transformer权重误差雪崩

OpenMRS Standalone 2.3.1 Java 8u261 部署避坑指南

今日更新

防幻觉智能体架构：电源硬件设计中的大模型编排与工具化实践

OpenGL面试核心考点：渲染管线、矩阵变换与性能排查全解析

全国计算机二级C语言真题PDF高效刷题指南

本周精选

防幻觉智能体架构：电源硬件设计中的大模型编排与工具化实践

OpenGL面试核心考点：渲染管线、矩阵变换与性能排查全解析

全国计算机二级C语言真题PDF高效刷题指南

本月热门

我发现了一个新思路：用 Remotion + Claude Code 像写代码一样自动化生成短视频

Windows下 Codex 中 Chrome 和 Computer Use 插件不可用问题排查及解决参考方式：TaoToken 统一 Key 配置与验证

2026 大模型集体涨价：用 Python 做企业 Token 成本测算与选型避坑（附配置）

读完文章，想聊聊您的网站？