ARTICLE DETAIL

资讯详情

深耕网站建设与运营推广的一线实战洞察。

极大似然法系统辨识:RML递推估计算法与工程实践

极大似然法系统辨识:RML递推估计算法与工程实践 简介极大似然法.zip面向系统辨识与参数估计学习者基于极大似然估计MLE和递归极大似然RML方法提供了从理论讲解到仿真验证的完整流程适用于需要理解模型辨识原理或完成相关实验的人群。压缩包共19个文件核心为6个MATLAB脚本.m用于实现RML算法、M序列生成与误差测试9张jpg图片展示白噪声、有色噪声等条件下的参数估计误差和输出误差对比2个mat数据文件保存仿真中间结果另有PPTX理论讲解与tmw仿真工程配置整体体积仅865KB。目前已有186人学习下载资源轻量但覆盖完整适合控制类专业学生或工程师快速上手。包内除RML_r.m、RML_1000.m等可直接运行的示例外还配有参数估计误差对比图、输出误差曲线和PRBS激励信号生成程序可直观观察不同噪声环境下的辨识效果。结合PPT中的步骤说明读者能够跟随代码掌握递归极大似然的参数更新逻辑进一步迁移到线性时不变系统的建模与控制场景中。1. 极大似然法做系统辨识为什么我放弃最小二乘改用 RML做系统辨识的工程师多半都经历过这个阶段先用最小二乘跑通流程然后发现有色噪声下估计结果总差一口气。极大似然法的核心思路很直接——把观测数据看作随机过程的一次实现找到让这些观测出现概率最大的那组参数。和最小二乘相比它显式建模噪声的统计特性在残差不是白噪声的场合也就是工程上最常见的场合能把参数估计的偏差压下去。这份《极大似然法.zip》里最值得看的就是 RML递推极大似然估计的完整实现。它解决的问题很具体带噪声的输入输出数据如何在线递推地辨识系统参数。适合正在做过程建模、控制器设计前需要拿到可靠模型、或者被最小二乘的有偏估计折腾过的人。2. 极大似然和 RML 的定位选型之前先搞清楚它在解决什么问题2.1 为什么最小二乘在有色噪声下会翻车最小二乘的估计式长这样theta_hat (Phi^T * Phi)^(-1) * Phi^T * Y它本身不假设噪声是白的还是色的形式上只是最小化残差平方和。问题出在一致性上当噪声是有色的也就是 e(k) 和过去的输入输出相关时Phi^T 和 e(k) 的期望不为零估计结果就不收敛到真值而是收敛到某个有偏的值。工程现场最常见的就是这种情况测量噪声经过传感器、滤波电路之后几乎不可能是严格白噪声。极大似然法换个角度切入。它先写出数据的联合概率密度然后最大化这个似然函数。在高斯噪声假设下极大似然估计等价于极小化带噪声模型加权的预报误差平方和。这个加权不是随便加的它由噪声模型的逆决定这就是为什么极大似然能把有色噪声的影响剥离掉。2.2 RML 在极大似然里的位置从离线到递推极大似然估计的离线版本很直接写出似然函数然后求极值就行但实际算起来要处理矩阵求逆和数值优化样本一长就慢。RML 解决的是在线问题每来一组新数据就更新一次参数估计不需要重算历史数据。RML 全称 Recursive Maximum Likelihood核心是把极大似然的梯度下降思想变成递推形式。它的标准做法是用预报误差的梯度构造修正方向用协方差矩阵控制步长同时把噪声模型的参数也放进估计向量一起辨识。这决定了 RML 的实际形态参数向量是扩展的不但包含过程模型参数还包含噪声模型参数。2.3 RML 和使用场景的匹配什么时候值得上 RML我在实际项目里对 RML 的定位就一句话在线辨识、有色噪声、需要白残差三个条件同时满足就值得上。如果只是离线建模数据量又不大直接用离线极大似然估计就行。如果噪声确实是白的那最小二乘本来就够用。RML 的典型使用场景包括自适应控制系统的对象在线辨识、过程模型随工况漂移时的跟踪估计、以及需要在辨识基础上做预测但残差必须满足白噪假设的场合。它最大的优点是把模型质量这个概念落到实处——不仅给参数还给你一组接近白噪声的残差后续模型的校验工作好做得多。3. 把 RML 跑起来核心递推公式与代码实现3.1 RML 的递推结构拆解RML 的递推结构是在递推最小二乘RLS的基础上加了一个噪声模型滤波器。最常用的形式是控制变量法也常被叫成增广矩阵法。核心思路把噪声模型的输出也当成输入数据的一部分然后构造增广的信息向量和参数向量。离散时间下假设对象模型写成A(z^(-1)) * y(k) B(z^(-1)) * u(k) C(z^(-1)) * e(k)其中 A、B、C 都是后移算子 z^(-1) 的多项式。RML 递推要同时估计 A、B、C 的系数。它的递推更新式在结构上和 RLS 一样但信息向量里包含的是经过滤波的变量而且是残差替代新息——这是因为真实的新息噪声无法直接测量只能用模型残差来近似。带遗忘因子的递推形式如下K(k) P(k-1) * phi_f(k) / (lambda phi_f(k)^T * P(k-1) * phi_f(k)) theta(k) theta(k-1) K(k) * (y(k) - phi_f(k)^T * theta(k-1)) P(k) (P(k-1) - K(k) * phi_f(k)^T * P(k-1)) / lambda其中 phi_f(k) 是滤波后的信息向量lambda 是遗忘因子P 是协方差矩阵。关键点在于 phi_f(k) 里有些分量是经过滤波器处理的不是直接用原始输入输出。3.2 一份可直接跑的 Python 参考实现下面给一份简化的 SISO 系统 RML 递推辨识代码模型阶次取 n_a 2、n_b 2、n_c 1。我习惯把它当作调试基线来用跑通之后再去加限制条件或者改成 C 代码。import numpy as np def rml_step(theta, P, phi_raw, yk, lam, n_a, n_b, n_c): 单步 RML 递推 theta: 当前参数向量 [a1, a2, b1, b2, c1] P: 协方差矩阵 phi_raw: 原始数据向量 [y(k-1), y(k-2), u(k-1), u(k-2), e(k-1)] yk: 当前输出 lam: 遗忘因子 # 预测残差用当前参数计算模型输出和真实输出的差 y_pred np.dot(theta, phi_raw) e_pred yk - y_pred # 增益向量 K phi_f phi_raw.copy() # 简化版用原始向量代替滤波结果完整实现见下文说明 K P.dot(phi_f) / (lam phi_f.T.dot(P).dot(phi_f)) # 更新参数和协方差 theta_new theta K * e_pred P_new (P - np.outer(K, phi_f.T.dot(P))) / lam return theta_new, P_new, e_pred这个函数是 RML 递推的骨架先算预测误差再算增益 K然后用 K 乘以误差修正参数。单看代码会发现phi_f直接用了原始向量这是简化版。理由很简单真实 RML 里信息向量各分量要用余差滤波器和辅助模型处理但第一版跑通时应该先把递推核心验证好再加滤波器逻辑。注意theta的初始化和P的初始化很关键。常见做法是 theta 初始化为零P 初始化为一个对角阵对角线取 100 到 10000 之间的值。太小收敛慢太大初期参数震荡剧烈。这个值本质上是你对初始参数不确定度的一种表达。3.3 完整的前向滤波真正可用的 RML 信息向量上一节的简化版没解决有色噪声和滤波的问题实际 RML 需要对信息向量做一步关键处理。RML 的信息向量由三个部分组成第一部分是过去的输出 y(k-1) 到 y(k-n_a)这部分直接测量得到第二部分是过去的输入 u(k-1) 到 u(k-n_b)同样直接测量得到第三部分是过去的残差 e(k-1) 到 e(k-n_c)它由 y_pred 减去 yk 算出来但要经过噪声模型的倒数滤波。具体的滤波关系是残差序列经过 1/C(z^(-1)) 滤波后得到的信号用于构造信息向量。它的工程含义是把残差变成经过系统动态加权的量让递推过程能区分开过程动态和噪声动态。完整代码在这个 zip 里给了测试完简化版再替换phi_f那一段就行。# 完整版中的信息向量构造 def build_phi_f(y_hist, u_hist, e_hist, c_coeff): # c_coeff 是 C 多项式系数例如 [c1, c2] # 对 e 序列做 1/C 滤波 n_c len(c_coeff) e_f np.zeros_like(e_hist) for k in range(1, len(e_hist)): e_f[k] e_hist[k] sum(c_coeff[j-1] * e_f[k-j] for j in range(1, min(k1, n_c1))) # 组装信息向量顺序与 theta 对齐 return np.concatenate([y_hist[-2:], u_hist[-2:], [e_f[-1]]])这段代码的逻辑先对残差序列做 1/C 滤波再把滤波后的值作为信息向量的最后一段。c_coeff是噪声模型参数在递推过程中会不断更新所以每次构造信息向量时都要用最新的theta里的 c 部分。顺序上有个坑theta 里 c 系数排在最后构造信息向量时也要对齐这个顺序。3.4 lambda 遗忘因子的选取从常数到自适应遗忘因子 lambda 决定模型跟踪时变系统的速度。lambda 取 1 表示不遗忘所有历史数据权重一样适合时不变系统。lambda 取 0.95 到 0.99 之间表示快速遗忘适合工况经常变化的系统。我踩过的坑是无脑取 0.98 跑时变系统结果系统其实没怎么变参数却因为遗忘过快而持续抖动。反过来说真正时变的系统取 0.999 也不合适跟不上变化速度。经验上先做一次离线辨识看残差大致稳定后再定 lambda。工程上还有一种做法是变遗忘因子残差大时减小 lambda 加大修正力度残差小时增大 lambda 保持稳定相当于给递推加了一个自适应步长。这个 zip 里带的 RML 主程序就实现了一套简单的变遗忘逻辑逻辑入口是残差能量的滑动窗口估计。4. RML 落地避坑指南参数设不对和数值不稳是两大根源4.1 初始协方差矩阵 P 设置不当导致的前几百步震荡现象递推刚开始时参数大幅度来回跳甚至符号都翻来覆去要等很久才慢慢收敛。原因P 初始值直接决定增益 K 的初始大小。P 太小相当于你认为初始参数非常准增益很低数据推不动参数P 太大则初始增益极大前几步的噪声和模型误差被严重放大。解决把 P 初始化为对角阵对角线取 1000 倍到 10000 倍的单位量级。简单做法是取 P0 1000 * I然后观察前 50 步的响应。如果震荡超过 200 步还在持续就减小到 100I如果收敛过慢增大到 5000I。多试几次就知道什么量级适合你的数据采样率。4.2 数据没有去趋势项或零均值化辨识结果整体偏置现象模型的直流增益和实际系统对不上误差主要存在于低频段残差序列里能看到明显的缓慢漂移成分。原因极大似然估计对均值结构很敏感。输入输出数据的直流分量和趋势项会被模型试图吸收到参数里导致参数估计被系统性带偏。这种情况在最小二乘里也会出现但 RML 对它的惩罚更狠因为噪声模型也在同步辨识。解决数据进 RML 之前先做预处理。每个信号减去自身均值有趋势就做一阶差分或者多项式拟合并减去。注意这里说的减均值是对整个批次做还是滑动窗口做——在线场景下要维护一组滑动均值估计不能等批次结束再处理否则实时性就没了。4.3 遗忘因子固定导致残差出现周期性的复发尖峰现象递推过程中参数已经稳定但每隔一段就出现一次明显的参数跳动然后慢慢回归原位残差里对应出现一个脉冲。原因固定遗忘因子下如果系统一段时间内激励不足比如说输入基本不变信息量变少协方差降不下去此时来一个较大的扰动就会被过度放大。这类问题在自适应控制里非常经典叫协方差爆发。解决换成变遗忘因子或者对协方差矩阵加一个下限约束避免 P 无限膨胀。具体的做法是每次更新后检查 P 的迹超过阈值就统一缩放。常见做法是检查迹的上限超过上限时对角元素乘个 0.1 的缩放系数。4.4 噪声模型阶次定错直接导致参数估计不一致现象过程参数估计结果看着还行但对残差做白噪检验不通过残差的自相关函数显著不为零。原因噪声模型 C(z^(-1)) 的阶次 n_c 低于实际的噪声动态阶次。C 阶次不够时残差里还有未被解释的相关结构直接表现为残差不白。反过来说如果 C 阶次过高参数增多反而加大了过拟合风险。解决做一组阶次扫描实验从 n_c 0 到 n_c 3 分别跑一遍辨识然后对比残差白噪检验的 p 值。选 p 值显著大于 0.05 的最小阶次。注意 n_c 0 时就退化成了递推增广最小二乘这也是 RML 的一个退化边界。4.5 数值溢出和协方差失去正定性现象递推到后半段P 矩阵不是对称正定的了有时候对角线直接变成负数参数更新开始发散。原因递推式的除法累计算误差加上输入数据持续很小导致 P 更新中减法项主导。更常见的原因是有限字长下P 矩阵的对称性被破坏后误差逐步累积最终失去正定性。解决强制对称化。每次更新 P 后执行 P (P P.T) / 2。并且每步检查最小特征值如果低于 1e-12 就直接重置 P 为初始值。这是工程里血泪经验最多的地方数值稳定性问题往往不是理论公式错了而是在有限精度下细节失控。另外建议全程使用 64 位浮点不要用 32 位。5. 验证 RML 辨识效果的三个硬指标残差白噪检验、变工况预测和对比实验拿到 RML 估计的参数后最重要的事不是看参数本身多接近真值而是验证模型能不能用。推荐三步走残差白噪检验、变工况预测、与最小二乘基线对比。这三个做完了模型的可靠性基本就有数了。残差白噪检验的做法把递推过程里保存的残差序列拿出来计算自相关函数 r(i) E[e(k) * e(k-i)]主要看前 20 个滞后。白噪声的自相关除了零滞后外都接近零。工程上的判定标准是95% 置信区间内的自相关点不超过总点数的 10% 就算合格否则就要回头调噪声模型阶次或者检查数据预处理。变工况预测要做的就是换一段没有参与辨识的数据用辨识出来的模型做多步预测比较预测误差。这一步的价值是区分拟合好和泛化好。我之前做温控系统的辨识时遇到过这种情况递推残差挺白但模型拿到另一条工况完全没辨识过的数据上一测预测误差直接翻了倍。后来查下来是激励信号不够丰富频段只覆盖了一个窄带模型在带外完全没约束。对比实验的做法就简单了同一份数据同时用 RLS 和 RML 跑一遍辨识对比两者参数估计值和真实系统的阶跃响应曲线。有色噪声场景下 RLS 的参数偏差通常在 5% 到 20%RML 能把它压到 2% 以内。这个对比不是用来证明 RML 更高级而是用来确认你当前数据的噪声特性确实值得上 RML。从那以后我每次做完 RML 辨识都强制走一遍这三步验证顺序也是固定的先白噪检验确认噪声模型选对了再变工况确认模型泛化了最后和基线对比确认投入产出比值了。这套流程帮我挡掉过三次准备直接交给控制组的模型也算是花钱买来的习惯。希望帮到你。本文还有配套的精品资源点击获取
返回列表