递推最小二乘参数辨识MATLAB程序详解与渐消记忆改进

简介最小二乘类参数辨识MATLAB程序包面向自动化、控制理论与信号处理方向的学习者与科研人员解决系统参数在线估计与辨识的实现问题。程序基于递推最小二乘与渐消记忆最小二乘算法配合500点M序列输入和白噪声扰动在稳定系统模型上演示了参数初值设定、增益矩阵K与方差P的递推更新及结果展示。程序包为1个docx文档容量78KB含完整代码、关键注释和辨识结果分析并解释了渐消因子0.998下K与P公式的调整方式便于对比两种算法的精度与收敛差异。文档还介绍白噪声、增益矩阵、估计方差等概念提供参数初值选取经验便于理解数学基础。已有2158人学习使用适合掌握最小二乘辨识原理并开展MATLAB编程实践也可用于课程设计或论文仿真。1. 最小二乘参数辨识这份MATLAB程序能帮你省掉哪些推导做控制、信号处理的人迟早会碰到参数辨识系统模型有了但里面的系数不知道要用输入输出数据把它估出来。我最初自己推RLS公式时被增益矩阵和协方差更新绕得晕头转向直到拿到一份能直接跑的最小二乘类MATLAB程序才真正理解递推过程。这份程序用500点M序列作为激励对二阶系统做递推最小二乘并演示了渐消记忆法修正稳态偏差。它适合刚刚接触辨识理论、想用MATLAB验证算法、或者在课程设计中需要完整代码的同学。下面我把这段代码拆开讲清楚包括参数怎么设、哪里容易踩坑。2. 从系统模型到递推最小二乘为什么用M序列和500个数据点2.1 系统方程与辨识目标先把对象写清楚程序里设计的系统是 Z(k2)1.8Z(k1)-0.9Z(k)0.75u(k1)0.5u(k)v(k)其中v(k)是白噪声。这个方程把未来输出表示成过去两拍输出、过去两拍输入和噪声的组合。实际生成观测数据时程序用的是z(i)1.8z(i-1)-0.9z(i-2)0.75M(i-1)0.5M(i-2)v(i-2)两者只是下标平移了两个周期本质相同。把方程改写成回归形式z(k) [-z(k-1); -z(k-2); u(k-1); u(k-2)] * theta v(k-2)其中theta [-1.8; 0.9; 0.75; 0.5]。注意h向量里带负号待估参数是-1.8而不是1.8。程序最后画图时拿Theta(1,:)和-1.8比别搞反。系统稳定性可以验证特征方程是λ^2 -1.8λ 0.90判别式为3.24-3.6-0.36两个共轭复根模长sqrt(0.9)≈0.9487位于单位圆内所以系统稳定。稳定是递推辨识能收敛的前提之一。如果极点靠近单位圆甚至在外输出会发散或者需要很长的数据才能保证持续激励。这里极点模长0.9487已经比较接近1实际上收敛速度不算快早期瞬态会影响初期的参数估计所以用500点数据是有道理的。参数辨识的目标就是用输入M(i)和输出z(i)序列把theta里的四个数估出来。程序里没有直接给出真值而是通过初始Theta[1;1;1;1]和递推公式逐步逼近。这个“从错误初值逼近真值”的过程恰好就是RLS最直观的教学演示。如果你把初值改成[0;0;0;0]再跑一遍会发现收敛过程相似只是路径不同。2.2 M序列与白噪声激励信号怎么选系统辨识需要输入信号能激发系统的所有动态模态。单位阶跃或单一正弦往往只能激励部分频率参数估计会不稳定。程序用的是M序列即伪随机二进制序列。M序列的频谱接近白噪声且幅值恒定属于持续激励信号能保证参数可辨识。生成代码如下x[0 1 0 1 1 0 1 1 1]; %寄存器初值不能全零 n503; %总共产生503个点 N500; %实际使用前500点 for i1:n tempxor(x(4),x(9)); M(i)x(9); %输出M序列 for j9:-1:2 x(j)x(j-1); end x(1)temp; endx是9位移位寄存器反馈抽头在第4位和第9位异或后送回第1位。这个配置对应的周期是2^9-1511取前500点。注意循环产生了503个点比N多3个这是给后续生成观测序列时留的余量避免下标越界。白噪声用randn(1,N)生成均值为0、方差为1的标准正态随机数恰好对应标准RLS递推公式分母里的那个“1”。如果实际噪声方差不是1比如你把randn的结果乘以sqrt(2)那K公式里的1就要改成2否则增益计算不匹配。很多初学者直接照搬公式换了噪声却不改分母结果收敛慢或偏差大这是第一个容易踩的坑。为什么选500点二阶系统有4个参数理论上每个参数至少需要几个数据点的信息量但考虑到噪声方差为1、系统极点模长0.9487接近单位圆早期瞬态需要一段时间衰减。500点覆盖了M序列的一个完整周期既保证激励充分又不会让计算量过大。如果换成100点你会看到参数曲线还没稳定就结束了。2.3 递推最小二乘的四个核心变量Theta、P、K、h标准RLS代码段如下P100*eye(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四个核心变量说清楚Theta是参数估计向量这里用4×(N1)矩阵存每一拍的结果。Theta(:,i-1)就是第i个数据点对应的估计值初值故意设成[1;1;1;1]离真值[-1.8;0.9;0.75;0.5]很远这样才能看出收敛过程。P是误差协方差矩阵初值100*eye(4)表示对初始估值非常不确信。P越大第一拍增益就越大修正幅度也大。随着递推进行P会逐渐变小增益也随之变小。K是增益向量由P和h计算得到决定当前新息z(i)-h*Theta的权重。新息越大修正越多。h是回归向量由上一拍输出、上两拍输出、上一拍输入、上两拍输入组成作用是把当前数据映射到参数空间。递推逻辑只有一句话用预测误差乘以增益修正参数再更新P。注意K公式里的inv在小维度里没问题但换成除法KP*h/(hPh1)会更稳。后面的避坑章节我会细说。Pstore单独存对角线是为了观察协方差收敛情况。如果Pstore下降很快并趋于平稳说明数据信息在积累如果一直震荡多半是输入激励不足或模型阶次不对。另外Pstore定义成4×(N1)但循环里实际只填到第N列最后一列保持0不影响画图。3. 跑通标准RLS程序Matlab代码逐行拆解3.1 生成M序列的移位寄存器实现M序列生成代码在上一章已经出现但有几个细节值得单独说。寄存器初值x不能全零。如果全部是0异或结果永远是0移位后序列也全是0输入就没有任何激励。这里用[0 1 0 1 1 0 1 1 1]非零初值且包含足够的0/1混合。抽头位置选择第4和第9这是基于本原多项式。抽头不同序列周期可能会缩短。如果周期缩短伪随机性变差频谱会出现较强的离散线谱激励不够平均。我一般会在仿真前先检查M序列的周期把M序列重复两次用find(diff(M)1)之类的办法找重复但最简单的是直接看前100点是否有明显规律。如果看到长串连续0或连续1说明寄存器配置不好。M序列的幅值这里是0和1相当于有直流分量。严格来说持续激励信号要求频谱丰富带直流也可以但直流分量会让输出z也带上偏置在有色噪声或非零均值环境下会影响参数估计。后面渐消记忆章节你会看到原程序最终结果仍有一点偏差输入有直流也是原因之一。把0改成-1可以消除直流这是常见做法。3.2 生成观测序列z(k)的递推式观测数据生成段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); endz(1)和z(2)是初始状态分别设为-1和0。由于系统稳定初始状态的影响会随时间衰减但前几个数据点仍会受影响。这相当于真实系统在开始采集前已有未知初始状态辨识算法需要靠后期数据慢慢纠正。如果你希望避免初始瞬态干扰可以先让系统空跑50拍再从第51拍开始记录。不过在这个例子里500点足够长末尾的收敛结果基本不受影响。这里有一个容易被忽视的下标映射z数组从1开始M数组从1开始。当i3时用M(2)和M(1)对应u(k-1)和u(k-2)。也就是说第1个可用输出是z(3)它依赖前两拍输出和M的前两个点。因此z的长度是N2多出两个元素给初值。3.3 递推循环增益矩阵与协方差更新把递推循环单独拿出来加一个残差变量便于调试for i3:N2 h[-z(i-1); -z(i-2); M(i-1); M(i-2)]; KP*h*inv(h*P*h1); e z(i) - h*Theta(:,i-2); Theta(:,i-1)Theta(:,i-2)K*e; P(eye(4)-K*h)*P; Pstore(:,i-1)[P(1,1),P(2,2),P(3,3),P(4,4)]; end这段代码是RLS的核心。hTheta(:,i-2)就是用上一拍参数对当前输出做预测e是预测误差。Ke就是对参数的修正量。P更新是标准公式没有遗忘因子所以P会随着数据增加不断减小。如果你打印Pstore的最后一列会发现P对角线已经非常小比如1e-3以下。这意味着参数估计的置信度很高但也意味着后续新数据的修正能力很弱。如果此时发现最终估计偏了想再拉回来几乎不可能。这正是第4章要解决的问题。调试时建议在循环里加一个计数器隔50拍打印一次Theta的当前值观察是否单调靠近真值。如果中途出现跳跃可能是h向量构造错了比如下标偏移一位。我在做模拟项目X时就曾把M(i-2)误写成M(i-1)结果前200拍参数完全乱跳后来逐行对公式才找到问题。3.4 画图与收敛性判断程序最后画图i1:N1; figure(1); plot(i,Theta(1,:),i,Theta(2,:),i,Theta(3,:),i,Theta(4,:),... i,-1.8*ones(1,N1),i,0.9*ones(1,N1),i,0.75*ones(1,N1),i,0.5*ones(1,N1)); grid on; title(参数估计过程);Theta是4×(N1)因为循环从i3到N2列号从1到N。画图时i取1:N1正好覆盖所有列。真值线用ones(1,N1)乘对应常数画成水平直线。注意Theta(1,:)对应的是-1.8Theta(2,:)对应0.9Theta(3,:)对应0.75Theta(4,:)对应0.5。判断收敛不能只看最后一拍还要看曲线是否在末尾平稳。如果末尾还在明显波动说明数据长度不够或遗忘因子太大。也可以把Pstore画出来figure(2); plot(i,Pstore(1,:),i,Pstore(2,:),i,Pstore(3,:),i,Pstore(4,:));Pstore下降速度反映参数信息积累速度。如果Pstore在几拍内突然掉到接近0说明P初值可能设太大了如果Pstore几乎不变说明输入激励不够。4. 渐消记忆最小二乘把遗忘因子0.998塞进哪里4.1 为什么普通RLS会有稳态偏差原程序注释说“辨识结果会与其精确值始终存在偏差”。这不是程序bug而是标准RLS在这个场景下的固有表现。原因有两层。第一数据长度有限。500点虽然是完整M序列周期但实际噪声是randn生成的一次性样本有限样本噪声均值不一定为0方差也不完全等于1这些随机波动会残留在参数估计里。第二标准RLS的P矩阵随递推收敛到很小增益K趋近于零。当P到后期特别小时参数修正量K*e几乎为零参数曲线变成一条水平线。如果前期的数据有较大的噪声干扰后期就没机会再纠正了。你可以做一个小实验把P初值从100改成1会发现最终偏差变大因为P太小导致早期修正不足。反过来把P初值改成10000会发现初阶段参数振荡剧烈但最终偏差不一定更小。这说明P初值也是一个需要配合数据长度调节的参数。4.2 两行语句的改动遗忘因子的作用位置渐消记忆最小二乘只改两行非常优雅KP*h*inv(h*P*h0.998); P(eye(4)-K*h)*P/0.998;对比标准RLSKP*h*inv(h*P*h1); P(eye(4)-K*h)*P;改动后的P更新多除了一个0.998。因为0.998小于1所以P每次更新后不是缩小而是被放大到约1.002倍。P不再趋近于零增益K也保持在一个合适范围。这样即使到了后期新数据到来时仍能对参数进行修正把早期残留的偏差逐步拉回真值。为什么K公式里的1也改成0.998这一项对应观测噪声方差在渐消记忆版本里遗忘因子λ以类似的形式参与分母。完整的渐消记忆RLS通常写成lambda 0.998; K P*h / (lambda h*P*h); P (eye(4)-K*h)*P / lambda;记住lambda0.998所以分母里那个数就是lambdaP更新除以lambda。0.998非常接近1基本保留了全部历史数据但又没有完全让P衰减到零属于“既维持原有记忆又保留修正能力”。如果lambda等于1就退化成标准RLS。4.3 遗忘因子选多大0.95~0.999的权衡遗忘因子决定了旧数据被淡忘的速度。lambda越接近1遗忘越慢稳态估计方差越小但对参数变化的跟踪越迟钝lambda越小跟踪越快但参数曲线抖动越大。遗忘因子稳态波动跟踪能力典型场景0.999极小几乎不跟踪参数恒定只想微调偏差0.998很小弱跟踪本例参数基本不变0.97中等较快跟踪参数缓慢时变0.95较大快速跟踪参数突变或强时变我做某跨平台系统仿真时曾把lambda设成0.98测试一个参数缓慢漂移的系统效果比0.998好很多但代价是稳态波动明显增大。在真实应用里选lambda的经验是先设0.99看收敛曲线再根据实际方差微调。如果滤波后曲线平稳但跟踪太慢减小lambda如果曲线毛刺太多增大lambda。一个常用技巧是让lambda先取0.98跑前100点之后切回0.998这样可以兼顾初期的快速收敛和后期的低波动。把这段切换逻辑写进循环里并不复杂if i 100 lambda 0.98; else lambda 0.998; end K P*h / (lambda h*P*h); P (eye(4)-K*h)*P / lambda;这样做之后你会发现前100点参数曲线迅速接近真值后面平稳微调整体效果比单一lambda更好。这个变体在一些教科书里叫“变遗忘因子RLS”很适合在课程设计里展示。5. 参数辨识避坑指南初值、噪声、激励和收敛判断5.1 收敛失败先查输入激励和初值现象Theta曲线发散或长时间在真值附近振荡不收敛。原因最常见的是输入信号不满足持续激励。如果输入是常值或正弦单频h矩阵的行向量之间线性相关信息矩阵不满秩参数不可辨识。另外P初值取得过小也会导致早期修正不足。比如Peye(4)而不是100*eye(4)第一拍增益就小从远初值靠近真值会很慢。解决确保输入用M序列或叠加多个不同频率的信号。M序列要检查周期如果抽头配置不对序列周期变短频谱出现线谱相当于激励不充分。初值方面P初值建议至少十倍于参数的量级。比如参数真值是1.8量级P初值取100就合理。我把P初值降为1试过参数曲线需要近200拍才靠近真值明显变慢。5.2 数值爆炸inv、NaN和P奇异现象递推几个循环后Theta突然变成NaN或Inf或者P矩阵对角线急剧增加。原因维度不高时inv一般没问题但如果P初值极大、输入又极小K可能会非常大导致P更新出现数值振荡。另外如果数据重复导致hPh1接近0inv的值爆炸。还有如果系统本身不稳定输出z发散RLS也会跟着发散。解决把inv改成右除写成KP*h/(hPh1)数值稳定性更好。在分母里加一个很小的正数作为正则项比如hPh1e-6能防止矩阵奇异。还要检查观测序列是不是真的稳定画一下z的曲线如果振幅随时间增长那就不是辨识算法的问题而是系统设置或数据生成代码下标错了。我见过有人把z(i-1)误写成z(i1)输出立刻发散。5.3 稳态偏差直流、有色噪声和渐消记忆现象参数曲线很平滑但不是停在真值上而是停在真值附近的一个固定偏差点。原因标准RLS后期增益太小早期偏差无法修正或者M序列是0/1带直流分量使得输出包含偏置估计参数会吸收这个偏置或者噪声不是纯白噪声而是有色噪声最小二乘估计会有偏。解决对于固定偏差用渐消记忆RLS如第4章。对于直流分量把M序列的0改成-11仍为1这样输入均值为0。对于有色噪声可以考虑广义最小二乘或增广最小二乘但那是另一个主题。有一个快速验证方法把randn噪声改成零序列vzeros(1,N)再跑一遍标准RLS如果偏差消失说明偏差主要来自噪声样本的有限均值如果偏差还在说明是输入直流或模型结构问题。5.4 数据长度不足500点到底够不够现象把N从500改成100后参数曲线不收敛末尾仍明显波动。原因100点不足以覆盖M序列完整周期511激励频率不完整且瞬态影响占比更大。另外系统极点模长0.9487时间常数约为1/(1-0.9487)≈19.5拍100点只包含约5个时间常数勉强够但如果有噪声统计精度就不足。解决N至少取一个M序列周期最好取2个周期。本例500接近511已经足够。如果你改周期或系统可以用经验公式N至少是待估参数个数乘10到20倍对于4个参数最少40到80点但那是理想情况有噪声时建议翻倍再翻倍。5.5 收敛判断别只看最后一点现象参数曲线在最后200拍已经平了但最后一拍值仍与真值有偏差有人误以为算法失败。原因RLS的收敛是一个统计过程最终估计是在噪声影响下的一个概率结果不会严格等于真值。曲线平了只代表参数不再大幅更新不代表就是真值。解决把最后50拍的Theta取平均作为最终估计值能有效降低单点随机波动。也可以看残差序列的均方根值如果接近噪声标准差1说明辨识质量可以接受。我一般会这样写theta_final mean(Theta(:,end-49:end),2);这会稍微改善估计精度尤其是当稳态偏差不大时。6. 验证与进阶用残差分析和改变输入信号确认辨识结果6.1 残差白噪声检验参数辨识出来之后不能只看曲线是否贴近真值。一个更严格的做法是检验残差序列是否接近白噪声。残差是最终模型对每个观测输出的预测误差。写一段代码theta_hat mean(Theta(:,end-49:end),2); % 用最后50拍均值 res zeros(1,N); for i3:N2 h[-z(i-1); -z(i-2); M(i-1); M(i-2)]; res(i-2)z(i)-h*theta_hat; end figure; subplot(2,1,1); plot(res); title(残差序列); subplot(2,1,2); autocorr(res); grid on;如果残差自相关函数在非零延迟处基本落在95%置信区间内说明模型已经捕获了系统的动态剩下的只是白噪声。如果自相关明显超界比如延迟为1或2时超出说明模型阶次偏低或者噪声本身有色。对于本例残差应该接近白噪声。如果残差均值明显不为零多半是输入直流或模型偏差。6.2 不同遗忘因子对比实验把渐消记忆RLS封装成一个循环分别取lambda1、0.998、0.97、0.95记录每种情况最后的估计误差。你几乎肯定能看到lambda1时误差最大lambda0.998时误差减小lambda0.95时误差可能略微增大但跟踪能力最强。这个对比实验是理解遗忘因子最好的方式。我建议把四条Theta(1,:)曲线画在同一张图里你能清楚看到哪个收敛最快、哪个波动最大。实际跑下来lambda0.95的曲线在末尾会像毛刺一样跳动而lambda0.999的曲线几乎是一条平线但位置偏高。6.3 扩展到更高阶系统与变体这套程序框架非常容易扩展。把h向量从4维改成2n维Theta初值也改成相同长度就能辨识一个n阶系统。递推公式本身不关心阶次只关心h与Theta的维度匹配。我曾经把一个二阶模型扩展到四阶只改了z生成代码里的系数、h的构造和真值线常量其他都没动。扩展时最容易翻车的是下标偏移。高阶差分方程往往涉及z(i-1)到z(i-n)输入m(i-1)到m(i-n)写循环时务必画一个时间轴验证i的起点和z的初值个数。我的习惯是先用零噪声跑一遍参数能精确回到真值再上噪声。这样能把代码bug和噪声影响分开。还有一个进阶技巧在标准RLS基础上加入输出权值限制可以防止参数突变。具体做法是每次更新后检查Theta的每个分量是否超出合理范围超了就按上界方向拉回。这个技巧在工业现场特别管用因为传感器偶尔会有野值一旦野值进入递推参数可能瞬间跑飞。从那以后我每次做辨识都强制走一遍“零噪声验证、残差白噪声检验、遗忘因子对比”这三步再也不敢只看估计曲线就收工。希望帮到你。本文还有配套的精品资源点击获取