Laboratory 3 - Estimation, filtering and system identification - Prof. M. Taragna

Exercise 1: Parametric estimation of a resistance from current-voltage data stored in resistor_data_1.mat

Contents

Introduction

The program code may be splitted in sections using the characters "%%". Each section can run separately with the command "Run Section" (in the Editor toolbar, just to the right of the "Run" button). You can do the same thing by highlighting the code you want to run and by using the button function 9 (F9). This way, you can run only the desired section of your code, saving your time. This script can be considered as a reference example.

clear all, close all, clc

Procedure

  1. Load the file resistor_data_1.mat containing the current-voltage data
  2. Plot the voltage measurements versus the current measurements
  3. Estimate the resistance using the Sample Mean algorithm
  4. Plot the estimated approximation versus the experimental data
  5. Estimate the resistance using the Least Squares algorithm
  6. Plot the estimated approximation versus the experimental data
  7. Estimate the resistance using the Gauss-Markov algorithm
  8. Plot the estimated approximation versus the experimental data
  9. Estimate the resistance using the Bayesian method
  10. Plot the estimated approximation versus the experimental data

Problem setup

% Step 1: load of data

load resistor_data_1.mat
% i = current measured in Ampere
% V = voltage measured in Volt

N = length(i); % N = number of data

% Step 2: plot of data

figure, plot(i,V,'ok'),
axis([min(i),max(i),0,35]), grid, title('Measured data and approximations'),
xlabel('Current{\it i} in A'), ylabel('Voltage{\it V} in V'),
legend('Measured data',4)

Estimation of the resistance using the Sample Mean algorithm

% Step 3: computation of the parameter estimate as mean of the R values
% derived from the Ohm's law for each couple of current-voltage data

R_SM = mean(V./i)   % Form #1: using the "mean" operator

R_SM_ = V'*(1./i)/N % Form #2: using the Sample Mean definition

% Step 4: graphical comparison of the results

V_SM = R_SM*i;
figure, plot(i,V,'ok', i,V_SM,'-b'),
axis([min(i),max(i),0,35]), grid, title('Measured data and approximations'),
xlabel('Current{\it i} in A'), ylabel('Voltage{\it V} in V'),
legend('Measured data','Sample Mean estimate',4)

R_SM =
    3.2324
R_SM_ =
    3.2324

Estimation of the resistance using the Least Squares algorithm

% Step 5: computation of the parameter estimate

Phi=i;
R_LS=Phi\V   % Form #1: using the "\" operator (more realiable)

A=pinv(Phi); % more realiable than inv(Phi'*Phi)*Phi'
R_LS_=A*V    % Form #2: using the pseudoinverse matrix
sigma_e = 5;
Sigma_ee = sigma_e^2*eye(N);
Sigma_R_LS = pinv(Phi)*Sigma_ee*pinv(Phi)'

% Step 6: graphical comparison of the results

V_LS = R_LS*i;
figure, plot(i,V,'ok', i,V_SM,'-b', i,V_LS,'-r'),
axis([min(i),max(i),0,35]), grid, title('Measured data and approximations'),
xlabel('Current{\it i} in A'), ylabel('Voltage{\it V} in V'),
legend('Measured data','Sample Mean estimate','Least Squares estimate',4)

R_LS =
    2.6979
R_LS_ =
    2.6979
Sigma_R_LS =
    0.0331

Estimation of the resistance using the Gauss-Markov algorithm

% Step 7: computation of the parameter estimate
% Note that the Gauss-Markov estimate coincides with the Least Squares
% estimate when the same weights are used for all the measurements

R_GM = inv(Phi'*inv(Sigma_ee)*Phi)*Phi'*inv(Sigma_ee)*V
Sigma_R_GM = inv(Phi'*inv(Sigma_ee)*Phi)

% Step 8: graphical comparison of the results

V_GM = R_GM*i;
figure, plot(i,V,'ok', i,V_SM,'-b', i,V_LS,'-r', i,V_GM,'--k'),
axis([min(i),max(i),0,35]), grid, title('Measured data and approximations'),
xlabel('Current{\it i} in A'), ylabel('Voltage{\it V} in V'),
legend('Measured data','Sample Mean estimate','Least Squares estimate',...
       'Gauss-Markov estimate',4)

R_GM =
    2.6979
Sigma_R_GM =
    0.0331

Estimation of the resistance using the Bayesian method

% Step 9: computation of the parameter estimate, assuming:
% - R_bar = a priori mean value of R = 3 (to be assumed early)
% - Sigma_RR = a priori variance of R = 10, 1, 0.1, 0.01 (to be assumed early)
% - V_bar = mean value of V = R_bar*i
% - Sigma_VV = Sigma_RR*i*i' + Sigma_ee
% - Sigma_RV = Sigma_RR*i'

R_bar = 3; % Try also with R_GM
Sigma_RR = [10, 1, 0.1, 0.01];
V_bar = R_bar*i;
for k=1:length(Sigma_RR),
    Sigma_VV = Sigma_RR(k)* i*i' + Sigma_ee;
    Sigma_RV = Sigma_RR(k)*i';
    R_B(k) = R_bar + Sigma_RV * inv(Sigma_VV) * (V-V_bar);
    Sigma_R_B(k) = Sigma_RV * inv(Sigma_VV) * Sigma_RV'; % estimate variance
    Sigma_R_minus_R_B(k) = Sigma_RR(k) - Sigma_R_B(k); % a posteriori variance
    V_B(:,k) = R_B(k)*i;
end
Sigma_RR, R_B, Sigma_R_B, Sigma_R_minus_R_B

% Step 10: graphical comparison of the results

figure, plot(i,V,'ok', i,V_SM,'-b', i,V_LS,'-r', i,V_GM,'--k', ...
             i,V_B(:,1),'-m', i,V_B(:,2),'--m', ...
             i,V_B(:,3),'-c', i,V_B(:,4),'--c'),
axis([min(i),max(i),0,35]), grid, title('Measured data and approximations'),
xlabel('Current{\it i} in A'), ylabel('Voltage{\it V} in V'),
legend('Measured data','Sample Mean estimate','Least Squares estimate',...
       'Gauss-Markov estimate', ...
       ['Bayesian estimate (\sigma_R^2=', num2str(Sigma_RR(1)),')'], ...
       ['Bayesian estimate (\sigma_R^2=', num2str(Sigma_RR(2)),')'], ...
       ['Bayesian estimate (\sigma_R^2=', num2str(Sigma_RR(3)),')'], ...
       ['Bayesian estimate (\sigma_R^2=', num2str(Sigma_RR(4)),')'],4)

Sigma_RR =
   10.0000    1.0000    0.1000    0.0100
R_B =
    2.6989    2.7076    2.7731    2.9300
Sigma_R_B =
    9.9670    0.9679    0.0751    0.0023
Sigma_R_minus_R_B =
    0.0330    0.0321    0.0249    0.0077


Published with MATLAB® R2014a