dimension error that prevents estimation
1 view (last 30 days)
Show older comments
Gostaria de realizar estimativas, em que eu necessito fazer um comparativo entre as minhas saídas y e minhas saídas estimadas e somente consigo a da primeira interação. Alguém pode ajudar? Segue o meu código até então. As I do not know how to program very well, forgive me any grotesque mistake in the code and questions I make that seem simple!
ps: I edited because I solved the problem and found another one that...
clc; close all; clear;dbstop if error;
R = 0.3456; %---->K
T = 0.9746;
Z = -1.6557;
% a_0 = [0 0 0]';
Q = 0.876; %----> cov e(k)
H =0;
y = [ 1 2 3 4 5 6 7 8 9 10 11 34 23 55 6;
2 3 8 1 1 3 5 6 7 4 12 8 5 6 3;
6 5 2 12 2 4 5 6 8 9 1 2 1 3 9];
[p,q] = size(y);
%inicicalização dados
A = eye(size(T,1));
P_oo = A*A';
P_finito_inicial = zeros(size(T,1));
P_finito = P_finito_inicial;
a_conhecido = zeros(size(T,1),1);
a_inicial = a_conhecido;
a_o = a_inicial;
k=1;
a_predicao(:,:,k) = zeros(q,size(T,1));
var_erro_finito = zeros(q,size(T,1));
inovacao_predicao=zeros(q,size(y,2));
inovacao_predicao(:,1) = y(1,:)'-Z*a_o(:,:,k);
saida_predicao(:,:,k) = zeros(p,q);
saida_predicao(1,:) = Z*a_predicao(1,:,:);
% Inicializacao Exata
for k =1:q
% Testa F_oo
F_oo = Z*P_oo*Z';
if all(F_oo(:)== 0) % testa para zero
%F e M
F_finito = Z*P_finito*Z' + H;
F_ini(:,:,k) = [zeros(size(F_finito )) zeros(size(F_finito ))];
M_finito = P_finito*Z';
K_o = (T*M_finito)/F_finito;
L_o = T-K_o*Z;
L_ini(:,:,k) = [L_o zeros(size(L_o))];
v_o = y(k,:)' - Z*a_o;
a_o = T*a_o +K_o*v_o; % nova estimativa
P_oo = T*P_oo*T';
P_finito(:,:,k+1) = T*P_finito*L_o' + R*Q*R'; % nova covariancia
elseif rank(F_oo)== min(size(F_oo,1),size(F_oo,2)) % testa nao singularidade
F_finito = Z*P_finito*Z' + H;
M_oo = P_oo*Z';
M_finito = P_finito*Z';
F_2 = -(F_oo\F_finito)/F_oo; %F_1 = inv(F_oo);
F_1 = eye(size(F_2))/F_oo;
F_ini(:,:,k) = [F_1 F_2];
K_o = (T*M_oo)*F_1;
K_1 = (T*M_finito)*F_1 +T*M_oo*F_2;
L_o = T-K_o*Z;
L_1 = -K_1*Z;
L_ini(:,:,k) = [L_o L_1];
v_o = y(k,:)' - Z*a_o;
a_o = T*a_o(:,:,k) +K_o*v_o; %erro
P_oo = T*P_oo*L_o';
else
error('Matriz F_oo pode ser singular')
end
%% Guarda valores
a_predicao(:,k) =a_o;
saida_predicao(k,:) = Z*a_predicao(:,k);
var_erro_finito(:,k) = diag(P_finito);
var_erro_oo(:,k) = diag(P_oo); % guarda var P_oo;
P_predicao(:,:,k) = [P_oo P_finito];
if all(P_oo<=1e-7)==1
d = k;
break;
end
end
%% FK preditor normal
a_estimado = a_predicao(:,d);
P = P_finito;
for k = 1:q
F = Z*P*Z' + H;
F_ref(:,:,k) = F; %guarda F para refinamento
K = T*P*Z'/F;
%L
L = T-K*Z;
% L_ref(:,:,k) = L;
%Inovacao
v= y(k,:)' - Z*a_estimado;
%% Atualizacao
a_estimado= T*a_estimado +K*v; % nova estimativa
P = T*P*L' + R*Q*R'; % nova covariancia
%% Guarda valores
a_predicao(:,k) =a_estimado;
saida_predicao(k,:) = Z*a_predicao(:,k);
inovacao_predicao(:,k) = v;
var_erro_finito(:,k) = diag(P);
P_predicao(:,:,k) = [P_oo P];
end
t = 0:1:q;
figure(1)
subplot(3,1,1)
plot(t,y(1,:),'r',t,saida_predicao(1,:)','-.b')
subplot(3,1,2)
plot(t,y(2,:),'r', t,saida_predicao(2,:)','-.b')
subplot(3,1,3)
plot(t,y(3,:),'r', t,saida_predicao(3,:)','-.b')
0 Comments
Answers (1)
Walter Roberson
on 9 Jul 2019
y = [ 1 2 3 4 5 6 7 8 9 10 11 34 23 55 6;
2 3 8 1 1 3 5 6 7 4 12 8 5 6 3;
6 5 2 12 2 4 5 6 8 9 1 2 1 3 9];
[p,q] = size(y);
so q is the number of columns in y
t = 0:1:q;
so t has q+1 entries. For example if q was 5 then t = [0 1 2 3 4 5], for a total of 6 entries.
plot(t,y(1,:),'r',t,saida_predicao(1,:)','-.b')
t has one more entry than y has columns; it also has one more entry than saida_predicao has columns.
5 Comments
Walter Roberson
on 9 Jul 2019
y(:,k)
is 3 x 1. y(:,k)' is 1 x 3. Z*a_estimado is 15 x 1. y(:,k)' - Z*a_estimado is an error in R2016a and earlier, but in R2016b it is like bsxfun(@minus, y(:,k).', Z*a_estimado) which would give you 15 x 3. That is not going to fit into a_predicao(:,k)
You need to rethink what you are doing. - Z*a_estimado is 15 x 1 so you probably want to the term before that to be 15 x 1 giving you a 15 x 1 result. But that would require indexing a row of y, and there are only 3 rows of y. Perhaps your for loop should be k=1:p instead of k=1:q
See Also
Categories
Find more on Matrix Indexing in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!