clear all;
clc;

% --- Charger les données du tableau Excel ---
% Remplacez 'data.xlsx' par le chemin vers votre fichier Excel
% Chargement des données depuis un fichier Excel
opts = detectImportOptions('test1_Log.xlsx');
opts.VariableNamingRule = 'preserve';
data = readtable('test1_Log.xlsx', opts);


% --- Données d'entrée ---
% Tableau OCV-SOC
soc_table = [100, 90, 80, 70, 60, 50, 40, 30, 20, 10, 0] / 100; % SOC en fraction
ocv_table = [4.1269, 4.0825, 4.0515, 3.9422, 3.8545, 3.7481, ...
             3.6186, 3.5414, 3.44, 3.2358, 2.7565]; % Tensions OCV (V)

% Paramètres du modèle 2RC
%params.R0 = 0.023;    
%params.R1 = 0.007296;    
%params.C1 = 3814.045;    
%params.R2 = 0.002747;   
%params.C2 = 1265.502;  

params.R0 = 0.023 ;
params.R1 = 0.05985;
params.C1 = 7345.99;
params.R2 = 0.0068;
params.C2 = 1073.25;
params.C_b = 18000;   % Capacité de la batterie en coulombs (Ah * 3600)

% Données d'échantillonnage
dt = 0.2; % Temps d'échantillonnage en secondes
soc_init = 0.872; % SOC initial (en fraction)

I = data.Current; % Courant en ampères
V_measured = data.Voltage; % Tension mesurée en volts

% --- Initialisations ---
n = length(I); % Nombre de points de données
SOC_coulomb = zeros(n, 1); % SOC estimé par comptage de Coulomb
SOC_corrige = zeros(n, 1); % SOC corrigé par la méthode modifiée
SOC_coulomb(1) = soc_init; % Initialisation SOC initial
SOC_corrige(1) = soc_init; % Initialisation SOC initial corrigé
V_OCV = zeros(n, 1); % Tension OCV estimée
V_R1 = 0; % Tension initiale sur le premier réseau RC
V_R2 = 0; % Tension initiale sur le deuxième réseau RC
SOC_OCV = zeros(n, 1); % SOC estimé par la tension OCV
I_c_real = zeros(n, 1); % Courant corrigé réel
errors_SOC = zeros(n, 1); % Erreur entre SOC_OCV et SOC_Coulomb


% --- Étape 1 : Simulation pas à pas pour OCV, SOC Coulomb et I_c_real ---
for k = 1:n
    % Mise à jour des tensions dynamiques (modèle 2RC)
    %V_R1 = V_R1 + (dt / (params.R1 * params.C1)) * (params.R1 * I(k) - V_R1);
    %V_R2 = V_R2 + (dt / (params.R2 * params.C2)) * (params.R2 * I(k) - V_R2);
    
    V_R1 = exp(-dt / (params.R1 * params.C1)) * V_R1 + params.R1 * (1 - exp(-dt / (params.R1 * params.C1))) * I(k); % Mise à jour de la tension V_RC1 en utilisant une équation exponentielle
    V_R2 = exp(-dt / (params.R2 * params.C2)) * V_R2 + params.R2 * (1 - exp(-dt / (params.R2 * params.C2))) * I(k); % Mise à jour de la tension V_RC2 en utilisant une équation exponentielle

    % Estimation de la tension OCV
    V_OCV(k) = V_measured(k) + params.R0 * I(k) + V_R1 + V_R2;

    % Estimation du SOC à partir de V_OCV (interpolation SOC-OCV)
    SOC_OCV(k) = interp1(ocv_table, soc_table, V_OCV(k), 'linear', 'extrap');

    % Calcul du SOC par comptage de Coulomb
    if k > 1
        SOC_coulomb(k) = SOC_coulomb(k-1) + (I(k) * dt) / params.C_b;
        % Limitation du SOC entre 0 % et 100 %
        if SOC_coulomb(k) > 1
            SOC_coulomb(k) = 1;
        elseif SOC_coulomb(k) < 0
            SOC_coulomb(k) = 0;
        end
        
    end

    
    % Calcul de l'erreur entre SOC_OCV et SOC_Coulomb
    errors_SOC(k) = SOC_OCV(k) - SOC_coulomb(k);

    % Calcul du courant corrigé réel I_c_real
    I_c_real(k) = I(k) +  errors_SOC(k)*2; % profil (charge et decharge ou bien seulement charge)
    %I_c_real(k) = I(k) +  errors_SOC(k)*4; % profil decharge
end

% --- Étape 2 : Interpolation pour déterminer les coefficients k2, k1, k0 ---
coeffs = polyfit(I, I_c_real, 2); % Ajustement quadratique entre I et I_c_real
k2 = coeffs(1);
k1 = coeffs(2);
k0 = coeffs(3);

% --- Étape 3 : Calcul du courant corrigé I_c avec les coefficients ---
I_c = k2 * I.^2 + k1 * I + k0;

% --- Étape 4 : Calcul du SOC corrigé avec la méthode modifiée ---
for k = 2:n
    SOC_corrige(k) = SOC_corrige(k-1) + (I_c(k) * dt) / params.C_b;
    % Limitation du SOC corrigé entre 0 % et 100 %
    if SOC_corrige(k) > 1
        SOC_corrige(k) = 1;
    elseif SOC_corrige(k) < 0
        SOC_corrige(k) = 0;
    end
     
end
% --- Calcul de l'erreur RMSE ---
% RMSE pour la méthode classique (Coulomb)
rmse_coulomb = sqrt(mean((SOC_OCV - SOC_coulomb).^2));

% RMSE pour la méthode modifiée (Coulomb modifié)
rmse_corrige = sqrt(mean((SOC_OCV - SOC_corrige).^2));


% Conversion en pourcentage des estimations de l'EKF (entre 0 et 100 %)
SOC_corrige = SOC_corrige * 100;  % Conversion en pourcentage

% Méthode Coulomb Counting
SOC_coulomb = SOC_coulomb * 100;
SOC_OCV = SOC_OCV * 100;
% Affichage des résultats complets avec SoC Coulomb après calcul
for k = 1:n
    fprintf('Étape %d: Courant = %.4f A, SoC estimé SOC_corrige = %.4f%%, SoC estimé Coulomb = %.4f%%\n', ...
        k, I(k), SOC_corrige(k), SOC_coulomb(k));
end



% --- Résultats ---
disp('Coefficients déterminés par interpolation quadratique :');
disp(['k2 = ', num2str(k2)]);
disp(['k1 = ', num2str(k1)]);
disp(['k0 = ', num2str(k0)]);

% Afficher les résultats
disp('Erreur RMSE :');
disp(['Méthode Coulomb classique : ', num2str(rmse_coulomb)]);
disp(['Méthode Coulomb modifiée : ', num2str(rmse_corrige)]);


% --- Tracés --

% SOC estimés
figure;
plot(1:n, SOC_coulomb, '--', 'LineWidth', 1.5);
hold on;
plot(1:n, SOC_corrige, ':', 'LineWidth', 1.5);
xlabel('Temps (échantillons)');
ylabel('SOC (%)');
title('Coulomb et Coulomb corrigée');
legend('SOC Coulomb', 'SOC Corrigé');
grid on;

% Courant corrigé vs mesuré
figure;
plot(1:n, I, '--', 'LineWidth', 1.5);
hold on;
plot(1:n, I_c, 'LineWidth', 1.5);
xlabel('Temps (échantillons)');
ylabel('Courant (A)');
title('Courant mesuré et corrigé');
legend('I mesuré', 'I corrigé');
grid on;