matlab-identify-linear-system

線形モデルの同定

MATLAB System Identification Toolbox を使用して、測定データから線形動的モデルを推定します。このスキルでは、System Identification Toolbox の専門家が用いる手法に従い、適切なモデルの種類を選択し、モデル次数を決定し、パラメーターを推定して、結果を検証します。

使用する場面

• I/O データから伝達関数、状態空間モデル、またはプロセスモデルを同定する
• 測定データからモデル次数を決定する
• シミュレーション、予測、または制御設計に用いるパラメトリックモデルを当てはめる
• 周波数応答データ（FRD）をパラメトリックモデルに変換する
• モデル構造を比較する（ARX と状態空間モデルと伝達関数）
• 時間領域データから周波数応答を求める
• PID 調整や制御設計のためにプラントモデルを決定する
• Simulink モデルの線形化が不可能または実用的でない場合に、データ駆動型の線形モデルを取得する
• 物理ベースモデル（グレーボックス）のパラメーターをデータで調整する
• 複数のモデルを比較して、データに最もよく適合するモデルを決定する
• 同定したモデルを使用してシステム応答をシミュレーションまたは予測する
• 高次システムや MIMO システムに対して部分空間同定を実行する。または固有値実現アルゴリズム（ERA）を使用する
• 周波数応答からモードパラメーター（固有振動数、減衰比、モード形状）を抽出する
• モデル構造を比較する（ARX と状態空間モデルと伝達関数）
• 入力信号と出力信号の相関を分析し、データにおけるフィードバックの可能性を調べる
• 入力信号の持続励振性を調べ、モデル同定に十分な情報をデータが含んでいることを確認する

使用しない場面

• システムが本質的に非線形で、線形モデルが適切でない場合
• 信頼性の高いモデル同定に必要なデータが不足している、または品質が低い場合
• 主な目的が非線形モデル（ニューラル状態空間モデル、NLARX、Hammerstein-Wiener など）の同定である場合
• 実験データを使用して Simulink モデルのパラメーターを推定する場合。代わりに Simulink Design Optimization Toolbox を使用してください
• コントローラーを設計する場合。プラントを同定した後、Control System Toolbox のスキルを使用してください
• 信号処理（モデル当てはめを伴わないフィルタリング、スペクトル解析）—信号処理のスキルを使用してください

実行戦略

エンドツーエンドの MATLAB スクリプトを 1 つ作成して実行してください。 フェーズごとにツールを 1 回ずつ呼び出して進めてはいけません。スクリプトでは次の処理を行います。

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} — 必ず 2 つ以上のモデルを 1 回の compare() 呼び出しに渡してください。モデルごとに compare() を個別に呼び出してはいけません
• 入力が滑らかなアナログ信号の場合、連続時間（CT）の推定前に data.InterSample = 'foh' を設定してください
• 複数入力の InterSample には、列方向のセル配列 {'zoh'; 'foh'} を使用してください（行方向のセル配列ではありません）
• 遅延を先に調べる：モデル推定の前に、必ず delayest を呼び出すかインパルス応答を調べてください（MIMO の場合も、遅延が小さく見える場合も同様です）
• 遅延に幅を持たせる：遅延推定値を 1 つだけ信用してはいけません。常に nk と nk±1 を試して適合度を比較し、最良のものを選んでください
• 次数範囲：次数を選択する場合は、単一の整数ではなく範囲（ベクトル）を使用してください — ssest(ze, 2:8, opt) であり、ssest(ze, 4, opt) ではありません

引数

ユーザーが指定するもの：$ARGUMENTS

以下を解析します：

• project_name（省略可）：projects/ 配下にある、SPEC.md を含むプロジェクト名
• data_source（省略可）：.mat ファイルへのパス、ワークスペース内の変数名、またはデータの説明文

どちらも指定されていない場合は、データソースを指定するか、同定したい問題について説明するようユーザーに尋ねます。

---

設計原則

1. まずはシンプルに始め、データが必要性を示す場合にのみ複雑さを加えます。 10～15の前に、2～4を試します。
2. 遅延を先に決めます。 遅延が間違っていると、高次にしても修正できず、致命的な問題になります。
3. Focus を正しく設定します。 最も見落とされやすいオプションです。デフォルトの 'prediction' は、シミュレーション用途では適切でない場合があります。
4. 必ず検証データを別に確保します。 学習データへの適合度を性能として報告してはいけません。
5. 高次化より正則化を優先します。 正則化した ARX(30) は、正則化していない ARX(5) を上回ることがよくあります。arxRegul または ssregest を使用します。
6. 状態空間モデルをデフォルトとします。 迷った場合は、MIMO、CT/DT に対応し、次数 n のみを必要とする ssest を使用します。
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 を参照してください。

非線形性の確認（STOP/GO ゲート）

線形同定を進める前に、線形モデルが適切であることを確認します。

方法	手順	解釈
振幅依存性	入力振幅の異なるデータセットからモデルを推定します	振幅によってゲイン／ダイナミクスが変化する場合 → 非線形
高調波解析	周期入力を適用し、出力スペクトルに偶数次高調波があるか確認します	偶数次高調波は非線形性を示します
モデル次数の増加	次数 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) を優先します
• パラメーターを1つ追加したときの改善が 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' で推定したモデルは、ダイナミクスが誤っていても、1ステップ先予測の適合度が非常に高くなる場合があります。

% 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) を追加します。ARX/ARMAX/BJ の nb を増やすか nk を調整します
ラグで有意	ラグ k で有意	プラントモデルとノイズモデルの両方が不十分です	モデル次数とノイズ次数の両方を増やし、遅延を確認します

適合度の解釈

適合度 %	評価	次の対応
> 95%	非常に良好	完了 — 結果を報告します
85-95%	良好	許容範囲です。確認のため別の方法を1つ試します
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. ボード線図 — 信頼境界を含む周波数応答
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. 図 — 比較、ボード線図、残差の 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 モデルに対して、名前と値のペアとして 'Ts',Value（Value>0）を ssest、tfest に渡さない	代わりに CT モデルが得られる	ssest(data, n, 'Ts', data.Ts, opt)
12	シミュレーションに誤った初期条件を使用する	検証データへのフィットが悪くなる	findstates' または data2state' を使用して、検証データへのフィットが最大になる初期条件を特定する
13	データのスケーリングを忘れる	条件の悪い同定問題となり、結果が悪くなる	特に数値最適化アルゴリズムを使用する場合は、入力、出力、時間単位が適切にスケーリングされていることを確認する。

---

Copyright 2026 The MathWorks, Inc.

---