First of all. This is not a programming question. I do not requesting help about programming. I'm requesting help about if I have made correct math steps here. Second, I got some issues with my web browser so I it can't display LaTeX for the moment. I'm working on it.
There is a method kalled Observer Kalman Filter Identification(OKID). The pourpose is that this method can estimate a kalman filter gain matrix from pure data. No need to use Algebraic Riccati Equations to compute the kalman gain matrix.
There is a paper by Jer-Nan Juang, NASA scientist. He is the author of Eigensystem Realization Theory(ERA). Created 1985. The ERA algorithm want to have a impulse response, e.g so called markov parameters. And the OKID can turn regular random input- and output data into a impulse response.
I have created an ERA MATLAB file to share with you. https://github.com/DanielMartensson/Mataveid/blob/master/sourcecode/era.m
The bad thing is that this ERA still need to have some work on the MIMO system option. And that is where OKID comes in handy. http://people.duke.edu/~hpgavin/SystemID/References/Juang+Phan+etal-JGCD-1993.pdf
or here at page 10: http://people.duke.edu/~hpgavin/SystemID/CourseNotes/EigensystemRealization.pdf
I have created a file named OKID3 in MATLAB and I want to share it with you, beacuse I need help to make it work.
function [Y] = okid3 (u, y)
% Get the dimensions first
q = size(y, 1); % Dimension of output
l = size(y, 2)/2; % Lenght of output and input and cut it to half
N = size(y, 2); % Total length
m = size(u, 1); % Dimension of input
% Create upper triangular V matrix and Ybar
V = zeros((q+m)*l + m, l);
% Begin with the first m rows for V
for k = 1:m
V(k, 1:l) = u(k, l+1:N);
end
% Now do the rest (q+m)*l rows for V. We want to implement v = [u;y] into V - This is equation 6
positionrow = 1;
for row = m:(q+m):(q+m)*l
% For u
for k = 1:m
V(row + k, 1:l) = u(k, l-positionrow+1:N-positionrow);
end
% For y
for k = 1:q
V(row + k + m, 1:l) = y(k, l-positionrow+1:N-positionrow);
end
positionrow = positionrow + 1;
end
% Now we have
% V = [u(k) u(k+1) ... u(N);
% v(k-1) v(k) ... v(N-1);
% v(k-2) v(k-1) ... v(N-2);
% . . . .
% . . . .
% . . . .
% . . . .
% v(0) v(1) ... v(N-k)]
% Important to have half of y
y_half = zeros(q, 1:l);
for k = 1:q
y_half(k, 1:l) = y(k, l+1:N);
end
% Solve for non-filtred markov parameters
[Q, R] = qr(V');
Ybar = linsolve(R, Q'*y_half')';
% Get D matrix
D = zeros(q, m);
for k = 1:q
D(k, 1:m) = Ybar(k, 1:m); % This is YBar_{-1}
end
% Remove D from Ybar
YbarNoD = zeros(q, (q+m)*l-m-m);
for i = 1:q
YbarNoD(i, 1:(q+m)*l-m-m) = Ybar(i, m+1:(q+m)*l-m);
end
% Time to filter markov parameters. Solve the system like this
% [I 0 0 0 0 [Y0 [Y1_0 + Y2_0*D
% Y0 I 0 0 0 Y1 Y1_1 + Y2_1*D
% Y1 Y0 I 0 0 * Y2 = Y1_2 + Y2_2*D
% Y2 Y1 Y0 I 0 Y3 Y1_3 + Y2_3*D
% ................ ... ..............
% Yn-1 Yn-2 Yn-3 Yn-4 ... Y0 I] Yn] Y1_n + Y2_n*D]
A = zeros(q*l-q, q*l-q);
% Create the diagonal
for i = 1:q*l-q
A(i, i) = 1;
end
% Insert into A - This is equation 20
for k = 0:(l-1)-2
[Ybar1_i, Ybar2_i] = getYbar1_2(YbarNoD, q, m, k);
% Insert into diagonal
for o = k:(l-1)-2
for i = 1:q
for j = 1:q
posX = (o+1)*q; % One step down to begin with, just to avoid the diagonals
posY = o*q-k*q; % One step to left and for every iteration, do less steps to posY
A(i+posX, j+posY) = -Ybar2_i(i, j);
end
end
end
end
% Create b matrix and b0 matrix - This is equation 20 and 28
b = zeros(q*l-q, q);
for k = 0:l-2
[Ybar1_i, Ybar2_i] = getYbar1_2(YbarNoD, q, m, k);
Ybar1Ybar2D = Ybar1_i + Ybar2_i*D;
for i = 1:q
for j = 1:q
b(i + k*q, j) = Ybar1Ybar2D(i, j);
end
end
end
% This is the impulse response
Y = linsolve(A, b);
end
% This will get [Ybar1, Ybar2]
function [Ybar1, Ybar2] = getYbar1_2(YbarNoD, q, m, k)
Ybar1 = zeros(q, m); % C(A+MC)^k(B+MD)
Ybar2 = zeros(q, q); % -C(A+MC)^kM
Ybar1Ybar2 = zeros(q, m+q);
% Do Ybar1Ybar2
for i = 1:q
Ybar1Ybar2(i, 1:m+q) = YbarNoD(i, 1+(q+m)*k: (q+m)*k + q+m);
end
% Split Ybar1Ybar2 to Ybar1 and Ybar2
for i = 1:q
Ybar1(i, 1:m) = Ybar1Ybar2(i, 1:m);
Ybar2(i, 1:q) = Ybar1Ybar2(i, m+1:q+q);
end
end
I will explain every part it do.
This part creates the $V$ matrix in the report:
% Now do the rest (q+m)*l rows for V. We want to implement v = [u;y] into V - This is equation 6
positionrow = 1;
for row = m:(q+m):(q+m)*l
% For u
for k = 1:m
V(row + k, 1:l) = u(k, l-positionrow+1:N-positionrow);
end
% For y
for k = 1:q
V(row + k + m, 1:l) = y(k, l-positionrow+1:N-positionrow);
end
positionrow = positionrow + 1;
end
This part create the lower triangular toeplitz matrix.
% Create the diagonal
for i = 1:q*l-q
A(i, i) = 1;
end
% Insert into A - This is equation 20
for k = 0:(l-1)-2
[Ybar1_i, Ybar2_i] = getYbar1_2(YbarNoD, q, m, k);
% Insert into diagonal
for o = k:(l-1)-2
for i = 1:q
for j = 1:q
posX = (o+1)*q; % One step down to begin with, just to avoid the diagonals
posY = o*q-k*q; % One step to left and for every iteration, do less steps to posY
A(i+posX, j+posY) = -Ybar2_i(i, j);
end
end
end
And this code part will create the right part of equation 20 in the paper.
for k = 0:l-2
[Ybar1_i, Ybar2_i] = getYbar1_2(YbarNoD, q, m, k);
Ybar1Ybar2D = Ybar1_i + Ybar2_i*D;
for i = 1:q
for j = 1:q
b(i + k*q, j) = Ybar1Ybar2D(i, j);
end
end
end
I get the markov parameters by splitting a combination of markov parameters.
% This will get [Ybar1, Ybar2]
function [Ybar1, Ybar2] = getYbar1_2(YbarNoD, q, m, k)
Ybar1 = zeros(q, m); % C(A+MC)^k(B+MD)
Ybar2 = zeros(q, q); % -C(A+MC)^kM
Ybar1Ybar2 = zeros(q, m+q);
% Do Ybar1Ybar2
for i = 1:q
Ybar1Ybar2(i, 1:m+q) = YbarNoD(i, 1+(q+m)*k: (q+m)*k + q+m);
end
% Split Ybar1Ybar2 to Ybar1 and Ybar2
for i = 1:q
Ybar1(i, 1:m) = Ybar1Ybar2(i, 1:m);
Ybar2(i, 1:q) = Ybar1Ybar2(i, m+1:q+q);
end
end
Let's do a simulation:
When I use Y = okid3(u, y) then I get these markov parameters back. This is not what it shold look like. The markov parameters should be smooth and not jumpy. If you want to try this algorithm. Just insert input data and output data that have a length that can be divided with 2 e.g length 200 or 100, for 300. Not 201, 299 etc..
Example:
>> u = 1:12; y = 20:31;
>> Y = okid3([u;2*u], [y;2*y])
Question:
I have work with this OKID algorithm for many days from now and I still don't get it why the markov parameters jump. There fore I asking you if it's the correct way to use this equation:
that cause the jumping, and I should use this instead?
