How can i solve this "Matrix dimensions must agree" error?
2 次查看(过去 30 天)
显示 更早的评论
I am getting error while run this code.
error:
Error using +
Matrix dimensions must agree.
Error in Untitled (line 77)
omega(:, :, i+1) = omega(:, :, i)+10 * x_hat2(:, i) * conj(e(i))';
clear all; close all; clc;
% Isometric model Kalman filter simulation Kalman Filter
%% Initial condition
ts = 1; % Sampling interval
t = [0:ts:100]; % Sampling moment
T = length(t); % Sampling points
%% Initial state
x = [0 40 0 20]';
%% Process Noise Covariance Process Noise Co - point Patrograph
q = 5; % Process noise covariance
Q1 = q*eye(4); % Process noise covariance matrix
%% Measurement Noise Covariance Observation Noise Co-Tariation Matrix
r = 5; % Observation noise covariance
Q2 = r*eye(2); % Observation noise covariance matrix
%% Process and measurement noise
v1 = sqrt(Q1)*randn(4,T); % Process noise sequence dimension:4 x T
v2 = sqrt(Q2)*randn(2,T); % Observing noise sequence dimension:2 x T
%% Estimate error covariance initialization
p = 5; % Estimation error covariance matrix
P(:,:,1) = p*eye(4); % Estimated covariance matrix
%==========================================================================
%% Continuous time status space model
A = [0 1 0 0;
0 0 0 0;
0 0 0 1;
0 0 0 0];
B = [0 0;
1 0;
0 0;
0 1];
C = [1 0 0 0;
0 0 1 0];
D = [1 0;
0 1];
%% Discrete time status space model
% F:State transfer matrix
% H:Observation matrix
sysc = ss(A,B,C,D); % Generate continuous time status space model
sysd = c2d(sysc, ts, 'zoh'); % Continuous time state space model to discrete time status space model
[F G C I] = ssdata(sysd); % Remove the discrete time status space model; process equation F:State transfer matrix; observation equationH:Observation matrix
%% Practice state of target
for i = 1:T-1
x(:,i+1) = F*x(:,i);
end
x = x + v1; % Estimation vector sequence under noise
y = C*x + v2; % Observation vector sequence under noise
%==========================================================================
%% Kalman Filter
% bookP179 algorithm5.4.1 Kalman Adaptive Filter Algorithm
K = P;
x_hat = [0 0 0 0].';
% x_hat = x;
for i = 1:T-1
G = F * K * C'* inv(C * K * C' + Q2); % Book P179 formula5.4.24Defined algorithm5.4.1Present a calculation method
alpha = y(:, i) - C * x_hat(:, i);
x_hat(:, i + 1) = F * x_hat(:, i) + G * alpha; % Calculate the current state vector x_hat
P = K - inv(F) * G * C * F;
K = F * P * F' + Q1; % Update Kalman Gain
end
%% LMS adaptive filter
% bookP183 algorithm5.5.1 LMSAdaptive filtering and basic deformation
omega(:, :, 1) = zeros(4, 4);
x_hat2(:, 1) = x(:, 1);
for i = 1:T - 1
x_hat2(:, i + 1) = omega(:, :, i)' * x_hat2(:, i); % Get the current estimate
e(:, i) = x(:, i + 1) - x_hat2(:, i + 1); % Calculate the error of the estimated value and the actual value
alpha = 1 / (0.5 + x(:, i)' * x(:, i));
omega(:, :, i + 1) = omega(:, :, i) + 10 * x_hat2(:, i) * conj(e(i))'; % Update weights
end
%==========================================================================
%% Estimate error
x_error = x-x_hat;
%% Graph 1 practical and tracking position
figure(1)
plot(x(1,:),x(3,:),'r');
hold on;
plot(x_hat(1,:),x_hat(3,:),'g.');
title('2D Target Position')
legend('Practical Position','Tracking Position')
xlabel('X axis [m]')
ylabel('Y axis [m]')
hold off;
%% Graph 2
figure(2)
plot(t,x(1,:)),grid on;
hold on;
plot(t,x_hat(1,:),'r'),grid on;
plot(t, x_hat2(1, :), 'k'), grid on
title('Practical and Tracking Position on X axis')
legend('Practical Position', 'Kalman Tracking Position', 'LMS Tracking Position')
xlabel('Time [sec]')
ylabel('Position [m]')
hold off;
%% Graph 3
figure(3)
plot(t,x_error(1,:)),grid on;
title('Position Error on X axis')
xlabel('Time [sec]')
ylabel('Position RMSE [m]')
hold off;
%% Graph 4
figure(4)
plot(t,x(2,:)),grid on;
hold on;
plot(t,x_hat(2,:),'r'),grid on;
plot(t,x_hat2(2,:),'k'),grid on;
title('Practical and Tracking Velocity on X axis')
legend('Practical Velocity','Tracking Velocity')
xlabel('Time [sec]')
ylabel('Velocity [m/sec]')
hold off;
%% Graph 5
figure(5)
plot(t,x_error(2,:)),grid on;
title('Velocity Error on X axis')
xlabel('Time [sec]')
ylabel('Velocity RMSE [m/sec]')
hold off;
How can i solve this error?
回答(1 个)
Jan
2021-7-6
编辑:Jan
2021-7-6
omega(:, :, i + 1) = omega(:, :, i) + 10 * x_hat2(:, i) * conj(e(i))';
This works since Matlab R2016b due to the auto-expanding. Do you use an older version? Then you need bsxfun:
omega(:, :, i + 1) = bsxfun(@plus, omega(:, :, i), 10 * x_hat2(:, i) * conj(e(i))';
You calculate alpha one line above, but do not use it anywhere.
What is the purpose of conj(e(i))' ? The matrix e(:, i) is calculated in the loop also. So do you really want to use e(i) and not e(:, i)? The conjugate transpose operator ' is strange also for a scalar. In addition you conjugate twice by conj and ' , so conj(e(i))' is the same as e(i). Remember that the transposition is .' including a dot.The quote ' is the conjugate transposition already. Maybe you want (bold guessing and not tested):
omega(:, :, i + 1) = omega(:, :, i) + 10 * x_hat2(:, i) * e(:, i)';
0 个评论
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!