---

name: matlab-identify-linear-system description: > Identify a linear dynamic model from input-output or time-series data using MATLAB System Identification Toolbox. Use when estimating transfer function, state-space, ARX, ARMAX, BJ, OE polynomial or process models from measurement data. license: https://www.mathworks.com/content/dam/mathworks/license/pmrl/license.md metadata: author: MathWorks version: "1.0"

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.

| Method | How | Interpretation | |--------|-----|----------------| | Amplitude dependence | Estimate models from datasets at different input amplitudes | If gain/dynamics change with amplitude -> nonlinear | | Harmonic analysis | Apply periodic input, check for even harmonics in output spectrum | Even harmonics indicate nonlinearity | | Model order escalation | Fit orders 2, 4, 8, 12, 16 — plot fit vs. order | Plateauing at low fit despite high order -> nonlinearity | | Residual structure | Inspect residuals vs. amplitude of u or y | Systematic patterns -> nonlinear | | Split-data test | Estimate on first half, validate on second half AND vice versa (use low model order, e.g. 2-4, to avoid false positives from estimation variance) | Asymmetric fits -> non-stationary or nonlinear | | ISNLARX | Use the isnlarx method on iddata to assess severity of nonlinearity |

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