I'm using subspace identification to identify a black-box state space model. To identify I follow these steps:
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$.
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
- 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
- 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']
- 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
- 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, :);
- 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;
- 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?