clear all;
clc;


% --- Charger les données du tableau Excel ---
opts = detectImportOptions('Log1.xlsx');
opts.VariableNamingRule = 'preserve';
data = readtable('Log1.xlsx', opts);

% --- Initialisation des données ---
dt = 0.2; % Intervalle de temps (secondes)
I = data.Current; % Courant (A)
V_meas = data.Voltage; % Tension mesurée (V)

% Tableau OCV-SOC
soc_table = ([100, 90, 80, 70, 60, 50, 40, 30, 20, 10, 0]) / 100;
ocv_table = [4.1269, 4.0825, 4.0515, 3.9422, 3.8545, 3.7481, ...
             3.6186, 3.5414, 3.44, 3.2358, 2.7565];

% Paramètres de la batterie
Q_nom = 18000; % Capacité nominale (Coulombs)
soc_initial = 0.6466; % SOC initial (100 %)
soc = soc_initial + cumsum(I) * dt / Q_nom; % Calcul du SOC
soc(soc < 0) = 0; % Limiter SoC à 0 %
soc(soc > 1) = 1; % Limiter SoC à 100 %

% Interpolation pour V_OC
Voc = interp1(soc_table, ocv_table, soc, 'linear', 'extrap');

% --- Initialisation des paramètres ---
% Initialisation des paramètres
params = [0.023, 0.05985, 7345.99, 0.0068, 1073.25]; % [R0, R1, C1, R2, C2]
lb = [0.025, 0.05, 500, 0.005, 500];% Limites inférieures
ub = [0.1, 0.06, 7350, 0.007, 1100];% Limites supérieures


t = (0:length(I)-1)' * dt; % Vecteur temps (colonne)

% --- Fonction pour le modèle 2RC avec mise à jour exponentielle ---
model_2RC = @(params, Voc, I, dt) ...
    calculate_vsim(params, Voc, I, dt);

% --- Boucle d'identification récursive ---
max_iterations = 10; % Limite de sécurité pour convergence
tolerance = 1e-2;
learning_rate = 0.01; % Taux d'apprentissage pour la correction des paramètres
prev_params = params;
error_history = zeros(max_iterations, 1); % Pour tracer les erreurs

for iter = 1:max_iterations
    % Étape 1 : Simulation de la tension prédite
    V_sim = model_2RC(prev_params, Voc, I, dt);
    
    % Étape 2 : Calcul de l'erreur de prédiction
    error = V_meas - V_sim; % \Delta Z


    error_norm = norm(error);
    error_history(iter) = error_norm;
    
    % Étape 3 : Ajustement des paramètres
    params_opt = lsqcurvefit(@(params, t) model_2RC(params, Voc, I, dt), ...
                             prev_params, t, V_meas, lb, ub);

    % Étape 4 : Correction des paramètres en fonction de l'erreur
    params_opt(1) = params_opt(1) + learning_rate * mean(error) * sign(params_opt(1)); % Mise à jour de R0
    params_opt(2) = params_opt(2) + learning_rate * mean(error) * sign(params_opt(2)); % Mise à jour de R1
    params_opt(3) = params_opt(3) + learning_rate * mean(error) * sign(params_opt(3)); % Mise à jour de C1
    params_opt(4) = params_opt(4) + learning_rate * mean(error) * sign(params_opt(4)); % Mise à jour de R2
    params_opt(5) = params_opt(5) + learning_rate * mean(error) * sign(params_opt(5)); % Mise à jour de C2
    
    % Vérification de convergence
    if norm(params_opt - prev_params) < tolerance
        fprintf('Convergence atteinte à l''itération %d\n', iter);
        break;
    end
    
    % Ajout de la vérification sur la variation de l'erreur
    if iter > 1 && abs(error_history(iter) - error_history(iter-1)) < 1e-1
        fprintf('Convergence atteinte après %d itérations.\n', iter);
        break;
    end

    % Mise à jour pour l'itération suivante
    prev_params = params_opt;
    
    % Affichage des résultats intermédiaires
    fprintf('Itération %d : R0 = %.4f, R1 = %.4f, C1 = %.4f, R2 = %.4f, C2 = %.4f\n', ...
        iter, params_opt(1), params_opt(2), params_opt(3), params_opt(4), params_opt(5));
end

% --- Résultats finaux ---
R0 = params_opt(1);
R1 = params_opt(2);
C1 = params_opt(3);
R2 = params_opt(4);
C2 = params_opt(5);

% Calcul des constantes de temps
tau1 = R1 * C1; % Constante de temps RC1
tau2 = R2 * C2; % Constante de temps RC2

fprintf('Paramètres finaux estimés :\n');
fprintf('R0 = %.4f Ω\n', R0);
fprintf('R1 = %.4f Ω\n', R1);
fprintf('C1 = %.4f F\n', C1);
fprintf('R2 = %.4f Ω\n', R2);
fprintf('C2 = %.4f F\n', C2);
fprintf('Tau1 = %.4f s, Tau2 = %.4f s\n', tau1, tau2);

% --- Visualisation des résultats ---
figure;
plot(t, V_meas, 'b', 'LineWidth', 1.5); hold on;
plot(t, V_sim, 'r--', 'LineWidth', 1.5);
xlabel('Temps (s)');
ylabel('Tension (V)');
legend('Tension mesurée', 'Tension simulée');
title('Comparaison des tensions mesurée et simulée après identification');
grid on;

figure;
plot(t, V_meas, 'b', 'LineWidth', 1.5); hold on;
plot(t, Voc, 'g--', 'LineWidth', 1.5);
plot(t, V_sim, 'r:', 'LineWidth', 1.5);
xlabel('Temps (s)');
ylabel('Tension (V)');
legend('Tension mesurée', 'Tension OCV', 'Tension simulée');
title('Décomposition des tensions simulées');
grid on;

% Calcul et affichage de l'erreur au fil du temps
figure;
plot(t, error, 'm', 'LineWidth', 1.5);
xlabel('Temps (s)');
ylabel('Erreur (V)');
title('Erreur entre la tension mesurée et simulée');
grid on;


% --- Fonction de calcul de la tension simulée ---
function V_sim = calculate_vsim(params, Voc, I, dt)
    % Initialisation des tensions dynamiques
    V_RC1 = 0; % Tension initiale sur RC1
    V_RC2 = 0; % Tension initiale sur RC2
    V_sim = zeros(size(I)); % Tension simulée

    for k = 2:length(I)
        % Mise à jour exponentielle pour RC1 et RC2
        V_RC1 = params(2) * (1 - exp(-dt / (params(2) * params(3)))) * I(k);
        V_RC2 = params(4) * (1 - exp(-dt / (params(4) * params(5)))) * I(k);
        
        % Tension simulée
        V_sim(k) = Voc(k) - params(1) * I(k) - V_RC1 - V_RC2;
    end
end