Non-parametric black-box identification of I/O models

Question 1

Possible representations of a discrete-time, linear, dynamical system

Answer 1

1. state space representation
2. transfer function
3. impulse response

Question 2

transformation from state-space to transfer function

Answer 2

state space: F,G,H

tf: W(z) = H * (zI - F)^-1 *G

Question 3

transformation from transfer function to state-space

Answer 3

• infinite realizations of a tf into a ss model
• control realization:

W(z) = (b0z^n-1 + b1z^n-2 + … + bn-1) / (z^n + a1*z^n-1 + … + an)

F = [0 1 0 … 0
0 0 1 … 0
…
-an -an-1 … -a1]

G = [0 0 0 … 1]’

H = [bn-1 bn-2 … b0]

D = 0

Question 4

transformation from t.f. to impulse response

Answer 4

• long division (tf in negative powers)
  or
• geometric series

Question 5

transformation from i.r. to t.f.

Answer 5

Z transform

not applicable in practice because an infinite number of samples of the i.r. should be used

Question 6

transformation from s.s. to i.r.

Answer 6

ω(t) = 0 if t=0 (and D=0)

H F^t-1 G if t > 0

Question 7

definition of full observability

Answer 7

A system F,G,H is fully observable if its observability matrix
O = [H
        HF
        HF^2
        …
        HF^n-1]
is full rank

Question 8

definition of reachability

Answer 8

A system F,G,H is fully reachable if its reachability matrix
R = [G FG F^2G … F^n-1G]
is full rank

Question 9

definition of Hankel matrix of order n

Answer 9

Given an impulse response
{ω(0), ω(1), …, ω(N)}

Hn = [ω(1) ω(2) ω(3) … ω(n)
         ω(2) ω(3) ω(4) … ω(n+1)
         ω(3) ω(4) ω(5) … ω(n+2)
          …
          ω(n) ω(n+1)      … ω(2n-1)] = 
= O * R

Question 10

factorization of the Hankel matrix

Answer 10

Hn = O * R

Henkel matrix of order n can be factorized into the product of the observability and reachability matrices of the system

Question 11

algorithm to estimate F,G,H from a noise-free impulse response

Answer 11

Step 1.

• build the Hankel matrix in increasing the order
• each time check its order
• if Hn and Hn+1 have the same rank, stop
• n is the estimated order of the system

step 2.
- factorize Hn+1 into two rectangular matrices: Hn+1 = On+1 * Rn+1,
where On+1 and Rn+1 are the extended observability and reachability matrices

step 3:
H^ = On+1(1;:) first row of On+1
G^ = Rn+1(:;1) first column of Rn+1
O1 = On+1(1:n;:)
O2 = On+1(2:n+1;:)
F^ = O1^-1 * O2

Question 12

observations on the method for noise-free estimation

Answer 12

• simple constructive method
• it uses only 2n-1 samples
• it is not parametric
• it is useless, because when noise is present, step 1 never stops, and even if n is known, the results are wrong

Question 13

4 SID procedure using noisy measurement of the i.r.

Answer 13

step 1: build the Hankel matrix using all the dataset
H~qd = [ω(1) ω(2) ω(3) … ω(d)
ω(2) ω(3) ω(4) … ω(d+1)
ω(3) ω(4) ω(5) … ω(d+2)
…
ω(q) ω(q+1) … ω(q+d-1)]

step 2: Singular Value Decomposition of H~qd:
H~qd = U~ S~ V~’
U~ (qq) and V~ (dd) are unitary matrices;
S~ (q*d) is a rectangular matrix with the singular values of H~qd in the diagonal, stored in decreasing order
σ1 > = σ2 > = σ3 > = … > = σq

step 3: definition of the order n of the system and of the matrices U^ S^ V^’
- from the plot of the singular values:
- in an ideal case there is a clear jump after n values: the last index before the jump is the estimated order of the system
- in a real case, there is not a clear jump but a knee; the value of n should be selected inside that range
- once the value of n has been decided,
U^: first n columns of U~
S^: n x n square diagonal matrix extracted from the left high corner of S~ (it contains the first n singular values)
V^’ = first n rows of V~’
- > H~qd = H^qd + Hresqd
where H^qd = U^ S^ V^’ is a q x d matrix with rank n (reduced from d)
Hresqd has rank q

step 4: estimation of H^, G^, F^
- H^qd = U^ S^ V^’ = U^ S^ ^1/2 S^ ^1/2 V^’
- Definition of the extended observability and reachability matrices
O^ = U^ S^ ^1/2
R^ = S^ ^1/2 V^’

- H^ = O^(1;:) first row of O^
G^ = R^(:;1) first column of R^
F^ estimated using the shift invariance property:
O^1 = O^(1:q-1 ; :)
O^2 = O^(2:q ; :)
O^1 F = O^2
since O^1 is not squared, the linear system is over-determined
- > least squares solution:
F^ = (O^1' O^1)^-1 O^1' O^2

Question 14

On the choice of q and d

Answer 14

q < d
d+d-1 = N to use the entire data-set

• q ~ d : best accuracy of the method, high computational effort
• q < < d : lowest computational effort
• if q > 1/2 d the sensitivity of the choice is small
• > q = 1/2 d as a rule of thumb

Question 15

def unitary matrix

Answer 15

a square matrix M is unitary if
det(M) = 1
M^-1 = M^T

Question 16

def singular values

Answer 16

real, positive numbers
they are a sort of “eigenvalues” for a rectangular matrix
if M is rectangular:
SV(M) = sqrt(EIG(MM^T)) = sqrt(EIG(M^TM))

Question 17

on the optimality of the procedure

Answer 17

SVD performs the optimal rank reduction, in the sense that the residual matrix is minimum in the Frobenius norm

Question 18

def Frobenius norm

Answer 18

Hresqd |F = sqrt( sum(i,j) Hresqd(i,j)^2 )

Question 19

on the choice of n

Answer 19

choice of n is a supervised step, but it can be automatized using a cross validation approach