MATLAB最小二乘参数辨识：批量与递推算法实现及避坑指南

简介这是一份面向控制、信号处理等领域学习者与科研人员的MATLAB参数辨识程序文档聚焦最小二乘类算法在系统参数在线估计与辨识中的应用。资源以单个docx文档形式提供体积仅78KB内含完整可运行的递推最小二乘与渐消记忆最小二乘代码实现并配有参数估计过程可视化与结果对比说明。文档基于一个极点位于单位圆内的稳定系统生成500点M序列作为输入并叠加白噪声扰动逐步展示增益矩阵、估计方差及参数更新公式的推导与编程细节同时对比了两种算法在相同条件下的辨识精度差异。该资源目前已有2158人学习下载适合正在学习系统辨识、自适应控制或MATLAB数值仿真的读者可将其作为理解最小二乘参数辨识原理、快速上手编程实践的直接参考资料。1. 最小二乘类的参数辨识MATLAB程序从“凑模型”到“算参数”的第一步做过系统辨识的人都有过这种经历手里握着几十秒的输入输出数据面对一个“大概是二阶、带纯延迟”的工业过程想把传递函数里的四个系数定下来。有人用频域法画Bode图再手动拟合费眼费神遇到噪声大一点的采样序列曲线毛刺能让人怀疑人生。我自己的做法是直接上最小二乘类的参数辨识MATLAB程序——把问题变成“构造回归矩阵解一次超定方程”几十行代码就能拿到参数估计值。这套方法不挑对象线性模型基本都能用既能离线批量算也能改成递推形式做在线辨识。适合正在做课程设计、写论文需要辨识模型的学生也适合在现场调试时想快速拿到一个可用的线性模型、后续再做控制的工程师。读完这篇你能跑通批量最小二乘和递推最小二乘两套程序知道噪声环境下怎么补救也清楚那些“同一份数据不同人辨识结果不一样”的坑到底出在哪。2. 为什么参数辨识首选最小二乘模型假设、解析解与适用边界2.1 先把要辨识的模型写出来ARX结构与回归形式参数辨识的第一步不是写代码而是把“要辨识的对象”用数学结构固定下来。工业场景里最常见的做法是把对象写成ARX模型也就是自回归各态外部输入模型离散形式如下y(k) a1·y(k-1) a2·y(k-2) … ana·y(k-na) b1·u(k-1) b2·u(k-2) … bnb·u(k-nb) e(k)这里y(k)是当前时刻输出u是输入e(k)是噪声项na和nb分别是输出和输入的阶次。这个模型为什么好用因为它把“过去时刻的输出”和“过去时刻的输入”当成已知量把未知参数a和b当成待求系数于是整个问题就变成了一个线性回归问题。把上式改写成矩阵形式会更直观定义回归向量φ(k) [-y(k-1), -y(k-2), …, u(k-1), u(k-2), …]ᵀ参数向量θ [a1, a2, …, b1, b2, …]ᵀ那么每个时刻的观测都能写成一个线性方程y(k) φ(k)ᵀ·θ e(k)。这个形式是后面所有代码的地基。无论是批量最小二乘、递推最小二乘还是辅助变量法本质上都是在处理“从一堆这样的线性方程里把θ解出来”这件事。所以动手写程序前先想清楚你的对象能不能用ARX结构描述。大多数线性定常系统、带稳定工作点的工业过程这层近似是合理的。如果你的对象是强非线性、时变或闭环运行的系统ARX只是起点后面需要更复杂的结构但最小二乘类方法仍然是入口。2.2 解析解与数值实现为什么不要用inv()对于批量处理的数据假设我们有N组观测把所有方程堆起来就得到矩阵形式Y Φ·θ E。其中Φ的第k行就是φ(k)ᵀ是一个N×(nanb)的矩阵Y是N×1的输出向量。最小二乘的准则是最小化残差平方和J(θ) ‖Y - Φ·θ‖²。对θ求导并令梯度为零能得到著名的正规方程θ_hat (Φᵀ·Φ)⁻¹·Φᵀ·Y理论形式很简洁但落到MATLAB里我有一个很明确的习惯永远不要直接用inv(Phi*Phi)去算。原因有两个一是Φᵀ·Φ的条件数是原矩阵条件数的平方数据稍微有些共线这个矩阵就可能病态到求逆结果完全失真二是MATLAB的反斜杠运算符有更聪明的求解路径。对超定方程直接写theta_hat Phi \ YMATLAB会自动选择QR分解路径来求解最小二乘问题数值稳定性比凑出正规方程再求逆好得多。有一个经验值供参考当cond(Phi*Phi)超过1e6时inv的解法大概率已经开始丢失有效数字。尤其现场采集的数据传感器噪声、量化误差、相邻采样点的相关性都会抬高条件数。你可以在程序里用cond()函数确认一下但顺手写成Phi\Y就能绕开这个雷区。这一点不夸张地说能救回不少“明明公式对但结果不对”的翻车现场。2.3 适用边界白噪声、充分激励与数据长度最小二乘不是万能药它的理论保证建立在几个前提上。第一噪声项e(k)必须是白噪声也就是与回归向量φ(k)不相关。这个条件在实际中很难严格满足——比如输出测量噪声经过传感器滤波后天然带相关性又比如闭环运行的系统里噪声会通过反馈回路污染输入信号。一旦噪声有色批量最小二乘的参数估计就是有偏的偏差大小跟噪声强度和相关性有关。第二输入信号必须持续激励。阶跃信号虽然工程上最容易做但它的频谱能量集中在低频对二阶以上的模型往往激不出足够的动态信息导致参数不可辨识或方差极大。第三数据长度要够。理论上N大于参数个数就能算但工程上我一般要求N至少是参数个数的20倍以上否则估计方差会大得没法看。这三个边界条件决定了你选什么工具数据离线、噪声白、激励充分用普通批量最小二乘即可数据在线到达、参数慢变用递推最小二乘加遗忘因子噪声明显有色就要上辅助变量法或增广最小二乘。下面两张表可以帮你快速选型。场景推荐方法关键参数离线数据白噪声信号充分激励批量最小二乘阶次na、nb在线数据参数缓慢变化递推最小二乘 遗忘因子遗忘因子λ噪声有色普通LS有偏辅助变量法IV工具变量构造方式怀疑阶次不对先用低阶LS再配合AIC/BIC定阶阶次搜索范围前提条件不满足时的典型症状补救方向白噪声假设参数估计有偏残差自相关明显IV、增广LS持续激励部分参数辨识不出方差大PRBS输入、叠加随机信号数据长度足够参数震荡换一段数据结果差异大加长采集时间、降阶次3. 批量最小二乘MATLAB程序从仿真数据到参数估计的一条龙实现3.1 先造一份能验证结果的仿真数据学辨识最怕的就是拿真机数据练手因为真实参数未知算出来的结果没法判断对错。所以我每次教学和自测都是先造一份“答案已知”的仿真数据跑通流程后再上真数据。这个习惯能帮你快速定位是程序写错了还是数据本身的问题。下面这段代码生成一个二阶ARX系统的输入输出序列真实参数设为a1-0.8、a20.15、b10.8、b20.2噪声取方差0.01的高斯白噪声。% 生成仿真数据y(k)a1*y(k-1)a2*y(k-2)b1*u(k-1)b2*u(k-2)e(k) rng(42); % 固定随机种子保证可复现 N 1000; % 数据长度 u randn(N, 1); % 白噪声激励信号 e 0.01 * randn(N, 1); % 过程噪声方差为 0.01^2 a1 -0.8; a2 0.15; % 真实参数输出侧 b1 0.8; b2 0.2; % 真实参数输入侧 y zeros(N, 1); for k 3:N y(k) -a1*y(k-1) - a2*y(k-2) b1*u(k-1) b2*u(k-2) e(k); end这段代码里有一个细节值得说明差分方程用的系数是“模型里带正号”的a1、a2而递推计算时负号被显式写了出来。这样做的意图是让真实参数和后面辨识结果直接对照不用来回变符号。随机种子rng(42)的作用是让每次运行生成相同的数据序列方便复现和调试。如果你希望更接近工程环境可以把输入u从randn换成PRBS伪随机二进制序列它的频谱更平直、幅值可控在后文避坑部分我还会专门提。3.2 批量辨识主程序构造回归矩阵、求解与残差分析数据有了接下来就是核心的辨识程序。这一步最关键的活是把数据矩阵Φ构造对——行对应时刻列对应回归变量顺序必须和参数向量θ一一对应。我习惯的写法是先用阶次na、nb声明维度再循环逐行填充这样逻辑清楚阶次改动时不容易漏改。% 批量最小二乘辨识 na 2; nb 2; % 模型阶次 L max(na, nb); % 构造回归矩阵所需的最小起始点 Phi zeros(N - L, na nb); % 回归矩阵每行对应一个时刻 Y y(L1:N); % 输出向量与Phi行数对齐 for k L1:N Phi(k-L, :) [-y(k-1) -y(k-2) u(k-1) u(k-2)]; % 注意列顺序与theta一致 end theta_hat Phi \ Y; % 核心一步反斜杠求解最小二乘问题 a1_hat theta_hat(1); a2_hat theta_hat(2); b1_hat theta_hat(3); b2_hat theta_hat(4); % 残差分析 resid Y - Phi * theta_hat; % 预测残差 V resid * resid / length(Y); % 残差方差 fprintf(估计参数: a1%.4f, a2%.4f, b1%.4f, b2%.4f\n, ... a1_hat, a2_hat, b1_hat, b2_hat); fprintf(真实参数: a1%.4f, a2%.4f, b1%.4f, b2%.4f\n, a1, a2, b1, b2); fprintf(残差方差: %.6f\n, V);这里有几个关键点必须说透。第一回归矩阵从第L1行开始构造是因为前L个时刻缺少足够的y延迟项硬凑只会引入无效行。第二Phi的列顺序必须和参数向量的位置严格一致我先把y的延迟排前面、u的延迟排后面这样做的好处是参数解释直观。第三theta_hat Phi \ Y这一行是整个程序的核心反斜杠在超定系统下自动走QR分解路径比(Phi*Phi)(Phi*Y)多一层数值保护。跑完这段程序你会看到估计参数与真实参数在小数点后两到三位吻合残差方差非常接近0.01²这就说明程序和模型结构都对。3.3 结果验证拟合优度与残差自相关拿到参数不等于完事。我的习惯是立刻做两个检验一是拟合优度二是残差自相关。拟合优度的指标看R²定义如下% 拟合优度检验 SST sum((Y - mean(Y)).^2); % 总平方和 SSE sum(resid.^2); % 残差平方和 R2 1 - SSE / SST; % R² 越接近1越好 fprintf(R² %.6f\n, R2); % 残差自相关检验粗略版 maxLag 20; R_auto zeros(maxLag, 1); for lag 1:maxLag R_auto(lag) resid(lag1:end) * resid(1:end-lag) / (resid * resid); end figure; stem(R_auto); title(残差自相关函数); xlabel(滞后);R²接近1说明模型解释了绝大部分输出变化但在辨识这件事上R²好不代表模型好因为ARX模型本身就包含过去的输出项哪怕参数全错只要数据平滑R²也可能很高。所以我更看重残差自相关如果残差在滞后1到20阶都接近0说明残差是白噪声模型已经把动态信息榨干了如果某个滞后处明显超出±1.96/sqrt(N)的置信带说明还有动态没建模进去或者噪声本身有色。上面这段程序没有调用任何工具箱函数全部用矩阵运算手写方便你移植到自己的工程里。4. 递推最小二乘在线辨识程序与遗忘因子调参实战4.1 从批量到递推RLS更新公式的工程含义批量最小二乘要求数据全部到位才能算但很多场景数据是一点一点来的在线监测、自适应控制、设备老化跟踪。这时批量法有两个硬伤一是每来一个新数据就要重新求解一次矩阵方程计算量随N线性增长二是老数据和新数据权重相同参数如果已经缓慢漂移老数据会拖住估计值不放。递推最小二乘RLS就是为了解决这两个问题而存在的。RLS的核心是一条递推更新链每来一个新观测(k时刻的φ(k)和y(k))就按下面的公式更新一次参数增益K(k) P(k-1)·φ(k) / (λ φ(k)ᵀ·P(k-1)·φ(k)) 参数θ(k) θ(k-1) K(k)·(y(k) - φ(k)ᵀ·θ(k-1)) 协方差P(k) (I - K(k)·φ(k)ᵀ)·P(k-1) / λ这里的λ是遗忘因子取值范围通常在0.95到1之间。λ1表示不遗忘所有历史数据权重相同RLS退化为普通最小二乘的递推形式λ越小旧数据被遗忘得越快参数跟踪能力越强但噪声也会跟着放大。工程直觉可以这样记λ决定“记忆长度”0.98大致相当于记得最近约50个点的信息0.95大约20个点。这个对应关系不是精确的但它能帮你把抽象的λ映射到实际需要的时间尺度上。4.2 带遗忘因子的RLS的MATLAB实现下面这段代码是完整的RLS辨识实现直接接着第3章的仿真数据运行每来一个数据点就更新一次参数并把参数轨迹记录下来方便绘制收敛曲线。% 递推最小二乘辨识带遗忘因子 n na nb; % 待辨识参数个数 theta_rls zeros(n, 1); % 参数初值先给零向量 P 1e6 * eye(n); % 协方差初值大数代表初始不确定度高 lambda 0.98; % 遗忘因子0.95~0.99之间调 theta_hist zeros(n, N); % 记录每步的参数用于画收敛曲线 for k L1:N phi_k [-y(k-1) -y(k-2) u(k-1) u(k-2)]; % 当前回归向量列向量 y_k y(k); % 当前输出 K P * phi_k / (lambda phi_k * P * phi_k); % 增益向量 e_k y_k - phi_k * theta_rls; % 预测误差 theta_rls theta_rls K * e_k; % 参数更新 P (eye(n) - K * phi_k) * P / lambda; % 协方差更新 theta_hist(:, k) theta_rls; % 记录 end % 绘制参数收敛曲线 figure; plot(theta_hist(1,:), b, LineWidth, 1.5); hold on; plot(theta_hist(2,:), r, LineWidth, 1.5); plot(theta_hist(3,:), g, LineWidth, 1.5); plot(theta_hist(4,:), m, LineWidth, 1.5); yline(a1, --b); yline(a2, --r); yline(b1, --g); yline(b2, --m); legend(a1,a2,b1,b2,真实值,Location,best); xlabel(采样时刻 k); ylabel(参数估计值); title(RLS参数收敛曲线);这段代码里有几个必须解释清楚的点。第一P的初值取1e6·I物理含义是“初始时对参数一无所知置信度极低”这样第一次更新时增益K会很大参数迅速从零跳到大概范围。如果初值取1或更小参数收敛会慢很多甚至会一路缓爬。第二分母里lambda加上一个标量φᵏᵀ·P·φₖ这个标量本质是“当前回归向量在新息方向上的不确定性”数值上永远是正数所以不存在除零风险。第三P的更新里除以λ是为了放大协方差、保留对新数据的开放度这也是遗忘因子的数学作用。运行后你会看到参数曲线在前200个点快速收敛之后在真实值附近小幅波动。波动幅度取决于λ和噪声方差λ越接近1波动越小但参数慢变时跟不上λ越小波动越大但能跟上缓变。这个矛盾就是调参的实质。4.3 遗忘因子和初值的调参经验实际调试RLS的时候我一般按下面的顺序做。先固定λ1跑一遍看参数能不能收敛到稳定值这一步验证模型结构和数据质量确认没问题后再把λ往小调观察参数轨迹的跟踪速度和波动幅度。一个常见做法是让λ0.98起步如果参数曲线抖得厉害就提到0.99如果发现参数明显滞后于实际变化就降到0.95。这个“步进搜索”的方式比直接猜一个值高效得多。参数取值范围调大时调小时遗忘因子λ0.95~1.0收敛更稳、跟踪变慢跟踪变快、波动变大P初值1e3~1e8初始收敛快、初期抖动大收敛慢、曲线平滑参数初值θ(0)零向量或先验值有先验时收敛更快无先验时从零学起还有一个容易忽略的细节RLS的收敛轨迹和输入信号的激励强度强相关。如果某段时间输入几乎不变比如阀门停在某个开度φ(k)的信息量不足P会慢慢膨胀或停滞此时参数会出现漂移。工程上我见过有人把λ调成0.9想加快跟踪结果输入一停参数立刻乱跳。正确做法是在线监测预测误差e_k如果e_k长期很小说明系统处于稳定工况此时应该暂停更新或把λ临时调回0.999防止参数被噪声牵着走。5. 最小二乘参数辨识避坑专题5个我反复踩的坑5.1 数据没去均值直流分量把辨识结果整体带偏现象辨识出来的模型稳态增益和真实值差一大截但残差看起来不大R²还挺高。原因采集的y和u往往带有直流分量也就是工作点不为零。ARX模型里如果没有显式的常数项c回归方程Y Φ·θ根本无法描述直流偏移最小二乘只能拿“倾斜”的参数去硬凑这个偏移结果就是动态参数都被污染。我在某温控系统辨识时就翻过这个车气体阀门开度信号有2V偏置温度信号还有1.5V偏置直接辨识出的增益比真值低了将近40%。解决辨识前先做去均值处理或者干脆用detrend函数去掉线性趋势。最稳妥的做法是分别计算u和y的均值并在构造回归矩阵前减掉。注意减均值要在同一时间段内计算别用整段数据的全局均值去减前面的数据那样会引入额外的斜坡失真。5.2 用inv()解正规方程病态矩阵下的“翻车”现场现象程序没报错但算出来的参数数值大得离谱正负号都在乱跳比如a1估计出-3.7而真实值是-0.8。原因直接用inv(Phi*Phi)时ΦᵀΦ的条件数是Φ的条件数的平方。现场数据各列之间天然有相关性比如y(k-1)和y(k-2)在采样周期短时几乎是线性相关的条件数轻易超过1e8矩阵求逆的结果已经丢失了所有有效数字。这不是算法错了是数值方法选错了。解决第一选择是phi Phi \ Y反斜杠自动处理病态问题。第二选择是pinv(Y)或对Phi做QR分解再回代。诊断手段是算一下cond(Phi*Phi)如果发现超过1e6首先要怀疑的不是算法而是数据——要么激励不够要么阶次选高了要么数据段里有大段重复模式。修好数据问题后条件数自然会降下来。5.3 激励信号不丰富阶跃信号辨识不出二阶参数现象用阶跃响应数据辨识二阶模型b2估计值接近0而且换一段数据结果差异极大参数方差大得没法看。原因阶跃信号只有两个幅值状态切换一次之后系统就进入稳态之后的数据几乎不携带动态信息。从频域看阶跃的能量集中在低频二阶模型的极点信息大多在中高频段激不出来。从数学上看输入序列接近常值向量回归矩阵的相关性极高信息矩阵接近奇异参数自然不可辨识。这是最小二乘最经典的失败模式之一。解决换激励信号。标准做法是用PRBS伪随机二进制序列它的频谱近似白噪声、幅值可控对系统冲击小。MATLAB里如果没有系统辨识工具箱可以手写一个简单的m序列生成器或者退一步用randn生成随机幅值序列直接激励。注意激励信号的幅值要根据对象的线性范围来定曾经有一个案例某同事把PRBS幅值取得过大系统进入了非线性区辨识出的模型在后续控制器设计时完全不匹配。% 简单的PRBS生成备选方案 nBits 10; N_prbs 2^nBits - 1; % m序列长度 reg ones(1, nBits); % 移位寄存器初值 prbs zeros(N_prbs, 1); for p 1:N_prbs prbs(p) reg(1); fb xor(reg(10), reg(1)); % 反馈抽头不同抽头对应不同序列 reg [fb, reg(1:end-1)]; end u (prbs - 0.5) * 2; % 转成±1幅值序列5.4 有色噪声使估计有偏残差检验与补救方向现象参数估计和真值明显不一致残差方差很小但自相关函数在低滞后处超出置信带。原因最小二乘的无偏性依赖“回归向量与噪声不相关”这个白噪声假设。如果测量噪声经过传感器低通滤波或者系统有未建模的高频动态残差就会带上相关性此时θ̂的期望值不再是θ真值。这个问题在开环数据上也可能出现闭环系统里则几乎是必现因为反馈回路会把噪声耦合回输入。解决先做残差自相关检验来确认问题存在然后换用辅助变量法IV或者增广最小二乘。IV的基本思路是找到一组与噪声不相关、但与信号强相关的工具变量来替代回归矩阵。后文第6章会给出具体代码。这里想强调一点残差小不是模型好的充分条件一定要做残差自相关检验这是区分“模型解释得不错”和“模型统计上正确”的分水岭。5.5 遗忘因子调太小递推参数震荡的排查现象RLS刚开始收敛得挺好跑到某一时刻参数突然剧烈震荡甚至发散把控制器输出都带抖了。原因遗忘因子设太小比如0.9以下模型只记住最近十几个点的信息此时如果输入信号恰好处于低激励时段信息矩阵接近奇异参数估计的方差会急剧放大。这是个放大器效应遗忘因子越大对噪声越不敏感但跟踪能力越差太小则反过来。我见过某装置的在线辨识迫降旁路阀抖动排查出来就是λ从0.98被运维人员改成了0.92。解决不要单独调λ要同时观察预测误差e_k和输入信号的活跃程度。工程上我通常的做法是做一个“条件更新”只有当输入变化超过某个阈值时才执行递推更新否则只保持参数不变。这样既保留了遗忘因子的跟踪能力又避开了低激励时段的参数漂移。此外可以把λ设成0.97到0.995之间别低于0.95。6. 再进一步辅助变量法与阶次辨识让最小二乘类工具更耐打6.1 用独立数据段验证辨识结果习惯比算法更重要做完辨识之后我建议你养成一个固定习惯把数据从中间切开用前60%做辨识后40%做验证。验证的做法很简单把辨识出的参数代入模型用后40%的输入重建输出序列再与真实输出对比计算验证集上的预测误差。这个方法成本极低却能暴露两类问题一是过拟合训练集上R²很高、验证集上一塌糊涂说明阶次偏高或数据量不够二是时变如果前段辨识的参数在后段明显失配说明系统参数随时间漂移了这时就该考虑RLS加遗忘因子而不是抱着批量解不放。% 用后40%数据做验证 N_val round(N * 0.4); Y_val y(N - N_val 1 : N); U_val u(N - N_val 1 : N); y_pred zeros(N_val, 1); for k 3:N_val y_pred(k) -a1_hat*y_pred(k-1) - a2_hat*y_pred(k-2) ... b1_hat*U_val(k-1) b2_hat*U_val(k-2); end val_err Y_val - y_pred; fprintf(验证集残差标准差: %.4f\n, std(val_err));注意这里验证时用的是预测输出y_pred自身的过去值而不是验证集里真实的y过去值这样才能真正检验模型的动态预测能力。如果只用真实y代入那验证的就是一步预测能力两者差距很大。这一步的代码很短但我认为它是整套辨识流程里最能拦下废模型的关口。6.2 辅助变量法有色噪声下的5行代码补救当残差检验确认噪声有色后最快的补救手段是辅助变量法IV。核心思路是用一组与噪声不相关的“工具变量”Z替代回归矩阵Φ做投影。常见的工具变量选择是输入的延迟序列因为输入u与外生白噪声不相关而u的延迟与Φ中的信号成分强相关。实现代码非常短构造一个Z矩阵再算一次投影即可。% 辅助变量法使用输入延迟作为工具变量 Z zeros(N - L, na nb); for k L1:N Z(k-L, :) [u(k-1) u(k-2) u(k-3) u(k-4)]; % 工具变量 end theta_iv (Z * Phi) \ (Z * Y); % IV估计这段代码能工作的前提是Z与Φ列空间有足够的重合度同时Z与噪声不相关。用u的延迟做工具变量在开环数据下这两个条件都容易满足。替换掉原来的Phi \ Y之后有偏问题会明显缓解。代价是方差比普通最小二乘略大这是用偏差换方差的经典权衡。如果你有系统辨识工具箱直接用辅助变量法命令会更方便但手写这段代码能让你真正理解它在干什么。6.3 阶次选择在最小二乘框架里用信息准则定阶模型阶次na、nb选错是辨识失败的重要原因但选多少其实可以由数据告诉你。信息准则的做法是对一组候选阶次分别做批量最小二乘计算残差方差V再按下面的公式打分取最小值对应的阶次。% 阶次辨识AIC / BIC 信息准则 maxOrder 4; AIC_mat zeros(maxOrder, maxOrder); BIC_mat zeros(maxOrder, maxOrder); for na_t 1:maxOrder for nb_t 1:maxOrder L_t max(na_t, nb_t); Phi_t zeros(N - L_t, na_t nb_t); for k L_t1:N Phi_t(k-L_t, :) [-y(k-1:-1:k-na_t) u(k-1:-1:k-nb_t)]; end Y_t y(L_t1:N); theta_t Phi_t \ Y_t; resid_t Y_t - Phi_t * theta_t; V_t resid_t * resid_t / length(Y_t); npar na_t nb_t; AIC_mat(na_t, nb_t) length(Y_t) * log(V_t) 2 * npar; BIC_mat(na_t, nb_t) length(Y_t) * log(V_t) npar * log(length(Y_t)); end end [minAIC, idxAIC] min(AIC_mat(:)); [minBIC, idxBIC] min(BIC_mat(:)); [na_AIC, nb_AIC] ind2sub(size(AIC_mat), idxAIC); [na_BIC, nb_BIC] ind2sub(size(BIC_mat), idxBIC);这里有个使用前提噪声必须是白噪声否则V的估计本身有偏信息准则会偏好高阶模型。所以在定阶之前先确认你已经处理好噪声相关性问题。另一个值得注意的坑是仿真数据里噪声设成0时V会小到让log(V)变成负无穷信息准则直接失效。我自己的调试经历里不止一次被这个问题坑过所以仿真时噪声方差不要低于1e-6。把这一段跑通后你手头的最小二乘类参数辨识工具链就齐了批量LS打底、RLS应对在线、IV处理有色噪声、信息准则定阶再加上交叉验证习惯这套组合能覆盖绝大多数线性系统辨识需求。希望你也能少走几步弯路一次就把参数辨识到位。本文还有配套的精品资源点击获取