3
$\begingroup$

I'm using subspace identification to identify a black-box state space model. To identify I follow these steps:

  1. We have measured a vector of inputs $u_k \in \Re^{m}$ and outputs $y_k \in \Re^{l}$. For SISO models, $m = l = 1$.

  2. Create hankel matrices for $U_p, U_f, U_p^+, U_f^-, Y_p, Y_f, Y_p^+, Y_f^-$ where $p$ means past and $f$ means future.

$$U_p = \begin{bmatrix} u_0 & u_1 & \dots & u_{j-1}\\ u_1 & u_2 & \dots & u_{j}\\ \vdots & \vdots& \ddots & \vdots\\ u_{i-1} & u_i & \dots & u_{i+j-2} \end{bmatrix}$$

$$U_p = \begin{bmatrix} u_i & u_{i+1} & \dots & u_{i+j-1}\\ u_{i+1} & u_{i+2} & \dots & u_{i+j}\\ \vdots & \vdots& \ddots & \vdots\\ u_{i+h-1} & u_{i+h} & \dots & u_{i+h+j-2} \end{bmatrix}$$

$$U_p^+ = \begin{bmatrix} u_i & u_{i+1} & \dots & u_{i+j-1}\\ u_{i+1} & u_{i+2} & \dots & u_{i+j}\\ \vdots & \vdots& \ddots & \vdots\\ u_{i+h-1} & u_{i+h} & \dots & u_{i+h+j-2} \\ u_{i} & u_{i+1} & \dots & u_{i+j-1} \end{bmatrix}$$

$$U_f^- = \begin{bmatrix} u_{i+1} & u_{i+2} & \dots & u_{i+j}\\ \vdots & \vdots& \ddots & \vdots\\ u_{i+h-1} & u_{i+h} & \dots & u_{i+h+j-2} \end{bmatrix}$$

Do exactly the same with $Y_p, Y_f, Y_p^+, Y_f^-$. Just replace $u_k$ with $y_k$. Matlab/Octave code:

  ny = size(y, 1); % Number of outputs
  nu = size(u, 1); % Number of inputs

  % Create hankel matrix for input
  Up = hank(u, i, j, 1); % h = 1
  Uf = hank(u, i, j, 2); % h = 2
  Upplus = [Up; Uf(nu,:)]; % Past plus 
  Ufminus = Uf((1+nu):end,:); % Future minus

  % Create hankel matrix for input
  Yp = hank(y, i, j, 1); % h = 1
  Yf = hank(y, i, j, 2); % h = 2
  Ypplus = [Yp; Yf(ny,:)]; % Past plus 
  Yfminus = Yf((1+ny):end,:); % Future minus

-

  function [H] = hank(g, i, j, h)
  % Create hankel matrix
  H = cell(i, j); 
  for i = 1:i
    for j = 1:j
      H{i,j} = g(:,h+i+j-2);
    end
  end

  % Cell to matrix
  H = cell2mat(H);
 end
  1. Create oblique projection matrix $\mathcal O$

$$\mathcal O_i = Y_f/_{U_f}W_p$$

$$\mathcal O_{i-1} = Y_f^-/_{U_f^-}W_p^+$$

Where $W_p = \begin{bmatrix} U_p\\ Y_p \end{bmatrix}, W_p^+ = \begin{bmatrix} U_p^+\\ Y_p^+ \end{bmatrix}$

Example:

$$ \underset{B}{A/C} = A \begin{bmatrix} C^\top & B^\top \end{bmatrix} \left( \begin{bmatrix} C \\ B \end{bmatrix} \begin{bmatrix} C^\top & B^\top \end{bmatrix} \right)^\dagger \begin{bmatrix} C \\ 0 \end{bmatrix} $$

The oblique projection matrix can be defined as: $$\mathcal O_i = \begin{bmatrix} C\\ CA\\ CA^2\\ \vdots\\ CA^{i-1} \end{bmatrix}\begin{bmatrix} B & AB & A^2B & \dots & A^{j-1}B \end{bmatrix}$$

Matlab/Octave code:

  % Create W matrix
  Wp = [Up; Yp];
  Wpplus = [Upplus; Ypplus];

  % Create oblique projection matrix
  O = obliqueprojectionmatrix(Yf, Uf, Wp)
  Ominus = obliqueprojectionmatrix(Yfminus, Ufminus, Wpplus)

-

function [O] = obliqueprojectionmatrix(A, B, C)
  O = A*[C' B']*pinv([C*C' C*B'; B*C' B*B'])*[C; 0*B]; % 0*B is added so we don't need to use "first r coloums"-metod
end
  1. Do singular value decomposition.

$$[U, S, V^T] = svd(\mathcal O_i)$$ $$\mathcal O_i = \begin{bmatrix} U_1 & U_2 \end{bmatrix}\begin{bmatrix} S_1 &0 \\ 0 & S_2 \end{bmatrix}\begin{bmatrix} V_1^T\\ V_2^T \end{bmatrix}$$

In literature you can find $$[U, S, V^T] = svd(W_1\mathcal O_i W_2)$$

Where $W_1, W_2$ are weighting matrices. But in this case, we don't need $W_1, W_2$ right now.

Matlab/Octave code:

[U, S, V] = svd(O, 'econ'); % [U1 U2]*[S1 0; 0 S2]*[V1'; V2']
  1. Do model reduction. Choose the state vector dimension $n_x$ depening on signular values.

$$U_1 = U(:, 1:n_x) \\ S_1 = S(1:n_x, 1:n_x)$$

Matlab/Octave code:

[U1, S1, V1, nx] = modelReduction(U, S, V);

-

function [U1, S1, V1, nx] = modelReduction(U, S, V)
  % Plot singular values 
  stem(1:length(S), 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!
  nx = inputdlg('Choose the state dimension by looking at hankel singular values: ');
  nx = str2num(cell2mat(nx));

  % Choose the dimension nx
  U1 = U(:, 1:nx);
  S1 = S(1:nx, 1:nx);
  V1 = V(:, 1:nx);
end
  1. Create the extended observability matrices $\Gamma_i, \Gamma_{i-1}$

$$\Gamma_i = U_1 S_1^{1/2}$$ $$\Gamma_{i-1} = \underline{\Gamma_i}$$

$\Gamma_{i-1}$ denotes the matrix $\Gamma_i$ without the last $l$ rows(number of outputs from $y_k \in \Re^{l}$)

Matlab/Octave code:

  % Create observability matrix
  Gamma = U1*sqrtm(S1);
  ny = size(y, 1); % Number of outputs
  Gammaminus = Gamma(1:end-ny, :);
  1. Create the state vector $x_i$ and $x_{i+1}$

$X_i$ can be defined as $\begin{bmatrix} B & AB & A^2B & \dots & A^{j-1}B \end{bmatrix}$

$$X_i = \Gamma_i^{\dagger}\mathcal O_i$$ $$X_{i+1} =\Gamma_{i-1}^{\dagger}\mathcal O_{i-1}$$

See step 3.

Matlab/Octave code:

  % Create state vectors
  X = pinv(Gamma)*O;
  Xplus = pinv(Gammaminus)*Ominus; 
  1. Now solve the matrices $A, B, C, D$ by using this least square $$\theta = YX^T(XX^T)^{-1}$$

$$\begin{bmatrix} X_{j+1}\\ Y_k \end{bmatrix}\begin{bmatrix} A & B\\ C & D \end{bmatrix}\begin{bmatrix} X_j\\ U_k \end{bmatrix}$$

$$\begin{bmatrix} A & B\\ C & D \end{bmatrix} = \begin{bmatrix} X_{j+1}\\ Y_k \end{bmatrix}\begin{bmatrix} X_{j}\\ U_k \end{bmatrix}^T(\begin{bmatrix} X_{j}\\ U_k \end{bmatrix}\begin{bmatrix} X_{j}\\ U_k \end{bmatrix}^T)^{-1} $$

Question:

The problem is that the length of $x_j$ and $u_k$ has not the same length and $x_{j+1}$ and $y_k$ has not the same length too. Because

$$U_k = {u_1, u_2, u_3, u_4,...,u_k}$$ $$Y_k = {y_1, y_2, y_3, y_4,...,y_k}$$ $$X_j = {x_1, x_2, x_3, x_4,...,x_j}$$ $$X_{j+1} = {x_2, x_3, x_4, x_5,...,x_{j+1}}$$

I know that is wrong to say $X_{j+1}, X_j$ because I haven't declare them. But if you look at step 3 and step 7. You will understand that the length of $X_i, X_{i+1}$ depends on $j$ only.

This is the problem. Mabey it's right, but I don't know how to find $A, B, C, D$ if the length differ.

Here is a demonstration how I use. Notice the length of X, Xplus, y, u. They are not equal.

Question -> Do them need to be equal?

>> u = 0.1:0.1:1;
>> y = sin(u)./u;
>> size(u)
ans =

    1   10

>> i = 5; j = 5;
>> dop(u, y, i, j, 0.1) % sampletime = 0.1
X =

  -1.0027e+00  -9.8765e-01  -9.6938e-01  -9.4805e-01  -9.2378e-01
   1.2002e-01   5.8200e-02  -3.2191e-03  -6.3907e-02  -1.2354e-01
  -2.0715e-03   8.1594e-04   2.2347e-03   1.4099e-03  -2.4157e-03

Xplus =

  -0.983448  -0.964376  -0.942236  -0.917154  -0.889272
   0.011996  -0.060384  -0.132032  -0.202623  -0.271841
   0.074765   0.018390  -0.040224  -0.101844  -0.167216

y =

   0.99833   0.99335   0.98507   0.97355   0.95885   0.94107   0.92031   0.89670   0.87036   0.84147

u =

    0.10000    0.20000    0.30000    0.40000    0.50000    0.60000    0.70000    0.80000    0.90000    1.00000

The whole function can be found here. Just copy this.

function [A, B, C, D] = dop(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 i - rows
  if(length(varargin) >= 3)
    i = varargin{3};
  else
    error('Missing amout of rows')
  end

  % Get j - columns
  if(length(varargin) >= 4)
    j = varargin{4};
  else
    error('Missing amout of columns')
  end

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

  % Get the delay
  if(length(varargin) >= 6)
    delay = varargin{6};
  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


  ny = size(y, 1); % Number of outputs
  nu = size(u, 1); % Number of inputs

  % Create hankel matrix for input
  Up = hank(u, i, j, 1); % h = 1
  Uf = hank(u, i, j, 2); % h = 2
  Upplus = [Up; Uf(nu,:)]; % Past plus 
  Ufminus = Uf((1+nu):end,:); % Future minus

  % Create hankel matrix for input
  Yp = hank(y, i, j, 1); % h = 1
  Yf = hank(y, i, j, 2); % h = 2
  Ypplus = [Yp; Yf(ny,:)]; % Past plus 
  Yfminus = Yf((1+ny):end,:); % Future minus

  % Create W matrix
  Wp = [Up; Yp];
  Wpplus = [Upplus; Ypplus];

  % Create oblique projection matrix
  O = obliqueprojectionmatrix(Yf, Uf, Wp);
  Ominus = obliqueprojectionmatrix(Yfminus, Ufminus, Wpplus);

  % Do SVD on CTRB matrix - nx is A matrix dimension
  [U, S, V] = svd(O, 'econ'); % [U1 U2]*[S1 0; 0 S2]*[V1'; V2']

  % Do model reduction
  [U1, S1, V1, nx] = modelReduction(U, S, V);

  % Create observability matrix
  Gamma = U1*sqrtm(S1);
  Gammaminus = Gamma(1:end-ny, :);

  % Create state vectors
  X = pinv(Gamma)*O
  Xplus = pinv(Gammaminus)*Ominus
  y
  u
  % Check if X, Xplus is equal to y and u


  % Use least square to find [A B; C D]
  ABCD = [X; u]\[Xplus; y]

end

function [H] = hank(g, i, j, h)
  % Create hankel matrix
  H = cell(i, j); 

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

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

function [O] = obliqueprojectionmatrix(A, B, C)
  O = A*[C' B']*pinv([C*C' C*B'; B*C' B*B'])*[C; 0*B]; % 0*B is added so we don't need to use "first r coloums"-metod
end

function [U1, S1, V1, nx] = modelReduction(U, S, V)
  % Plot singular values 
  stem(1:length(S), 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!
  nx = inputdlg('Choose the state dimension by looking at hankel singular values: ');
  nx = str2num(cell2mat(nx));

  % Choose the dimension nx
  U1 = U(:, 1:nx);
  S1 = S(1:nx, 1:nx);
  V1 = V(:, 1:nx);
end

The information about this subspace identification can be found from:

Subspace identification article

Or from this book "Overschee and De Moor (1996)" from : Dr. Reza Shadmehr's home page Book about system identification at Dr.Reza Shadmehr's home page

See page 55 at that book. What does $U_{i|i}$ mean?

$\endgroup$

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.