Chapter 02

Subspace Identification — the 4SID Algorithm

Given nothing but a measured impulse response, recover a state-space model. The Hankel factorisation makes this constructive rather than an optimisation: grow Hankel matrices until the rank stops climbing, split the result into observability and reachability factors, and read the matrices off. Then do it again with noise, where the rank test collapses and the SVD takes over.

Reading: ~35 min Interactive: 1 widgets Source: Polimi MIDA2 2025/26 — Lecture 4 (Ch. 2.1 — the 4SID algorithm) · Polimi MIDA2 2025/26 — Lecture 5, slides 1–8 (optimality of the SVD truncation; worked example) · Polimi MIDA2 2025/26 — Exercise session 1 (SS identification, part 1) · Polimi MIDA2 2025/26 — Exercise session 2 (SS identification, part 2) · Polimi MIDA2 2025/26 — Exercise session 3 (SS identification, part 3 — MATLAB)

01 · Motivation

From data back to a model

Chapter 01 ended on an asymmetry. Going from a state-space model to an impulse response is a formula. Going the other way — the direction an engineer with a data logger actually needs — was marked with a question mark on the transformation triangle. This chapter removes it.

The setup is deliberately spare. You may excite the system and record its response; you may not open it. What you obtain is a truncated impulse response

The dataset

$\{\omega(0),\,\omega(1),\,\omega(2),\dots,\omega(N)\}$

and what you want is a quadruple $\{\hat F,\hat G,\hat H,\hat D\}$ reproducing it. The method is called 4SID — subspace-based state-space system identification — and its engine is the factorisation from the end of Chapter 01:

$\mathcal H_n=\mathcal O\cdot\mathcal R.$

Measured numbers on the left; the unknown model on the right. Every step below is a consequence.

why

Constructive, not optimised — and why that is unusual

MIDA1 identified models by optimisation: pick a parametric family $M(\vartheta)$ , define a performance index

$J(\vartheta)=\frac1N\sum_{t=1}^{N}\big(y(t)-W(z;\vartheta)u(t)\big)^{2},$

and minimise it. There is a model class chosen up front, a cost function, and a search.

4SID has none of these. It is non-parametric and constructive: no family is postulated, no cost is minimised, no iteration runs. The order is discovered from a rank, and the matrices fall out of a matrix factorisation. This is the conceptual hinge between the two modules, and it is worth being able to state in one sentence — the theory questions ask for exactly this contrast.

tip

Why the method sat unused for thirty years

4SID for noise-free data was known in the 1960s and dismissed as useless, because it is not robust to noise: the rank test that carries step 1 disintegrates the moment measurements are imperfect (§05). It was rediscovered in the mid-1990s, when numerical linear algebra had matured enough to deliver a robust singular value decomposition — the tool that replaces the rank test. The method did not change; the numerics underneath it did.

02 · The method

4SID in three steps

The single hottest topic in the course

4SID carries 10 questions worth 56 marks across 9 of the 12 sittings in the archive — 14% of every mark available, the largest single block in the bank. It is almost always question 1, it is almost always numerical, and it is almost always the same three steps applied to a list of impulse coefficients. Marks here are the most predictable in the course, and the procedure below is worth drilling until it is mechanical.

Assume for now that the impulse response is measured noise-free. The whole algorithm is three steps.

Memorise this The 4SID procedure

Step 1 — find the order from the Hankel rank

Build Hankel matrices of increasing size from $\omega(1)$ onwards and watch the rank:

$\mathcal H_1=[\omega(1)],\quad \mathcal H_2=\begin{bmatrix}\omega(1)&\omega(2)\\ \omega(2)&\omega(3)\end{bmatrix},\quad \mathcal H_3=\begin{bmatrix}\omega(1)&\omega(2)&\omega(3)\\ \omega(2)&\omega(3)&\omega(4)\\ \omega(3)&\omega(4)&\omega(5)\end{bmatrix},\ \dots$

The ranks climb $1,2,3,\dots$ and then stop. The value at which the rank stops climbing is the order $n$ . Section §03 is entirely about the ways this step can mislead you.

Step 2 — factorise $\mathcal H_{n+1}$

Take the $(n+1)\times(n+1)$ Hankel matrix — one size past the order — which has rank $n$ , and split it into a tall factor and a wide one:

$\mathcal H_{n+1}=\mathcal O_{n+1}\cdot\mathcal R_{n+1},\qquad \mathcal O_{n+1}:(n+1)\times n,\quad \mathcal R_{n+1}:n\times(n+1)$

These are the extended observability and extended reachability matrices:

$\mathcal O_{n+1}=\begin{bmatrix}H\\HF\\ \vdots\\HF^{\,n}\end{bmatrix},\qquad \mathcal R_{n+1}=\begin{bmatrix}G&FG&\cdots&F^{\,n}G\end{bmatrix}$

Step 3 — read off $\hat H$, $\hat G$, $\hat D$, and shift for $\hat F$

Three of the four matrices are simply extracted:

$\hat H=\mathcal O_{n+1}(1,:)\ \ \text{(first row)},\qquad \hat G=\mathcal R_{n+1}(:,1)\ \ \text{(first column)},\qquad \hat D=\omega(0)$

$\hat F$ needs one idea. Split $\mathcal O_{n+1}$ into its top $n$ rows and its bottom $n$ rows:

$\mathcal O_1=\mathcal O_{n+1}(1{:}n,:),\qquad \mathcal O_2=\mathcal O_{n+1}(2{:}n{+}1,:)$

Because each row of $\mathcal O$ is the previous one times $F$ , these two square matrices satisfy the shift invariance property $\mathcal O_2=\mathcal O_1F$ , so

$\hat F=\mathcal O_1^{-1}\,\mathcal O_2.$

Choosing the factorisation

Step 2 says “split it” without saying how — and a rank- $n$ matrix factorises in infinitely many ways. Any of them works. The standard recipe makes the arithmetic almost disappear:

✓

The general factorisation rule

Take $\mathcal R_{n+1}=\mathcal H_{n+1}(1{:}n,:)$ — the first $n$ rows of the Hankel matrix. They are guaranteed independent, because $\mathcal H_n$ was full rank.
Then $\mathcal O_{n+1}$ must reproduce those rows unchanged and express the last one as a combination of them, so it has the shape

$\mathcal O_{n+1}=\begin{bmatrix}I_n\\ \square\end{bmatrix}$

and only the last row has to be computed — solve for the coefficients writing row $n+1$ of $\mathcal H_{n+1}$ in terms of the first $n$ .

The payoff: with $\mathcal O_1=I_n$ , step 3’s inverse vanishes and

$\hat F=\mathcal O_1^{-1}\mathcal O_2=\mathcal O_{n+1}(2{:}n{+}1,:)$

— $\hat F$ is the bottom $n$ rows of $\mathcal O_{n+1}$ , read straight off. No matrix inversion anywhere in the noise-free algorithm.

Do not expect the true matrices back

Chapter 01 §08 established that identification recovers a realization, never the realization. 4SID returns whichever minimal realization its factorisation happened to pick, related to the true one by some unknown change of basis. Comparing $\hat F$ entry-by-entry against a “true” $F$ is a category error; compare the invariants — eigenvalues, transfer function, impulse response.

The same warning has a sharper form: 4SID recovers only the reachable-and-observable subsystem. Feed it data from a second-order plant with a cancellation and it returns a first-order model — and that is the correct answer, not a failure. The data never contained the hidden mode.

tip

An unstable plant in, a stable model out — correctly

Take $W(z)=\dfrac{z+2}{\left(z+\tfrac12\right)(z+2)}$ . Realised in controllable canonical form it is second order with

$F=\begin{bmatrix}0&1\\ -1&-\tfrac52\end{bmatrix},\quad G=\begin{bmatrix}0\\1\end{bmatrix},\quad H=\begin{bmatrix}2&1\end{bmatrix},$

which is fully reachable ( $\det\mathcal R_2=-1$ ) but not observable ( $\det\mathcal O_2=0$ ) — the cancellation at $z=-2$ showing up exactly as a rank loss, per Chapter 01 §07. Simulating the impulse response by hand makes the consequence vivid:

$t$	$x_1(t)$	$x_2(t)$	$y(t)$
0	$0$	$0$	$0$
1	$0$	$1$	$1$
2	$1$	$-\tfrac52$	$-\tfrac12$
3	$-\tfrac52$	$\tfrac{21}{4}$	$\tfrac14$
4	$\tfrac{21}4$	$-\tfrac{85}{8}$	$-\tfrac18$

The output halves every step while the state roughly doubles. Now run 4SID on that output. The ranks give $\operatorname{rank}\mathcal H_1=1$ and $\operatorname{rank}\mathcal H_2=1$ , so $n=1$ , and the algorithm returns $\hat F=-\tfrac12$ , $\hat G=1$ , $\hat H=1$ , $\hat D=0$ — a first-order, asymptotically stable model of a second-order unstable plant.

That is not an error. 4SID works from input/output data, and the diverging mode leaves no trace in the data. The model is a perfect description of everything that was measurable.

2023-07-q12023Q014SIDmedium8 pts

Given the following impulse-response coefficients: $$\omega(0)=1,\quad \omega(1)=0,\quad \omega(2)=-\tfrac12,\quad \omega(3)=0,\quad \omega(4)=\tfrac14,\quad \omega(5)=0,\quad \omega(6)=-\tfrac18,\quad \omega(7)=0$$ a. Identify the system matrices using the 4SID method. What is the system order? b. Compute the transfer function of the identified system. c. Identify the system matrices from the transfer function. Did you retrieve the same matrices as in (a)? d. Provide the analytical expression of the impulse response $\omega(t)$ for the identified system.

03 · The trap

Reading the order off the rank — and the plateau that lies

Where the marks are lost

Order determination is examined in its own right — 2 questions, 12 marks, 2 sittings — and it is step 1 of all ten 4SID questions, so an error here voids the entire answer. Both standalone questions were built around exactly the failure mode below, and both are graded hard. This is the highest mark-per-line section of the chapter.

Step 1 was stated as “the rank stops climbing”. Read carelessly, that becomes “the first time the rank fails to increase” — and those are not the same sentence.

The rank can pause and then resume

Take $W(z)=z^{-4}+z^{-1}$ , whose impulse response is $\omega(1)=1$ , $\omega(4)=1$ , and zero everywhere else. The Hankel ranks run

$\operatorname{rank}\mathcal H_1=1,\quad \operatorname{rank}\mathcal H_2=1,\quad \operatorname{rank}\mathcal H_3=3,\quad \operatorname{rank}\mathcal H_4=4,\quad \operatorname{rank}\mathcal H_5=4.$

A student who stops at the first repeat reports $n=1$ and identifies $\tilde W(z)=z^{-1}$ — a completely different system. The true order is $4$ .

The rule is “stops increasing definitively”. Always extend at least one matrix past the first repeat before committing, and prefer two when the response contains zeros.

Run that example to the end and the identified $\hat F$ comes out as a pure shift — ones on the superdiagonal, zeros elsewhere — which is nilpotent: $\hat F^{4}=0$ . That is exactly right for a finite impulse response, since a state matrix that annihilates itself after $n$ steps produces a response that terminates. A nilpotent $\hat F$ is a useful sanity signal: it says the data came from an FIR system, and it is impossible for any plant whose response decays geometrically forever.

Two structures produce these false plateaus, and both are common in exam data:

Pure delay

A system that does nothing for $k$ steps has $\omega(1)=\dots=\omega(k-1)=0$ , so the early Hankel matrices are rank-deficient for a reason that has nothing to do with the order. In a worked lecture example $\mathcal H_1=[0]$ outright — rank $0$ — while the true order is $2$ .

Interleaved zeros

A response alternating between zero and non-zero entries (typical of systems whose poles sit on the imaginary axis) makes ranks climb in jumps of two, so a single-step check lands inside a gap.

✓

A reliable procedure

Build one more Hankel matrix than you think you need, every time. Concretely: find the first $i$ with $\operatorname{rank}\mathcal H_{i+1}=\operatorname{rank}\mathcal H_i$ , then check $\mathcal H_{i+2}$ as well. If the rank is still unchanged, the order is that common rank. If it has climbed, keep going. The cost is one extra determinant; the benefit is not answering the wrong question.

Note that the order is the rank at the plateau, not the index of the matrix where you stopped — these coincide in the clean case and diverge in exactly the cases above.

2025-09-q1b2025Q01Hankel rank and system orderhard6 pts

Given the following impulse-response coefficients: $$\omega(0)=1;\;\omega(1)=0;\;\omega(2)=0;\;\omega(3)=1;\;\omega(4)=\tfrac12;\;\omega(5)=\tfrac14;\;\omega(6)=\tfrac18;\;\omega(7)=\tfrac1{16}$$ Identify the system matrices using the 4SID method. What is the system order?

Building Hankel matrices from an impulse response you obtain ranks 1, 2, 2, 3, 3, 3. What is the system order?

04 · Before you start

Turning a real test into an impulse response

The algorithm wants $\{\omega(t)\}$ : the response to a unit impulse applied at $t=0$ . Real experiments rarely oblige, and exam questions almost never do. Two corrections cover nearly every case.

Amplitude. If the input is $u(0)=A$ rather than $1$ , linearity means the recorded output is $A\,\omega(t)$ . You can either divide the whole response by $A$ before starting, or — equivalently, and the way MATLAB code in the exercise sessions does it — run the algorithm on the raw data and divide the result:

Correcting for amplitude

$\hat G=\frac{\mathcal R_{n+1}(:,1)}{A}.$

Only $\hat G$ needs it. The rank, and therefore the order, is unaffected — scaling every entry of a matrix cannot change its rank — and $\hat F$ and $\hat H$ come out of ratios that cancel the factor.

Delay. If the impulse arrives at $t=t_0$ rather than $t=0$ , shift the time axis: the sample recorded at $t_0+k$ is $\omega(k)$ . Re-index first, then build the Hankel matrix.

Two index errors that cost whole questions

Putting $\omega(0)$ into the Hankel matrix. It starts at $\omega(1)$ . Including the feedthrough changes every rank and destroys the factorisation. $\omega(0)$ has exactly one job: it is $\hat D$ .
Forgetting the amplitude divisor. The identified $\hat G$ comes out $A$ times too large, so the model’s gain is wrong by $A$ while its poles are perfectly right — a failure mode that survives every sanity check that only looks at eigenvalues.

✓

A non-zero $\omega(0)$ is harmless

A system that is not strictly proper has $\omega(0)\ne 0$ , and students often assume 4SID needs adjusting. It does not. $\omega(0)$ enters the algorithm in exactly one place — the assignment $\hat D=\omega(0)$ — and never touches the Hankel matrix at all. Contrast this with the controllable canonical form of Chapter 01 §09, where a non-zero feedthrough must be split off before you can do anything: same plant, opposite handling.

2024-07-q1a2024Q014SIDmedium5 pts

A system is subject to the following input: $$u(t)=\begin{cases}3 & \text{if } t=2\\ 0 & \text{else}\end{cases}$$ The response of the system is $$y(0)=0,\; y(1)=0,\; y(2)=6,\; y(3)=0,\; y(4)=\tfrac38,\; y(5)=0,\; y(6)=\tfrac3{32},\; y(7)=0$$ a. Identify the system order. b. Identify the system matrices using the 4SID method. c. Provide the MATLAB code to compute the eigenvalues of the matrix $F$ found in question (b), then write the system of (b) as a discrete-time state-space MATLAB object.

F = [0 1; 1/4 0];
G = [0; 1/8];
H = [1 0];
D = 2;
eig(F)
sys = ss(F, G, H, D, -1);

05 · Real data

4SID with noise — one big Hankel and an SVD

Everything above assumed exact measurements. Add noise and step 1 fails — not gracefully, but completely.

The measured response is now

Noisy measurement

$\tilde\omega(t)=\underbrace{\omega(t)}_{\text{true}}+\underbrace{\eta(t)}_{\text{measurement noise}},\qquad t=1,\dots,N$

with $N$ typically 500–1000 samples. And here is the problem: a matrix of noisy numbers is essentially always full rank. In a lecture case study on a second-order pendulum, the noiseless ranks ran $1,2,2$ — order 2, correct. With noise on the same system they ran $1,2,3,4$ : climbing forever, never plateauing. The rank test does not degrade under noise; it stops working.

why

Rank is a yes/no question about an approximate object

Rank asks whether a row is exactly a combination of others. Noise makes the answer “no” every time, however nearly true it is. What is needed is a graded version — a measure of how close the matrix is to having low rank — and that is exactly what singular values provide.

Step 1 (noisy) — one big rectangular Hankel matrix

Instead of a sequence of square matrices, build a single rectangular one from all the data:

Rectangular Hankel

$\tilde{\mathcal H}_{q,d}=\begin{bmatrix}\tilde\omega(1)&\tilde\omega(2)&\cdots&\tilde\omega(d)\\ \tilde\omega(2)&\tilde\omega(3)&\cdots&\tilde\omega(d+1)\\ \vdots& & &\vdots\\ \tilde\omega(q)&\tilde\omega(q+1)&\cdots&\tilde\omega(q+d-1)\end{bmatrix},\qquad q+d-1=N.$

The constraint $q+d-1=N$ is arithmetic, not advice: the bottom-right entry must be the last sample you actually recorded. Within it there is a genuine trade-off, taking $d>q$ :

choice	model quality	computational cost
$q\approx d$	best	worst
$q \ll d$	worst	best

The rule of thumb is $q>\tfrac12 d$ . For $N=1000$ , the choices $q=350$ , $400$ and $450$ all give very similar results.

Step 2 — the singular value decomposition

SVD

$\tilde{\mathcal H}_{qd}=\tilde U\,\tilde S\,\tilde V^{\top},\qquad \tilde U:q\times q,\quad \tilde S:q\times d,\quad \tilde V^{\top}:d\times d$

with $\tilde U$ and $\tilde V$ unitary ( $M^{-1}=M^{\top}$ ), and $\tilde S$ rectangular-diagonal carrying the singular values, which are real, positive, and sorted:

$\sigma_1\ge\sigma_2\ge\cdots\ge\sigma_q.$

The SVD is a diagonalisation for rectangular matrices, and singular values are the corresponding stand-in for eigenvalues — formally $\mathrm{SV}(M)=\sqrt{\mathrm{EIG}(MM^{\top})}$ .

Step 3 — cut the spectrum into “system” and “noise”

Plot the singular values in order. In the ideal case there is a clear jump after a few of them: the large ones are system, the flat tail is noise, and the index at the jump is the order $n$ . In real data there is usually no jump but a knee, and you choose the last singular value you judge to belong to the system.

Truncating to rank $n$ is the denoising step: everything below the cut is discarded as the residual.

This step is a judgement call, and the exam knows it

There is no formula for where to cut. Typical lecture figures use $N=1000$ , $q=400$ , and a chosen $n=8$ out of 400 singular values. Choosing too low gives a model that cannot represent the plant; too high fits noise as though it were dynamics. When a question asks you to “determine the order from the singular values”, the expected answer names the jump and acknowledges that the knee is read, not computed.

Having chosen $n$ , truncate: keep the first $n$ columns of $\tilde U$ , the leading $n\times n$ block of $\tilde S$ , and the first $n$ rows of $\tilde V^{\top}$ , giving

Rank-reduced reconstruction

$\underbrace{\tilde{\mathcal H}_{qd}}_{\operatorname{rank}\,q}=\underbrace{\hat{\mathcal H}_{qd}}_{\operatorname{rank}\,n}+\underbrace{\mathcal H_{\text{RES}}}_{\text{discarded}},\qquad \hat{\mathcal H}_{qd}=\hat U\hat S\hat V^{\top}.$

This is where the denoising actually happens — a rank reduction from $400$ to $8$ in the figures above, with the residual thrown away.

Step 4 — extract the matrices, with one change

Split the square root of $\hat S$ symmetrically between the two factors:

Balanced factorisation

$\hat{\mathcal H}_{qd}=\big(\hat U\hat S^{1/2}\big)\big(\hat S^{1/2}\hat V^{\top}\big)=\hat{\mathcal O}\cdot\hat{\mathcal R},\qquad \hat S^{1/2}=\operatorname{diag}\big(\sqrt{\sigma_1},\dots,\sqrt{\sigma_n}\big)$

Then $\hat H=\hat{\mathcal O}(1,:)$ and $\hat G=\hat{\mathcal R}(:,1)$ exactly as before. The change is in $\hat F$ :

$\hat{\mathcal O}_1$ is no longer square — `inv` becomes `pinv`

In the noise-free route $\mathcal O_1$ was $n\times n$ and invertible. Here $\hat{\mathcal O}_1=\hat{\mathcal O}(1{:}q-1,:)$ is $(q-1)\times n$ with $q-1$ enormously larger than $n$ . Shift invariance still holds — $\hat{\mathcal O}_1F=\hat{\mathcal O}_2$ — but this is now an over-determined system with no exact solution, so it is solved in the least-squares sense:

$\hat F=\big(\hat{\mathcal O}_1^{\top}\hat{\mathcal O}_1\big)^{-1}\hat{\mathcal O}_1^{\top}\,\hat{\mathcal O}_2$

The matrix $A^{+}=(A^{\top}A)^{-1}A^{\top}$ is the pseudo-inverse — MATLAB’s pinv. Using inv or \ here is a shape error, not a stylistic one.

Deep dive Why over-determined systems need least squares at all

Solving $Ax=B$ with $A$ of size $h\times n$ has three regimes:

case	shape	verdict
more unknowns	$h \lt n$	under-determined — infinitely many solutions
square	$h=n$	one solution, $x=A^{-1}B$ (if $A$ is regular)
more equations	$h \gt n$	over-determined — generally no solution

The third case is ours, with $h=q-1$ in the hundreds and $n$ in the single digits. There is no $x$ satisfying every equation, so one asks instead for the $x$ minimising the total squared residual. Differentiating and setting to zero gives the normal equations

$A^{\top}Ax=A^{\top}B\ \Longrightarrow\ \hat x=\big(A^{\top}A\big)^{-1}A^{\top}B,$

with $A^{\top}A$ square and invertible whenever $A$ has full column rank — which $\hat{\mathcal O}_1$ does, by construction, since we truncated at rank $n$ .

Worth noticing: having far more equations than unknowns is not a nuisance here, it is the whole point. Each surplus row is another noisy measurement of the same shift relation, and least squares averages them. The over-determination is where the noise rejection comes from.

06 · Why it works

Why the truncated SVD is the right answer

Step 3 discarded a residual matrix. That should prompt a question: of the infinitely many ways to write a rank- $q$ matrix as “rank $n$ plus a remainder”, why is this one privileged?

A rank reduction can be done arbitrarily — paste in an identity block and you have split a rank-3 matrix into a rank-2 piece and a remainder, having learned nothing. What makes a reduction good is that the discarded part is as small as possible. Measuring size entrywise with the Frobenius norm

Frobenius norm

$\big\|\mathcal H_{\text{RES}}\big\|_F=\sqrt{\sum_{i,j}\big(\mathcal H^{\,ij}_{\text{RES}}\big)^{2}},$

the truncated SVD is provably the minimiser — no rank- $n$ matrix approximates $\tilde{\mathcal H}_{qd}$ better. The optimality comes precisely from the singular values being sorted: cutting after $\sigma_n$ discards the smallest possible amount of the matrix.

This is why the method waited for the SVD rather than for a better idea. The same result explains why the SVD turned up everywhere in machine learning at the same time — order reduction, separating signal from noise, and compression are one theorem wearing three hats.

07 · In practice

4SID in MATLAB

MATLAB is now on nearly every paper

The course states plainly that the exam contains simple MATLAB questions, and the archive bears it out: MATLAB is absent from all three sittings before June 2024 and present in 8 of the 9 papers since — every one except February 2026. Three of those questions are graded under identification specifically (7 marks, 3 sittings). The examinable function list is closed and short; the questions ask for three to five lines. These are the cheapest marks on the paper.

The examinable function set, as stated in the exercise sessions:

sys = ss(F, G, H, D, Ts);        % or sys = ss(W);  -- omit Ts and you get a CONTINUOUS-time system
eigenvalues = eig(sys);
O = obsv(sys);                   % or obsv(F, H)
R = ctrb(sys);                   % or ctrb(F, G)
obs_rank = rank(O);
W = tf(sys);                     % or tf(num, den, Ts)
W = zpk(W);                      % zero-pole-gain form
zeros_vec = zero(W);   pole_vec = pole(W);   sys_gain = dcgain(W);
bode(sys)   impulse(sys)   step(sys)   lsim(sys, u, t)
sys = n4sid(data);               % 4SID
[U, S, V] = svd(Hqd);            % 4SID

The two tests you are most often asked to write out, in idiom form:

nx = size(sys.A, 1);
all(abs(eig(sys)) < 1)                  % asymptotically stable  (<= 1 for simply stable)
rank(obsv(sys)) == nx                   % fully observable
rank(ctrb(sys)) == nx                   % fully reachable

Identifying a model from a general input/output experiment is three lines:

data = iddata(y_vec, u_vec, Ts);   % iddata(OUTPUT, INPUT, sampling TIME)
sys  = n4sid(data, n);             % order n; omit n and n4sid picks it from the singular values
y_sim = lsim(sys, u_vec, t_vec);

Four errors that appear in real answer keys

Argument order in iddata. The signature is iddata(y, u, Ts) — output first. Reversing it identifies the inverse system and everything downstream is wrong. One official solution in the archive gets this backwards; a later paper’s key has it right.
Sampling frequency instead of sampling time. At 1000 Hz, Ts is 1/1000, not 1000.
impulse or step when the question supplies an input. Those generate their own inputs; to drive a model with given data you need lsim.
Omitting Ts in ss(F,G,H,D). Without it MATLAB builds a continuous-time system and every stability conclusion inverts. Use Ts = -1 if the period is genuinely unspecified but the system is discrete.

MATLAB indexes from 1; the impulse response indexes from 0

Stated more times than anything else in the exercise sessions. The mapping is

$\omega(0)=\texttt{y(1)},\quad \omega(1)=\texttt{y(2)},\quad \omega(2)=\texttt{y(3)},\ \dots$

so ir_data(1) is the feedthrough $\hat D$ , and the Hankel matrix starts at ir_data(2). This off-by-one compounds with the $\omega(0)$ -versus- $\omega(1)$ rule of §04, and the two together are the most reliable source of wrong numerical answers in the course.

✓

The experiment must excite the system

When identifying from a general input rather than an impulse, the input has to be rich enough to reveal the dynamics — persistent excitation. The lecture case study uses filtered white Gaussian noise for exactly this reason. A single sinusoid excites one frequency and identifies almost nothing; a step is barely better. If a question asks why a particular input was chosen, this is the answer.

Two further practical notes: with n4sid the order is usually chosen a priori (too low gives an incorrect model, too high an over-parametrised one that fits noise), and n4sidOptions('EnforceStability', true) constrains the result to a stable model when you know the plant is stable.

2025-09-q1c2025Q01MATLAB for identificationeasy3 pts

You are given a generic dataset containing the input and output timeseries signals of a system: `u_vec` and `y_vec` with time vector `t_vec`. The sampling frequency of the vectors is 1000 Hz. Write the MATLAB code to identify the system `sys` using the 4SID method, assuming that the system has order 5. Write the code to compute the response of the identified system `sys` to the input vector `u_vec`.

data  = iddata(y_vec, u_vec, 1/1000);   % iddata(OUTPUT, INPUT, Ts)
sys   = n4sid(data, 5);                 % 4SID, order 5
y_sim = lsim(sys, u_vec, t_vec);

08 · End to end

A complete worked identification

One pass through everything: start from known matrices, generate the impulse response they produce, throw the matrices away, and recover a model from the numbers alone.

Full circuit From matrices, to data, and back

The data-generating system

$F=\begin{bmatrix}\tfrac12&0\\1&\tfrac14\end{bmatrix},\quad G=\begin{bmatrix}1\\0\end{bmatrix},\quad H=\begin{bmatrix}0&1\end{bmatrix},\quad D=0$

Order $n=2$ , SISO, eigenvalues $\{\tfrac12,\tfrac14\}$ — both inside the unit circle, so asymptotically stable.

Check it is fully observable and reachable

$\mathcal O=\begin{bmatrix}H\\HF\end{bmatrix}=\begin{bmatrix}0&1\\1&\tfrac14\end{bmatrix},\quad \operatorname{rank}=2;\qquad \mathcal R=\begin{bmatrix}G&FG\end{bmatrix}=\begin{bmatrix}1&\tfrac12\\0&1\end{bmatrix},\quad \operatorname{rank}=2$

Both full rank, so the realization is minimal and identification has a chance of recovering all of it. (Had either failed, 4SID would return a lower-order model — correctly.)

Transfer function, by direct manipulation

$(z-\tfrac12)x_1=u\Rightarrow x_1=\frac{u}{z-\tfrac12};\qquad (z-\tfrac14)x_2=x_1\Rightarrow x_2=\frac{u}{(z-\tfrac12)(z-\tfrac14)}$

and since $y=x_2$ ,

$W(z)=\frac{1}{(z-\tfrac12)(z-\tfrac14)}.$

Its poles are $\tfrac12$ and $\tfrac14$ — identical to the eigenvalues of $F$ , with no cancellation, exactly as full observability and reachability guarantee.

Impulse response by long division

In negative powers, $W(z)=\dfrac{z^{-2}}{1-\tfrac34z^{-1}+\tfrac18z^{-2}}$ , and dividing out:

$\omega(0)=0,\quad \omega(1)=0,\quad \omega(2)=1,\quad \omega(3)=\tfrac34,\quad \omega(4)=\tfrac7{16},\quad \omega(5)=\tfrac{15}{64}$

Note $\omega(0)=0$ (strictly proper) and $\omega(1)=0$ — a two-step delay. That second zero is about to matter.

Three other routes reach the same numbers, and any of them is acceptable in an exam: the formula $\omega(t)=HF^{\,t-1}G$ ; the geometric-series expansion when the denominator is simple enough; or simply rolling out the recursion — set $u(0)=1$ , $x(0)=0$ , and iterate $x(t+1)=Fx(t)+Gu(t)$ , $y(t)=Hx(t)+Du(t)$ by hand. The last is slowest per sample but needs nothing memorised, and it is the one to fall back on under pressure.

Now forget the matrices — find the order

$\mathcal H_1=[0],\qquad \mathcal H_2=\begin{bmatrix}0&1\\1&\tfrac34\end{bmatrix}\ \operatorname{rank}=2,\qquad \mathcal H_3=\begin{bmatrix}0&1&\tfrac34\\ 1&\tfrac34&\tfrac7{16}\\ \tfrac34&\tfrac7{16}&\tfrac{15}{64}\end{bmatrix}\ \operatorname{rank}=2$

$\mathcal H_1$ has rank zero — the delay again. Applying “stop at the first matrix that is not full rank” mechanically would report $n=1$ and be wrong. Reading where the rank saturates gives $n=2$ , correct.

Factorise and read off

$\mathcal H_3=\mathcal O_3\mathcal R_3$ with the extended matrices

$\mathcal O_3=\begin{bmatrix}0&1\\1&\tfrac14\\ \tfrac34&\tfrac1{16}\end{bmatrix},\qquad \mathcal R_3=\begin{bmatrix}1&\tfrac12&\tfrac14\\ 0&1&\tfrac34\end{bmatrix}$

from which $\hat H=\mathcal O_3(1,:)$ , $\hat G=\mathcal R_3(:,1)$ , $\hat D=\omega(0)=0$ , and $\hat F=\mathcal O_3(1{:}2,:)^{-1}\mathcal O_3(2{:}3,:)$ — recovering a model with the same transfer function, the same eigenvalues, and the same impulse response as the system we started from.

✓

What to check when you finish

An identification you cannot validate is a guess. Three checks, in increasing strength:

Eigenvalues of $\hat F$ should be plausible — inside the unit circle for a stable plant, and matching any decay rate visible in the data.
Reconstruct the impulse response from $\{\hat F,\hat G,\hat H,\hat D\}$ and compare with the measurements. It should agree to within the noise.
Simulate a different input — a step is standard — and compare against a held-out experiment. This is the only check that tests generalisation rather than fit, and it is what the MATLAB case studies do with lsim.

Load-bearing ideas

The three steps: grow Hankel matrices until the rank stops climbing (that rank is $n$ ); factorise $\mathcal H_{n+1}=\mathcal O_{n+1}\mathcal R_{n+1}$ ; read off $\hat H=\mathcal O_{n+1}(1,:)$ , $\hat G=\mathcal R_{n+1}(:,1)$ , $\hat D=\omega(0)$ , and $\hat F=\mathcal O_1^{-1}\mathcal O_2$ by shift invariance.
The factorisation shortcut: take $\mathcal R_{n+1}$ = first $n$ rows of $\mathcal H_{n+1}$ , so $\mathcal O_{n+1}=[I_n;\,\square]$ and $\hat F$ is just its bottom $n$ rows — no inverse needed.
“Stops increasing” means definitively. The rank can pause — a pure delay makes early Hankel matrices rank-deficient for reasons unrelated to the order. Always check one matrix past the first repeat.
Normalisation: divide $\hat G$ by the impulse amplitude; re-index if the impulse is delayed; never put $\omega(0)$ in the Hankel matrix; a non-zero $\omega(0)$ is harmless and only sets $\hat D$ .
With noise the rank test dies. Build one $q\times d$ Hankel with $q+d-1=N$ , take the SVD, cut the singular values at the jump or knee, and rebuild at rank $n$ . Truncating the sorted SVD is the best possible rank- $n$ approximation in Frobenius norm.
$\hat{\mathcal O}_1$ becomes rectangular, so the shift step is over-determined and solved by least squares — the pseudo-inverse $(A^{\top}A)^{-1}A^{\top}$ , i.e. pinv, not inv.
4SID recovers a realization, not the realization — and only the reachable-and-observable part. Compare invariants, never entries.

Exam radar. The heaviest chapter in the course: 4SID alone is 10 questions and 14% of all marks, graded in 9 of 12 sittings, and it is almost always question 1. Drill §02’s three steps until they need no thought, then drill §03 until a pausing rank sequence makes you build one more matrix automatically. Finish with §07 — three lines of iddata/n4sid/lsim have been worth marks on 8 of the last 9 papers, and getting iddata(y, u, Ts) the right way round is the whole question.