简介：这份MATLAB程序包面向自动化、控制工程及信号处理方向的学习者与研究人员，用于系统辨识课程的实验与算法验证。包内共7个文件，全部为.m脚本，压缩包约7KB，涵盖普通最小二乘、递归最小二乘、扩展递归最小二乘、广义最小二乘以及基于神经网络的参数辨识等多类实现，可对照比较批处理与在线学习在参数估计上的差异。其中递归类脚本侧重实时数据流下的参数更新与跟踪，广义最小二乘部分则涉及异方差或相关误差下的加权处理，神经网络脚本展示了用最小二乘优化权值以拟合非线性系统的思路。目前已有2948人学习下载，适合希望从基础公式推导过渡到代码落地、并借助现成脚本快速复现实验的读者参考。
1. 最小二乘法系统辨识：从一堆采样点到能用的模型，中间差了什么
如果你手上有两组数据——一组是输入信号，一组是系统响应——想反推出这个系统的传递函数或者差分方程，最小二乘法系统辨识就是最直接的路径。它不依赖复杂的优化器，不需要神经网络调参那套玄学，核心就是解一个矩阵方程。很多人第一次接触系统辨识，都是从最小二乘法 matlab 程序开始的：几行代码就能跑出一个模型，但跑出来的结果能不能用、参数有没有物理意义、阶次选得对不对，这些才是真正区分“跑通”和“能用”的地方。
这篇文章面向的是需要在 MATLAB 里完成系统辨识的工程师和学生。你可能有实测的输入输出数据，也可能是在做仿真验证，目标都是拿到一个可以预测、可以控制的模型。我会从最小二乘法的数学形式讲起，然后落到 MATLAB 的具体实现，包括数据预处理、阶次选择、参数求解、模型验证，最后讲几个我踩过的坑。读完你应该能自己写出一套可复用的辨识脚本，而不是每次都在网上找一段代码改改就用。
2. 最小二乘法辨识的数学形式与 MATLAB 落地路径
2.1 从差分方程到回归向量：最小二乘到底在解什么
系统辨识的最小二乘法，本质上是把一个差分方程整理成线性回归的形式。假设系统可以用 ARX 模型描述：
y(k) + a₁·y(k-1) + … + aₙ·y(k-n) = b₁·u(k-1) + … + bₘ·u(k-m) + e(k)
把未知参数移到一边，写成向量形式：
y(k) = φᵀ(k)·θ + e(k)
其中 θ = [-a₁, …, -aₙ, b₁, …, bₘ]ᵀ 是待估参数，φ(k) = [y(k-1), …, y(k-n), u(k-1), …, u(k-m)]ᵀ 是回归向量。对于 N 个采样点，把所有方程堆起来：
Y = Φ·θ + E
最小二乘的解就是让残差平方和最小的 θ：
θ̂ = (ΦᵀΦ)⁻¹·ΦᵀY
这个公式看起来简单，但直接求逆在数值上不稳定。MATLAB 里实际用的是左除运算符或者 QR 分解，避免显式求逆。我一般会优先用\，因为 MATLAB 会根据矩阵结构自动选择最稳的解法。
理解了这个形式，你就知道为什么数据长度 N 必须大于参数个数 n+m，否则 ΦᵀΦ 奇异，解不出来。另外，输入信号必须持续激励，否则某些参数不可辨识——这是后面避坑章节要展开的。
2.2 用 MATLAB 构造回归矩阵：一个可复用的函数
下面这个函数把输入输出数据转成回归矩阵和目标向量。我把它写成独立函数，方便在不同项目里复用。
function [Phi, Y] = build_regressor(u, y, na, nb) % build_regressor 构造ARX模型的回归矩阵 % 输入： % u - 输入信号向量（列向量） % y - 输出信号向量（列向量） % na - 输出滞后阶数 % nb - 输入滞后阶数 % 输出： % Phi - 回归矩阵，每行对应一个采样时刻 % Y - 目标向量，对应Phi每行的输出值 u = u(:); y = y(:); N = length(y); start_idx = max(na, nb) + 1; % 保证回归向量完整的最早时刻 Phi = zeros(N - start_idx + 1, na + nb); Y = zeros(N - start_idx + 1, 1); for k = start_idx:N row = k - start_idx + 1; % 输出滞后项：y(k-1) ... y(k-na) for i = 1:na Phi(row, i) = -y(k - i); end % 输入滞后项：u(k-1) ... u(k-nb) for j = 1:nb Phi(row, na + j) = u(k - j); end Y(row) = y(k); end end逻辑说明：这个函数把每个采样时刻的滞后项拼成一行。注意输出滞后项前面加了负号，因为 ARX 模型的标准形式是 y(k) + a₁y(k-1) + … = b₁u(k-1) + …，移项后 y(k) = -a₁y(k-1) - … + b₁u(k-1) + …，所以回归向量里 y 的滞后项带负号，这样解出来的 θ 直接就是 [-a₁, …, b₁, …]。
参数说明：na和nb决定了模型阶次。na是输出自回归阶数，nb是输入延迟阶数。起始索引start_idx取max(na, nb) + 1，保证第一个回归向量的所有滞后项都有数据。如果你的系统有纯延迟 d，可以把输入项改成u(k-d-j+1)，这个后面讲。
2.3 参数求解与模型验证：别只看拟合曲线
拿到 Phi 和 Y 之后，求解就一行：
theta = Phi \ Y;但解出来只是开始。我一般会做三件事来验证模型质量。
第一，看残差。残差 e = Y - Phi*theta，如果残差有自相关，说明模型结构不对或者阶次不够。用autocorr(e)看一眼，残差应该在零附近随机波动。
第二，看参数置信区间。最小二乘的参数协方差估计是：
sigma2 = (e' * e) / (length(Y) - size(Phi,2)); cov_theta = sigma2 * inv(Phi' * Phi); std_theta = sqrt(diag(cov_theta));如果某个参数的绝对值小于两倍标准差，这个参数可能不显著，考虑降阶。
第三，做交叉验证。把数据分成两段，前一段辨识，后一段预测。预测输出和实际输出的拟合度用fitpercent或者归一化均方误差衡量。我习惯算 NRMSE：
y_pred = Phi * theta; nrmse = norm(Y - y_pred) / norm(Y - mean(Y)) * 100;如果 NRMSE 低于 80%，模型基本可用；低于 60% 就要检查数据质量或者模型结构了。
3. 阶次选择、数据预处理与递推最小二乘的实现细节
3.1 阶次选不对，后面全白费：AIC 与试凑的结合
阶次选择是系统辨识里最容易被忽视的一步。很多人直接设 na=2, nb=2 就跑，结果要么欠拟合要么过拟合。我一般用 AIC 准则先粗选，再结合物理意义微调。
AIC 的定义是：
AIC = N·ln(σ²) + 2·(na + nb)
其中 σ² 是残差方差，N 是样本数。阶次越高，σ² 越小，但第二项惩罚越大。把 na 和 nb 从 1 到 6 遍历，找 AIC 最小的组合。
best_aic = inf; best_na = 1; best_nb = 1; for na = 1:6 for nb = 1:6 [Phi, Y] = build_regressor(u, y, na, nb); theta = Phi \ Y; e = Y - Phi * theta; sigma2 = (e' * e) / length(Y); aic = length(Y) * log(sigma2) + 2 * (na + nb); if aic < best_aic best_aic = aic; best_na = na; best_nb = nb; end end end但 AIC 只是参考。如果被辨识对象是一个二阶机械系统，na=2 就有物理意义，即使 AIC 建议 na=3，我也会优先选 na=2。过拟合的模型在控制里会引入不必要的振荡。
3.2 数据预处理：去均值、去趋势、滤波的先后顺序
实测数据往往有直流偏置和低频漂移。不处理直接辨识，会得到一个带积分特性的模型，参数物理意义全乱。我一般按这个顺序处理：
第一步，去均值。输入输出都减去各自的均值，消除直流分量。
u = u - mean(u); y = y - mean(y);第二步，去趋势。如果数据有明显的线性漂移，用detrend去掉。
y = detrend(y);第三步，滤波。如果噪声带宽远高于系统带宽，可以低通滤波。但要注意滤波会引入相位延迟，如果输入输出用同一个滤波器，延迟可以抵消。我一般用零相位滤波filtfilt。
[b_f, a_f] = butter(2, 0.1); % 截止频率根据实际调整 u = filtfilt(b_f, a_f, u); y = filtfilt(b_f, a_f, y);顺序很重要：先去均值再去趋势，最后滤波。如果先滤波，直流分量会影响滤波器初始状态，引入额外瞬态。
3.3 递推最小二乘：在线辨识的 MATLAB 实现
离线最小二乘适合批量数据，但如果要做自适应控制或者在线监测，就需要递推形式。递推最小二乘（RLS）的公式如下：
θ(k) = θ(k-1) + K(k)·[y(k) - φᵀ(k)·θ(k-1)] K(k) = P(k-1)·φ(k) / [λ + φᵀ(k)·P(k-1)·φ(k)] P(k) = [I - K(k)·φᵀ(k)]·P(k-1) / λ
其中 λ 是遗忘因子，通常取 0.95 到 1.0。λ=1 时所有历史数据等权，适合时不变系统；λ<1 时旧数据权重衰减，适合时变系统。
function [theta_hist, y_pred] = rls_identify(u, y, na, nb, lambda, P0) % rls_identify 递推最小二乘辨识 % 输入： % u, y - 输入输出数据 % na, nb - 阶次 % lambda - 遗忘因子，0.95~1.0 % P0 - 初始协方差矩阵，通常取 1e4*eye(na+nb) % 输出： % theta_hist - 参数历史，每行一个时刻 % y_pred - 预测输出 u = u(:); y = y(:); N = length(y); n_params = na + nb; theta = zeros(n_params, 1); P = P0; theta_hist = zeros(N, n_params); y_pred = zeros(N, 1); start_idx = max(na, nb) + 1; for k = start_idx:N phi = zeros(n_params, 1); for i = 1:na phi(i) = -y(k - i); end for j = 1:nb phi(na + j) = u(k - j); end y_hat = phi' * theta; e = y(k) - y_hat; K = P * phi / (lambda + phi' * P * phi); theta = theta + K * e; P = (P - K * phi' * P) / lambda; theta_hist(k, :) = theta'; y_pred(k) = y_hat; end end参数说明：lambda越接近 1，参数收敛越稳但跟踪时变能力越弱；越接近 0.9，跟踪快但参数波动大。我一般先用 1.0 跑一遍看收敛，如果系统有时变特性再降到 0.98 左右。P0取大一点表示初始参数不确定性大，让算法快速调整；取小一点则初始收敛慢但更稳。
4. 避坑与排查：最小二乘辨识里那些让人翻车的细节
4.1 输入信号持续激励不足，参数解出来是“假”的
现象：辨识出来的模型拟合训练数据很好，但换一组数据预测就完全不对。参数值看起来也合理，但物理意义对不上。
原因：输入信号频谱太窄，比如只用一个频率的正弦波做输入，只能激励系统在该频率附近的动态。其他频段的参数没有被数据约束，最小二乘会给出一个在训练集上拟合好但泛化差的解。
解决：用伪随机二进制序列（PRBS）或者扫频信号作为输入。PRBS 的频谱覆盖宽，实现也简单。MATLAB 里可以用idinput生成：
u = idinput(1023, 'prbs', [0 0.5], [-1 1]);如果只能用实际运行数据，检查输入的自功率谱，确保在系统带宽范围内都有能量。
4.2 数据长度不够或者矩阵病态，解出来的参数大得离谱
现象：theta里某些参数是 10⁶ 量级，或者 MATLAB 警告“矩阵接近奇异”。
原因：数据长度 N 小于参数个数，或者输入输出高度相关导致 ΦᵀΦ 条件数极大。另一个常见原因是数据没有去均值，直流分量让矩阵第一列和常数项共线。
解决：先检查cond(Phi'*Phi)，如果大于 10⁸ 就要警惕。增加数据长度、去均值、降低阶次都是办法。如果实在无法增加数据，可以用正则化最小二乘（岭回归）：
lambda_reg = 1e-6; theta = (Phi' * Phi + lambda_reg * eye(size(Phi,2))) \ (Phi' * Y);正则化会引入偏差，但能换来数值稳定。lambda_reg从小往大试，看参数是否稳定。
4.3 纯延迟没处理，模型阶次虚高
现象：AIC 建议的 nb 比实际物理延迟大很多，辨识出来的 b 参数前面几个接近零。
原因：系统有纯延迟 d，但回归向量里输入项从 u(k-1) 开始，模型需要用额外的 b 参数来“凑”出这个延迟。这会导致阶次虚高，参数冗余。
解决：先估计延迟。可以用互相关函数找输入输出的最大相关滞后：
[c, lags] = xcorr(y, u, 'normalized'); [~, idx] = max(abs(c)); delay_est = lags(idx);然后在构造回归矩阵时，输入项从u(k-delay_est)开始。这样 nb 可以设得小一些，参数更简洁。
4.4 残差自相关严重，说明模型结构不对
现象：残差e的autocorr图在滞后几阶处有明显尖峰，超出置信区间。
原因：ARX 模型假设残差是白噪声。如果实际系统有未建模的动态，或者噪声是有色的，残差就会自相关。这时候最小二乘估计是有偏的。
解决：换模型结构。如果残差自相关，可以试试 ARMAX 或者输出误差（OE）模型。MATLAB 的 System Identification Toolbox 里有armax和oe函数。如果不想用工具箱，可以把残差也作为回归项加入，变成 ARARX 结构，但实现复杂度会上升。
4.5 训练集和验证集划分不当，过拟合被掩盖
现象：训练集 NRMSE 95%，验证集 NRMSE 只有 50%。
原因：数据划分时没有考虑系统的工作点变化。如果训练集只覆盖低负载工况，验证集包含高负载工况，模型外推能力差是必然的。
解决：按工况分层划分，确保训练集和验证集覆盖相似的工作范围。如果数据本身就不覆盖全工况，那模型只在该工况附近有效，不要指望外推。我一般会留出最后 20% 的数据做验证，并且检查两段数据的统计特性是否一致。
5. 从离线到在线：最小二乘辨识的进阶用法与验证习惯
5.1 用频域验证代替只看时域拟合
时域拟合好不代表模型频域特性对。我习惯在辨识完成后，用etfe或者spa估计非参数频率响应，再和辨识模型的bode图对比。如果在中频段有超过 3dB 的偏差，模型用于控制器设计时可能出问题。
% 非参数频率响应估计 [G_etfe, w] = etfe([y u], [], 1024); % 辨识模型的频率响应 sys_arx = tf(theta(nb+1:end)', [1 theta(1:na)'], 1); [G_arx, ~] = freqresp(sys_arx, w); % 对比幅频特性 semilogx(w, 20*log10(abs(G_etfe)), 'b', w, 20*log10(abs(squeeze(G_arx))), 'r--'); legend('ETFE', 'ARX模型'); xlabel('频率 (rad/s)'); ylabel('幅值 (dB)');这个对比能暴露很多时域看不出来的问题。比如高频段模型增益偏高，可能是噪声建模不当；低频段偏差大，可能是去趋势不彻底。
5.2 参数不确定性与蒙特卡洛验证
最小二乘给出的参数是一个点估计，但实际参数是有不确定性的。我一般会做蒙特卡洛验证：在残差里加不同随机种子，重新辨识 100 次，看参数分布。
n_runs = 100; theta_samples = zeros(n_runs, length(theta)); for run = 1:n_runs e_sim = e .* randn(size(e)); % 用残差尺度生成噪声 y_sim = Phi * theta + e_sim; theta_samples(run, :) = Phi \ y_sim; end % 看参数标准差 std_theta_mc = std(theta_samples);如果某个参数的蒙特卡洛标准差和解析协方差算出来的差很多，说明残差不是白噪声，解析公式不可靠。
5.3 一个我坚持了多年的习惯：先仿真验证再上实测数据
每次写新的辨识脚本，我一定先用仿真数据跑通。构造一个已知传递函数的系统，生成输入输出，加噪声，然后看辨识结果能不能还原参数。这个习惯帮我省了很多时间——如果仿真数据都辨识不准，实测数据更不可能。
% 仿真验证：已知系统 sys_true = tf([1 2], [1 3 2]); % 传递函数 u_sim = idinput(2000, 'prbs', [0 0.5], [-1 1]); y_sim = lsim(sys_true, u_sim, (0:length(u_sim)-1)'); y_sim = y_sim + 0.01 * randn(size(y_sim)); % 加噪声 % 辨识 [Phi_sim, Y_sim] = build_regressor(u_sim, y_sim, 2, 2); theta_sim = Phi_sim \ Y_sim; % 对比真实参数 % 离散化后的真实参数需要根据采样时间计算仿真验证能确认脚本逻辑正确、阶次选择合理、预处理没有引入偏差。这一步花十分钟，能避免后面几小时的排查。
5.4 模型交付前必做的三个检查
第一，检查参数是否在物理合理范围内。比如时间常数应该是正的，增益符号应该和实际系统一致。如果出现负的时间常数，模型肯定有问题。
第二，检查模型在训练数据之外的预测能力。用验证集算 NRMSE，低于 70% 就要谨慎。
第三，检查模型是否可稳定用于控制。如果辨识出的离散模型有单位圆外的极点，闭环控制可能不稳定。用abs(roots([1 theta(1:na)']))检查极点模值。
这三个检查做完，模型才能交给控制组或者用于后续设计。我见过太多人跳过这些步骤，直接把辨识结果拿去调 PID，最后振荡得没法收拾。希望帮到你。
本文还有配套的精品资源，点击获取
