|
|
|
clear all;
|
|
clc;
|
|
|
|
|
|
% Chargement des données depuis un fichier Excel
|
|
% Chargement des données depuis un fichier Excel
|
|
opts = detectImportOptions('Log6.xlsx');
|
|
opts.VariableNamingRule = 'preserve';
|
|
data = readtable('Log6.xlsx', opts);
|
|
|
|
params.R0 = 0.023 ;
|
|
params.R1 = 0.05985;
|
|
params.C1 = 7345.99;
|
|
params.R2 = 0.0068;
|
|
params.C2 = 1073.25;
|
|
|
|
% Définition des paramètres de la batterie
|
|
%params.R0 = 0.0232;
|
|
%params.R1 = 0.0482;
|
|
%params.C1 = 7399.9982;
|
|
%params.R2 = 0.0032;
|
|
%params.C2 = 1265.502;
|
|
params.C_b = 18000; %capacité de la batterie
|
|
dt = 0.2; %temps d'échantillonage
|
|
soc_init = 0.872;
|
|
|
|
% Initialisation de l'état et de la covariance
|
|
x = [soc_init; 0; 0]; % SoC initial, V_RC1 et V_RC2 initialisés à 0
|
|
P = eye(3)*1e-6; % Matrice de covariance initiale avec régularisation
|
|
|
|
% Covariance des bruits
|
|
Q = eye(3)*1e-9; % Bruit de processus ajusté
|
|
R = 0.1; % Bruit de mesure ajusté
|
|
|
|
% Simulation du système avec données de courant et de tension observées
|
|
current = data.Current; %courant de décharge
|
|
voltage_measured = data.Voltage; %tension de décharge
|
|
soc_real = data.SOC; % SoC réel depuis le tableau Logpoly3, soc de la methode de notre client
|
|
n_steps = length(current);
|
|
|
|
% Initialisation des matrices pour stocker les SoC estimés par EKF et Coulomb Counting
|
|
SoC_EKF_estimates = zeros(n_steps, 1);
|
|
P_EKF_estimates = zeros(size(P, 1), size(P, 2), n_steps);
|
|
|
|
for k = 1:n_steps
|
|
% Prédiction avec le modèle d'état
|
|
[x_pred, F] = f_battery_ekf(x, current(k), params, dt);
|
|
P_pred = F * P * F' + Q;
|
|
|
|
% Mise à jour avec le modèle de mesure
|
|
[z_pred, H] = h_battery_ekf(x_pred, current(k), params);
|
|
y = voltage_measured(k) - z_pred; % Innovation (erreur de mesure)
|
|
S = H * P_pred * H' + R; % Covariance de l'innovation
|
|
K = P_pred * H' / S; % Gain de Kalman
|
|
|
|
x = x_pred + K * y; % Mise à jour de l'état
|
|
P = P_pred - K * S *K'; % Mise à jour de la covariance
|
|
|
|
% Contrainte pour s'assurer que le SoC reste entre 0 et 1
|
|
%x(1) = max(0, min(1, x(1)));
|
|
|
|
% Stockage des résultats
|
|
SoC_EKF_estimates(k) = x(1);
|
|
P_EKF_estimates(:, :, k) = P;
|
|
end
|
|
|
|
% Conversion en pourcentage des estimations de l'EKF (entre 0 et 100 %)
|
|
SoC_EKF_estimates = SoC_EKF_estimates * 100; % Conversion en pourcentage
|
|
|
|
% Méthode Coulomb Counting
|
|
SoC_Coulomb = Coulomb_Counting(current, params.C_b, soc_init*100);
|
|
|
|
% Affichage des résultats complets avec SoC Coulomb après calcul
|
|
for k = 1:n_steps
|
|
fprintf('Étape %d: Courant = %.4f A, SoC estimé EKF = %.4f%%, SoC estimé Coulomb = %.4f%%\n', ...
|
|
k, current(k), SoC_EKF_estimates(k), SoC_Coulomb(k));
|
|
end
|
|
|
|
% Calcul du RMSE entre l'EKF et Coulomb Counting
|
|
RMSE = sqrt(mean((SoC_EKF_estimates - SoC_Coulomb').^2));
|
|
fprintf('RMSE entre EKF et Coulomb Counting: %.6f\n', RMSE);
|
|
|
|
% Calcul du RMSE entre le SoC estimé par EKF et le SoC réel
|
|
RMSE_EKF_SOC = sqrt(mean((SoC_EKF_estimates - soc_real).^2));
|
|
fprintf('RMSE entre SoC estimé par EKF et SoC réel calculé par la méthode de notre client: %.6f\n', RMSE_EKF_SOC);
|
|
|
|
|
|
|
|
% Tracé des résultats
|
|
figure;
|
|
hold on;
|
|
plot(SoC_EKF_estimates, 'b', 'LineWidth', 1.5);
|
|
plot(SoC_Coulomb, 'r--', 'LineWidth', 1.5);
|
|
ylim([0 100]); % Limite les valeurs entre 0 et 100 %
|
|
xlabel('Time Step');
|
|
ylabel('SoC Estimate (%)');
|
|
title('SoC Estimation: EKF vs Coulomb Counting');
|
|
legend('EKF Estimation', 'Coulomb Counting');
|
|
grid on;
|
|
hold off;
|
|
|
|
figure;
|
|
hold on;
|
|
plot(SoC_EKF_estimates, 'b', 'LineWidth', 1.5);
|
|
plot(soc_real, 'g--', 'LineWidth', 1.5);
|
|
ylim([0 100]); % Limite les valeurs entre 0 et 100 %
|
|
xlabel('Time Step');
|
|
ylabel('SoC Estimate (%)');
|
|
title('SoC Estimation: EKF vs Real SoC');
|
|
legend('EKF Estimation', 'Real SoC');
|
|
grid on;
|
|
hold off;
|
|
|
|
%%%%%%%%%%%%%%%%%%% Fonctions pour EKF %%%%%%%%%%%%%%%%%%%
|
|
% Modèle de transition d'état
|
|
function [x_next, F] = f_battery_ekf(x, I, params, dt)
|
|
SoC = x(1);
|
|
V_RC1 = x(2);
|
|
V_RC2 = x(3);
|
|
|
|
SoC_next = SoC + (I / params.C_b) * dt;
|
|
SoC_next = max(0, min(1, SoC_next));
|
|
|
|
exp_R1C1 = exp(-dt / (params.R1 * params.C1));
|
|
exp_R2C2 = exp(-dt / (params.R2 * params.C2));
|
|
|
|
V_RC1_next = exp_R1C1 * V_RC1 + params.R1 * (1 - exp_R1C1) * I;
|
|
V_RC2_next = exp_R2C2 * V_RC2 + params.R2 * (1 - exp_R2C2) * I;
|
|
|
|
x_next = [SoC_next; V_RC1_next; V_RC2_next];
|
|
|
|
% Calcul de la jacobienne F
|
|
dV_RC1_dV_RC1 = exp_R1C1;
|
|
dV_RC2_dV_RC2 = exp_R2C2;
|
|
F = [1, 0, 0; 0, dV_RC1_dV_RC1, 0; 0, 0, dV_RC2_dV_RC2];
|
|
end
|
|
|
|
%%%%%%%%%%%%%%%%%%% Coulomb Counting Method %%%%%%%%%%%%%%%%%%%
|
|
function SOC = Coulomb_Counting(courant_decharge, C_N, SoC0)
|
|
% Number of points
|
|
N_points = length(courant_decharge);
|
|
delta_t = 0.2; % Time step in
|
|
% Initialize the SOC array
|
|
SOC = zeros(1, N_points); % SoC at each time step
|
|
SOC(1) = SoC0 ; % Initial SoC in percentage
|
|
|
|
% Calculate the state of charge (SOC) for each point with the Coulomb counting method
|
|
for i = 2:N_points
|
|
% Calculate the consumed or added charge in Ah (current * time in seconds)
|
|
charge_consomme = courant_decharge(i) * delta_t; % Current in A, delta_t in seconds
|
|
|
|
% Update the SOC according to the Coulomb counting method
|
|
SOC(i) = SOC(i-1) + ((charge_consomme / C_N) * 100); % Discharge is negative
|
|
if SOC(i) > 100
|
|
SOC(i) = 100;
|
|
elseif SOC(i) < 0
|
|
SOC(i) = 0;
|
|
end
|
|
|
|
|
|
end
|
|
end
|
|
|
|
|
|
|
|
% Modèle de mesure
|
|
function [z, H] = h_battery_ekf(x, I, params)
|
|
SoC = x(1);
|
|
V_RC1 = x(2);
|
|
V_RC2 = x(3);
|
|
|
|
OCV = OCV_from_SOC(SoC);
|
|
z = OCV - V_RC1 - V_RC2 - I * params.R0;
|
|
|
|
% Calcul de la jacobienne H
|
|
dOCV_dSoC = dOCV_from_SOC(SoC);
|
|
H = [dOCV_dSoC, -1, -1];
|
|
end
|
|
|
|
% Fonction pour calculer l'OCV en fonction du SoC
|
|
function OCV = OCV_from_SOC(SoC)
|
|
% Courbe de SoC vs V_OC (données fournies)
|
|
soc_values = [0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100] / 100; % SoC en fraction
|
|
ocv_values = [2.7668, 3.2943, 3.4666, 3.5866, 3.6586, 3.7473, 3.854, 3.9416, 4.051, 4.0824, 4.1257]; % Tensions en Volts
|
|
|
|
% Interpolation pour calculer l'OCV en fonction du SoC
|
|
OCV = interp1(soc_values, ocv_values, SoC, 'linear', 'extrap');
|
|
end
|
|
|
|
|
|
% Dérivée de l'OCV par rapport au SoC
|
|
function dOCV = dOCV_from_SOC(SoC)
|
|
soc_values = [0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100] / 100;
|
|
ocv_values = [2.7668, 3.2943, 3.4666, 3.5866, 3.6586, 3.7473, 3.854, 3.9416, 4.051, 4.0824, 4.1257];
|
|
p = polyfit(soc_values, ocv_values, 2);
|
|
dOCV = polyval(polyder(p), SoC);
|
|
end
|