% Simple algorithm to estimate R0 of the battery based on Gregory 
% Plett book Battery management systems: volume II

clear all;
clc;

%% Load data to analyze

load("Time_test_val.mat");
load("Vbat_test_val.mat");
load("Ibat_test_val.mat");

% Here is the plot of the current and voltage
figure(); plot(Time_test_val/3600, Ibat_test_val);
title('Current during the test'); grid on;
xlabel('Time(hours)'); ylabel('Current(A)');
figure(); plot(Time_test_val/3600, Vbat_test_val);
title('Voltage during the test'); grid on;
xlabel('Time(hours)'); ylabel('Voltage(v)');

%% Initialization of the real value of R0

R0 = 5e-3; % This is the real value of R0 obtained from the ident. process
R0_vect = R0*ones(1, length(Vbat_test_val));

%% Estimation of R0 without any filter

% You have to define a threshold for the calculation

threshold = 5;

for k=2:length(Vbat_test_val)
    num = Vbat_test_val(k)-Vbat_test_val(k-1);
    den = Ibat_test_val(k-1)-Ibat_test_val(k);
    R0_hat(k) = abs(num/den);
    if (abs(den)<threshold)
        R0_hat(k) = R0_hat(k-1);
    end
end

% Here is the plot of the unfilteres version
figure(); plot(Time_test_val/3600,1000*R0_vect,Time_test_val/3600, R0_hat*1000);
title('Unfiltered R0 estimate'); grid on;
xlabel('Time(hours)'); ylabel('R0(ohm)');
legend('Actual value', 'Estimate','location','southeast');

%% Estimation of R0 with a filter

% Here we define an alpha value for the filter
alpha = 0.999;
R0_hat_filt = R0_hat;

for k=2: length(Vbat_test_val)
    R0_hat_filt(k) = alpha*R0_hat_filt(k-1) + (1-alpha)*R0_hat(k);
end

% Here we plot the filtered R0 estimate
figure(); plot(Time_test_val/3600,1000*R0_vect,Time_test_val/3600, R0_hat_filt*1000);
title('Filtered R0 estimate'); grid on;
xlabel('Time(hours)'); ylabel('R0(ohm)');
legend('Actual value', 'Estimate','location','southeast');