2
$\begingroup$

This article is from 1989 and it contains an algorithm how to estimate a state space model from arbitrary input and output.

The algorithm begin with to create two hankel matrices.

$$H_1 = \begin{bmatrix} u [ k ] & u [ k + 1 ] & \dots & \dots & \dots & u[ k + j - 1 ] \\ y [ k ] & y [ k + 1 ] & \dots & \dots & \dots & y[ k + j - 1 ] \\ u [ k + 1 ] & u [ k + 2 ] & \dots & \dots & \dots & u[ k + j ] \\ y [ k + 1 ] & y [ k + 2 ] & \dots & \dots & \dots & y[ k + j ] \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ u [ k + i -1] & u [ k + i ] & \dots & \dots & \dots & u[ k + j + i -2 ] \\ y [ k + i -1 ] & y [ k + i ] & \dots & \dots & \dots & y[ k + j + i -2 ] \\ \end{bmatrix}$$

$$H_2 = \begin{bmatrix} u [ k + i ] & u [ k + i + 1 ] & \dots & \dots & \dots & u[ k +i + j - 1 ] \\ y [ k + i ] & y [ k + i + 1 ] & \dots & \dots & \dots & y[ k +i + j - 1 ] \\ u [ k + i + 1 ] & u [ k +i + 2 ] & \dots & \dots & \dots & u[ k +i + j ] \\ y [ k + i + 1 ] & y [ k +i + 2 ] & \dots & \dots & \dots & y[ k +i + j ] \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ u [ k + 2i -1] & u [ k + 2i ] & \dots & \dots & \dots & u[ k + j + 2i -2 ] \\ y [ k + 2i -1 ] & y [ k + 2i ] & \dots & \dots & \dots & y[ k + j + 2i -2 ] \\ \end{bmatrix}$$

$u[k]$ is the input and has the size $m x p$ where $m$ is the dimension (number of inputs) and $p$ is the measured data length.

$y[k]$ is the output and has the size $l x p$ where $l$ is the dimension (number of outputs) and $p$ is the measured data length.

So the hankel matrices got the size $H \in \Re^{mi x lj}$ where $lj$ is the column length of the matrix. If $l = 1$ and $m = 1$ then the row and length of the matrix will be $mxj$.

Here is the matlab code for computing the hankel matriecs.

% Get hankel matrices
i = 2;
H1 = hank([u;y], i - 1); % Indexing starting at 1.
H2 = hank([u;y], i);

I would allays set $i = j$ because then the hankel matrix row will be exactly twice as large as the hankel matrix column.

function [H] = hank(g, k)
  % Create hankel matrix
  H = cell(length(g)/2,length(g)/2); 

  for i = 1:length(g)/2
     for j = 1:length(g)/2
        H{i,j} = g(:,k+i+j-2);
     end
  end

  % Cell to matrix
  H = cell2mat(H);
end

Next step is to compute the singular value decomposition of those hankel matrices.

$$svd([H_1; H_2]) = \begin{bmatrix} U_{11} & U_{12} \\ U_{21} & U_{22} \end{bmatrix} \begin{bmatrix} S_{11} & 0\\ 0 & 0 \end{bmatrix}V^T$$

Matlab code

[U,S,V] = svd([H1;H2], 'econ');

Those matrices got the dimension $$dim ( U_{11} ) = ( mi + li ) x ( 2mi + n ) \\ dim ( U_{12} ) = ( mi + li ) x ( 2li - n ) \\ dim ( U_{21} ) = ( mi + li ) x ( 2mi + n ) \\ dim ( U_{22} ) = ( mi + li ) x ( 2li - n ) \\ dim ( S_{11} ) = ( 2mi + n ) x ( 2mi + n ) \\$$

Where $n$ is the state dimension $x$.

I have solved this problem by using this:

% Do model reduction
[U11, U12, U21, U22, S11, n] = modelReduction(U, S, V); % 

This function splits up all the $U$ matrices and $S$ matrecies and then scale down those to a $nxn$ dimension depending on user choise.

function [U11, U12, U21, U22, S11, n] = modelReduction(U, S, V)
  % Plot singular values 
  stem(diag(S));
  title('Hankel Singular values');
  xlabel('Amount of singular values');
  ylabel('Value');

  % Choose system dimension n - Remember that you can use modred.m to reduce some states too!
  n = inputdlg('Choose the state dimension by looking at hankel singular values: ');
  n = str2num(cell2mat(n));

  % Find the size
  [um, un] = size(U);
  [sm, sn] = size(S);

  % Split them up
  U11 = U(1:um/2, 1:un/2);
  U12 = U(1:um/2, (1 + un/2):un);
  U21 = U((1 + um/2):um, 1:un/2);
  U22 = U((1 + um/2):um, (1 + un/2):un);
  S11 = S(1:sm/2, 1:sn/2);

  % Do reduction - No need to do reduction for Uq later
  U11 = U11(1:n, 1:n);
  U12 = U12(1:n, 1:n);
  U21 = U21(1:n, 1:n);
  U22 = U22(1:n, 1:n);
  S11 = S11(1:n, 1:n);
end

The next step is to compute this:

$$U_{12}^TU_{11}S_{11} = \begin{bmatrix} U_q & U_q^{\perp} \end{bmatrix}\begin{bmatrix} S_q & 0\\ 0& 0 \end{bmatrix}\begin{bmatrix} V_q^T\\ Vq^{\perp t} \end{bmatrix}$$

My matlab code is:

 % Compute Uq
 [Uq,Sq,Vq] = svd(U12'*U11*S11, 'econ');

I don't think we need to do any reduction for $U_q$ because I all ready done it.

Next step we going to find matrecies $A, B, C, D$.

$$\begin{bmatrix} x [ k + i + 1 ]& \dots &x[ k + i + j - 1 ] \\ y [ k + i ] & \dots & y [ k + i + j - 2 ] \end{bmatrix} = \begin{bmatrix} A & B \\ C & D \end{bmatrix} \begin{bmatrix} x [ k + i ]& \dots &x[ k + i + j - 2 ] \\ u [ k + i ] & \dots & y [ k + i + j - 2 ] \end{bmatrix}$$

Matlab code for this is:

  i = n/(l+m); % Important! Or else the algorithm won't work - My own modification :)
  dxk = Uq'*U12'*U((m + l + 1 ):( i + 1 )*( m + l ), :)*S;
  yk  = U((m*i + l*i + m + 1):(m + l)*(i + 1), :)*S;
  xk  = Uq'*U12'*U(1:(m*i + l*i ), :)*S;
  uk  = U((m*i + l*i + 1):(m*i + l*i + m), :)*S;

  % Solve matrices by using least square
  ABCD = [dxk; yk]/[xk; uk];

  % Extract the system matrices
  Ad = ABCD(1:n, 1:n);
  Bd = ABCD(1:n, (n + 1):(n + m));
  Cd = ABCD((n + 1):(n + l), 1:n);
  Dd = ABCD((n + 1):(n + l), (n + 1):(n + m));

For the testing. Assume that we have a step response which look like this:

enter image description here

Matlab code

>> x =  0:0.1:15;
>> y = 30*(1-exp(-x).*cos(x));
>> u = linspace(1, 1, length(x));

Then we try the code. I choose the state vector dimension $n = 6$.

enter image description here

Ok. It works great! Lets say that we find $y2$ from

>> y2 = 30*(1-exp(-x).*cos(x.^2));

enter image description here

Now I use both $y$ and $y2$ into the algorithm. I choose the state vector dimension $n = 6$.

enter image description here

The problem is that this is the best results I can have. I don't know why. So now to my question:

Is it the right method to split up the matrecies $U$:

  U11 = U(1:um/2, 1:un/2);
  U12 = U(1:um/2, (1 + un/2):un);
  U21 = U((1 + um/2):um, 1:un/2);
  U22 = U((1 + um/2):um, (1 + un/2):un);
  S11 = S(1:sm/2, 1:sn/2);

Or do I need to use this?

$$dim ( U_{11} ) = ( mi + li ) x ( 2mi + n ) \\ dim ( U_{12} ) = ( mi + li ) x ( 2li - n ) \\ dim ( U_{21} ) = ( mi + li ) x ( 2mi + n ) \\ dim ( U_{22} ) = ( mi + li ) x ( 2li - n ) \\ dim ( S_{11} ) = ( 2mi + n ) x ( 2mi + n ) \\$$

Because I got some trobbles to choose the right dimension. That's why I asking you if it's the right method I have been choose?

You can find the complete algorithm here

function [Ad, Bd, Cd, Dd] = mdvv(varargin)
  % Check if there is any input
  if(isempty(varargin))
    error('Missing imputs')
  end

  % Get input 
  if(length(varargin) >= 1)
    u = varargin{1};
  else
    error('Missing input')
  end

  % Get output 
  if(length(varargin) >= 2)
    y = varargin{2};
  else
    error('Missing output')
  end

  % Get the sample time
  if(length(varargin) >= 3)
    sampleTime = varargin{3};
  else
    error('Missing sample time');
  end

  % Get the delay
  if(length(varargin) >= 4)
    delay = varargin{4};
  else
    delay = 0; % If no delay was given
  end

  % Check if u and y has the same length
  if(length(u) ~= length(y))
    error('Input(u) and output(y) has not the same length')
  end

  % Get number of inputs (m) and outputs (l)
  m = size(u, 1); 
  l = size(y, 1);

  % Get hankel matrices
  i = 2;
  H1 = hank([u;y], i - 1); % Indexing starting at 1.
  H2 = hank([u;y], i);

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

  % Do model reduction
  [U11, U12, U21, U22, S11, n] = modelReduction(U, S, V); % 

  % Compute Uq
  [Uq,Sq,Vq] = svd(U12'*U11*S11, 'econ');

  % Find the states and input/outputs vectors
  i = n/(l+m); % Important! Or else the algorithm won't work - My own modification :)
  dxk = Uq'*U12'*U((m + l + 1 ):( i + 1 )*( m + l ), :)*S;
  yk  = U((m*i + l*i + m + 1):(m + l)*(i + 1), :)*S;
  xk  = Uq'*U12'*U(1:(m*i + l*i ), :)*S;
  uk  = U((m*i + l*i + 1):(m*i + l*i + m), :)*S;

  % Solve matrices by using least square
  ABCD = [dxk; yk]/[xk; uk];

  % Extract the system matrices
  Ad = ABCD(1:n, 1:n);
  Bd = ABCD(1:n, (n + 1):(n + m));
  Cd = ABCD((n + 1):(n + l), 1:n);
  Dd = ABCD((n + 1):(n + l), (n + 1):(n + m));

end


function [U11, U12, U21, U22, S11, n] = modelReduction(U, S, V)
  % Plot singular values 
  stem(diag(S));
  title('Hankel Singular values');
  xlabel('Amount of singular values');
  ylabel('Value');

  % Choose system dimension n - Remember that you can use modred.m to reduce some states too!
  n = inputdlg('Choose the state dimension by looking at hankel singular values: ');
  n = str2num(cell2mat(n));

  % Find the size
  [um, un] = size(U);
  [sm, sn] = size(S);

  % Split them up
  U11 = U(1:um/2, 1:un/2);
  U12 = U(1:um/2, (1 + un/2):un);
  U21 = U((1 + um/2):um, 1:un/2);
  U22 = U((1 + um/2):um, (1 + un/2):un);
  S11 = S(1:sm/2, 1:sn/2);

  % Do reduction - No need to do reduction for Uq later
  U11 = U11(1:n, 1:n);
  U12 = U12(1:n, 1:n);
  U21 = U21(1:n, 1:n);
  U22 = U22(1:n, 1:n);
  S11 = S11(1:n, 1:n);
end


function [H] = hank(g, k)
  % Create hankel matrix
  H = cell(length(g)/2,length(g)/2); 

  for i = 1:length(g)/2
    for j = 1:length(g)/2
      H{i,j} = g(:,k+i+j-2);
    end
  end

  % Cell to matrix
  H = cell2mat(H);
end
$\endgroup$
  • $\begingroup$ I do not think that $cos(t^2)$ can be a part of the step response of any LTI system. So it is probably impossible to identify a system which has an identical step response. $\endgroup$ – Kwin van der Veen Dec 27 '17 at 8:09
  • $\begingroup$ So this algorithm is only for dermististic systems and not stochastic systems ? $\endgroup$ – Daniel Mårtensson Dec 27 '17 at 10:17
  • $\begingroup$ By the way. I have other algorithms which can estimate a system which looks like this. I think there is some problems with my algorithm above. Something with choosing matrecies U11, U12, U21 etc. $\endgroup$ – Daniel Mårtensson Dec 27 '17 at 10:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.