# Identify a state space model from measured inputs and outputs

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


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