MSI Ex2 Matlab code

Dateiinhalt

%% Exercise 2 - Least Squares estimation
%%

close all
clear all

N = 20; % number of measurements
M = 1000; %number of experiments
R = 1; % true value of the resistance
E = 2; % true value of the voltage source
sigma = 0.2; % standard deviation used in the generation of the normal noise

%% Data generation - Case 1

i_min = 5; % minimum value of the input current
i_max = 10; % maximum value of the input current
n_u = randn(N,M)*sigma; % measurement noise
i = linspace(i_min,i_max,N)'; % input current
u = E+R*i*ones(1,M)+n_u; % voltage measurements

%% Least Squares estimator - Case 1
R_est1 = zeros(1,M);
E_est1 = zeros(1,M);

for ind = 1:M
    K = [i, ones(N,1)]; % coefficient matrix
    theta = (K'*K)\(K'*u(:,ind)); % solution of the LS problem
    R_est1(ind) = theta(1); % estimated value of R
    E_est1(ind) = theta(2); % estimated value of E
end

figure
hold on

plot(E_est1, R_est1, 'r+');
title('Resistance and voltage estimation')
xlabel('E');
ylabel('R');

%% Data generation - Case 2

i_min = -10; % minimum value of the input current
i_max = 10; % maximum value of the input current
n_u = randn(N,M)*sigma; % measurement noise
i = linspace(i_min,i_max,N)'; % input current
u = E+R*i*ones(1,M)+n_u; % voltage measurement
n_u = randn(N,M)*sigma; % measurement noise
i = linspace(i_min,i_max,N)'; % input current
u = E+R*i*ones(1,M)+n_u; % output voltage

%% Least Squares estimator - Case 2
R_est2 = zeros(1,M);
E_est2 = zeros(1,M);

for ind = 1:M
    K = [i, ones(N,1)]; % coefficient matrix
    theta = (K'*K)\(K'*u(:,ind)); % solution of the LS problem
    R_est2(ind) = theta(1); % estimated value of R
    E_est2(ind) = theta(2); % estimated value of E
end

plot(E_est2, R_est2, 'b+');

%% Covariance matrices computation and confidence ellipses

R_mean1 = mean(R_est1); % mean value of R for case 1
E_mean1 = mean(E_est1); % mean value of E for case 1
Sigma1 = cov([E_est1' R_est1']); % covariance matrix for case 1
[V1,D1] = eig(Sigma1); % eigenvalues and eigenvector computation

R_mean2 = mean(R_est2); % mean value of R for case 2
E_mean2 = mean(E_est2); % mean value of E for case 2
Sigma2 = cov([E_est2' R_est2']); % covariance matrix for case 2
[V2,D2] = eig(Sigma2); % eigenvalues and eigenvector computation

% generate the coordinates of 50 points on a unit circle
xy=[cos(linspace(0,2*pi,50)); sin(linspace(0,2*pi,50))];

% generate the points of the confidence ellipse
xy_ellipse1 = [E_mean1; R_mean1]*ones(1,50) + V1*sqrt(D1)*xy;
xy_ellipse2 = [E_mean2; R_mean2]*ones(1,50) + V2*sqrt(D2)*xy;

plot(xy_ellipse1(1,:), xy_ellipse1(2,:), 'k');
plot(xy_ellipse2(1,:), xy_ellipse2(2,:), 'g');