Question 3 of 6: Least-squares identification of a second-order discrete model

Question 3: Least-squares identification of a second-order discrete model (25 marks)

Question text not reproduced: the examination questions are © Engineers and Geoscientists BC. Open the official past paper (linked at the top of this page) to read the question, then follow the worked solution below.

Given. A sampled input–output record from an unknown plant and a fixed model structure with four free coefficients (five written coefficients, one of which is absorbed by normalisation).

Find. (a) the regression that turns the model into a linear least-squares problem and its closed-form solution; (b) the conditions under which that estimate converges to the true parameters; (c) the predicted steady-state output for a constant input of 2.

Approach. Convert the transfer function into a difference equation, recognise that the unknown coefficients enter it linearly, stack one equation per usable sample and solve the resulting over-determined system in the least-squares sense; then evaluate the identified model at $z=1$.

Part (a) — the method, in prose before the algebra. The structure is an ARX (auto-regressive with exogenous input) model, and the whole point of the exercise is that although $P(z)$ is a nonlinear function of the data, the difference equation it implies is linear in the unknown coefficients. That single observation is what makes an ordinary least-squares fit possible instead of an iterative nonlinear search. One normalisation is mandatory first: multiplying numerator and denominator by any constant leaves $P(z)$ unchanged, so the five written coefficients are only four independent parameters and the fit would otherwise collapse onto the trivial all-zero solution. The standard choice is $a_{2}=1$, which corresponds to writing the difference equation with a unit coefficient on the most advanced output sample.

1. Part (a) — convert to a difference equation. Normalise $a_{2}=1$ and cross-multiply $Y(z)\left(z^{2}+a_{1}z+a_{0}\right)=U(z)\left(b_{1}z+b_{0}\right)$. Reading $z^{n}$ as an $n$-step advance,$$y(k+2)+a_{1}y(k+1)+a_{0}y(k)=b_{1}u(k+1)+b_{0}u(k)$$so that$$y(k+2)=-a_{1}y(k+1)-a_{0}y(k)+b_{1}u(k+1)+b_{0}u(k)$$
2. Part (a) — put it in regression form. Define the regressor row and the parameter vector$$\varphi(k)^{\mathsf{T}}=\begin{pmatrix}-y(k+1)&-y(k)&u(k+1)&u(k)\end{pmatrix},\qquad \theta=\begin{pmatrix}a_{1}&a_{0}&b_{1}&b_{0}\end{pmatrix}^{\mathsf{T}}$$so that each sample contributes the scalar equation $y(k+2)=\varphi(k)^{\mathsf{T}}\theta$. The unknowns appear linearly, which is the structural fact the whole method rests on.
3. Part (a) — stack the equations. A record of $N$ samples supplies usable equations only for $k=0,\dots,N-3$, because each equation needs $y(k+2)$, $y(k+1)$, $y(k)$, $u(k+1)$ and $u(k)$ simultaneously. Hence$$Y=\Phi\theta,\qquad Y=\begin{pmatrix}y(2)\\ \vdots\\ y(N-1)\end{pmatrix},\qquad \Phi=\begin{pmatrix}\varphi(0)^{\mathsf{T}}\\ \vdots\\ \varphi(N-3)^{\mathsf{T}}\end{pmatrix}$$with $\Phi$ of size $(N-2)\times4$. At least six samples are therefore needed for a determined problem and more for a genuinely over-determined one — counting $N$ equations instead of $N-2$ is the single commonest slip in this question.
4. Part (a) — minimise the equation error. Real data will not satisfy $Y=\Phi\theta$ exactly, so choose $\theta$ to minimise the sum of squared residuals$$J(\theta)=\left\|Y-\Phi\theta\right\|^{2}=\left(Y-\Phi\theta\right)^{\mathsf{T}}\left(Y-\Phi\theta\right)$$Setting $\partial J/\partial\theta=0$ gives $-2\Phi^{\mathsf{T}}\left(Y-\Phi\theta\right)=0$, i.e. the normal equations, whose solution is$$\boxed{\hat{\theta}=\left(\Phi^{\mathsf{T}}\Phi\right)^{-1}\Phi^{\mathsf{T}}Y}$$$J$ is a positive-semidefinite quadratic in $\theta$, so this stationary point is the global minimum whenever $\Phi^{\mathsf{T}}\Phi$ is invertible. In practice one solves the normal equations by QR or SVD factorisation of $\Phi$ rather than forming $\Phi^{\mathsf{T}}\Phi$ explicitly, because squaring the matrix squares its condition number.
5. Part (a) — a worked eight-sample illustration. The exam supplies no data table, so take the record below (a pseudo-random binary input, which is the standard identification test signal). Every entry is quoted to four decimals, as a real instrument would report it.

6. Part (a) — solve the illustration. Six regressor rows are available ($k=0$ to $5$). Solving the normal equations gives$$\hat{\theta}=\begin{pmatrix}\hat{a}_{1}\\ \hat{a}_{0}\\ \hat{b}_{1}\\ \hat{b}_{0}\end{pmatrix}=\begin{pmatrix}-0.900\\ 0.200\\ 0.300\\ 0.200\end{pmatrix},\qquad \boxed{\hat{P}(z)=\frac{0.300z+0.200}{z^{2}-0.900z+0.200}}$$with a largest residual of $3.9\times10^{-5}$, which is the rounding of the tabulated data and nothing else. The identified denominator factors as $(z-0.5)(z-0.4)$, so the model poles are $0.5$ and $0.4$ — both inside the unit circle, which part (c) will need.

Part (b) — conditions for convergence. Four requirements must hold together, and it is worth naming them separately because the exam asks for conditions, not for a single formula.

Identifiability of the data. $\Phi^{\mathsf{T}}\Phi$ must be non-singular, and for convergence as $N\to\infty$ the stronger requirement is that $\left(1/N\right)\Phi^{\mathsf{T}}\Phi$ tend to a positive-definite limit. That is the statement that the input is persistently exciting of order $n_{a}+n_{b}=4$: its spectrum must be non-zero at at least four distinct frequencies (two sinusoids, a PRBS, or filtered white noise all qualify; a step or a constant does not, and a single sinusoid can identify at most a two-parameter model).

Noise structure. The residual must be zero-mean and uncorrelated with the regressors, $\mathbb{E}\left[\varphi(k)\,e(k)\right]=0$. For the ARX structure this holds when the disturbance enters as white equation error, i.e. as $e(k)$ added directly to the difference equation. If instead the measurement noise is coloured, or is added to the output and then fed back through the $-y$ entries of $\varphi$, the estimate is biased however long the record: the regressors then contain the noise, $\mathbb{E}\left[\varphi e\right]\neq0$, and one must move to instrumental variables, an ARMAX/output-error structure, or prediction-error minimisation.

Adequate model set. The true plant must be representable by the assumed second-order, one-zero structure. If the real system is third order, or carries an unmodelled transport delay, least squares still converges — but to the best four-parameter approximation in the equation-error sense, not to the truth.

Open-loop, bounded data. The record must be taken in open loop (closed-loop data correlates $u$ with the noise through the controller and biases the estimate) and all signals must remain bounded, which for an unstable plant means the experiment has to be run under a stabilising controller and handled with a closed-loop identification method. Under all four conditions $\hat{\theta}\to\theta$ with probability one, and the covariance of the estimate decays as $\sigma^{2}\left(\Phi^{\mathsf{T}}\Phi\right)^{-1}$, i.e. as $1/N$.

7. Part (c) — the steady-state output is the DC gain times the input. A constant input is the $z$-domain point $z=1$ (the discrete equivalent of $s=0$), so provided the identified model is stable the final-value theorem gives$$y_{ss}=u_{\infty}\,\hat{P}(1)=u_{\infty}\,\frac{b_{1}+b_{0}}{a_{2}+a_{1}+a_{0}}$$This is the general answer to part (c) and is what should be written first, because it is valid whatever the fitted numbers turn out to be.
8. Part (c) — check stability, then substitute. The identified poles are $0.5$ and $0.4$, both of modulus less than one, so the constant-input response does converge and the final-value theorem is legitimate. With $u(k)=2$,$$\hat{P}(1)=\frac{0.300+0.200}{1-0.900+0.200}=\frac{0.500}{0.300}=1.667$$$$\boxed{y_{ss}=2\times1.667=3.333}$$Iterating the difference equation forward from rest with $u\equiv2$ reaches $3.3334$ within about twenty samples, which confirms the algebra and, more usefully, confirms the sign convention on $a_{1}$ — a sign slip there would give $\hat{P}(1)=0.5/2.1=0.238$ and an obviously wrong answer.