企业官网建设流程全解析

简介：面向控制系统、信号处理与机器学习方向的学习者与工程人员，这份资源提供一套最小二乘类参数辨识的MATLAB实现，用于对系统参数进行在线估计与辨识。包内共1个docx文档，压缩包约78KB，以文字与代码片段形式呈现完整程序思路，便于对照理解算法流程。文档以极点位于单位圆内的稳定系统为例，生成500点M序列输入并叠加白噪声扰动，给出递推最小二乘的参数估计值Theta、估计方差P与增益矩阵K的更新过程，并借助MATLAB绘图展示参数估计曲线；随后改用渐消记忆法，将增益与方差更新语句中的系数调整为0.998，得到更准确的辨识结果。已有2161人学习，适合希望掌握递推最小二乘与渐消记忆最小二乘编程实现、并理解参数估计偏差来源的读者参考。

1. 从一条“参数永远差一点”的曲线说起

如果你在 MATLAB 里跑过系统辨识，大概率见过这种图：四条参数估计曲线一路爬升，最后稳稳停在真值旁边，但就是差那么一点点，怎么加数据都收不拢。我第一次跑递推最小二乘（RLS）时也卡在这儿，盯着Theta矩阵怀疑自己公式抄错了。后来才明白，这不是代码写错，而是算法本身的“记忆”在作祟——老数据一直占着权重，新数据压不过它。

这份最小二乘类的参数辨识 MATLAB 程序，干的就是把这件事讲透：先用一个二阶差分方程造出带白噪声的观测数据，再用 RLS 在线估参数，最后换成渐消记忆法把偏差压下去。它适合正在学系统辨识、控制工程或者要做在线参数估计的开发者，代码不长，但把“为什么偏”和“怎么纠”这两步都摆在了明面上。下面我按自己拆项目的顺序，从数据怎么造、参数怎么递推、偏差怎么来的，一路讲到怎么改成能复用的辨识脚本。

2. 系统模型与数据生成：M 序列和白噪声是怎么造出来的

2.1 被辨识系统的差分方程长什么样

整个程序的起点是一个二阶线性系统，写成差分方程是：

Z(k+2) = 1.8*Z(k+1) - 0.9*Z(k) + 0.75*u(k+1) + 0.5*u(k) + v(k)

这里Z是系统输出，u是输入，v是白噪声扰动。待辨识的参数一共四个：两个自回归系数1.8、-0.9，两个输入增益0.75、0.5。特征方程z^2 - 1.8z + 0.9 = 0的极点模长约为 0.949，落在单位圆内，所以系统稳定——这一点很关键，不稳定系统的输出会发散，递推过程直接炸掉，新手拿到的数据就没法用。

程序里把z(1)设成-1、z(2)设成0，相当于给系统一个非零初始状态，避免前几个采样点全是零导致回归向量退化。这个细节原文没多解释，但实际跑的时候如果初始值全给 0，前几步的h向量接近零向量，增益矩阵K会算出奇怪的值，曲线开头会抖一下。

2.2 M 序列：给系统一个“够激励”的输入

辨识参数不能拿常数输入去喂，否则回归向量线性相关，参数根本解不唯一。程序用 9 级移位寄存器生成 M 序列（伪随机二进制序列），周期是2^9 - 1 = 511，取前 500 点当输入：

% 9 级移位寄存器生成 M 序列 x = [0 1 0 1 1 0 1 1 1]; % 寄存器初值，不能全 0 n = 503; % 多跑几步让序列进入稳定周期 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

这段代码的逻辑是标准的 Fibonacci 型线性反馈移位寄存器：每拍把第 4 位和第 9 位做异或，结果塞回最高位，其余位依次右移，输出取最低位。抽头选(4, 9)对应一个本原多项式，保证周期跑满 511 而不是提前循环。n取 503 比N大 3，是为了让序列先“热身”几拍，避免初值影响前几个输入点。M 序列的好处是它的自相关接近冲激函数，频谱够宽，能把系统的动态特性充分激励出来——这是参数可辨识的前提，换成阶跃或者正弦，四个参数不一定能同时估准。

2.3 白噪声与观测序列的组装

输入有了，扰动用randn(1, N)生成方差为 1 的高斯白噪声。然后按差分方程逐点算出z：

v = randn(1, N); % 方差为 1 的白噪声 z = zeros(N+2, 1); z(1) = -1; z(2) = 0; % 非零初始状态 for i = 3:N+2 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从第 3 点开始递推，对应方程里的k+2；输入M(i-1)、M(i-2)对应u(k+1)、u(k)；噪声v(i-2)对应v(k)。这种“输出比输入多两拍”的错位是差分方程编程里最容易翻车的地方，索引差一位，辨识出来的参数就会整体偏移。我一般会在生成数据后先画一段z和M的对比图，确认输入输出之间的延迟关系符合预期，再往下走。

3. 递推最小二乘核心：一个 for 循环里的四行更新

3.1 回归向量与参数向量的构造

RLS 的本质是把差分方程改写成线性回归形式z(k) = h' * Theta + v(k)，其中回归向量h由历史输出和历史输入拼成：

h = [-z(i-1); -z(i-2); M(i-1); M(i-2)];

负号是因为把1.8*z(i-1)移到等式左边后系数变成-1.8，这样Theta的第一个分量估的就是-1.8，和程序最后画图时参考线画在-1.8对得上。四个分量依次对应a1、a2、b1、b2。这个符号约定一定要和你的模型写法统一，否则估出来的数和真值符号相反，排查半天以为是算法问题，其实是移项时漏了负号。

3.2 增益矩阵、参数更新与协方差递推

核心循环只有四行，但每一行都有讲究：

P = 100*eye(4); % 协方差矩阵初值，表示初始“很不确定” Theta = zeros(4, N+1); Theta(:,1) = [1; 1; 1; 1]; % 参数初值，随便给，靠递推收敛 for i = 3:N+2 h = [-z(i-1); -z(i-2); M(i-1); M(i-2)]; K = P*h*inv(h'*P*h + 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

P初值取100*eye(4)是常见做法，数值大表示对初始参数估计“不信任”，让前几步的修正幅度大一些，收敛更快。Theta初值给[1;1;1;1]纯属随手，RLS 对初值不敏感，只要P够大，几十步内就能拉到真值附近。K的计算里h'*P*h + 1的分母那个1是噪声方差，这里假设噪声方差为 1，和前面randn的设定一致；如果你的实际噪声方差不是 1，这个位置要换成对应的sigma^2，否则增益尺度不对，收敛速度会变。

参数更新那行Theta + K*(z(i) - h'*Theta)里的z(i) - h'*Theta是预测误差，也就是新息。整个 RLS 就是在用新息不断修正旧估计，K决定每次修正多少。协方差更新P = (I - K*h')*P让P逐步缩小，表示随着数据增多，参数估计越来越确定。Pstore把对角线存下来，后面可以画收敛过程，看哪个参数收敛慢。

3.3 结果可视化与偏差的直观呈现

循环结束后画四条参数曲线，再叠上四条真值参考线：

i = 1:N+1; 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('参数估计过程');

跑出来会看到曲线快速上升然后走平，但走平的位置和参考线之间始终留着一道缝。这道缝不是随机误差，而是 RLS 的“记忆”造成的：标准 RLS 等价于对所有历史数据等权最小二乘，早期数据里的噪声被永久保留在估计里，新数据只能按1/k的权重慢慢稀释它，所以估计值收敛到一个有偏的极限。数据越长偏差越小，但永远不会归零。理解这一点，才能明白下一步为什么要改算法。

4. 渐消记忆法改造：两行代码把偏差压下去

4.1 遗忘因子加在哪里

原文给的改法非常直接：把增益公式分母里的1换成0.998，协方差更新里的P除以0.998：

lambda = 0.998; % 遗忘因子 K = P*h*inv(h'*P*h + lambda); % 分母由 1 改为 lambda Theta(:,i-1) = Theta(:,i-2) + K*(z(i) - h'*Theta(:,i-2)); P = (eye(4) - K*h')*P / lambda; % 协方差整体放大 1/lambda

lambda叫遗忘因子，取值在 0 到 1 之间，越接近 1 记忆越长、越接近 0 忘得越快。它的作用是在每一步把旧协方差放大一点，相当于给历史数据乘一个衰减权重，让新数据的话语权更大。0.998对应的时间常数大约是1/(1-0.998) = 500步，和 500 点的数据长度匹配，所以在这个例子里效果刚好。

4.2 为什么这样改就准了

标准 RLS 的协方差更新会让P单调减小到零，增益K也随之趋零，算法“停止学习”，早期噪声造成的偏差就固化下来。除以lambda之后，P每步被放大1/0.998 ≈ 1.002，和更新带来的缩小相互抵消一部分，P稳定在一个非零值上，K不会归零，算法保持持续修正能力。直观地说，标准 RLS 是“越老越信”，渐消记忆是“只信最近一段”。

这里有个容易忽略的点：分母里的lambda和P除的lambda必须是同一个值。有人只改一处，结果增益和协方差不同步，曲线要么震荡要么收敛更慢。我一般把lambda提成变量放在循环外面，两处都引用它，改参数时只动一个地方，避免这种低级错误。

4.3 遗忘因子的取值与验证方法

lambda不是越小越好。太小（比如 0.9）会让算法对噪声过度敏感，参数曲线抖动明显，稳态误差反而变大；太大（比如 0.9999）又退化成标准 RLS，偏差压不下去。工程上的常见做法是先按数据长度估一个量级，再在附近扫几个值对比：

验证方法很朴素：把Theta最后 100 点的均值打印出来，和真值[-1.8; 0.9; 0.75; 0.5]逐项对比，算相对误差。标准 RLS 的相对误差通常在百分之几量级，换成0.998后能压到千分之几。如果换了lambda误差反而变大，先检查是不是数据长度太短——500 点配0.998是合适的，如果只有 100 点还用0.998，等效记忆长度超过数据长度，等于没忘。

5. 避坑与排查：跑不出正确曲线时先看这几处

5.1 曲线发散或出现 NaN

现象：Theta数值越来越大，或者某一步之后全是NaN。原因多半是P初值给得过大又碰上回归向量接近零，inv里的矩阵接近奇异；也可能是系统本身不稳定，输出指数增长把数值撑爆。解决：先确认极点模长小于 1，再把P初值从100*eye(4)降到10*eye(4)试试，同时检查h向量前几步是否全零。如果用的是真实采集数据，先做去均值，别让直流分量把回归矩阵搞病态。

5.2 参数符号和真值相反

现象：估出来的Theta(1)收敛到+1.8而不是-1.8。原因：回归向量里z(i-1)前面的负号漏了，或者模型移项时符号约定和画图参考线不一致。解决：把h的定义和差分方程写在纸上对一遍，确认Theta每个分量对应的到底是原系数还是移项后的系数。这个坑我踩过不止一次，尤其是换系统阶数的时候，符号特别容易乱。

5.3 改了遗忘因子但结果没变化

现象：把1改成0.998后曲线和之前几乎一样。原因：只改了K的分母，忘了改P的除法；或者lambda定义在循环里每次被重置。解决：确认lambda在循环外定义，两处引用同一个变量，改完在循环里加一句disp(lambda)确认值真的传进去了。另外检查是不是把0.998写成了0.998的倒数，方向搞反。

5.4 前几十步曲线剧烈抖动

现象：开头一段参数估计值上下跳，之后才平稳。原因：P初值大、初值Theta又离真值远，前几步增益很大，属于正常现象；但如果抖动幅度大到离谱，可能是噪声方差和分母里的1不匹配。解决：把Theta初值设得更接近真值（比如[0;0;0;0]或直接给真值附近），或者把P初值调小。真实项目里如果前几十点不能丢，可以加一段“预热”不计入结果。

5.5 换数据后完全不收敛

现象：用自己的输入输出数据跑，参数一直乱跳。原因：输入激励不够，回归向量线性相关，参数不可辨识；或者数据里输出和输入的延迟对不上。解决：先算输入的自相关，确认不是常数或单一频率；再画输入输出互相关，找延迟峰值，把h里的索引对齐。M 序列之所以好用，就是因为它天然满足持续激励条件，换成其他输入要额外验证。

6. 二次开发：把脚本改成能复用的辨识函数

原程序是一个从头跑到尾的脚本，参数、数据、算法全揉在一起。真要拿去做项目，第一步是把它拆成函数，让系统阶数、遗忘因子、数据长度都变成入参。我一般会拆成三个文件：gen_data.m负责造 M 序列和观测数据，rls_ident.m负责递推，main.m负责调用和画图。这样换一个系统只需要改gen_data里的差分方程系数，算法部分不用动。

rls_ident.m的骨架大概是这样：

function [Theta, Pstore] = rls_ident(z, u, na, nb, lambda, P0) % z: 输出序列, u: 输入序列 % na: 自回归阶数, nb: 输入阶数 % lambda: 遗忘因子, 1 表示标准 RLS % P0: 协方差初值 n = na + nb; N = length(z) - max(na, nb); Theta = zeros(n, N+1); P = P0 * eye(n); Pstore = zeros(n, N+1); for k = 1:N idx = k + max(na, nb); h = [-z(idx-1:-1:idx-na); u(idx-1:-1:idx-nb)]; K = P*h / (h'*P*h + lambda); Theta(:,k+1) = Theta(:,k) + K*(z(idx) - h'*Theta(:,k)); P = (eye(n) - K*h')*P / lambda; Pstore(:,k+1) = diag(P); end end

几个改造点值得说。第一，阶数na、nb变成参数后，h的拼接用切片z(idx-1:-1:idx-na)自动适配，不用手写四行。第二，增益公式从inv改成左除/，数值上更稳，也更快。第三，lambda传 1 就退化成标准 RLS，一套代码覆盖两种算法，对比实验时不用维护两份。第四，Pstore存对角线而不是整个矩阵，省内存，画收敛曲线够用。

验证改造是否正确的办法：拿原脚本的数据跑一遍新函数，把Theta和原脚本的Theta逐点做差，最大绝对误差应该在1e-10量级，说明只是重构没有改逻辑。然后再换一组阶数（比如三阶系统）测试，确认切片索引没越界。

还有一个实用技巧：把遗忘因子做成随步数变化的。前 100 步用lambda=1让算法快速收敛，之后切到0.998保持跟踪能力，这样开头不抖、稳态又准。实现上就是在循环里加一句if k > 100, lambda = 0.998; end。这个思路在时变系统辨识里很常见，参数缓慢漂移时比固定遗忘因子更稳。

从那以后我每次改辨识脚本，都强制先用已知真值的仿真数据跑一遍，把估计值和真值做差画出来，确认误差曲线收敛到零附近，再往真实数据上套。仿真这关过不了，真实数据只会更玄学。希望这套拆法帮到你，拿到程序后先跑原版看偏差，再改遗忘因子看收敛，最后按函数拆一遍，基本就吃透了。

本文还有配套的精品资源，点击获取