★ FEATURED ARTICLE 简介这是一套面向系统辨识学习者的极大似然法实践材料聚焦极大似然估计与递归极大似然RML在动态系统参数估计中的应用适合已具备概率统计与控制理论基础、希望理解在线辨识流程的读者用于课程实验或算法复现。压缩包共19个文件以MATLAB的.m脚本、jpg实验曲线、mat数据、pptx讲稿及tmw工程文件为主整体约865KB便于按理论、代码与结果三条线对照查阅。已有189人学习下载。内容围绕RML的参数递推更新、白噪声与有色噪声下的估计误差对比、PRBS与M序列激励实验展开配有传递函数或状态空间模型参数估计结果及误差分析图可帮助读者把极大似然原理落到可运行的代码和可验证的实验结论上并理解非高斯噪声、非线性系统等场景下的适应性与局限。1. 结论先行带噪数据里最小二乘为什么“准不了”极大似然又是怎么扳回来的同一组实测输入输出数据最小二乘能把参数算出来可换一段激励信号重跑零点位置就飘出十几倍——这是系统辨识入门阶段最容易踩坑的一次。问题不在解方程在于最小二乘默认了方程误差是白的而真实扰动几乎永远是有色的传感器漂移、量化台阶、机械谐振全都在残差里排着队。极大似然法Maximum Likelihood走的是另一条路把噪声的谱形状也当成待估对象一起建模再让观测数据在这套模型下出现的概率最大。RMLRecursive Maximum Likelihood是它的递推形态每采到一个新样本就把参数和梯度往正确方向推一格适合在线、缓变、算力有限的场合。搞懂它你就搞懂了辨识里那个绕不开的取舍为什么多估一组噪声参数反而能把系统参数估得更准。2. 似然函数是怎么从噪声模型里长出来的多估一组噪声参数为什么反而更准很多人第一次看到 ARMAX 结构时会本能地反感本来只要估 A 和 B 两组系数为什么还要硬塞一个 C 多项式进去参数变多了方差不应该更大吗这个直觉在无偏性面前是失效的——多出来的那组参数不负责描述系统它负责描述“残差长什么样”而只有把残差描述对了A 和 B 的无偏性才成立。2.1 噪声假设才是偏差的源头方程误差与输出误差的分野辨识里有两套完全不同的误差注入方式混淆它们是偏差的第一来源。方程误差结构写成 y(t) φᵀ(t)θ e(t)回归向量 φ 只含 t-1 及更早的 y 和 u。只要 e(t) 是白的φ(t) 与 e(t) 不相关最小二乘就是无偏的这是教科书里那条漂亮结论。可真实系统更像 A(q)y(t) B(q)u(t) v(t)其中 v 是有色扰动。此时残差不再是白噪声而 v(t) 里含有过去扰动的成分过去扰动又已经渗进了 y(t-1)、y(t-2) 里于是回归向量和残差相关最小二乘的偏差不会随样本数趋于无穷而消失。输出误差结构写成 y(t) G(q)u(t) v(t)噪声直接叠在输出端。这时最小二乘的回归向量里含着带噪的 y属于典型的变量误差errors-in-variables问题结果是一个系统性的衰减偏置噪声越大估计出的增益越小频率响应被整体“压扁”。这类偏差最坑人的地方是它看起来很稳定——多跑几组数据、换个初值估计值都一致地偏在同一个位置你甚至会以为那就是真值。极大似然法的处理方式是不假装噪声不存在而是给它一个参数化模型 H(q,θ)然后定义白化后的预测误差ε(t,θ) H⁻¹(q,θ)·[ y(t) − G(q,θ)u(t) ]只要 G 和 H 的结构与真实系统匹配ε 就是白的与过去数据渐近不相关。这一步“白化”就是全部秘密所在它把有色扰动摊平让后续的二次型准则重新变得无偏。2.2 从联合概率密度到可算的代价函数σ² 被“浓缩”掉了假设 e(t) 是零均值、方差 σ² 的高斯白噪声且独立同分布。N 个样本的联合密度是p(Y|θ,σ²) (2πσ²)^(−N/2)·exp[ −(1/(2σ²))·Σ ε²(t,θ) ]取负对数再除以 N得到V(θ,σ²) ½·ln(2πσ²) (1/(2Nσ²))·Σ ε²(t,θ)对 σ² 求偏导置零解出 σ̂² (1/N)·Σ ε²(t,θ)。把它代回去σ² 这个讨厌的未知量就消失了剩下所谓的浓缩似然concentrated likelihoodV(θ) ½·ln[ (1/N)·Σ ε²(t,θ) ] 常数这个式子值得盯着看几秒。它说明对单输入单输出系统而言极大似然估计和“让预测误差平方和最小”的预测误差法在最优点上是同一件事——因为 ln 是单调函数最小化 ln(Σε²) 和最小化 Σε² 给出同一个 θ̂。所以极大似然的价值不在于换了一个更花哨的准则函数而在于换了一个 ε 的定义ε 被 H⁻¹ 滤过而不是原始方程误差。注意这个等价性只在“σ² 自由、单输出”时成立。多输出时准则里出现的是 ln det ΛΛ 为残差协方差矩阵ln det 与 trace 不是单调关系通道之间会产生加权竞争这时极大似然才真正区别于普通最小二乘。做多通道辨识的人如果直接套单通道结论很容易在一组弱激励通道上翻车。2.3 频域视角似然最大化等价于按噪声谱的倒数加权把上面那套推导搬到频域会得到一个更有工程直觉的结论当数据量趋于无穷、模型阶次随数据增长时极大似然估计的渐近代价近似为∫ |G(e^{jω}) − G₀(e^{jω})|² · Φu(ω) / Φv(ω) dω这里 G₀ 是真系统Φu 是输入谱Φv 是噪声谱。它的含义很直白拟合误差在信噪比高的频段被重罚在噪声淹没的频段被轻罚。换句话说极大似然自动帮你做了一个频率加权权重正好是噪声谱的倒数。这也解释了工程上常说的“预滤波”为什么有用。给输入输出数据先过一个 L 滤波器其中 Φv L·L*本质上就是在手工近似这个权重。RML 里的那个 1/Ĉ 滤波做的就是同一件事只不过它不用你猜而是根据当前参数估计自动更新。3. RML 的递推内核每一步到底在算什么把这套准则搬到在线场景就得回答三个问题梯度怎么算、步长怎么定、历史信息怎么衰减。RML 的答案分别是一步预测误差的滤波、协方差矩阵递推、遗忘因子。3.1 梯度 ψ 而不是回归向量 φRML 与 ELS 的分水岭以 ARMAX 为例A(q)y(t) B(q)u(t) C(q)e(t)三个多项式都取首一形式。预测误差满足C(q)·ε(t) A(q)·y(t) − B(q)·u(t)对参数逐项求偏导把 C 视为定值得到∂ε/∂aᵢ y(t−i)/C ∂ε/∂bᵢ −u(t−i)/C ∂ε/∂cᵢ −ε(t−i)/C按 ψ −∂ε/∂θ 定义梯度方向于是ψ(t) [ −y_f(t−1..t−na), u_f(t−1..t−nb), ε_f(t−1..t−nc) ]ᵀ下标 f 表示该信号被 1/C(q) 滤过。把它和增广最小二乘ELS的回归向量φ(t) [ −y(t−1..t−na), u(t−1..t−nb), ε(t−1..t−nc) ]ᵀ摆在一起对比差别只有一处RML 多了一个 1/C 白化滤波。就这一处决定了 ELS 在有色噪声下渐近有偏、而 RML结构正确时渐近无偏。原因在 3.1 的开头已经埋好了φ 里的 ε(t−i) 本身就是被噪声污染的与当前残差相关ψ 里的 ε_f 是白化后的量渐近与 ε(t) 不相关更新方向不再带系统性偏差。提示ψ 的符号约定必须和更新式保持一致。常见写法是 θ̂(t) θ̂(t−1) P(t)ψ(t)ε(t)也有代码写成 φ/C 而不带负号靠 R 和 ε 的符号互相抵消。两种都能跑但混用会出现“参数朝反方向爬、代价函数不降”的诡异现象排查时优先核对这里。3.2 协方差递推与遗忘因子R 矩阵到底承担了什么Gauss-Newton 形式的更新是 θ̂(t) θ̂(t−1) R⁻¹(t)ψ(t)ε(t)其中 R(t) Σ λ^(t−k)·ψ(k)ψᵀ(k)λ 是遗忘因子。用矩阵求逆引理把它改写成 P R⁻¹ 的递推形式就得到能在线跑的版本K(t) P(t−1)ψ(t) / (λ ψᵀ(t)P(t−1)ψ(t)) θ̂(t) θ̂(t−1) K(t)·ε(t) P(t) [ P(t−1) − K(t)ψᵀ(t)P(t−1) ] / λ遗忘因子 λ 的物理意义是有效记忆长度约等于 1/(1−λ)。λ 0.99 对应大约 100 个样本的记忆λ 0.999 对应约 1000 个。系统时不变就取 λ 1此时 P 单调收缩参数收敛后基本静止系统缓慢漂移才取小于 1用方差换跟踪能力。P 矩阵承担的角色常被误解。它不是“学习率”而是一个方向加权的度量激励充分的参数组合方向上 P 迅速收缩、步长变小激励不足的方向上 P 迟迟不降、步长保持很大。这既是他能快速收敛的原因也是数值不稳定的根源——见 3.4。3.3 一份可跑的最小实现下面这段把上面的公式逐行落地故意不做任何工程化包装方便对照公式阅读。function [theta, P, epsv] rml_armax(u, y, na, nb, nc, lambda) % 递推极大似然ARMAX 结构 A(q)y B(q)u C(q)e % theta [a1..ana, b1..bnb, c1..cnc] % lambda 为遗忘因子时不变系统取 1缓变系统取 0.98~0.999 N numel(y); n na nb nc; theta zeros(n, 1); P 1e4 * eye(n); % 大初值 对初值几乎不信任先大步走 L max([na nb nc]) 1; yf zeros(L,1); uf zeros(L,1); ef zeros(L,1); % 1/C 滤波后的历史 ebuf zeros(nc,1); % 未滤波的过去预测误差 epsv zeros(N,1); for t max([na nb nc]) 1 : N c theta(nanb1 : end); % (1) 用 t-1 时刻的参数算先验预测误差 phi [-y(t-1:-1:t-na); u(t-1:-1:t-nb); ebuf]; e y(t) - phi * theta; epsv(t) e; % (2) 用当前的 C 把 y、u、eps 各过一遍 1/C得到滤波信号 yf [y(t) - c*yf(1:nc); yf(1:L-1)]; uf [u(t) - c*uf(1:nc); uf(1:L-1)]; ef [e - c*ef(1:nc); ef(1:L-1)]; % (3) 梯度 滤波后的回归向量注意跳过 yf(1) 这类当前时刻项 psi [-yf(2:na1); uf(2:nb1); ef(2:nc1)]; % (4) 带遗忘因子的 Gauss-Newton 更新 K P*psi / (lambda psi*P*psi); theta theta K * e; P (P - K*psi*P) / lambda; % (5) 滚动缓存过去误差 ebuf [e; ebuf(1:nc-1)]; end end几个参数的作用需要说清楚。na、nb、nc决定模型结构nb里已经含了纯延迟如果你的系统从 q⁻² 才开始有响应要把那一项补零而不是把nb减一否则 φ 与 ψ 的时间对齐会整体错位代价函数降不下去。lambda取 1 时 P 每步都不会被放大收敛后参数静止取小于 1 时记得关注 P 的迹它长期不降说明激励不足。1e4 * eye(n)这个初值不是玄学它对应“参数方差先验极大”让前几十步的步长足够大代价是先期会有几次明显超调。第 (2) 步是整段代码最容易被质疑的地方C 是每步都在变的那么“1/C 滤波”严格来说是一个时变滤波算出来的 ψ 只是准梯度而不是精确梯度。这是所有递推预测误差类算法的共同妥协——只要 Ĉ 收敛得比 θ̂ 慢实际上它们同步收敛这个近似就不影响渐近性质。3.4 数值稳定性P 为什么会在没有新信息的方向上爆掉P 的递推式里有一个除以 λ 的操作。当 λ 1 时每一步 P 都被整体放大 1/λ只有在 ψ 张成的方向上K ψᵀ P 这一项才会把它压回去。如果某个方向长期没有激励ψ 在该方向的分量为零压低项不起作用而放大项照旧生效P 就在这个零空间里指数增长。跑上几千点之后P 的量级能冲到 1e12 以上此时浮点舍入误差足以让 P 失去对称正定性接着参数开始乱跳、发散看起来很像是“算法崩了”其实只是共线性在算数上发作了。工程上有几条常见的加固手段一是改用 UD 分解或平方根滤波把 P 写成 U·D·Uᵀ只递推 D 的对角元和 U 的乘子正定性由结构本身保证不靠数值运气二是给更新加一个死区当 |ε(t)| 小于某个阈值时干脆跳过这一步的更新避免纯噪声驱动参数三是监控 P 的迹超过阈值就重置为初始大矩阵等于承认“过去的信息已经不可信了”。提示初次上手时建议先设 λ 1 跑一遍看看参数是否收敛到和批处理 PEM 一致的位置再切到 λ 1。直接上遗忘因子会同时引入跟踪能力和数值风险两个变量一起动出问题根本分不清是谁的责任。4. 收敛性似然面为什么不是一碗水梯度下降类算法在凸问题上收敛是理所当然的可似然函数一般不是凸的。参数化之后的预测误差对参数是非线性的因为 ψ 依赖 θ似然面可能有多座山峰也可能沿着某些方向完全平坦。搞不清这两件事调参就全靠试。4.1 持续激励与可辨识性R 什么时候会病态ARMAX 结构下的可辨识性对输入信号有硬要求u 必须持续激励到足够阶次常见的结论是阶次至少要覆盖 na nc。道理很直观——要让 A 和 C 的多项式系数被分开辨识输入必须先把系统的动态模态都激发起来如果输入是阶跃或者几条正弦某些模态从未被激励A 和 C 的系数就可以互相补偿而输出完全不变似然面在这些方向上是一条平谷。闭环数据是另一个高频雷区。闭环里 u 是 y 的反馈与噪声天然相关如果控制器里还带积分作用u 的低频分量几乎完全由噪声积分决定ψ 的某些列近似线性相关R 的条件数会非常难看。标准解法是加 dither 信号小幅伪随机扰动打破相关性或者改用间接法、两步法先估闭环传递函数再反解。判断激励是否足够最直接的办法是看 P 的特征值谱。收敛后最大特征值与最小特征值之比超过 1e8基本可以判定存在不可辨识方向继续跑下去得到的参数没有统计意义即便输出拟合得不错。4.2 初值与参数化为什么先用批处理跑一段再切递推递推算法对初值敏感因为它是局部搜索。默认全零初值意味着从一个“无动态”的模型出发如果真系统有靠近单位圆的极点前几百步的梯度方向可能把参数推向一个局部极小之后再也出不来——表现是代价函数早早变平但残差自相关明显P 也降不下去。一个稳妥的做法是先用批处理 PEM 在全量数据上做一次 Gauss-Newton哪怕只迭代十几轮把结果作为递推的初值同时把批处理估出的渐近协方差当作 P 的初值而不是一律 1e4。这样递推阶段基本只做微调不容易走岔。过参数化是另一个陷阱。na 取得比真实阶次大A 和 C 之间会出现近似零极点对消似然面沿着对消方向极其平坦参数可以在这条线上任意滑动而输出不变。这时候你会看到 a 系数和 c 系数此消彼长、数值都很大模型本身却拟合得很好。这不是算法的问题是模型结构冗余靠 AIC、FPE 或者交叉验证降阶才能解决。4.3 收敛的三种验证手段不能只看“参数不怎么动了”就宣布收敛那可能只是 P 变小导致步长变小。下面三类检查一起看判断才靠谱。第一是代价函数的单调性。带遗忘因子时准则为 V(t) Σ λ^(t−k) ε²(k)正常情况下它应该总体下降并在收敛后小幅波动。如果它稳步上升说明 ψ 的符号或者时间对齐错了。第二是梯度范数。收敛点满足 E[ψ(t)ε(t)] 0实操里可以监控 ψᵀε 的滑动平均它应该围绕零附近抖动而不是恒为正或恒为负。第三是残差白性检验这是最能说明问题的一项。import numpy as np def residual_check(eps, u, maxlag20): 残差白性 输入残差互相关检验95% 置信带 ±1.96/sqrt(N) N len(eps) eps eps - eps.mean() u u - u.mean() # 归一化自相关 acf np.correlate(eps, eps, full)[N-1:] / (N * eps.var()) # 输入与残差的互相关越界说明有未建模动态或反馈回路 ccf np.correlate(eps, u, full)[N-1:] / (N * eps.std() * u.std()) bound 1.96 / np.sqrt(N) print(fN{N}, 95% 置信带 ±{bound:.4f}) for k in range(1, maxlag 1): tag_a 越界 if abs(acf[k]) bound else ok tag_c 越界 if abs(ccf[k]) bound else ok print(flag {k:2d} acf{acf[k]:.3f} {tag_a} ccf{ccf[k]:.3f} {tag_c}) return acf, ccf, bound三个参数的用法maxlag取到 20 左右就够看趋势取太大置信带又不变越界概率自然上升容易误判bound是白噪声假设下的 95% 带N 越大带越窄所以长数据下偶尔一两个 lag 越界是正常的ccf那一列如果整体呈现缓慢衰减的包络几乎可以确定存在未建模动态或者闭环反馈此时该回去加阶次或加 dither而不是继续调 λ。5. 避坑与排查RML 上线后最常翻的五个车5.1 现象参数不收敛一直在缓慢漂移原因通常是遗忘因子设得太小。λ 0.9 对应有效记忆只有十个样本噪声对参数的随机推动远大于真实漂移的信号量估计值就变成了一场随机游走。另一种可能是系统本来就时不变而你误开了遗忘。解决把 λ 提到 0.98 以上重新观察 P 的迹确认系统确实时变后再考虑自适应遗忘因子——令 λ(t) 1 − 1/(t·κ)前期快速遗忘、后期自动趋近 1兼顾跟踪和收敛。5.2 现象输出拟合度很高但残差自相关明显越界这是最迷惑人的一种。拟合度高说明模型对输出的预测能力强残差有自相关说明预测误差里还留着可被利用的结构也就是模型阶次不够或者 C 的阶次错了。此时估出来的 A、B 系数是无偏的吗通常不是——残差里的结构会渗进梯度产生偏差。解决先看残差 ACF 在哪个 lag 上有峰那个位置对应缺失的动态再试 nc 加一观察代价函数是否显著下降。如果加了 nc 还是越界考虑换 Box-Jenkins 结构把 G 和 H 完全独立参数化或者检查是否存在非线性。5.3 现象跑几千点后 P 突然失去正定性参数爆炸原因是 3.4 里说的共线性。激励不足的方向上 P 被 1/λ 反复放大浮点误差累积到破坏正定性。某些实现里 P 的更新用了简化式丢掉了对称化步骤误差累积会更快。解决每步强制 P (P P)/2 保证对称改用 UD 分解递推或加一步 P ← P δI 的正则化δ 取 1e-10 量级就够。彻底一点的办法是监控 P 的最大特征值超限就重置。5.4 现象闭环数据下参数永远偏一个固定量偏置方向和幅度都稳定重跑几次都一样说明不是随机误差。原因是闭环里 u 与噪声相关ψ 与 ε 不再渐近不相关一致性条件被破坏如果控制器含积分偏置会更明显。解决在参考信号上叠加小幅 dither人为恢复激励或者放弃直接辨识改用间接法先估闭环传递函数。判断是否属于这一类看 4.3 里的 ccf 那一列整体有缓慢包络就基本可以确认。5.5 现象换个采样率辨识结果全变参数本身没有可比性离散模型的系数本来就随采样率变但更严重的是物理含义变了同样的 C 多项式在采样率改变之后描述的噪声谱形状完全不同。采样太快噪声在带内近似白C 退化成一阶采样太慢高频噪声全被折叠进低频C 被迫去拟合混叠后的谱阶次需要跟着提高。解决固定采样率与抗混叠滤波器的带宽一起调参跨采样率比较时只看连续域指标极点位置、频率响应、物理参数不要比离散系数。这个问题在从仿真数据迁到实测数据时几乎必踩一次。6. 设计取舍批处理 PEM 与 RML 到底该怎么选把两种做法摆在一起看取舍其实相当清楚。维度批处理 PEMRML / 递推形式每步计算量一次迭代 O(N·n²)每样本 O(n²)内存占用需要保留全量数据只留状态与 P局部极小换多个初值重跑即可一旦跑偏很难回头时变跟踪需要分段滑窗遗忘因子天然支持收敛诊断可看整条迭代轨迹只看滑动统计量真正的取舍点不在计算量在“可回退性”。批处理 PEM 随时可以换个初值重新开始RML 没有后悔药——一旦参数被推到一个平坦的局部极小附近P 已经缩得很小后续数据再怎么进来也拉不回去。所以工程上的常见做法是混合先用批处理跑出初值和协方差再切递推到在线模式并且保留一个周期性重跑的守护任务。另一个容易被忽略的取舍是模型结构。RML 的一致性严格依赖“结构正确”这个前提而结构是否正确只能靠残差检验事后确认。这意味着 RML 天然是一个需要人在环路上的算法而不是一个设好参数就能撒手不管的黑匣子。凡是宣称“全自动在线辨识”的方案都要追问一句它的结构是怎么定的。从那以后我每次切到递推模式前都强制先跑一遍批处理 PEM 对齐初值、看一眼残差的 ACF 和 CCF再放手让它在线跑。这一步多花十分钟能省掉后面一整天的排查。希望帮到你。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
