Finding discrete matrix $A$ from time continuous data - System identification

Question

I'm seeking a way to find discrete $A$ from $P_k = [CA^kB, CA^kM]$ inside a hankel matrix.

I can show you an example to find $A$ if we only have $CA^kB$ for $i = 1, 2, 3, 4, 5, \dots, n$

Assume that we have the data points $Y_k = CA^kB$ for $i = 1, 2, 3, 4, 5, 6, \dots, n$

A = [0  1;
    -1 -3];

C = [1 0];

B = [0;
     1];

%% CAB     
k = 1:50;
Yk = zeros(1, 50);
for i = k
  Yk(i) = C*A^i*B;
end

And our goal is to find discrete $A$ from the data $Y_k$.

Then we can set up the hankel matrices

$$H_j = \begin{bmatrix} Y_1 & Y_2 & Y_3 & Y_3 & \dots & Y_{n/2} \\ Y_2 & Y_3 & \vdots & \vdots & \vdots & Y_{n/2+1}\\ Y_3 & Y_4 & \vdots & \vdots & \vdots & Y_{n/2+2} \\ Y_4 & \vdots & \vdots & \vdots & \vdots & Y_{n/2+3}\\ \vdots & \vdots & \vdots & & & \vdots \\ Y_{n/2} & Y_{n/2 +1} & Y_{n/2+2} & Y_{n/2+3} & \dots & Y_{n} \end{bmatrix}$$

%% Hankel matrix
H0 = hank(Yk, 1);
H1 = hank(Yk, 2);

% Create the half square hankel matrix
function [H] = hank(g, k)
  % We got markov parameters g = [g0 g1 g2 g2 g3 ... gl]; with size m*m. g0 = D
  m = size(g, 1);
  n = size(g, 2);
  l = length(g)/(m*2);
  H = zeros(l*m, l*m);
  for i = 1:l
    if(and(i == l, k == 2))
      % This is a special case when g is out of index, just add zeros instead!
      row = g(:, 1 + (k+i-1)*m:(k+i-2)*m + l*m);
      H(1 + (i-1)*m:(i-1)*m + m, 1:l*m) = [row zeros(m, m)]; 
   else
      row = g(:, 1 + (k+i-1)*m:(k+i-1)*m + l*m);
      H(1 + (i-1)*m:(i-1)*m + m, 1:l*m) = row;
    end
  end
end

Notice that if we have noise, we can reduce it a little bit.

%% Reduce noise
R0 = H0*H0';
R1 = H1*H0';

Now we using SVD on $R_0$ and cut the matrices with dimension $nx$ with is the same dimension of square matrix $A$.

$$[U, S, V^T] = svd(R_0)$$

% Do SVD on R0
[U,S,V] = svd(R0, 'econ');

%% Cut U, S, V with dimension nx = 2
nx = 2;
Un = U(:, 1:nx);
Sn = S(1:nx, 1:nx);
Vn = V(:, 1:nx);

And now we can find our dicrete matrix $A$

$$A_d = S_n^{-1/2}U_n^TR_1V_nS_n^{-1/2}$$

%% Find A-matrix
Ad = Sn^(-1/2)*Un'*R1*Vn*Sn^(-1/2)

Also we can find $C$ and $B$ as well, by using pseudo inverse

$$P_a = U_nS_n^{1/2}$$ $$X = P_a^{+}H_0$$ $$B_d = X(1:nx, 1:nx)$$ $$C_d = P_a(1:nx, 1:nx)$$

  Pa = Un*Sn^(1/2);
  X = pinv(Pa)*H0; 
  Bd = X(1:nx, 1:2:nu);
  K = X(1:nx, 2:2:nu);
  % From Pa we can find Cd
  Cd = Pa(1:nu, 1:nx);

Question:

If we have $P_k = [CA^kB, CA^kM]$ instead of $Y_k = CA^kB$. How would we find discrete $A$ from $P_k$?

Literature:

See equation 29 here and see equation 18 here

Hint:

Read function [sysd, K] = eradcokid(g, sampleTime, delay, systemorder) from here. It will exlain that how you can at least find $B_d$ and $C_d$. But not $A_d$.

Answer 1

I solved it!

Have a look at okid.m file here