Matlab Identify Linear System

Linear Model Identification

Estimate a linear dynamic model from measurement data using MATLAB System Identification Toolbox. This skill selects the right model type, determines model order, estimates parameters, and validates results — following the methodology a System Identification Toolbox expert would use.

When to Use

• Identify a transfer function, state-space, or process model from I/O data

• Determine model order from measurement data

• Fit a parametric model for simulation, prediction, or control design

• Convert frequency response data (FRD) to a parametric model

• Compare model structures (ARX vs state-space vs transfer function)

• Determine frequency response from time-domain data

• Determine a plant model for PID tuning or control design

• Obtain a data-driven linear model when linearization of a Simulink model is not possible or practical

• Tune parameters of a physics-based model (grey-box) using data

• Compare multiple models to determine which best fits the data

• Simulate or predict system response using the identified model

• Perform subspace identification for high-order systems or MIMO systems, or use Eigenvalue Realization Algorithm (ERA)

• Extract modal parameters (natural frequencies, damping ratios, mode shapes) from frequency response

• Compare model structures (ARX vs state-space vs transfer function)

• Study the possibility of feedback in data by analyzing the correlation between input and output signals

• Study persistence of excitation in the input signals to ensure that the data is informative enough for model identification

When NOT to Use

• When the system is inherently nonlinear and a linear model is not appropriate

• When the available data is insufficient or of poor quality for reliable model identification

• When the primary goal is to identify a nonlinear model (e.g., neural state-space, NLARX, Hammerstein-Wiener)

• When estimating the parameters of a Simulink model using experimental data; use Simulink Design Optimization Toolbox instead

• When designing a controller; use Control System Toolbox skills after identifying the plant

• Signal processing (filtering, spectral analysis without model fitting) — use signal processing skills

Execution Strategy

Write a single end-to-end MATLAB script and run it. Do NOT step through phases one tool call at a time. The script should:

1. Create/load data + split into estimation/validation

2. Estimate delay, select structure, estimate model(s)

3. Validate on held-out data by simulation

4. Print results

Only break into multiple steps if the first script fails or produces poor results (fit < 70%).

Critical rules for every script (MANDATORY — violating any of these is a bug):

• InteractiveOrderSelection=false when using order vectors — this is a HIDDEN property (not visible in disp() or tab-complete) on BOTH ssestOptions AND n4sidOptions. It MUST be set explicitly or a GUI popup HALTS execution

• EstimateCovariance=false during ANY search loop (order scan, delay scan) — covariance for discarded models wastes time

• Focus='simulation' for simulation/control use on ssestOptions, n4sidOptions, procestOptions (NOT available on tfestOptions — tfest has no Focus)

• Multi-model compare returns CELL: [~, fits] = compare(zv, m1, m2); fits{1}, fits{2} — ALWAYS pass 2+ models to ONE compare() call, NEVER call compare() separately per model

• data.InterSample = 'foh' BEFORE CT estimation if input is smooth analog

• For multi-input InterSample: use column cell {'zoh'; 'foh'} (NOT row cell)

• Delay-first: ALWAYS call delayest or inspect impulse response BEFORE any model estimation (even for MIMO, even when delay seems small)

• Hedge delays: NEVER trust a single delay estimate — always try nk AND nk±1, compare fits, pick best

• Order range: When selecting order, use a RANGE (vector) not a single integer — ssest(ze, 2:8, opt) not ssest(ze, 4, opt)

Arguments

The user provides: $ARGUMENTS

Parse:

• project_name (optional): name of a project under projects/ that has a SPEC.md

• data_source (optional): path to a .mat file, variable name in workspace, or inline description of the data

If neither is provided, ask the user to specify a data source or describe the identification problem.

Design Principles

1. Start simple, add complexity only when data justifies it. Try order 2-4 before 10-15.

2. Delay first. A wrong delay cannot be fixed by higher order — it's catastrophic.

3. Set Focus correctly. The #1 missed option. Default 'prediction' is sometimes wrong for simulation use.

4. Always hold out validation data. Never report training fit as performance.

5. Regularization > high order. A regularized ARX(30) often outperforms unregularized ARX(5). Use arxRegul, or ssregest.

6. State-space is the default. When unsure, ssest handles MIMO, CT/DT, needs only order n.

7. Compare 2-3 structures. The first model is rarely the best.

8. Check residuals. A high fit with correlated residuals means the model is missing dynamics.

9. Know when to stop. >90% fit with white residuals on validation data is success.

Fast Path — Use When Problem Is Clear

If the problem maps directly to one of these patterns, write a single script immediately:

Step/impulse response → process model (do NOT split single-transient data):

% 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 time-domain → transfer function:

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 → state-space:

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 / large periodic data → frequency-domain path:

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

If the fast path gives fit > 85%, you're done. Report results and move on.

Deep Path — For Ambiguous or Failed First Attempts

Use this structured investigation when the fast path fails (fit < 70%), the problem is ambiguous, or the user asks for deeper analysis.

Problem Characterization

Determine: SISO/MIMO, time/frequency domain, intended use (simulation/prediction/control), known constraints. See references/data-preparation.md for preprocessing details.

Nonlinearity Check (STOP/GO Gate)

Before committing to linear identification, verify that a linear model is appropriate.

Decision

• If nonlinearity is mild (gain varies <20%), proceed with linear ID but note limitations

• If strong nonlinearity detected, recommend: Hammerstein-Wiener, NLARX, or neural state-space

• If non-stationary (time-varying), consider segmented estimation or recursive methods

Data Preparation

See references/data-preparation.md for the full preprocessing workflow including:

• Loading and inspection (advice, plot)

• Preprocessing checklist (offsets, missing data, outliers, non-uniform sampling)

• Prefiltering (band-pass, frequency weighting)

• InterSample behavior for CT models

• Train/validation split

• Frequency-domain conversion

• Data quality red flags

• Probability of output feedback in the data (checkFeedback)

• Persistence of excitation check (pexcit)

Model Type Selection

Apply this decision tree. The FIRST matching branch is the recommendation:

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

See references/model-structures.md for detailed guidance on process models, polynomial models, and when to use each.

Order Determination

See references/order-determination.md for methods:

• Delay estimation (delayest, impulse response)

• ARX structure search (arxstruc, selstruc)

• Subspace order selection (n4sid with order range)

• Iterative complexity (transfer function ladder)

• Process model ladder

• Frequency-domain path (AAA initialization)

Rules of thumb:

• Max useful order: n_max ~ min(N/20, 30)

• MIMO state-space: start with n = max(ny, nu) * 2 up to 5

• If ARX(10) and ssest(4) give similar fits, prefer ssest(4)

• Stop increasing order when improvement < 2% per additional parameter

Estimation

See references/estimation.md for the full estimation workflow including:

• Critical options (Focus, InitialState, Stability, Regularization, OutputWeight, WeightingFilter, ErrorThreshold)

• Initialization methods for ssest and tfest

• Why ssest outperforms n4sid

• Frequency-domain estimation path (AAA, lsrf)

• Estimation commands for all model types

• Regularized ARX to state-space (ssregest)

• Structured estimation with parameter bounds

• Multi-model strategy

Key reminders:

% 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'));

Validation (MANDATORY)

Compare on Validation Data (MANDATORY)

CRITICAL: Always validate by simulation (infinite prediction horizon). A model estimated with Focus='prediction' can show excellent 1-step-ahead fits even when the dynamics are wrong.

% 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

Residual Analysis

% 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));

Residual interpretation guide:

Fit Interpretation

When model is inadequate — re-estimation recipe:

% 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{:});

See references/validation-and-diagnostics.md for uncertainty analysis, stability assessment, diagnostic checklist, and initial conditions guidance.

Results and Artifacts

Generate Outputs

1. Model summary — type, order, fit percentage, key parameters

2. Comparison plot — measured vs. simulated on validation data

3. Bode plot — frequency response with confidence bounds

4. Residual plot — autocorrelation and cross-correlation

% 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);

Results Summary

================================================================
  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>
================================================================

Save Artifacts

1. MATLAB script — identify_model.m containing the full reproducible workflow

2. Model MAT file — identified_model.mat with the final model object

3. Figures — PNG files for comparison, Bode, residuals

4. Results summary — printed to console

Special Topics

For advanced scenarios, see:

• references/special-topics.md — closed-loop ID, control-oriented ID, MIMO, continuous-time, physics-informed constraints, frequency-domain limitations, recursive estimation, tfest vs. tfestimate

Common Gotchas — Quick Reference

Copyright 2026 The MathWorks, Inc.