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

Question 3: Least-squares identification of a first-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 discrete model $Y(z) = P(z)U(z)$ with $P(z) = \dfrac{b}{z-a}$, and four sampled input–output pairs:

Find. The least-squares estimates of $a$ and $b$, and the steady-state output the identified model predicts for a constant input $u(k) = 2$.

Approach. Turn the transfer function into the difference equation it represents, write one equation per usable data triple, recognise the result as an over-determined linear system that is linear in the unknowns, and solve the normal equations.

1. Part (a) — convert the model to a difference equation. Cross multiplying $Y(z)(z-a) = bU(z)$ and using the forward-shift property $zY(z) \leftrightarrow y(k+1)$ for zero initial conditions, $$y(k+1) = a\,y(k) + b\,u(k).$$ This is the key structural observation of the whole question: although $P(z)$ is a ratio, the equation it induces is linear in the unknown parameters $a$ and $b$, which is exactly what makes ordinary least squares applicable without any iteration.
2. Write one equation per usable data point. The table supplies four samples, and each equation needs $y(k)$, $u(k)$ and $y(k+1)$; the last row therefore contributes only as a left-hand side. Three equations result: $$\begin{aligned} k = 0:&\quad 300\,a + 75\,b = 387\\ k = 1:&\quad 387\,a + 30\,b = 375\\ k = 2:&\quad 375\,a + 15\,b = 324 \end{aligned}$$ Three equations in two unknowns: the system is over-determined and, because the data carry measurement noise, inconsistent. Solving any two of them exactly (the first pair gives $a = 0.8247$, $b = 1.861$) leaves the third badly violated, which is precisely why a least-squares solution is demanded rather than an exact one.
3. Put the equations in regressor form. Collecting the unknowns into $\theta = [\,a\ \ b\,]^{T}$, $$Y = \Phi\,\theta, \qquad \Phi = \begin{bmatrix} 300 & 75\\ 387 & 30\\ 375 & 15\end{bmatrix}, \qquad Y = \begin{bmatrix} 387\\ 375\\ 324\end{bmatrix}.$$ The least-squares estimate minimises $J(\theta) = \lVert \Phi\theta - Y\rVert^2$, whose gradient vanishes at the solution of the normal equations $\Phi^{T}\Phi\,\hat{\theta} = \Phi^{T}Y$.
4. Form and solve the normal equations. Accumulating the products, $$\Phi^{T}\Phi = \begin{bmatrix} 380\,394 & 39\,735\\ 39\,735 & 6\,750\end{bmatrix}, \qquad \Phi^{T}Y = \begin{bmatrix} 382\,725\\ 45\,135\end{bmatrix}.$$ The determinant is $380\,394 \times 6\,750 - 39\,735^2 = 9.881\times10^{8} \ne 0$, so the data are persistently exciting enough to identify both parameters. Solving the two-by-two system, $$\boxed{\;\hat{a} = 0.7989, \qquad \hat{b} = 1.9837\;}$$ so the identified model is $P(z) = 1.9837/(z - 0.7989)$.
5. Check the fit before trusting it. Substituting the estimates back gives one-step prediction residuals of $+1.45$, $-6.31$ and $+5.35$, with a sum of squares of $70.5$ against outputs of order 350 — residuals under 2 % of the signal, and they alternate in sign rather than drifting, which is what an unbiased fit to noisy data looks like. The estimated pole $z = 0.799$ lies inside the unit circle, so the identified model is stable and the steady-state calculation that follows is meaningful.
6. Part (b) — evaluate the DC gain. For a constant input the steady-state output of a discrete system is its transfer function evaluated at $z = 1$ (the $z$-plane image of DC), which for this first-order model is $$P(1) = \frac{\hat{b}}{1 - \hat{a}} = \frac{1.9837}{1 - 0.7989} = 9.865 .$$ The same number follows from setting $y(k+1) = y(k) = y_{ss}$ in the difference equation, which is the more physical route: $y_{ss}(1-\hat{a}) = \hat{b}\,u$.
7. Predict the steady-state output. With $u(k) = 2$ held constant, $$\boxed{\;y_{ss} = \frac{\hat{b}}{1-\hat{a}}\,u = 9.865 \times 2 = 19.73\;}$$ This is far below the 300–390 range of the measured data, which is consistent rather than alarming: the recorded run was a decay from a large initial condition under a rapidly falling input, so the model is being asked to extrapolate to a much smaller sustained drive. The prediction should be reported with that caveat attached — identification data should bracket the operating point at which the model is later used.