matlab-identify-linear-system

线性模型辨识

使用 MATLAB System Identification Toolbox 根据测量数据估计线性动态模型。此技能会选择合适的模型类型、确定模型阶次、估计参数并验证结果，遵循 System Identification Toolbox 专家的方法。

适用场景

• 根据输入/输出数据辨识传递函数、状态空间模型或过程模型
• 根据测量数据确定模型阶次
• 拟合参数化模型，用于仿真、预测或控制设计
• 将频率响应数据（FRD）转换为参数化模型
• 比较模型结构（ARX、状态空间与传递函数）
• 根据时域数据确定频率响应
• 确定用于 PID 调优或控制设计的被控对象模型
• 当无法或不便对 Simulink 模型进行线性化时，获取数据驱动的线性模型
• 使用数据调整基于物理原理的模型（灰箱模型）参数
• 比较多个模型，以确定哪个模型最符合数据
• 使用辨识出的模型仿真或预测系统响应
• 对高阶系统或 MIMO 系统执行子空间辨识，或使用 Eigenvalue Realization Algorithm（ERA）
• 从频率响应中提取模态参数（固有频率、阻尼比、振型）
• 比较模型结构（ARX、状态空间与传递函数）
• 通过分析输入与输出信号之间的相关性，研究数据中是否存在反馈
• 研究输入信号的持续激励性，确保数据足以支持模型辨识

不适用场景

• 系统本质上是非线性的，不适合使用线性模型时
• 可用数据不足或质量较差，无法可靠地辨识模型时
• 主要目标是辨识非线性模型（例如神经状态空间、NLARX、Hammerstein-Wiener）时
• 使用实验数据估计 Simulink 模型参数时；请改用 Simulink Design Optimization Toolbox
• 设计控制器时；辨识被控对象后，请使用 Control System Toolbox 技能
• 信号处理（滤波、未拟合模型的频谱分析）——请使用信号处理技能

执行策略

**编写并运行一个端到端的 MATLAB 脚本。**不要一次只通过一次工具调用逐步执行各阶段。脚本应：

1. 创建/加载数据，并将其拆分为估计集和验证集
2. 估计延迟、选择结构并估计模型
3. 使用留出数据通过仿真进行验证
4. 打印结果

只有在首个脚本运行失败或结果较差（拟合度 < 70%）时，才拆分为多个步骤。

每个脚本的关键规则（强制要求——违反任意规则都属于错误）：

• 使用阶次向量时，设置 InteractiveOrderSelection=false —— 这是 ssestOptions 和 n4sidOptions 上的隐藏属性（不会显示在 disp() 或选项卡补全中）。必须显式设置，否则 GUI 弹窗会导致执行停滞
• 在任何搜索循环（阶次扫描、延迟扫描）期间设置 EstimateCovariance=false —— 对已舍弃模型计算协方差会浪费时间
• 用于仿真/控制时，在 ssestOptions、n4sidOptions 和 procestOptions 上设置 Focus='simulation'（tfestOptions 不支持此选项——tfest 没有 Focus）
• 多模型比较会返回 CELL：[~, fits] = compare(zv, m1, m2); fits{1}, fits{2} —— 始终在一次 compare() 调用中传入 2 个或更多模型，绝不针对每个模型分别调用 compare()
• 如果输入是平滑模拟信号，则在连续时间（CT）估计之前设置 data.InterSample = 'foh'
• 对于多输入 InterSample，使用列单元格 {'zoh'; 'foh'}（不要使用行单元格）
• 延迟优先：任何模型估计之前，都必须调用 delayest 或检查脉冲响应（MIMO 也一样，即使延迟看似很小）
• 延迟试探：绝不要只相信单次延迟估计——始终尝试 nk 及 nk±1，比较拟合度并选出最佳值
• 阶次范围：选择阶次时，使用范围（向量）而非单个整数——使用 ssest(ze, 2:8, opt)，不要使用 ssest(ze, 4, opt)

参数

用户提供：$ARGUMENTS

解析：

• project_name（可选）：projects/ 下包含 SPEC.md 的项目名称
• data_source（可选）：.mat 文件路径、工作区中的变量名，或对数据的内联描述

如果两者都未提供，请让用户指定数据源或描述辨识问题。

---

设计原则

1. 从简单开始，仅在数据证明有必要时增加复杂度。 优先尝试 2–4 阶，而不是 10–15 阶。
2. 先确定延迟。 错误的延迟无法通过提高阶数来修正，后果严重。
3. 正确设置 Focus。 这是最常被忽略的选项。默认的 'prediction' 有时不适用于仿真用途。
4. 始终留出验证数据。 绝不要将训练拟合度作为性能报告。
5. 正则化优于高阶。 正则化的 ARX(30) 往往优于未正则化的 ARX(5)。使用 arxRegul 或 ssregest。
6. 状态空间是默认选择。 不确定时，ssest 可处理 MIMO 和 CT/DT，且只需要阶数 n。
7. 比较 2–3 种结构。 第一个模型很少是最好的。
8. 检查残差。 拟合度高但残差相关，意味着模型遗漏了动态特性。
9. 知道何时停止。 验证数据上的拟合度 >90%，且残差为白噪声，即为成功。

---

快速路径 — 问题明确时使用

如果问题直接对应以下某种模式，立即编写单个脚本：

阶跃／脉冲响应 → 过程模型（不要拆分单次瞬态数据）：

% Step data is one transient — splitting creates IC discontinuity. Use full data.
opt = procestOptions('Focus', 'simulation');
m1 = procest(data, "P1D", opt); m2 = procest(data, "P2D", opt);
[~, fits] = compare(data, m1, m2); fprintf('P1D: %.1f%%, P2D: %.1f%%\n', fits{:});
fprintf('K=%.2f, Tp=%.1f, Td=%.1f\n', m1.Kp, m1.Tp1, m1.Td);

SISO 时域 → 传递函数：

ze = data(1:floor(end*0.7)); zv = data(floor(end*0.7)+1:end);
nk = delayest(ze);
% Hedge delay: try nk-1, nk, nk+1
delays = max(1, nk + (-1:1));
opt = ssestOptions('Focus', 'simulation', InteractiveOrderSelection=false, EstimateCovariance=false);
models = cell(1, numel(delays));
for i = 1:numel(delays)
    models{i} = tfest(ze, 3, 1, delays(i)*ze.Ts);
end
[~, fits] = compare(zv, models{:}); fprintf('Delay hedge fits: '); fprintf('%.1f%% ', fits{:}); fprintf('\n');
[~, best] = max(cell2mat(fits)); nk_best = delays(best);
% Final estimation with best delay
m1 = tfest(ze, 2, 0, nk_best*ze.Ts); m2 = tfest(ze, 3, 1, nk_best*ze.Ts);
m3 = ssest(ze, 2:6, opt);
[~, fits] = compare(zv, m1, m2, m3); fprintf('Fits: %.1f%%, %.1f%%, %.1f%%\n', fits{:});

MIMO → 状态空间：

ze = data(1:floor(end*0.7)); zv = data(floor(end*0.7)+1:end);
nk = delayest(ze);  % delay-first, even for MIMO
opt = ssestOptions('Focus', 'simulation', InteractiveOrderSelection=false, EstimateCovariance=false);
sys = ssest(ze, 1:10, opt);  % order RANGE, not single integer
[~, fit] = compare(zv, sys); fprintf('Fit: %.1f%%\n', fit);
% For MIMO bandwidth, use per-channel: bandwidth(sys(i,j))
for i = 1:size(sys,1), for j = 1:size(sys,2)
    fprintf('BW(%d,%d)=%.2f rad/s\n', i, j, bandwidth(sys(i,j)));
end, end

FRD／大型周期数据 → 频域路径：

opt = ssestOptions('InitializeMethod', 'AAA', 'Focus', 'simulation', ...
    InteractiveOrderSelection=false, EstimateCovariance=false);
opt.SearchOptions.MaxIterations = 0;
sys = ssest(Gfrd, 1:maxOrder, opt);

如果快速路径的拟合度 > 85%，即可完成。报告结果并继续。

---

深入路径 — 适用于情况不明确或首次尝试失败

快速路径失败（拟合度 < 70%）、问题不明确，或用户要求深入分析时，使用此结构化调查流程。

问题特征

确定：SISO／MIMO、时域／频域、预期用途（仿真／预测／控制）、已知约束。预处理详情请参见 references/data-preparation.md。

非线性检查（继续／停止关卡）

在确定采用线性辨识之前，先确认线性模型是否适用。

方法	操作	解释
幅值依赖性	根据输入幅值不同的数据集估计模型	如果增益／动态特性随幅值变化 → 非线性
谐波分析	输入周期信号，检查输出频谱中是否存在偶次谐波	偶次谐波表明存在非线性
模型阶数递增	拟合 2、4、8、12、16 阶模型 — 绘制拟合度与阶数的关系图	即使阶数很高，拟合度仍较低并趋于平稳 → 非线性
残差结构	检查残差与 u 或 y 的幅值之间的关系	存在系统性模式 → 非线性
数据拆分测试	在前半部分上估计，在后半部分上验证，反之亦然（使用较低模型阶数，例如 2–4 阶，避免估计方差造成误报）	拟合度不对称 → 非平稳或非线性
ISNLARX	在 iddata 上使用 isnlarx 方法评估非线性程度

决策

• 如果非线性较弱（增益变化 <20%），继续进行线性 ID，但注明局限性
• 如果检测到强非线性，建议使用：Hammerstein-Wiener、NLARX 或神经状态空间
• 如果系统非平稳（随时间变化），可考虑分段估计或递归方法

数据准备

完整的预处理工作流请参见 references/data-preparation.md，包括：

• 加载和检查（advice、plot）
• 预处理检查清单（偏移量、缺失数据、异常值、非均匀采样）
• 预滤波（带通、频率加权）
• CT 模型的 InterSample 行为
• 训练／验证集拆分
• 频域转换
• 数据质量风险信号
• 数据中输出反馈的可能性（checkFeedback）
• 持续激励检查（pexcit）

模型类型选择

应用此决策树。推荐采用第一个匹配的分支：

1. Physical structure known (ODEs with unknown parameters)?
   --> idgrey + greyest (outside this skill's scope)

2. Frequency-domain data (idfrd), very large dataset (N > 50k), periodic input, or high modal density?
   --> Frequency-domain path:
       - ssest with InitializeMethod='AAA' (SISO/SIMO/MISO/MIMO — only option for full MIMO FRD)
       - ssest with InitializeMethod='lsrf' (SISO/SIMO/MISO only — vector fitting)
       - tfest on idfrd/etfe/spa data (uses lsrf internally; SISO/SIMO/MISO only)

3. Low-order process (1-3 poles, <=1 zero, with gain+delay)?
   --> idproc + procest

4. SISO, continuous-time, moderate complexity (np <= 10)?
   --> idtf + tfest

5. MIMO, or high-order, or "just need a good model quickly"?
   --> idss + ssest (with n4sid for initialization)

6. Need explicit noise model (prediction/filtering application)?
   --> Polynomial models: ARX, IV4, ARMAX, OE, BJ

7. Time-series (no input, output only)?
   --> ar() for AR, or ssest with nu=0 for state-space

有关过程模型、多项式模型以及各自适用场景的详细指南，请参阅 references/model-structures.md。

阶次确定

有关以下方法，请参阅 references/order-determination.md：

• 延迟估计（delayest、脉冲响应）
• ARX 结构搜索（arxstruc、selstruc）
• 子空间阶次选择（使用阶次范围的 n4sid）
• 迭代复杂度（传递函数阶次递增）
• 过程模型阶次递增
• 频域路径（AAA 初始化）

经验法则：

• 最大有效阶次：n_max ~ min(N/20, 30)
• MIMO 状态空间：从 n = max(ny, nu) * 2 开始，最高到 5
• 如果 ARX(10) 和 ssest(4) 的拟合效果相近，优先选择 ssest(4)
• 每增加一个参数，改进幅度低于 2% 时停止增加阶次

估计

有关完整的估计工作流，请参阅 references/estimation.md，其中包括：

• 关键选项（Focus、InitialState、Stability、Regularization、OutputWeight、WeightingFilter、ErrorThreshold）
• ssest 和 tfest 的初始化方法
• ssest 优于 n4sid 的原因
• 频域估计路径（AAA、lsrf）
• 所有模型类型的估计命令
• 正则化 ARX 到状态空间（ssregest）
• 带参数边界的结构化估计
• 多模型策略

关键提醒：

% Focus — available on ssestOptions, n4sidOptions, arxOptions, etc. (NOT tfestOptions)
opt = ssestOptions('Focus', 'simulation');   % for simulation/control
opt = ssestOptions('Focus', 'prediction');   % for forecasting
% NOTE: tfest does NOT have a Focus option. For time-domain data, tfest always
% produces a stable model. Use WeightingFilter for frequency emphasis with tfest.

% Disable interactive order selection when using an order vector
opt = n4sidOptions(InteractiveOrderSelection=false);
% or: opt = ssestOptions(InteractiveOrderSelection=false);

% Continuous-time estimation — set 'Ts',0 for ssest; tfest is CT by default
opt = ssestOptions('Focus', 'simulation', InteractiveOrderSelection=false);
model_ct = ssest(data, n, 'Ts', 0, opt);
model_ct = tfest(data, np, nz);  % CT by default from sampled data
data.InterSample = 'foh';        % set BEFORE estimation for smooth analog inputs

% Regularization for high-order ARX
[Lambda, R] = arxRegul(data, orders, arxRegulOptions('RegularizationKernel', 'TC'));

---

验证（必须执行）

使用验证数据进行比较（必须执行）

重要：始终通过仿真（无限预测时域）进行验证。使用 Focus='prediction' 估计的模型，即使其动态特性有误，也可能表现出极佳的一步超前拟合效果。

% Validate by simulation (default of compare)
[yhat, fit] = compare(zv, model);
fprintf('Validation fit (simulation): %.1f%%\n', fit);

% Multi-model comparison — fit is a CELL ARRAY, not a numeric vector
[yhat, fits] = compare(zv, m1, m2, m3);
fprintf('Fits: %.1f%%, %.1f%%, %.1f%%\n', fits{:});
% Extract as numeric vector: cell2mat(fits)

% 1-step prediction fit (for forecasting models ONLY)
[yhat_pred, fit_pred] = compare(zv, model, 1);
% WARNING: fit_pred >> fit_sim means the noise model is doing the heavy lifting

残差分析

% Programmatic residual analysis (no plots — suitable for batch/agent mode)
[e, r] = resid(zv, model);
% e = residual iddata object
% r = 3D array [M x nz x nz] where nz = ny + nu, M = number of lags (26 default)
%   r(:,1:ny,1:ny)       = residual autocovariance (RAW, not normalized)
%   r(:,ny+1:end,1:ny)   = cross-covariance between input and residual
% Quick whiteness check (SISO: ny=1, nu=1, nz=2):
acf = r(:,1,1) / r(1,1,1);  % normalize by lag-0 to get correlation
N = size(zv.y, 1);
conf99 = 2.58 / sqrt(N);    % 99% confidence bound
is_white = all(abs(acf(2:end)) < conf99);
fprintf('Residuals white: %s (99%% bound = %.4f)\n', string(is_white), conf99);

% Cross-correlation: input-residual (SISO)
xcf = r(:,2,1) / sqrt(r(1,1,1) * r(1,2,2));  % normalized cross-covariance
is_uncorr = all(abs(xcf) < conf99);
fprintf('Residuals uncorrelated with input: %s\n', string(is_uncorr));

残差解读指南：

自相关（acf）	互相关（xcf）	诊断	操作
白噪声	不相关	模型足够	完成
在多个滞后处显著	不相关	噪声模型不足，但用于仿真的对象模型可能足够	增加噪声模型阶次（ARMAX/BJ 中的 C/D）；对象 G 仍可使用
白噪声	在滞后 k 处显著	缺少滞后 k 处的输入动态	添加回归项 u(t-k)：增加 nb 或调整 ARX/ARMAX/BJ 中的 nk
在多个滞后处显著	在滞后 k 处显著	对象模型和噪声模型都不够	增加对象模型阶次和噪声阶次；检查延迟

拟合效果解读

拟合率	结论	后续操作
> 95%	极佳	完成 — 报告结果
85-95%	良好	可以接受；尝试一种替代方案进行确认
70-85%	中等	增加阶次、尝试不同结构、检查数据
50-70%	较差	结构错误、缺少非线性，或数据质量差
< 50%	失败	重新评估基本假设（延迟？反馈？非线性？）

模型不够充分时的重新估计方案：

% Scan an order RANGE (never just guess one number)
opt = ssestOptions('Focus', 'simulation', InteractiveOrderSelection=false, EstimateCovariance=false);
model_new = ssest(ze, 2:8, opt);
% Compare old and new with multi-model compare (MANDATORY pattern)
[~, fits] = compare(zv, model_old, model_new);
fprintf('Old: %.1f%%, New: %.1f%%\n', fits{:});

有关不确定性分析、稳定性评估、诊断清单和初始条件指南，请参阅 references/validation-and-diagnostics.md。

---

结果和产物

生成输出

1. 模型摘要 — 类型、阶次、拟合百分比、关键参数
2. 对比图 — 验证数据上的测量值与仿真值
3. Bode 图 — 带置信界限的频率响应
4. 残差图 — 自相关和互相关

% Comparison plot
[yhat, fit] = compare(zv, model);
title(sprintf('Validation: %.1f%% fit', fit));

% Bode with confidence
h = bodeplot(model);
showConfidence(h, 3);

% Residuals
resid(zv, model);

结果摘要

================================================================
  Linear Model Identification — Results
================================================================
  Model type:     <idss / idtf / idpoly / idproc>
  Order:          <n / [na nb nk] / np poles, nz zeros>
  Delay:          <nk samples (X seconds)>
  Focus:          <simulation / prediction>

  Validation fit: XX.X% (NRMSE on held-out data)
  FPE:            <value>
  AIC:            <value>

  Key dynamics:
    Poles:        <dominant pole locations>
    Zeros:        <zero locations if few>
    DC gain:      <value>
    Bandwidth:    <-3dB frequency (SISO only; for MIMO use bandwidth(model(i,j)) per channel)>

  Status: FIT ACHIEVED / BELOW TARGET — <recommendation>
================================================================

保存产物

1. MATLAB 脚本 — identify_model.m，包含完整且可复现的工作流
2. 模型 MAT 文件 — identified_model.mat，包含最终模型对象
3. 图表 — 用于对比、Bode 图和残差的 PNG 文件
4. 结果摘要 — 输出到控制台

---

特殊主题

有关高级场景，请参阅：

• references/special-topics.md — 闭环 ID、面向控制的 ID、MIMO、连续时间、物理信息约束、频域局限性、递归估计、tfest 与 tfestimate

---

常见问题 — 快速参考

#	错误	后果	修正方法
1	在训练数据上验证	未能发现过拟合	始终留出验证集
2	使用预测拟合度判断仿真质量	对动态特性盲目自信	通过仿真验证（horizon=Inf）
3	InterSample 设置错误	CT 模型在高频处产生系统性偏差	对平滑模拟输入设置 'foh'
4	忽略 Focus 选项	模型优化了错误的准则	对仿真/控制用途设置 Focus='simulation'
5	未去除数据中的偏移趋势	DC 增益错误，整体拟合效果差	使用 detrend(data) 或偏移量
6	在闭环中使用 ARX	植物估计有偏	使用 iv4、BJ 或间接方法
7	假设了错误的延迟	后果严重——无论阶数多高都无济于事	先估计延迟，再尝试 +/-1
8	参数过多	训练拟合效果很好，但泛化能力差	正则化或降低阶数
9	忽略不确定性	对模型精度过于自信	始终检查 Bode 图上的 showConfidence
10	假设训练数据的初始条件适用于新数据	验证拟合效果差	对每个新数据集重新估计初始条件
11	对 DT 模型使用 ssest、tfest 时，未将 'Ts',Value（Value>0）作为名称-值对传入	得到 CT 模型	ssest(data, n, 'Ts', data.Ts, opt)
12	仿真时使用了错误的初始条件	验证数据拟合效果差	使用 findstates' 或 data2state' 确定能使验证数据拟合效果最佳的初始条件
13	忘记缩放数据	导致病态的辨识问题，结果不佳	确保输入、输出和时间单位都经过适当缩放，尤其是在使用数值优化算法时。

---

版权所有 2026 The MathWorks, Inc.

---