Exercise13

13.1 Model Reduction and Model Reduction Error

Find balanced state space realizations  for the following systems $G\left(s\right)\:=\:C\left(sI-A\right)^{-1}B$ and truncate them to find suitable approximations $\widehat{G}$ that make  $|G(i\omega) -\widehat{G}(i\omega)|$small for all $\omega$ (hint: balreal+modred). Plot Bode diagrams for $G$ and $\widehat{G}$ to illustrate the results for different model orders.

a) $A = \begin{pmatrix} -1 & 0 & 0 \\ 0 & -0.1 & 0 \\ 0 & 0 & -0.101 \end{pmatrix}, \quad B = \begin{pmatrix} 1 \\ 1 \\ -1 \end{pmatrix}, \quad C = \begin{pmatrix} 1 & 1 & 1 \end{pmatrix}$

b) $A = \begin{pmatrix} -1 & 0 & 0 \\ 0 & -0.1 & 0 \\ 0 & 0 & -0.101 \end{pmatrix}, \quad B = \begin{pmatrix} 1 \\ 1 \\ 1 \end{pmatrix}, \quad C = \begin{pmatrix} 1 & 1 & 1 \end{pmatrix}$

c) $G(s) = \dfrac{(s+2)(s+4)(s+6)(s+8)}{(s+1)(s+3)(s+5)(s+7)}$

d) It is possible to show that the truncation error lies between the following bounds

$\sigma_{r+1} \leq |G(i\omega)-\widehat{G}(i\omega)| \leq 2 (\sigma_{r+1} + \ldots \sigma_{n})$      (the "twice the tail-bound")

The bound is valid if $A$ is stable, and you are keeping the first $r$ states and have $\sigma_{r} > \sigma_{r+1}$.

Compare these bounds with the result you got in c).

Hint: Either use a Bode  plot of $G-\widehat{G}$, or use the function hinfnorm(H) which calculates $\sup_\omega|H(i\omega)|$.

13.2 Balanced Realization Transformation

Show that the coordinate transformation $x=T\bar x$ changes matrices the following way

a) $\bar A = T^{-1} A T$, $\quad \bar B = T^{-1}B$, $\quad \bar C = C T, \quad \bar D = D$ (state-space matrices)

b) $\bar P = T^{-1}PT^{-T}, \quad \bar Q = T^TQT$ (reachability and observability Gramians)

c) Show that the following procedure transforms the system to a form where $\bar P = \bar Q = \Sigma$

• Compute Gramians $P$ and $Q$ for the system given by $(A,B,C,D)$
• Compute a matrix $R$ so that $P = RR^T$ (this is called Cholesky factorization)
• Compute a SVD of $R^TQR$ and write it on the form $R^TQR = U\Sigma^2 U^T$
• Use the coordinate transformation $x=T\bar x$ with $T = RU\Sigma^{-1/2}$

13.3 Grey box Identification of a Heated Rod

The file greyrod.m on Lecture 12 does Grey Box Identification of the parameters kappa and htf. A compartment model of the system is produced by heatd.m.

Try different model orders 1,2,5 and 100 and compare the result of the grey box identification. Study both the Bode diagrams, and the estimated values of kappa and htf and comment on the results.

 (You will need greyrod.m Download greyrod.m and heatd.m Download heatd.m )

13.4 State Space Identification with measurable states

If the full state vector $x$ is known, then the state space model $(A,B,C,D)$ can be found by standard least squares regression in

$\begin{bmatrix} x_2 & \cdots & x_{N+1} \\ y_1 & \cdots & y_{N} \end{bmatrix} = \begin{bmatrix} A &B \\ C & D \end{bmatrix} \begin{bmatrix} x_1 & \cdots & x_{N} \\ u_1 & \cdots & u_{N} \end{bmatrix} + \textrm{noise}$

a) Study the file ex13_4.m Download ex13_4.m

b) Suggest a practical method to reduce the number of nonzero elements in $(A,B,C,D)$. Such a method can be useful if you suspect that many elements in theses matrices are zero. (You don't need to implement it).

-----------------------------------

Solutions:

Matlab files with solutions

ex 13_1: See the code ex13_1.m Download ex13_1.m

ex 13_2:  See this solution exercise-13solutions.pdf Download exercise-13solutions.pdf 

ex 13_3: See the code ex13_3.m Download ex13_3.m (you will need greyrod.m Download greyrod.m and heatd.m Download heatd.m )

ex 13_4: 

a) The code in the exercise implements solutions to the normal equations for a linear matrix equation of the form

Measurement-matrix = Theta-matrix * Regressor-matrix

in the code the  Measurement-matrix is called "lhs", for "left hand side", and Regressor-matrix is called "rhs".

We have talked about such a generalization of linear regression to the matrix case at an earlier exercise, and concluded then that the corresponding normal equations has similar structure as before. In the code the solution is obtained both from such a normal equation, and from matlab's internal implementation of solution to over and underdetermined linear equations (which often have better numerics)

b) One can use the "lasso" method from an early lecture in the course, adding an L_1 penalty on the sum of all coefficients in the A,B,C and D matrices. The simplest way to implement this would probably be to vectorize the matrix equation, i.e. rewrite the given linear equations system in the problem into the traditional form $y = \Phi \theta$, where $\theta$ would be a parameter vector with all the matrix elements of A,B,C and D.