diff --git a/Battery/EquivalentCircuitModel/calibration/FittingEIS.m b/Battery/EquivalentCircuitModel/calibration/FittingEIS.m new file mode 100644 index 000000000..c1812db7a --- /dev/null +++ b/Battery/EquivalentCircuitModel/calibration/FittingEIS.m @@ -0,0 +1,587 @@ +classdef FittingEIS + + properties + params0 + Z_re_exp + Z_im_exp + omega + scales + + % best_params_found = [] + % fitting_error = [] + + end + + + + methods + + + function feis = FittingEIS(params0, scales, Z_re_exp, Z_im_exp, omega) + + % feis.params0 = params0; + % feis.scales = scales; + % feis.Z_re_exp = Z_re_exp; + % feis.Z_im_exp = Z_im_exp; + % feis.omega = omega; + + + if istable(Z_re_exp), Z_re_exp = table2array(Z_re_exp); end + if istable(Z_im_exp), Z_im_exp = table2array(Z_im_exp); end + if istable(omega), omega = table2array(omega); end + if istable(params0), params0 = table2array(params0); end + if istable(scales), scales = table2array(scales); end + + % Everything has to be "double" and colon (:) + feis.params0 = double(params0(:)); + feis.scales = double(scales(:)); + feis.Z_re_exp= double(Z_re_exp(:)); + feis.Z_im_exp= double(Z_im_exp(:)); + feis.omega = double(omega(:)); + + %pb here because the way I take exp data changes a lot the + %result... + + end + + function p_norm = unscaled2scaled(feis, p) + pmin = feis.scales(1:5); + pmax = feis.scales(6:10); + p_norm = (p-pmin)./(pmax-pmin); + end + + function p = scaled2unscaled(feis, p_norm) + pmin = feis.scales(1:5); + pmax = feis.scales(6:10); + p = (pmax-pmin).*p_norm +pmin; + end + + %% Thevenin Model + + function [min_value, history, best_params, fitting_error] = optimizationBFGS(feis) + + + f_opt = @(p_norm) feis.optifunc(scaled2unscaled(feis, p_norm)); + + params0_norm = unscaled2scaled(feis, feis.params0); + % best_params_norm = lsqnonlin(deltagap, params0_norm, lb_norm, ub_norm, [0,0,-1,0,10], 0); + params0_norm = params0_norm(:); + + pmin = feis.scales(1:5); + pmax = feis.scales(6:10); + p03min = pmin(3); + p03max = pmax(3); + p05min = pmin(5); + p05max = pmax(5); + + A_custom = [0, 0, -1, 0, 2*(p05max - p05min)/(p03max- p03min)]; + b_custom = (p03min-2*p05min)/(p03max - p03min); + + max_iter = 300; + tol_obj = 1e-15; + tol_grad = 1e-5; + + [min_value, best_params_norm, history] = unitBoxBFGS(... + params0_norm, ... + f_opt, ... + 'maximize', false, ... + 'linIneq', struct('A', A_custom, 'b', b_custom), ... %A*u<=b + 'enforceFeasible', true, ... + 'maxIt', max_iter, ... + 'objChangeTol', tol_obj, ... + 'gradTol', tol_grad, ... + 'lineSearchMaxIt', 100 ... + ); + best_params = scaled2unscaled(feis, best_params_norm(:) ); + + % explications of why it stopped + fitting_error = feis.optifunc(best_params); + + it_count = length(history.val) - 1; + final_pg = history.pg(end); + + % Cost function variation + if it_count >= 1 + delta_v = abs(history.val(end) - history.val(end-1)); + else + delta_v = inf; + end + + if it_count >= max_iter + fprintf('Stopped because %d iterations reached.\n', max_iter); + elseif final_pg < tol_grad + fprintf('Stopped because reached gradTol = %e.\n', tol_grad); + elseif delta_v < tol_obj + fprintf('Stopped because reached objChangeTol = %e.\n', tol_obj); + else + fprintf('Unexpected error or relative tolerance reached.\n'); + end + + + end + + + + function [v, g_norm] = optifunc(feis, p) + + % Preventing physical values to get too close from 0 + p_safe = max(p, 1e-8); + + [Z_re_param, Z_im_param] = load_nyquist(p_safe, feis.omega); + + + if istable(Z_re_param), Z_re_param = table2array(Z_re_param); end + if istable(Z_im_param), Z_im_param = table2array(Z_im_param); end + + + Z_re_param = double(Z_re_param(:)); + Z_im_param = double(Z_im_param(:)); + + + if ~isnumeric(Z_re_param) || ~isreal(Z_re_param) || any(isnan(Z_re_param(:))) + v = 1e10; + g_norm = zeros(5,1); + return; + end + + Modulus_Z = sqrt(feis.Z_re_exp.^2 + feis.Z_im_exp.^2); + Modulus_Z = Modulus_Z(:); + + err_re = (feis.Z_re_exp(:) - Z_re_param(:)) ./ Modulus_Z(:); + err_im = (feis.Z_im_exp(:) - Z_im_param(:)) ./ Modulus_Z(:); + + v = sum(err_re.^2 + err_im.^2); + + + if nargout > 1 + + w = feis.omega(:); + R0 = p_safe(1); + R1 = p_safe(2); + C1 = p_safe(3); + R2 = p_safe(4); + C2 = p_safe(5); + + g_re_dR0 = -ones(size(w)); + g_re_dR1 = -(1-(R1*C1.*w).^2) ./ (1+(R1*C1.*w).^2).^2; + g_re_dC1 = (2*C1*R1^3.*w.^2) ./ (1+(R1*C1.*w).^2).^2; + g_re_dR2 = -(1-(R2*C2.*w).^2) ./ (1+(R2*C2.*w).^2).^2; + g_re_dC2 = (2*C2*R2^3.*w.^2) ./ (1+(R2*C2.*w).^2).^2; + + J_re = [g_re_dR0, g_re_dR1, g_re_dC1, g_re_dR2, g_re_dC2]; + + g_im_dR0 = zeros(size(w)); + g_im_dR1 = (2*R1*C1.*w) ./ (1+(R1*C1.*w).^2).^2; + g_im_dC1 = (R1^2.*w - C1^2*R1^4.*w.^3) ./ (1+(R1*C1.*w).^2).^2; + g_im_dR2 = (2*R2*C2.*w) ./ (1+(R2*C2.*w).^2).^2; + g_im_dC2 = (R2^2.*w - C2^2*R2^4.*w.^3) ./ (1+(R2*C2.*w).^2).^2; + + J_im = [g_im_dR0, g_im_dR1, g_im_dC1, g_im_dR2, g_im_dC2]; + derr_re_dp = J_re ./ Modulus_Z; + derr_im_dp = J_im ./ Modulus_Z; + + g_true = 2 * (derr_re_dp' * err_re) + 2 * (derr_im_dp' * err_im); + % Chain rule : + pmin = feis.scales(1:5); + pmax = feis.scales(6:10); + dp_dpnorm = (pmax - pmin); % Dérivée de p par rapport à p_norm + g_norm = g_true .* dp_dpnorm; + + + end + end + + % finite difference method + % if nargout > 1 + % g_norm = zeros(5,1); + % dp_norm = 1e-6; + % + % for i = 1:5 + % p_perturb = p; + % p_perturb(i) = p_perturb(i) + dp_norm * feis.scales(i); + % + % [Z_re_pert, Z_im_pert] = load_nyquist(p_perturb, feis.omega); + % + % err_re_pert = (feis.Z_re_exp(:) - double(Z_re_pert(:))) ./ Modulus_Z; + % err_im_pert = (feis.Z_im_exp(:) - double(Z_im_pert(:))) ./ Modulus_Z; + % + % v_perturb = sum(err_re_pert.^2 + err_im_pert.^2); + % + % g_norm = (v_perturb - v) / dp_norm; + % end + % end + + +%% Warburg Model + + function [min_value, history, best_params, fitting_error] = optimizationBFGS_warburg(feis) + + f_opt = @(p_norm) feis.optifunc_warburg(p_norm.* feis.scales) ; + + params0_norm = feis.params0 ./ feis.scales; + params0_norm = params0_norm(:); + + + + [min_value, best_params_norm, history] = unitBoxBFGS(... + params0_norm, ... + f_opt, ... + 'maximize', false, ... + 'enforceFeasible', true, ... + 'maxIt', 300, ... + 'objChangeTol', 1e-6, ... + 'gradTol', 1e-5, ... + 'lineSearchMaxIt', 200 ... + ); + best_params = best_params_norm(:) .* feis.scales(:); + + fitting_error = feis.optifunc_warburg(best_params); + + end + + function [v, g_norm] = optifunc_warburg(feis, p) + + % Preventing physical values to get too close from 0 + p_safe = max(p, 1e-7); + + + w = feis.omega(:); + R0 = p_safe(1); + R1 = p_safe(2); + Q1 = p_safe(3); + a1 = p_safe(4); + R2 = p_safe(5); + Q2 = p_safe(6); + a2 = p_safe(7); + Q = p_safe(8); + L = p_safe(9); + + s = 1i .* w; + s = s(:); + + sa1 = s.^a1; + sa2 = s.^a2; + + Z_param_cplx = R0 + s.*L + R1 ./ (1 + R1.*Q1.*sa1) + R2 ./ (1 + R2.*Q2.*sa2) + 1 ./ (Q.*(s.^0.5)); + + Z_re_param = real(Z_param_cplx); + Z_im_param = imag(Z_param_cplx); + + if istable(Z_re_param), Z_re_param = table2array(Z_re_param); end + if istable(Z_im_param), Z_im_param = table2array(Z_im_param); end + + + Z_re_param = double(Z_re_param(:)); + Z_im_param = double(Z_im_param(:)); + + Modulus_Z = sqrt(Z_re_param.^2 + Z_im_param.^2); + Modulus_Z = Modulus_Z(:); + + err_re = (feis.Z_re_exp(:) - Z_re_param(:)) ./ Modulus_Z(:); + err_im = (feis.Z_im_exp(:) - Z_im_param(:)) ./ Modulus_Z(:); + + v = sum(err_re.^2 + err_im.^2); + + + if nargout > 1 + + den1 = (1 + R1 .* Q1 .* sa1).^2; + den2 = (1 + R2 .* Q2 .* sa2).^2; + + + dZ_dR0 = ones(size(s)); % R0 + dZ_dR1 = 1 ./ den1; % R1 + dZ_dR2 = 1 ./ den2; % R2 + + dZ_dQ1 = - (R1^2 .* sa1) ./ den1; % Q1 + dZ_dQ2 = - (R2^2 .* sa2) ./ den2; % Q2 + + dZ_da1 = - (Q1 * R1^2 .* sa1 .* log(s)) ./ den1; % a1 + dZ_da2 = - (Q2 * R2^2 .* sa2 .* log(s)) ./ den2; % a2 + + dZ_dQ = - 1 ./ (Q^2 .* (s.^0.5)); % Q + dZ_dL = s; % L + + J = [dZ_dR0, dZ_dR1, dZ_dQ1, dZ_da1, dZ_dR2, dZ_dQ2, dZ_da2, dZ_dQ, dZ_dL]; + + J_im = imag(J); + J_re = real(J); + + derr_re_dp = -J_re ./ Modulus_Z; + derr_im_dp = -J_im ./ Modulus_Z; + + g_true = 2 * (derr_re_dp' * err_re) + 2 * (derr_im_dp' * err_im); + g_norm = g_true .* feis.scales(:); + + + if any(isinf(g_norm)) || any(isnan(g_norm)) || isinf(v) || isnan(v) v = 1e10; + g_norm = zeros(9,1); + return; + end + end + end + + + + %% lsq method + + function residuals = optifunc_lsq(feis, p) + % 1. On empêche physiquement les paramètres de valoir zéro + p = max(p, 1e-8); + + % 2. Chargement du modèle + [Z_re_param, Z_im_param] = load_nyquist(p, feis.omega); + if istable(Z_re_param), Z_re_param = table2array(Z_re_param); end + if istable(Z_im_param), Z_im_param = table2array(Z_im_param); end + + Z_re_param = double(Z_re_param(:)); + Z_im_param = double(Z_im_param(:)); + + % Si le modèle renvoie du NaN (physiquement impossible), on pénalise + if ~isnumeric(Z_re_param) || ~isreal(Z_re_param) || any(isnan(Z_re_param(:))) + residuals = 1e10 * ones(2 * length(feis.omega), 1); + return; + end + + % 3. Calcul du module pour l'erreur relative + Modulus_Z = sqrt(feis.Z_re_exp.^2 + feis.Z_im_exp.^2); + Modulus_Z = Modulus_Z(:); + + % 4. Vecteurs de résidus (Non mis au carré, lsqnonlin s'en charge) + err_re = (feis.Z_re_exp(:) - Z_re_param(:)) ./ Modulus_Z; + err_im = (feis.Z_im_exp(:) - Z_im_param(:)) ./ Modulus_Z; + + % 5. On empile les parties réelles et imaginaires en un seul vecteur colonne + residuals = [err_re; err_im]; + end + + function [min_value, best_params, resnorm] = optimizationLsqnonlin(feis) + % Fonction objectif : renvoie le vecteur de résidus + f_res = @(p_norm) feis.optifunc_lsq(p_norm(:) .* feis.scales(:)); + + % Paramètres initiaux normés + params0_norm = feis.params0(:) ./ feis.scales(:); + + % Bornes (Limites inférieures strictes > 0 pour éviter les crashs) + lb_norm = ones(5,1) * 1e-5; + ub_norm = ones(5,1) * 1e5./ feis.scales(:); + + % Options de lsqnonlin + options = optimoptions('lsqnonlin', ... + 'Display', 'iter-detailed', ... % Affichage très détaillé + 'Algorithm', 'trust-region-reflective', ... % Parfait pour l'EIS + 'MaxFunctionEvaluations', 3000, ... + 'MaxIterations', 500, ... + 'StepTolerance', 1e-8, ... + 'FunctionTolerance', 1e-6); + + % Lancement de l'optimiseur natif MATLAB + [best_params_norm, resnorm, ~, exitflag, output] = lsqnonlin(... + f_res, params0_norm, lb_norm, ub_norm, options); + + % Récupération des vraies valeurs (dénormalisation) + best_params = best_params_norm(:) .* feis.scales(:); + min_value = resnorm; % Erreur finale (somme des carrés des résidus) + + % Affichage du statut de fin + disp('=== RAISON DE FIN DE LSQNONLIN ==='); + disp(output.message); + end + +%% lsq method warburg + function v = optifunclsq_warburg(feis, p) + p_safe = max(p, 1e-15); + + R0 = p_safe(1); R1 = p_safe(2); Q1 = p_safe(3); a1 = p_safe(4); + R2 = p_safe(5); Q2 = p_safe(6); a2 = p_safe(7); Q = p_safe(8); L = p_safe(9); + + w = feis.omega(:); + s = 1i .* w; + + sa1 = s.^a1; + sa2 = s.^a2; + Z_param_cplx = R0 + s.*L + R1 ./ (1 + R1.*Q1.*sa1) + R2 ./ (1 + R2.*Q2.*sa2) + 1 ./ (Q.*(s.^0.5)); + + Z_re_param = real(Z_param_cplx); + Z_im_param = imag(Z_param_cplx); + + if any(isnan(Z_re_param)) || any(isnan(Z_im_param)) || any(isinf(Z_re_param)) || any(isinf(Z_im_param)) + v = 1e6; % Pénalité forte mais gérable + return; + end + + Modulus_Z = sqrt(Z_re_param.^2 + Z_im_param.^2); + Modulus_Z = Modulus_Z(:); + + err_re = (feis.Z_re_exp(:) - Z_re_param) ./ Modulus_Z; + err_im = (feis.Z_im_exp(:) - Z_im_param) ./ Modulus_Z; + + v = sum(err_re.^2 + err_im.^2); + + end + + + function [min_value, history, best_params, fitting_error] = optimizationlsq_warburg(feis) + f_opt = @(p_norm) feis.optifunclsq_warburg(p_norm .* feis.scales); + + params0_norm = feis.params0 ./ feis.scales; + params0_norm = params0_norm(:); + + if length(params0_norm) ~= 9 + error('feis.params0 should have 9 parameters but it has %d).', length(params0_norm)); + end + + lb_norm = [1e-10; 1e-10; 1e-10; 1e-3; 1e-10; 1e-10; 0.01; 1e-5; 1e-15] ./ feis.scales; + ub_norm = [1e5; 1e5; 1e6; 1.0; 1e5; 1e6; 1.0; 1e15; 1e-2 ] ./ feis.scales; + + options = optimoptions('fmincon', ... + 'Display', 'iter', ... + 'Algorithm', 'sqp', ... % Très bon pour ce type de fitting + 'SpecifyObjectiveGradient', false, ... + 'MaxFunctionEvaluations', 10000, ... + 'MaxIterations', 1000); + + [best_params_norm, min_value] = fmincon(f_opt, params0_norm, [], [], [], [], lb_norm, ub_norm, [], options); + + history = []; % Si vous n'en avez pas besoin + best_params = best_params_norm(:) .* feis.scales(:); + fitting_error = feis.optifunclsq_warburg(best_params); + end + + %% showing the results + + function plotresults_thevenin(feis, best_params, fitting_error) + + [Z_re_fit, Z_im_fit] = load_nyquist(best_params, feis.omega); + + figure; + subplot(3,1,1); + semilogx(feis.omega, feis.Z_re_exp, 'r', 'MarkerFaceColor', 'r'); + hold on; + semilogx(feis.omega, Z_re_fit, 'b'); + legend('experience', 'fitted model'); + title('Fitting Real Impedance'); + xlabel('Omega'); + ylabel('Z_{re} '); + + subplot(3,1,2); + semilogx(feis.omega, -feis.Z_im_exp, 'r', 'MarkerFaceColor', 'r'); + hold on; + semilogx(feis.omega, -Z_im_fit, 'b'); + legend('experience', 'fitted model'); + title('Fitting Imaginary Impedance'); + xlabel('Omega'); + ylabel('-Z_{im} '); + + subplot(3,1,3); + plot(feis.Z_re_exp, -feis.Z_im_exp, 'r', 'MarkerFaceColor', 'r'); + hold on; + plot(feis.Z_re_exp, -Z_im_fit, 'b'); + legend('experience', 'fitted model'); + title('Nyquist Diagram'); + xlabel('Z_{re}'); + ylabel('-Z_{im} '); + axis equal; + + grid on; + + text_error = sprintf('Fitting error : %.2e', fitting_error); + + subplot(3,1,3); + + text(0.05, 0.90, text_error, 'Units', 'normalized', ... + 'BackgroundColor', 'white', ... + 'EdgeColor', 'black', ... + 'FontSize', 11, ... + 'FontWeight', 'bold'); + end + + + function plotresults_warburg(feis, best_params, fitting_error) + + % [Z_re_fit, Z_im_fit] = load_nyquist(best_params, feis.omega); + w = feis.omega(:); + R0 = best_params(1); + R1 = best_params(2); + Q1 = best_params(3); + a1 = best_params(4); + R2 = best_params(5); + Q2 = best_params(6); + a2 = best_params(7); + Q = best_params(8); + L = best_params(9); + s = 1i .* w; + s = s(:); + + sa1 = s.^a1; + sa2 = s.^a2; + + Z_param_cplx = R0 + s.*L + R1 ./ (1 + R1.*Q1.*sa1) + R2 ./ (1 + R2.*Q2.*sa2) + 1 ./ (Q.*(s.^0.5)); + + Z_re_fit = real(Z_param_cplx); + Z_im_fit = imag(Z_param_cplx); + + + figure; + subplot(3,1,1); + semilogx(feis.omega, feis.Z_re_exp, 'r', 'MarkerFaceColor', 'r'); + hold on; + semilogx(feis.omega, Z_re_fit, 'b'); + legend('experience', 'fitted model'); + title('Fitting results'); + xlabel('Omega'); + ylabel('Z_{re} '); + + subplot(3,1,2); + semilogx(feis.omega, -feis.Z_im_exp, 'r', 'MarkerFaceColor', 'r'); + hold on; + semilogx(feis.omega, -Z_im_fit, 'b'); + legend('experience', 'fitted model'); + title('Fitting results'); + xlabel('Omega'); + ylabel('-Z_{im} '); + + subplot(3,1,3); + plot(feis.Z_re_exp, -feis.Z_im_exp, 'r', 'MarkerFaceColor', 'r'); + hold on; + plot(Z_re_fit, -Z_im_fit, 'b'); + legend('experience', 'fitted model'); + title('Nyquist'); + xlabel('Z_{re}'); + ylabel('-Z_{im} '); + axis equal; + + grid on; + text_error = sprintf('Fitting error : %.2e', fitting_error); + text(0.05, 0.85, text_error, 'Units', 'normalized', ... + 'BackgroundColor', 'white', ... + 'EdgeColor', 'black', ... + 'FontSize', 11, ... + 'FontWeight', 'bold'); + end + + function printResults(feis, best_params, fitting_error) + + fprintf('\n=== FITTING SCORE ===\n'); + fprintf('Error : %e\n', fitting_error); + + fprintf('\n=== PARAMETERS FOUND ===\n'); + fprintf('R0 = %.2e Ohms\n', best_params(1)); + fprintf('R1 = %.2e Ohms\n', best_params(2)); + fprintf('C1 = %.2e Farads\n', best_params(3)); + fprintf('R2 = %.2e Ohms\n', best_params(4)); + fprintf('C2 = %.2e Farads\n', best_params(5)); + + % fprintf('\n=== PARAMETERS FOUND ===\n'); + % fprintf('R0 = %.4e Ohms\n', best_params(1)); + % fprintf('R1 = %.4e Ohms\n', best_params(2)); + % fprintf('Q1 = %.4e s^a/Ohm\n', best_params(3)); + % fprintf('a1 = %.4f \n', best_params(4)); % %f suffit pour 'a' car il est entre 0 et 1 + % fprintf('R2 = %.4e Ohms\n', best_params(5)); + % fprintf('Q2 = %.4e s^a/Ohm\n', best_params(6)); + % fprintf('a2 = %.4f \n', best_params(7)); + % fprintf('Q = %.4e \n', best_params(8)); + % fprintf('L = %.4e Henrys\n', best_params(9)); + end + end +end diff --git a/Battery/EquivalentCircuitModel/calibration/plotting/plot_nyquist_ank.m b/Battery/EquivalentCircuitModel/calibration/plotting/plot_nyquist_ank.m new file mode 100644 index 000000000..2044245c4 --- /dev/null +++ b/Battery/EquivalentCircuitModel/calibration/plotting/plot_nyquist_ank.m @@ -0,0 +1,24 @@ + + + +load_experimental_data(); + +% --- 5. Tracé du diagramme de Nyquist --- +figure('Name', 'EIS', 'Color', 'w'); + +% On trace -Im(Z) en fonction de Re(Z) +plot(Z_real, Z_imag, 'o', ... + 'LineWidth', 1.5, ... + 'MarkerSize', 2, ... + 'MarkerFaceColor', [0 0.4470 0.7410], ... + 'MarkerEdgeColor', 'k'); +grid on; +axis equal; % Force la même échelle en X et Y + +% Ajout des labels +xlabel('Z_{re} ', 'FontSize', 12, 'FontWeight', 'bold'); +ylabel('-Z_{im}', 'FontSize', 12, 'FontWeight', 'bold'); +title('Nyquist diagram from data', 'FontSize', 14); + +% Amélioration des axes +set(gca, 'FontSize', 11, 'LineWidth', 1); \ No newline at end of file diff --git a/Battery/EquivalentCircuitModel/calibration/plotting/plot_nyquist_params.m b/Battery/EquivalentCircuitModel/calibration/plotting/plot_nyquist_params.m new file mode 100644 index 000000000..2afa33f60 --- /dev/null +++ b/Battery/EquivalentCircuitModel/calibration/plotting/plot_nyquist_params.m @@ -0,0 +1,26 @@ + +% params = [6.2261e-3, 0.3e-3, 10000.1264, 0.00097, 2e5]; % best params +% handmade + +params = [0.003, 0.00134, 5, 0.00330, 20000]; +[Z_re_exp, Z_im_exp, omega] = load_experimental_data(); +omega = logspace(-2, 4, 100); + +[Z_real, Z_imag] = load_nyquist(params, omega) + +figure; +plot(Z_re_exp, Z_im_exp, 'ro', 'MarkerFaceColor', 'r'); +hold on; +plot(Z_real, Z_imag, '-o', 'LineWidth', 1.5, 'MarkerFaceColor', 'b'); +axis equal; % Essentiel pour Nyquist +grid on; +xlabel('Z_{réel} (\Omega)', 'FontWeight', 'bold'); +ylabel('-Z_{imaginaire} (\Omega)', 'FontWeight', 'bold'); +title('Nyquist diagram simulated'); + + + + + + + diff --git a/Battery/EquivalentCircuitModel/calibration/pxdParamsImpact/activationEnergyOfReaction_impact.m b/Battery/EquivalentCircuitModel/calibration/pxdParamsImpact/activationEnergyOfReaction_impact.m new file mode 100644 index 000000000..6ff436588 --- /dev/null +++ b/Battery/EquivalentCircuitModel/calibration/pxdParamsImpact/activationEnergyOfReaction_impact.m @@ -0,0 +1,48 @@ + +mrstModule add ad-core mrst-gui mpfa agmg linearsolvers + +jsonstruct_material = parseBattmoJson(fullfile('ParameterData','ParameterSets','Chen2020','chen2020_lithium_ion_battery.json')); +jsonstruct_geometry = parseBattmoJson(fullfile('Examples', 'JsonDataFiles', 'geometryChen.json')); + +jsonstruct = mergeJsonStructs({jsonstruct_material, ... + jsonstruct_geometry}); + + + +includeDoubleLayer = false; + +if includeDoubleLayer + + jsonstruct.(ne).(co).(am).(itf).useDoubleLayerCapacity = true; + jsonstruct.(ne).(co).(am).(itf).doubleLayerCapacitance = 0.2; + +end + + +[model, inputparams, ~] = setupModelFromJson(jsonstruct); + +c_ne = 29.866*mol/litre; % initial concentration at negative electrode +c_pe = 17.038*mol/litre; % initial concentration at positive electrode + +initstate = initStateChen2020(model, c_ne, c_pe); + +b = linspace(1e4, 5e4, 20); %change here + +frequences = logspace(-2, 4, 30); +figure; +hold on; +for i = 1:length(b) + + inputparams.NegativeElectrode.Coating.ActiveMaterial.Interface.activationEnergyOfReaction = b(i); %change here + inputparams.PositiveElectrode.Coating.ActiveMaterial.Interface.activationEnergyOfReaction = b(i)/2; %change here + + impsolv = ImpedanceSolver(inputparams, 'initstate', initstate, 'computeSteadyState', false); + + Z = impsolv.computeImpedance(frequences); + Z_re = real(Z); + Z_im = imag(Z); + curve = sprintf('activation Energy Of Reaction = %.2e', b(i)); %change here + plot(Z_re, -Z_im, 'DisplayName', curve); +end + +legend('show'); \ No newline at end of file diff --git a/Battery/EquivalentCircuitModel/calibration/pxdParamsImpact/bruggeman_impact.m b/Battery/EquivalentCircuitModel/calibration/pxdParamsImpact/bruggeman_impact.m new file mode 100644 index 000000000..1412d6644 --- /dev/null +++ b/Battery/EquivalentCircuitModel/calibration/pxdParamsImpact/bruggeman_impact.m @@ -0,0 +1,46 @@ + +mrstModule add ad-core mrst-gui mpfa agmg linearsolvers + +jsonstruct_material = parseBattmoJson(fullfile('ParameterData','ParameterSets','Chen2020','chen2020_lithium_ion_battery.json')); +jsonstruct_geometry = parseBattmoJson(fullfile('Examples', 'JsonDataFiles', 'geometryChen.json')); + +jsonstruct = mergeJsonStructs({jsonstruct_material, ... + jsonstruct_geometry}); + + + +includeDoubleLayer = false; + +if includeDoubleLayer + + jsonstruct.(ne).(co).(am).(itf).useDoubleLayerCapacity = true; + jsonstruct.(ne).(co).(am).(itf).doubleLayerCapacitance = 0.2; + +end + + +[model, inputparams, ~] = setupModelFromJson(jsonstruct); + +c_ne = 29.866*mol/litre; % initial concentration at negative electrode +c_pe = 17.038*mol/litre; % initial concentration at positive electrode + +initstate = initStateChen2020(model, c_ne, c_pe); +b = linspace(0.5,2,10); + +frequences = logspace(-2, 4, 30); +figure; +hold on; +for i = 1:length(b) + + inputparams.NegativeElectrode.Coating.bruggemanCoefficient = b(i); + inputparams.PositiveElectrode.Coating.bruggemanCoefficient = b(i); + impsolv = ImpedanceSolver(inputparams, 'initstate', initstate, 'computeSteadyState', false); + + Z = impsolv.computeImpedance(frequences); + Z_re = real(Z); + Z_im = imag(Z); + curve = sprintf('Bruggeman coefficient = %.2f', b(i)); + plot(Z_re, -Z_im, 'DisplayName', curve); +end + +legend('show'); \ No newline at end of file diff --git a/Battery/EquivalentCircuitModel/calibration/pxdParamsImpact/doubleLayerCapacitance_impact.m b/Battery/EquivalentCircuitModel/calibration/pxdParamsImpact/doubleLayerCapacitance_impact.m new file mode 100644 index 000000000..df5714137 --- /dev/null +++ b/Battery/EquivalentCircuitModel/calibration/pxdParamsImpact/doubleLayerCapacitance_impact.m @@ -0,0 +1,61 @@ +mrstModule add ad-core mrst-gui mpfa agmg linearsolvers + +jsonstruct_material = parseBattmoJson(fullfile('ParameterData','ParameterSets','Chen2020','chen2020_lithium_ion_battery.json')); +jsonstruct_geometry = parseBattmoJson(fullfile('Examples', 'JsonDataFiles', 'geometryChen.json')); + +jsonstruct = mergeJsonStructs({jsonstruct_material, ... + jsonstruct_geometry}); + +ne = 'NegativeElectrode'; +co = 'Coating'; +am = 'ActiveMaterial'; +itf = 'Interface'; + +includeDoubleLayer = true; + +if includeDoubleLayer + + jsonstruct.(ne).(co).(am).(itf).useDoubleLayerCapacity = true; + jsonstruct.(ne).(co).(am).(itf).doubleLayerCapacitance = 0.2; + +end + + +[model, inputparams, ~] = setupModelFromJson(jsonstruct); + +c_ne = 29.866*mol/litre; % initial concentration at negative electrode +c_pe = 17.038*mol/litre; % initial concentration at positive electrode + +initstate = initStateChen2020(model, c_ne, c_pe); + +b = logspace(-2, 2, 20); %change here + +frequences = logspace(-2, 4, 30); + + +figure; +hold on; +for i = 1:length(b) + + fprintf('compute impedance for doubleLayerCapacitance = %.2e ... ', b(i)); %change here + + tic + jsonstruct.(ne).(co).(am).(itf).doubleLayerCapacitance = b(i); %change here + [model, inputparams, ~] = setupModelFromJson(jsonstruct); + + % inputparams.(ne).(co).(am).(itf).doubleLayerCapacitance = b(i); + impsolv = ImpedanceSolver(inputparams, 'initstate', initstate, 'computeSteadyState', false); + + Z = impsolv.computeImpedance(frequences); + + fprintf('done in %g s\n', toc); + + Z_re = real(Z); + Z_im = imag(Z); + curve = sprintf('doubleLayerCapacitance = %.2e', b(i)); %change here + + plot(Z_re, -Z_im, 'DisplayName', curve); + +end + +legend('show'); \ No newline at end of file diff --git a/Battery/EquivalentCircuitModel/calibration/pxdParamsImpact/electronicConductivity_impact.m b/Battery/EquivalentCircuitModel/calibration/pxdParamsImpact/electronicConductivity_impact.m new file mode 100644 index 000000000..df70d8b3f --- /dev/null +++ b/Battery/EquivalentCircuitModel/calibration/pxdParamsImpact/electronicConductivity_impact.m @@ -0,0 +1,46 @@ +mrstModule add ad-core mrst-gui mpfa agmg linearsolvers + +jsonstruct_material = parseBattmoJson(fullfile('ParameterData','ParameterSets','Chen2020','chen2020_lithium_ion_battery.json')); +jsonstruct_geometry = parseBattmoJson(fullfile('Examples', 'JsonDataFiles', 'geometryChen.json')); + +jsonstruct = mergeJsonStructs({jsonstruct_material, ... + jsonstruct_geometry}); + + + +includeDoubleLayer = false; + +if includeDoubleLayer + + jsonstruct.(ne).(co).(am).(itf).useDoubleLayerCapacity = true; + jsonstruct.(ne).(co).(am).(itf).doubleLayerCapacitance = 0.2; + +end + + +[model, inputparams, ~] = setupModelFromJson(jsonstruct); + +c_ne = 29.866*mol/litre; % initial concentration at negative electrode +c_pe = 17.038*mol/litre; % initial concentration at positive electrode + +initstate = initStateChen2020(model, c_ne, c_pe); +b = linspace(10, 1000, 20); %change here + +frequences = logspace(-2, 4, 30); +figure; +hold on; +for i = 1:length(b) + + inputparams.NegativeElectrode.Coating.ActiveMaterial.electronicConductivity = b(i); %change here + inputparams.PositiveElectrode.Coating.ActiveMaterial.electronicConductivity = b(i); %change here + + impsolv = ImpedanceSolver(inputparams, 'initstate', initstate, 'computeSteadyState', false); + + Z = impsolv.computeImpedance(frequences); + Z_re = real(Z); + Z_im = imag(Z); + curve = sprintf('Electronic Conductivity = %.2f', b(i)); %change here + plot(Z_re, -Z_im, 'DisplayName', curve); +end + +legend('show'); \ No newline at end of file diff --git a/Battery/EquivalentCircuitModel/calibration/pxdParamsImpact/param_dependency.m b/Battery/EquivalentCircuitModel/calibration/pxdParamsImpact/param_dependency.m new file mode 100644 index 000000000..8a44ae62d --- /dev/null +++ b/Battery/EquivalentCircuitModel/calibration/pxdParamsImpact/param_dependency.m @@ -0,0 +1,130 @@ +function param_dependency() + mrstModule add ad-core mrst-gui mpfa agmg linearsolvers + + % 1. Chargement des données (Conservez votre méthode de chargement ici) + jsonstruct_material = parseBattmoJson(fullfile('ParameterData','ParameterSets','Chen2020','chen2020_lithium_ion_battery.json')); + jsonstruct_geometry = parseBattmoJson(fullfile('Examples', 'JsonDataFiles', 'geometryChen.json')); + jsonstruct = mergeJsonStructs({jsonstruct_material, jsonstruct_geometry}); + ne = 'NegativeElectrode'; + co = 'Coating'; + am = 'ActiveMaterial'; + itf = 'Interface'; + + includeDoubleLayer = true; + + if includeDoubleLayer + + jsonstruct.(ne).(co).(am).(itf).useDoubleLayerCapacity = true; + jsonstruct.(ne).(co).(am).(itf).doubleLayerCapacitance = 0.2; + + end + + + + c_ne = 29.866*mol/litre; + c_pe = 17.038*mol/litre; + + % 2. Création de l'interface graphique (Fenêtre plus haute pour accueillir les curseurs) + fig = uifigure('Name', 'Impedance Explorer', 'Position', [100, 50, 600, 650]); + + % Graphique + ax = uiaxes(fig, 'Position', [50, 350, 500, 250]); + title(ax, 'Nyquist Diagram'); + xlabel(ax, 'Re(Z)'); + ylabel(ax, '-Im(Z)'); + grid(ax, 'on'); + + % Label de statut + lbl = uilabel(fig, 'Position', [50, 310, 500, 22], 'Text', 'Ready, use slider to start.'); + lbl.FontWeight = 'bold'; + + % --- CURSEUR 1 : Diffusion Solide (Anode) --- + % Modifie le transport du Lithium dans les particules de graphite. + uilabel(fig, 'Position', [50, 260, 500, 40], 'WordWrap', 'on', ... + 'Text', '1. Anode refDiffusionCoefficient * 10^x (initial value: 1.3135e-15)'); + + sld1 = uislider(fig, 'Position', [50, 240, 500, 3], ... + 'Limits', [-5, 5], ... % Définit le minimum et le maximum + 'MajorTicks', -5:0.5:5, ... % Ajoute une graduation tous les 1 + 'Value', 0); + + % --- CURSEUR 2 : Taux de réaction (Anode) --- + % Modifie la cinétique de transfert de charge. + uilabel(fig, 'Position', [50, 170, 500, 40], 'WordWrap', 'on', ... + 'Text', '2. Anode reactionRateConstant * 10^x (initial value: 5.031e-11) '); + sld2 = uislider(fig, 'Position', [50, 150, 500, 3],... + 'Limits', [-5, 5], ... % Définit le minimum et le maximum + 'MajorTicks', -5:0.5:5, ... % Ajoute une graduation tous les 1 + 'Value', 0); + + % --- CURSEUR 3 : Taux de réaction (Cathode) --- + % Modifie la cinétique de transfert de charge. + uilabel(fig, 'Position', [50, 80, 500, 40], 'WordWrap', 'on', ... + 'Text', '3. doubleLayerCapacitance * x (initial value: 0.2)'); + sld3 = uislider(fig, 'Position', [50, 60, 500, 3], ... + 'Limits', [-5, 5], ... % Définit le minimum et le maximum + 'MajorTicks', -5:0.5:5, ... % Ajoute une graduation tous les 1 + 'Value', 0); + + % 4. Connexion aux événements + % On crée une fonction anonyme qui passe TOUS les curseurs + callback_fcn = @(src, event) updatePlot(ax, lbl, jsonstruct, c_ne, c_pe, sld1, sld2, sld3); + + sld1.ValueChangedFcn = callback_fcn; + sld2.ValueChangedFcn = callback_fcn; + sld3.ValueChangedFcn = callback_fcn; + + % Optionnel : Tracer l'état initial + % updatePlot(ax, lbl, jsonstruct, c_ne, c_pe, sld1, sld2, sld3); +end + +% --- Fonction de mise à jour --- +function updatePlot(ax, lbl, jsonstruct, c_ne, c_pe, sld1, sld2, sld3) + + % Figer l'affichage + lbl.Text = 'Processing... '; + lbl.FontColor = '#D95319'; + drawnow; + + % 1. Récupération des valeurs de base depuis votre JSON + base_diff_anode = 1.3135e-15; + base_rate_anode = 5.031e-11; + base_doubleLayerCapacitance = 0.2; + + % 2. Application des multiplicateurs (10^slider_val) + val_diff_anode = base_diff_anode * (10^(sld1.Value)); + val_rate_anode = base_rate_anode * (10^(sld2.Value)); + val_doubleLayerCapacitance = base_doubleLayerCapacitance * 10^(sld3.Value); + + % 3. Injection dans la structure JSON AVANT de construire le modèle + jsonstruct.NegativeElectrode.Coating.ActiveMaterial.SolidDiffusion.referenceDiffusionCoefficient = val_diff_anode; + jsonstruct.NegativeElectrode.Coating.ActiveMaterial.Interface.reactionRateConstant = val_rate_anode; + jsonstruct.NegativeElectrode.Coating.ActiveMaterial.Interface.doubleLayerCapacitance = val_doubleLayerCapacitance; + + try + % 4. Exécution de BattMo + [model, inputparams, ~] = setupModelFromJson(jsonstruct); + initstate = initStateChen2020(model, c_ne, c_pe); + + impsolv = ImpedanceSolver(inputparams, 'initstate', initstate, 'computeSteadyState', false); + + frequences = logspace(-2, 3, 30); + Z = impsolv.computeImpedance(frequences); + + % 5. Tracé + plot(ax, real(Z), -imag(Z), '-o', 'LineWidth', 2, 'Color', '#0072BD'); + + % Retour à la normale + lbl.Text = 'Done'; + lbl.FontColor = '#77AC30'; + + catch ME + lbl.Text = 'Erreur lors du calcul (vérifiez la console).'; + lbl.FontColor = '#A2142F'; + disp('--- ERREUR DANS LE CALCUL ---'); + disp(ME.message); + for k = 1:length(ME.stack) + disp(['Ligne ', num2str(ME.stack(k).line), ' dans ', ME.stack(k).name]); + end + end +end \ No newline at end of file diff --git a/Battery/EquivalentCircuitModel/calibration/pxdParamsImpact/porosity_impact.m b/Battery/EquivalentCircuitModel/calibration/pxdParamsImpact/porosity_impact.m new file mode 100644 index 000000000..47d1a7cce --- /dev/null +++ b/Battery/EquivalentCircuitModel/calibration/pxdParamsImpact/porosity_impact.m @@ -0,0 +1,45 @@ +mrstModule add ad-core mrst-gui mpfa agmg linearsolvers + +jsonstruct_material = parseBattmoJson(fullfile('ParameterData','ParameterSets','Chen2020','chen2020_lithium_ion_battery.json')); +jsonstruct_geometry = parseBattmoJson(fullfile('Examples', 'JsonDataFiles', 'geometryChen.json')); + +jsonstruct = mergeJsonStructs({jsonstruct_material, ... + jsonstruct_geometry}); + + + +includeDoubleLayer = false; + +if includeDoubleLayer + + jsonstruct.(ne).(co).(am).(itf).useDoubleLayerCapacity = true; + jsonstruct.(ne).(co).(am).(itf).doubleLayerCapacitance = 0.2; + +end + + +[model, inputparams, ~] = setupModelFromJson(jsonstruct); + +c_ne = 29.866*mol/litre; % initial concentration at negative electrode +c_pe = 17.038*mol/litre; % initial concentration at positive electrode + +initstate = initStateChen2020(model, c_ne, c_pe); +b = linspace(0.1, 0.9, 20); %change here + +frequences = logspace(-2, 4, 30); +figure; +hold on; +for i = 1:length(b) + + inputparams.Separator.porosity = b(i); %change here + + impsolv = ImpedanceSolver(inputparams, 'initstate', initstate, 'computeSteadyState', false); + + Z = impsolv.computeImpedance(frequences); + Z_re = real(Z); + Z_im = imag(Z); + curve = sprintf('porosity = %.2f', b(i)); %change here + plot(Z_re, -Z_im, 'DisplayName', curve); +end + +legend('show'); \ No newline at end of file diff --git a/Battery/EquivalentCircuitModel/calibration/pxdParamsImpact/reactionRateConstant_impact.m b/Battery/EquivalentCircuitModel/calibration/pxdParamsImpact/reactionRateConstant_impact.m new file mode 100644 index 000000000..12baa2706 --- /dev/null +++ b/Battery/EquivalentCircuitModel/calibration/pxdParamsImpact/reactionRateConstant_impact.m @@ -0,0 +1,47 @@ +mrstModule add ad-core mrst-gui mpfa agmg linearsolvers + +jsonstruct_material = parseBattmoJson(fullfile('ParameterData','ParameterSets','Chen2020','chen2020_lithium_ion_battery.json')); +jsonstruct_geometry = parseBattmoJson(fullfile('Examples', 'JsonDataFiles', 'geometryChen.json')); + +jsonstruct = mergeJsonStructs({jsonstruct_material, ... + jsonstruct_geometry}); + + + +includeDoubleLayer = false; + +if includeDoubleLayer + + jsonstruct.(ne).(co).(am).(itf).useDoubleLayerCapacity = true; + jsonstruct.(ne).(co).(am).(itf).doubleLayerCapacitance = 0.2; + +end + + +[model, inputparams, ~] = setupModelFromJson(jsonstruct); + +c_ne = 29.866*mol/litre; % initial concentration at negative electrode +c_pe = 17.038*mol/litre; % initial concentration at positive electrode + +initstate = initStateChen2020(model, c_ne, c_pe); + +b = linspace(1e-11, 1e-12, 20); %change here + +frequences = logspace(-2, 4, 30); +figure; +hold on; +for i = 1:length(b) + + inputparams.NegativeElectrode.Coating.ActiveMaterial.Interface.reactionRateConstant = b(i); %change here + inputparams.PositiveElectrode.Coating.ActiveMaterial.Interface.reactionRateConstant = b(i)*5; %change here + + impsolv = ImpedanceSolver(inputparams, 'initstate', initstate, 'computeSteadyState', false); + + Z = impsolv.computeImpedance(frequences); + Z_re = real(Z); + Z_im = imag(Z); + curve = sprintf('reactionRateConstant = %.2e', b(i)); %change here + plot(Z_re, -Z_im, 'DisplayName', curve); +end + +legend('show'); \ No newline at end of file diff --git a/Battery/EquivalentCircuitModel/calibration/pxdParamsImpact/referenceDiffusionCoefficient_impact.m b/Battery/EquivalentCircuitModel/calibration/pxdParamsImpact/referenceDiffusionCoefficient_impact.m new file mode 100644 index 000000000..88aebd8d1 --- /dev/null +++ b/Battery/EquivalentCircuitModel/calibration/pxdParamsImpact/referenceDiffusionCoefficient_impact.m @@ -0,0 +1,46 @@ +mrstModule add ad-core mrst-gui mpfa agmg linearsolvers + +jsonstruct_material = parseBattmoJson(fullfile('ParameterData','ParameterSets','Chen2020','chen2020_lithium_ion_battery.json')); +jsonstruct_geometry = parseBattmoJson(fullfile('Examples', 'JsonDataFiles', 'geometryChen.json')); + +jsonstruct = mergeJsonStructs({jsonstruct_material, ... + jsonstruct_geometry}); + + + +includeDoubleLayer = false; + +if includeDoubleLayer + + jsonstruct.(ne).(co).(am).(itf).useDoubleLayerCapacity = true; + jsonstruct.(ne).(co).(am).(itf).doubleLayerCapacitance = 0.2; + +end + + +[model, inputparams, ~] = setupModelFromJson(jsonstruct); + +c_ne = 29.866*mol/litre; % initial concentration at negative electrode +c_pe = 17.038*mol/litre; % initial concentration at positive electrode + +initstate = initStateChen2020(model, c_ne, c_pe); +b = logspace(-15, -10, 20); %change here + +frequences = logspace(-2, 4, 30); +figure; +hold on; +for i = 1:length(b) + + inputparams.NegativeElectrode.Coating.ActiveMaterial.SolidDiffusion.referenceDiffusionCoefficient = b(i); %change here + inputparams.PositiveElectrode.Coating.ActiveMaterial.SolidDiffusion.referenceDiffusionCoefficient = b(i)/10; %change here + + impsolv = ImpedanceSolver(inputparams, 'initstate', initstate, 'computeSteadyState', false); + + Z = impsolv.computeImpedance(frequences); + Z_re = real(Z); + Z_im = imag(Z); + curve = sprintf('Diffusion Coefficient = %.2e', b(i)); %change here + plot(Z_re, -Z_im, 'DisplayName', curve); +end + +legend('show'); \ No newline at end of file diff --git a/Battery/EquivalentCircuitModel/calibration/pxdParamsImpact/saturationConcentration_impact.m b/Battery/EquivalentCircuitModel/calibration/pxdParamsImpact/saturationConcentration_impact.m new file mode 100644 index 000000000..ac40f3e89 --- /dev/null +++ b/Battery/EquivalentCircuitModel/calibration/pxdParamsImpact/saturationConcentration_impact.m @@ -0,0 +1,47 @@ +mrstModule add ad-core mrst-gui mpfa agmg linearsolvers + +jsonstruct_material = parseBattmoJson(fullfile('ParameterData','ParameterSets','Chen2020','chen2020_lithium_ion_battery.json')); +jsonstruct_geometry = parseBattmoJson(fullfile('Examples', 'JsonDataFiles', 'geometryChen.json')); + +jsonstruct = mergeJsonStructs({jsonstruct_material, ... + jsonstruct_geometry}); + + + +includeDoubleLayer = false; + +if includeDoubleLayer + + jsonstruct.(ne).(co).(am).(itf).useDoubleLayerCapacity = true; + jsonstruct.(ne).(co).(am).(itf).doubleLayerCapacitance = 0.2; + +end + + +[model, inputparams, ~] = setupModelFromJson(jsonstruct); + +c_ne = 29.866*mol/litre; % initial concentration at negative electrode +c_pe = 17.038*mol/litre; % initial concentration at positive electrode + +initstate = initStateChen2020(model, c_ne, c_pe); + +b = linspace(3e4, 7e4, 20); %change here + +frequences = logspace(-2, 4, 30); +figure; +hold on; +for i = 1:length(b) + + inputparams.NegativeElectrode.Coating.ActiveMaterial.Interface.saturationConcentration = b(i); %change here + inputparams.PositiveElectrode.Coating.ActiveMaterial.Interface.saturationConcentration = b(i)+2e4; %change here + + impsolv = ImpedanceSolver(inputparams, 'initstate', initstate, 'computeSteadyState', false); + + Z = impsolv.computeImpedance(frequences); + Z_re = real(Z); + Z_im = imag(Z); + curve = sprintf('saturation Concentration = %.2e', b(i)); %change here + plot(Z_re, -Z_im, 'DisplayName', curve); +end + +legend('show'); \ No newline at end of file diff --git a/Battery/EquivalentCircuitModel/calibration/pxdParamsImpact/volumefraction_impact.m b/Battery/EquivalentCircuitModel/calibration/pxdParamsImpact/volumefraction_impact.m new file mode 100644 index 000000000..3eff22ccb --- /dev/null +++ b/Battery/EquivalentCircuitModel/calibration/pxdParamsImpact/volumefraction_impact.m @@ -0,0 +1,46 @@ +mrstModule add ad-core mrst-gui mpfa agmg linearsolvers + +jsonstruct_material = parseBattmoJson(fullfile('ParameterData','ParameterSets','Chen2020','chen2020_lithium_ion_battery.json')); +jsonstruct_geometry = parseBattmoJson(fullfile('Examples', 'JsonDataFiles', 'geometryChen.json')); + +jsonstruct = mergeJsonStructs({jsonstruct_material, ... + jsonstruct_geometry}); + + + +includeDoubleLayer = false; + +if includeDoubleLayer + + jsonstruct.(ne).(co).(am).(itf).useDoubleLayerCapacity = true; + jsonstruct.(ne).(co).(am).(itf).doubleLayerCapacitance = 0.2; + +end + + +[model, inputparams, ~] = setupModelFromJson(jsonstruct); + +c_ne = 29.866*mol/litre; % initial concentration at negative electrode +c_pe = 17.038*mol/litre; % initial concentration at positive electrode + +initstate = initStateChen2020(model, c_ne, c_pe); +b = linspace(0.5,0.95,10); + +frequences = logspace(-2, 4, 30); +figure; +hold on; +for i = 1:length(b) + + inputparams.NegativeElectrode.Coating.volumeFraction = b(i); + inputparams.PositiveElectrode.Coating.volumeFraction = b(i); + + impsolv = ImpedanceSolver(inputparams, 'initstate', initstate, 'computeSteadyState', false); + + Z = impsolv.computeImpedance(frequences); + Z_re = real(Z); + Z_im = imag(Z); + curve = sprintf('Volume Fraction = %.2f', b(i)); + plot(Z_re, -Z_im, 'DisplayName', curve); +end + +legend('show'); \ No newline at end of file diff --git a/Battery/EquivalentCircuitModel/calibration/pxdParamsImpact/volumetricSurfaceArea_impact.m b/Battery/EquivalentCircuitModel/calibration/pxdParamsImpact/volumetricSurfaceArea_impact.m new file mode 100644 index 000000000..3e949d075 --- /dev/null +++ b/Battery/EquivalentCircuitModel/calibration/pxdParamsImpact/volumetricSurfaceArea_impact.m @@ -0,0 +1,47 @@ +mrstModule add ad-core mrst-gui mpfa agmg linearsolvers + +jsonstruct_material = parseBattmoJson(fullfile('ParameterData','ParameterSets','Chen2020','chen2020_lithium_ion_battery.json')); +jsonstruct_geometry = parseBattmoJson(fullfile('Examples', 'JsonDataFiles', 'geometryChen.json')); + +jsonstruct = mergeJsonStructs({jsonstruct_material, ... + jsonstruct_geometry}); + + + +includeDoubleLayer = false; + +if includeDoubleLayer + + jsonstruct.(ne).(co).(am).(itf).useDoubleLayerCapacity = true; + jsonstruct.(ne).(co).(am).(itf).doubleLayerCapacitance = 0.2; + +end + + +[model, inputparams, ~] = setupModelFromJson(jsonstruct); + +c_ne = 29.866*mol/litre; % initial concentration at negative electrode +c_pe = 17.038*mol/litre; % initial concentration at positive electrode + +initstate = initStateChen2020(model, c_ne, c_pe); + +b = linspace(1e5, 7e5, 20); %change here + +frequences = logspace(-2, 4, 30); +figure; +hold on; +for i = 1:length(b) + + inputparams.NegativeElectrode.Coating.ActiveMaterial.Interface.volumetricSurfaceArea = b(i); %change here + inputparams.PositiveElectrode.Coating.ActiveMaterial.Interface.volumetricSurfaceArea = b(i); %change here + + impsolv = ImpedanceSolver(inputparams, 'initstate', initstate, 'computeSteadyState', false); + + Z = impsolv.computeImpedance(frequences); + Z_re = real(Z); + Z_im = imag(Z); + curve = sprintf('volumetric Surface Area = %.2e', b(i)); %change here + plot(Z_re, -Z_im, 'DisplayName', curve); +end + +legend('show'); \ No newline at end of file diff --git a/Battery/EquivalentCircuitModel/calibration/runClassFitting.m b/Battery/EquivalentCircuitModel/calibration/runClassFitting.m new file mode 100644 index 000000000..59dbedd95 --- /dev/null +++ b/Battery/EquivalentCircuitModel/calibration/runClassFitting.m @@ -0,0 +1,80 @@ + +%% Main code: + +% filename = 'C:\Users\Alexandre Fichter\Documents\stage_3A\contenu stage\data_August\ank_data\Supplementary material\02_Electrical_characterization\EIS\131-828_EIS_01_MB_CD8.txt'; +% filename = '/home/xavier/Matlab/Projects/battmo/Data/131-828_EIS_01_MB_CD8.txt'; +% +% [Z_re_exp, Z_im_exp, omega] = load_experimental_data(filename); +% [Z_re_exp, Z_im_exp, omega] = load_santoni_data(); +[Z_re_exp, Z_im_exp, omega] = load_chen_data(); + +% omega = logspace(-4, 2, 50); +% params = [0.05052, 1.12673, 59119.9, 0.03155, 11054.0]; +% [Z_re_exp, Z_im_exp] = load_nyquist(params, omega); + +params0 = [0.05052, 1.12673, 59119.9, 0.03155, 11054.0]; % initial condition: C1>2*C2 + +a = 1000; + +pmin = params0 / a; +pmax = params0 * a; +scales = [pmin, pmax]; + +params0_w = [1e-1,1e-1,1e-1,1e-1,1e-2,1,1,1e2,1e-07]; +scales_w = params0_w*10; %to be changed between lsq and bfgs + +feis = FittingEIS(params0, scales, Z_re_exp, Z_im_exp, omega); + +[~, ~, best_params, fitting_error] = feis.optimizationBFGS(); + + +% [min_value, best_params, fitting_error] = feis.optimizationLsqnonlin(); +% % Lancement avec LSQNONLIN + + +%% Printing results +feis.plotresults_thevenin(best_params, fitting_error); + +feis.printResults(best_params, fitting_error); % to be changed btw Warburg and Thevenin + + + + +%% Robustness test +% for i = 1:5 + % params0 = pmin + (pmax - pmin) .* rand(1, 5); + % disp(['Tour ', num2str(i), ' - params0 :']); + % disp(params0); + % feis = FittingEIS(params0, scales, Z_re_exp, Z_im_exp, omega); + % + % [~, ~, best_params, fitting_error] = feis.optimizationBFGS(); + % feis.plotresults_thevenin(best_params, fitting_error); + % feis.printResults(best_params, fitting_error); + % + % drawnow; + +% end +%% + +% result with lsqnonlin: +% === FITTING SCORE === +% Error : 7.324871e-02 +% +% === PARAMETERS FOUND === +% R0 = 0.05052 Ohms +% R1 = 1.12673 Ohms +% C1 = 59119.9 Farads +% R2 = 0.03155 Ohms +% C2 = 11054.0 Farads + + +% results warburg chen +% 4.7197e-02 ... +% 1.0000e-2 ... +% 2.3392e-03 ... +% 0.0103 ... +% 6.8573e-04 ... +% 1.0000e-1 ... +% 0.0307 ... +% 7.0505e+02 ... +% 2.9134e-06 \ No newline at end of file diff --git a/Battery/EquivalentCircuitModel/calibration/runFittingEIS.m b/Battery/EquivalentCircuitModel/calibration/runFittingEIS.m new file mode 100644 index 000000000..33c8d940a --- /dev/null +++ b/Battery/EquivalentCircuitModel/calibration/runFittingEIS.m @@ -0,0 +1,172 @@ +%% Fitting Script + +% initial values [R0, R1, C1, R2, C2] +% params0 = [8.2e-3, 1.5e-3, 1, 1.5e-2, 200]; % handmade best ones + +params0 = [3e-3, 0.00260, 100, 0.00132, 1]; + +[Z_re_exp, Z_im_exp, omega] = load_experimental_data(); + +% lengthZ = length(Z_im_exp); + +% synthetic data +% opti_params = [0.00341, 0.00260, 5663.9, 0.00132, 0.9]; % params with +% lsqnonlin +% +% params with unitboxBFGS: +% R0 = 0.00318 Ohms +% R1 = 0.00247 Ohms +% C1 = 10654.4 Farads +% R2 = 0.00155 Ohms +% C2 = 0.4 Farads +% + +% [Z_re_exp, Z_im_exp] = load_nyquist(opti_params, omega); + +%% + +scales = [10e-2, 10e-3, 15000, 10e-3, 5]; + +params0_norm = params0 ./ scales; + +lb_true = [1e-5, 1e-5, 1e-5, 1e-5, 1e-5]; +lb_norm = lb_true ./ scales; +ub_norm = []; + +f_opt = @(p_norm) optifunc(p_norm, scales, Z_re_exp, Z_im_exp, omega) ; + +% best_params_norm = lsqnonlin(deltagap, params0_norm, lb_norm, ub_norm, [0,0,-1,0,10], 0); +params0_norm = params0_norm(:); + +[min_value, best_params_norm, history] = unitBoxBFGS(... + params0_norm, ... + f_opt, ... + 'maximize', false, ... + 'linIneq', struct('A', [0, 0, -1, 0, 2*scales(5)/scales(3)], 'b', 0), ... + 'enforceFeasible', true, ... + 'maxIt', 300, ... + 'objChangeTol', 1e-6, ... + 'gradTol', 1e-5 ... +); + + +best_params = best_params_norm(:) .* scales(:); + +[v_opt, g_opt] = optifunc(best_params_norm, scales, Z_re_exp, Z_im_exp, omega) + +fprintf('\n=== FITTING SCORE ===\n'); +fprintf('Error : %e\n', v_opt); + +fprintf('\n=== PARAMETERS FOUND ===\n'); +fprintf('R0 = %.5f Ohms\n', best_params(1)); +fprintf('R1 = %.5f Ohms\n', best_params(2)); +fprintf('C1 = %.1f Farads\n', best_params(3)); +fprintf('R2 = %.5f Ohms\n', best_params(4)); +fprintf('C2 = %.1f Farads\n', best_params(5)); + +%% Graphique + + +parameters.R0 = best_params(1); +parameters.R1 = best_params(2); +parameters.C1 = best_params(3); +parameters.R2 = best_params(4); +parameters.C2 = best_params(5); + +[Z_re_fit, Z_im_fit] = load_nyquist(best_params, omega); + + +%% +figure; +subplot(3,1,1); +semilogx(omega, Z_re_exp, 'ro', 'MarkerFaceColor', 'r'); +hold on; +semilogx(omega, Z_re_fit, 'bo', 'LineWidth', 2); +% legend('Expérimental', 'Modèle (Fitted)'); +title('Fitting results'); +xlabel('Omega'); +ylabel('Z_{re} '); + +subplot(3,1,2); +semilogx(omega, Z_im_exp, 'ro', 'MarkerFaceColor', 'r'); +hold on; +semilogx(omega, Z_im_fit, 'bo', 'LineWidth', 2); +% legend('Expérimental', 'Modèle (Fitted)'); +title('Fitting results'); +xlabel('Omega'); +ylabel('-Z_{im} '); + +subplot(3,1,3); +plot(Z_re_exp, Z_im_exp, 'ro', 'MarkerFaceColor', 'r'); +hold on; +plot(Z_re_exp, Z_im_fit, 'bo', 'LineWidth', 2); +% legend('Expérimental', 'Modèle (Fitted)'); +title('Nyquist'); +xlabel('Z_{re}'); +ylabel('-Z_{im} '); +axis equal; + +grid on; + +text_error = sprintf('Fitting error : %.2e', v_opt); + +subplot(3,1,3); + +text(0.05, 0.90, text_error, 'Units', 'normalized', ... + 'BackgroundColor', 'white', ... + 'EdgeColor', 'black', ... + 'FontSize', 11, ... + 'FontWeight', 'bold'); + +function [v, g_norm] = optifunc(p_norm, scales, Z_re_exp, Z_im_exp, omega) + + p = p_norm(:) .* scales(:); + + [Z_re_param, Z_Im_param] = load_nyquist(p, omega); + + Module_Z = sqrt(Z_re_exp.^2 + Z_im_exp.^2); + Module_Z = Module_Z(:); + + err_re = (Z_re_exp(:) - Z_re_param(:)) ./ Module_Z(:); + err_im = (Z_im_exp(:) - Z_Im_param(:)) ./ Module_Z(:); + + v = sum(err_re.^2 + err_im.^2); + + +if nargout > 1 + w = omega(:); + R0 = p(1); + R1 = p(2); + C1 = p(3); + R2 = p(4); + C2 = p(5); + + + g_re_dR0 = -ones(size(w)); + g_re_dR1 = -(1-(R1*C1.*w).^2) ./ (1+(R1*C1.*w).^2).^2; + g_re_dC1 = (2*C1*R1^3.*w.^2) ./ (1+(R1*C1.*w).^2).^2; + g_re_dR2 = -(1-(R2*C2.*w).^2) ./ (1+(R2*C2.*w).^2).^2; + g_re_dC2 = (2*C2*R2^3.*w.^2) ./ (1+(R2*C2.*w).^2).^2; + + J_re = [g_re_dR0, g_re_dR1, g_re_dC1, g_re_dR2, g_re_dC2]; + + + g_im_dR0 = zeros(size(w)); + g_im_dR1 = -(2*R1*C1.*w) ./ (1+(R1*C1.*w).^2).^2; + g_im_dC1 = -(R1^2.*w - C1^2*R1^4.*w.^3) ./ (1+(R1*C1.*w).^2).^2; + g_im_dR2 = -(2*R2*C2.*w) ./ (1+(R2*C2.*w).^2).^2; + g_im_dC2 = -(R2^2.*w - C2^2*R2^4.*w.^3) ./ (1+(R2*C2.*w).^2).^2; + + J_im = [g_im_dR0, g_im_dR1, g_im_dC1, g_im_dR2, g_im_dC2]; + derr_re_dp = J_re ./ Module_Z; + derr_im_dp = J_im ./ Module_Z; + + + g_true = 2 * (derr_re_dp' * err_re) + 2 * (derr_im_dp' * err_im); + g_norm = g_true .* scales(:); +end + +end + + + diff --git a/Battery/EquivalentCircuitModel/calibration/utils/costFunctionEIS.m b/Battery/EquivalentCircuitModel/calibration/utils/costFunctionEIS.m new file mode 100644 index 000000000..28a284906 --- /dev/null +++ b/Battery/EquivalentCircuitModel/calibration/utils/costFunctionEIS.m @@ -0,0 +1,15 @@ +function [total_error, Z_re_exp, Z_im_exp] = costFunctionEIS(best_params) + + + Z_re_exp = plot_nyquist_ank.Z_real; + Z_im_exp = plot_nyquist_ank.Z_imag; + + + omega = logspace(-2, 4, 100) + + Z_re_params = plot_nyquist_params.Z_real(best_params); + Z_im_params = plot_nyquist_params.Z_imag(best_params); + + totalerror = sum((Z_re_exp(omega)-Z_re_params(omega)).^2) + sum((Z_im_exp(omega)-Z_im_params(omega))^2) + +end diff --git a/Battery/EquivalentCircuitModel/calibration/utils/load_chen_data.m b/Battery/EquivalentCircuitModel/calibration/utils/load_chen_data.m new file mode 100644 index 000000000..a05b021e5 --- /dev/null +++ b/Battery/EquivalentCircuitModel/calibration/utils/load_chen_data.m @@ -0,0 +1,56 @@ +function [Z_real, Z_imag, omegas] = load_chen_data() + + + % We define some shorthand names for simplicity. + ne = 'NegativeElectrode'; + pe = 'PositiveElectrode'; + elyte = 'Electrolyte'; + thermal = 'ThermalModel'; + co = 'Coating'; + am = 'ActiveMaterial'; + itf = 'Interface'; + sd = 'SolidDiffusion'; + ctrl = 'Control'; + cc = 'CurrentCollector'; + + mrstModule add ad-core mrst-gui mpfa agmg linearsolvers + + jsonstruct_material = parseBattmoJson(fullfile('ParameterData','ParameterSets','Chen2020','chen2020_lithium_ion_battery.json')); + jsonstruct_geometry = parseBattmoJson(fullfile('Examples', 'JsonDataFiles', 'geometryChen.json')); + + jsonstruct = mergeJsonStructs({jsonstruct_material, ... + jsonstruct_geometry}); + + + + includeDoubleLayer = true; + + if includeDoubleLayer + + jsonstruct.(ne).(co).(am).(itf).useDoubleLayerCapacity = true; + jsonstruct.(ne).(co).(am).(itf).doubleLayerCapacitance = 0.2; + + end + + [model, inputparams, ~] = setupModelFromJson(jsonstruct); + + c_ne = 29.866*mol/litre; % initial concentration at negative electrode + c_pe = 17.038*mol/litre; % initial concentration at positive electrode + + initstate = initStateChen2020(model, c_ne, c_pe); + + options = []; + options.stateInitialization.initializationSetup = 'given state'; + options.stateInitialization.computeSteadyState = false; + + extrastructs = []; + extrastructs.initstate = initstate; + + impsolv = ImpedanceSolver(inputparams, options, extrastructs); + + omegas = logspace(-4, 2, 50); + Z = impsolv.computeImpedance(omegas); + + Z_real = real(Z); + Z_imag = imag(Z); +end diff --git a/Battery/EquivalentCircuitModel/calibration/utils/load_experimental_data.m b/Battery/EquivalentCircuitModel/calibration/utils/load_experimental_data.m new file mode 100644 index 000000000..b2e0c2643 --- /dev/null +++ b/Battery/EquivalentCircuitModel/calibration/utils/load_experimental_data.m @@ -0,0 +1,55 @@ +function [Z_real, Z_imag, omega] = load_experimental_data(filename) + + % --- 1. Configuration et Importation des données --- + % filename = 'C:\Users\Alexandre Fichter\Documents\stage_3A\contenu stage\data_August\ank_data\Supplementary material\02_Electrical_characterization\EIS\131-828_EIS_01_MB_CD8.txt'; + + % Le délimiteur est la tabulation ('\t') + opts = detectImportOptions(filename, 'Delimiter', '\t'); + + opts.VariableNamesLine = 1; + opts.DataLine = 2; + opts.VariableNamingRule = 'preserve'; + data = readtable(filename, opts); + colNames = data.Properties.VariableNames; + + + idx_real = contains(colNames, 'Re(Z)', 'IgnoreCase', true); + idx_imag = contains(colNames, 'Im(Z)', 'IgnoreCase', true); + idx_freq = contains(colNames, 'freq', 'IgnoreCase', true); + + % Extraction des données brutes (peu importe où elles se trouvent) + Z_real_raw = data{:, idx_real}; + Z_imag_raw = data{:, idx_imag}; + freq_raw = data{:, idx_freq}; + + if iscell(Z_real_raw) || isstring(Z_real_raw) + Z_real = str2double(strrep(Z_real_raw, ',', '.')); + Z_imag = str2double(strrep(Z_imag_raw, ',', '.')); + freq = str2double(strrep(freq_raw, ',', '.')); + else + Z_real = Z_real_raw; + Z_imag = Z_imag_raw; + freq = freq_raw; + end + + idx_valides = (Z_real ~= 0) & ~isnan(Z_real); + Z_real = Z_real(idx_valides); + Z_imag = -Z_imag(idx_valides); % in the file is given -Im + freq = freq(idx_valides); + + Z_real = Z_real(8:60); + Z_imag = Z_imag(8:60); + omega = 2*pi*freq(8:60); + + Z_real = Z_real(:); + Z_imag = Z_imag(:); + omega = omega(:); +%% Figure test + % figure; + % semilogx(omega, Z_real, 'ro', 'MarkerFaceColor', 'r'); + % % legend('Expérimental', 'Modèle (Fitted)'); + % title('Fitting results'); + % xlabel('Omega'); + % ylabel('Z_{re} '); + +end diff --git a/Battery/EquivalentCircuitModel/calibration/utils/load_nyquist.m b/Battery/EquivalentCircuitModel/calibration/utils/load_nyquist.m new file mode 100644 index 000000000..1cad56a49 --- /dev/null +++ b/Battery/EquivalentCircuitModel/calibration/utils/load_nyquist.m @@ -0,0 +1,15 @@ +function [Z_real, Z_imag] = load_nyquist(params, omega) + + R_0 = params(1); + R_1 = params(2); + C_1 = params(3); + R_2 = params(4); + C_2 = params(5); + + + + Z_real = R_0 + R_1 ./ (1+(R_1* C_1.*omega).^2) + R_2 ./ (1+(R_2* C_2.*omega).^2); + Z_imag = -(R_1*R_1*C_1.*omega ./ (1+(R_1* C_1.*omega).^2)) ... + - (R_2*R_2*C_2.*omega ./ (1+(R_2* C_2.*omega).^2)); + +end \ No newline at end of file diff --git a/Battery/EquivalentCircuitModel/calibration/utils/load_santoni_data.m b/Battery/EquivalentCircuitModel/calibration/utils/load_santoni_data.m new file mode 100644 index 000000000..5a84bcf63 --- /dev/null +++ b/Battery/EquivalentCircuitModel/calibration/utils/load_santoni_data.m @@ -0,0 +1,12 @@ +function [Z_real, Z_imag, omegas] = load_santoni_data() + + path_data = 'C:\Users\Alexandre Fichter\Documents\stage_3A\contenu stage\papers\LiPo Battery LP-503562-IS-3 EIS, Capacity, ECM Data\LiPO_1\EIS_Charge_discharge\EIS_45\1_EIS.csv'; + data_exp = readtable(path_data, 'Delimiter', '\t'); + + + freq = table2array(data_exp(:, 1)); + Z_real = table2array(data_exp(:, 2)); + Z_imag = table2array(data_exp(:, 3)); + + omegas = 2 * pi .* freq; +end \ No newline at end of file diff --git a/Battery/EquivalentCircuitModel/examples/exampleECM.m b/Battery/EquivalentCircuitModel/examples/exampleECM.m new file mode 100644 index 000000000..e555eb2c0 --- /dev/null +++ b/Battery/EquivalentCircuitModel/examples/exampleECM.m @@ -0,0 +1,38 @@ +%% Equivalent Circuit Model Simulation + +% Load data set + +params = createParametersECM(); +inputparams = EquivalentCircuitModelInputParams(params); + +% Setup model +model = EquivalentCircuitModel(inputparams); + +% Run Simulation +[t, U, I, SOC] = model.solve(); + + +% Plotting +figure(Name='Highlighting the difference in resistance depending on the state of charge'); +tiledlayout(3, 1); + +nexttile +plot(t, U, 'LineWidth', 2) +title('Battery Voltage') +xlabel('Time /s') +ylabel('Voltage /V') +grid on + +nexttile +plot(t, I, 'LineWidth', 2, 'Color', 'r') +title('Applied Current') +xlabel('Time /s') +ylabel('Current /A') +grid on + +nexttile +plot(t, SOC, 'LineWidth', 2, 'Color', 'g') +title('SOC') +xlabel('Time /s') +ylabel('SOC') +grid on diff --git a/Battery/EquivalentCircuitModel/examples/notebooks/calibrationECM.mlx b/Battery/EquivalentCircuitModel/examples/notebooks/calibrationECM.mlx new file mode 100644 index 000000000..f158497b4 Binary files /dev/null and b/Battery/EquivalentCircuitModel/examples/notebooks/calibrationECM.mlx differ diff --git a/Battery/EquivalentCircuitModel/examples/notebooks/exampleECM.ipynb b/Battery/EquivalentCircuitModel/examples/notebooks/exampleECM.ipynb new file mode 100644 index 000000000..6dafa2118 --- /dev/null +++ b/Battery/EquivalentCircuitModel/examples/notebooks/exampleECM.ipynb @@ -0,0 +1,104 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "# Equivalent Circuit Model\n", + "\n", + "## Load data set" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "source": [ + "jsonstruct = parseBattmoJson('exampleECM.json');\n", + "inputparams = EquivalentCircuitModelInputParams(jsonstruct);" + ], + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Setup model" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "source": [ + "model = EquivalentCircuitModel(inputparams);" + ], + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## run Simulation" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "source": [ + "[t, U, I] = model.solve();" + ], + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Plotting" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "source": [ + "tiledlayout(1, 2);\n", + "nexttile\n", + "plot(t, U)\n", + "nexttile\n", + "plot(t, I)" + ], + "outputs": [ + { + "data": { + "text/html": [ + "
\"figure_0.png\"
" + ] + }, + "metadata": {}, + "execution_count": 4, + "output_type": "execute_result" + } + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "MATLAB (matlabkernel)", + "language": "matlab", + "name": "matlab" + }, + "language_info": { + "file_extension": ".m", + "mimetype": "text/matlab", + "name": "matlab", + "nbconvert_exporter": "matlab", + "pygments_lexer": "matlab", + "version": "24.1.0.2689473" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} \ No newline at end of file diff --git a/Battery/EquivalentCircuitModel/examples/notebooks/exampleECM.mlx b/Battery/EquivalentCircuitModel/examples/notebooks/exampleECM.mlx new file mode 100644 index 000000000..3274a0afa Binary files /dev/null and b/Battery/EquivalentCircuitModel/examples/notebooks/exampleECM.mlx differ diff --git a/Battery/EquivalentCircuitModel/examples/runECMfitting.m b/Battery/EquivalentCircuitModel/examples/runECMfitting.m new file mode 100644 index 000000000..7f908d14d --- /dev/null +++ b/Battery/EquivalentCircuitModel/examples/runECMfitting.m @@ -0,0 +1,48 @@ +%% Fitting Script + +% Charge les données expérimentales +t_exp = [0:1:1000]'; % Temps en secondes +V_exp = 4.1 * ones(1001, 1) - 0.1 * (t_exp/1000); % Fausse courbe de tension +I_exp = ones(1001, 1); % Courant expérimental + + +parameters = createParametersECM(); + +% Remplacer le courant du modèle par le courant expérimental +parameters.I.dataX = t_exp; +parameters.I.dataY = I_exp; + +% Valeurs initiales [R0, R1, C1, R2, C2] +initial_values = [0.010, 0.010, 3000, 0.010, 100000]; + +options = optimset('Display', 'iter', 'TolX', 1e-4); +disp('beginning of optimization'); + +% Optimization with fminsearch +best_params = fminsearch(@(p) costFunctionECM(p, parameters, initial_values, t_exp, V_exp), initial_values, options); + +fprintf('\n=== PARAMETERS FOUND ===\n'); +fprintf('R0 = %.5f Ohms\n', best_params(1)); +fprintf('R1 = %.5f Ohms\n', best_params(2)); +fprintf('C1 = %.1f Farads\n', best_params(3)); +fprintf('R2 = %.5f Ohms\n', best_params(4)); +fprintf('C2 = %.1f Farads\n', best_params(5)); + +%% Vérification Graphique + + +parameters.R0 = best_params(1); +parameters.R1 = best_params(2); +parameters.C1 = best_params(3); +parameters.R2 = best_params(4); +parameters.C2 = best_params(5); + +inputparams = EquivalentCircuitModelInputParams(parameters); +model = EquivalentCircuitModel(inputparams); +[t_sim, V_sim, ~] = model.solve(); + +figure; +plot(t_exp, V_exp, 'k', 'LineWidth', 2); hold on; +plot(t_sim, V_sim, 'r--', 'LineWidth', 2); +legend('Expérimental', 'Modèle (Fitted)'); +title('Résultat du Fitting'); diff --git a/Battery/EquivalentCircuitModel/models/EquivalentCircuitModel.m b/Battery/EquivalentCircuitModel/models/EquivalentCircuitModel.m new file mode 100644 index 000000000..ff38ffcc0 --- /dev/null +++ b/Battery/EquivalentCircuitModel/models/EquivalentCircuitModel.m @@ -0,0 +1,231 @@ +classdef EquivalentCircuitModel < BaseModel + + properties + + nominalcellcapacity + + OCP + R0 + R1 + C1 + R2 + C2 + initSOC + initOverpotential2 + lowerVoltageCutoff + + I + totalTime + + %% helpers + OCPfunc + Ifunc + R0func + R1func + C1func + R2func + C2func + end + + methods + + function model = EquivalentCircuitModel(inputparams) + + model = model@BaseModel(); + + %% Setup the model using the input parameters + + fdnames = {'nominalcellcapacity', ... + 'OCP' , ... + 'R0' , ... + 'R1' , ... + 'C1' , ... + 'R2' , ... + 'C2' , ... + 'initSOC' , ... + 'initOverpotential2' , ... + 'lowerVoltageCutoff' , ... + 'I' , ... + 'totalTime'}; + + model = dispatchParams(model, inputparams, fdnames); + + model.OCPfunc = setupFunction(model.OCP); + model.Ifunc = setupFunction(model.I); + model.R0func = setupFunction(model.R0); + model.R1func = setupFunction(model.R1); + model.C1func = setupFunction(model.C1); + model.R2func = setupFunction(model.R2); + model.C2func = setupFunction(model.C2); + + end + + function model = registerVarAndPropfuncNames(model) + + model = registerVarAndPropfuncNames@BaseModel(model); + + varnames = {}; + % Current + varnames{end + 1} = 'I'; + % + varnames{end + 1} = 'OCP'; + % + varnames{end + 1} = 'U'; + + varnames{end + 1} = 'UR'; + varnames{end + 1} = 'U1'; + varnames{end + 1} = 'U2'; + varnames{end + 1} = 'SOC'; + + varnames{end + 1} = 'dU1dt'; + varnames{end + 1} = 'dU2dt'; + varnames{end + 1} = 'dSOCdt'; + + model = model.registerVarNames(varnames); + + fn = @EquivalentCircuitModel.updateUR; + inputnames = {'I'}; + model = model.registerPropFunction({'UR', fn, inputnames}); + + fn = @EquivalentCircuitModel.updatedU1dt; + inputnames = {'U1', 'I'}; + model = model.registerPropFunction({'dU1dt', fn, inputnames}); + + fn = @EquivalentCircuitModel.updatedU2dt; + inputnames = {'U2', 'I'}; + model = model.registerPropFunction({'dU2dt', fn, inputnames}); + + fn = @EquivalentCircuitModel.updatedSOCdt; + inputnames = {'I'}; + model = model.registerPropFunction({'dSOCdt', fn, inputnames}); + + fn = @EquivalentCircuitModel.updateU; + inputnames = {'SOC', 'U1', 'U2', 'UR'}; + model = model.registerPropFunction({'U', fn, inputnames}); + + end + + function state = updateUR(model, state) + + state.UR = model.R0func(state.SOC)*state.I; + + end + + function state = updatedU2dt(model, state) + + state.dU2dt = model.updatedUdt(state.U2, model.R2func(state.SOC), model.C2func(state.SOC), state.I); + + end + + function state = updatedU1dt(model, state) + + state.dU1dt = model.updatedUdt(state.U1, model.R1func(state.SOC), model.C1func(state.SOC), state.I); + + end + + function state = updatedSOCdt(model, state) + + Qmax = model.nominalcellcapacity; + + state.dSOCdt = -(1 / (Qmax*3600)) * state.I; + + end + + function state = updateU(model, state) + + OCP = model.OCPfunc(state.SOC); + state.U = OCP - (state.U1 + state.U2 + state.UR); + + end + + function state0 = setupInitialCondition(model, SOC0) + + state0.U1 = 0; + state0.U2 = 0; + state0.SOC = SOC0; + + end + + function state = setupStateFromY(model, y) + + state.U1 = y(1); + state.U2 = y(2); + state.SOC = y(3); + + end + + function y = setupYfromState(model, state) + + y(1) = state.U1; + y(2) = state.U2; + y(3) = state.SOC; + + end + + function y = setupFfromState(model, state) + + y(1) = state.dU1dt; + y(2) = state.dU2dt; + y(3) = state.dSOCdt; + y = y'; + + end + + function f = ode(model, t, y) + + state = model.setupStateFromY(y); + + state.I = model.Ifunc(t); + + state = model.evalVarName(state, 'dU1dt'); + state = model.evalVarName(state, 'dU2dt'); + state = model.evalVarName(state, 'dSOCdt'); + + f = model.setupFfromState(state); + + end + + function [t, U, I, SOC] = solve(model) + + if isempty(model.computationalGraph) + model = setupComputationalGraph(model); + end + + % setup initial condition + state0 = model.setupInitialCondition(model.initSOC); + y0 = model.setupYfromState(state0); + + % setup time span + tspan = [0, model.totalTime]; + + % solve ode + [t_out, y_out] = ode45(@(t, y) model.ode(t, y), tspan, y0); + + SOC = zeros(length(t_out), 1); + U = zeros(length(t_out), 1); + I = zeros(length(t_out), 1); %gagne du temps en allouant déjà la mémoire + for i = 1:length(t_out) + + state = model.setupStateFromY(y_out(i, :)); + I(i) = model.Ifunc(t_out(i)); + state.I = I(i); + state = model.evalVarName(state, 'U'); + U(i) = state.U; + state = model.evalVarName(state, 'SOC'); + SOC(i) = state.SOC; + end + + t = t_out; + end + + end + + methods (Static) + + function dUdt = updatedUdt(U, R, C, I) + dUdt = -(1 / (R * C)) * U + (1 / C) * I; + end + + end + +end diff --git a/Battery/EquivalentCircuitModel/models/EquivalentCircuitModelInputParams.m b/Battery/EquivalentCircuitModel/models/EquivalentCircuitModelInputParams.m new file mode 100644 index 000000000..10b6a0e12 --- /dev/null +++ b/Battery/EquivalentCircuitModel/models/EquivalentCircuitModelInputParams.m @@ -0,0 +1,55 @@ +classdef EquivalentCircuitModelInputParams < InputParams +% +% Input parameter class for the :code:`Battery` model. +% + + properties + + nominalcellcapacity + OCP + R0 + R1 + C1 + R2 + C2 + initSOC + initOverpotential2 + lowerVoltageCutoff + I + totalTime + + end + + methods + + function inputparams = EquivalentCircuitModelInputParams(jsonstruct) + + inputparams = inputparams@InputParams(jsonstruct); + + end + + end + +end + + + +%{ +Copyright 2021-2024 SINTEF Industry, Sustainable Energy Technology +and SINTEF Digital, Mathematics & Cybernetics. + +This file is part of The Battery Modeling Toolbox BattMo + +BattMo is free software: you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +BattMo is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with BattMo. If not, see . +%} diff --git a/Battery/EquivalentCircuitModel/utils/costFunctionECM.m b/Battery/EquivalentCircuitModel/utils/costFunctionECM.m new file mode 100644 index 000000000..0c12089f5 --- /dev/null +++ b/Battery/EquivalentCircuitModel/utils/costFunctionECM.m @@ -0,0 +1,27 @@ +function total_error = costFunctionECM(params_test, parameters, initial_values, t_exp, V_exp) + + + parameters.R0 = abs(params_test(1)); + parameters.R1 = abs(params_test(2)); + parameters.C1 = max(abs(params_test(3)), 0.1); + parameters.R2 = abs(params_test(4)); + parameters.C2 = max(abs(params_test(5)), 0.1); + + inputparams = EquivalentCircuitModelInputParams(parameters); + model = EquivalentCircuitModel(inputparams); + [t_sim, V_sim, ~] = model.solve(); + + % Time synchronization + V_sim_aligne = interp1(t_sim, V_sim, t_exp, 'linear', 'extrap'); + + voltage_error = sum((V_exp - V_sim_aligne).^2); + + lambda = 0.5; + + rel_error = sum((params_test-initial_values)./initial_values); + + tikhonov_error = lambda * rel_error.^2; + + total_error = voltage_error + tikhonov_error; + +end diff --git a/Battery/EquivalentCircuitModel/utils/createParametersECM.m b/Battery/EquivalentCircuitModel/utils/createParametersECM.m new file mode 100644 index 000000000..4e2c35f7d --- /dev/null +++ b/Battery/EquivalentCircuitModel/utils/createParametersECM.m @@ -0,0 +1,59 @@ +function p = createParametersECM() + + p.nominalcellcapacity = 63.7; + p.initSOC = 1.0; + p.lowerVoltageCutoff = 3.0; + p.initOverpotential2 = 0; + p.totalTime = 4200; + + p.I.functionFormat = 'tabulated'; + p.I.argumentList = {'time'}; + p.I.dataX = [0, 100, 100.1, 120, 120.1, 300, 300.1, 3800, 3800.1, 4000, 4000.1, 4020, 4020.1, 4200]; + p.I.dataY = [0, 0, 60, 60, 0, 0, 60, 60, 0, 0, 60, 60, 0, 0]; + + p.OCP.functionFormat = 'tabulated'; + p.OCP.argumentList = {'SOC'}; + p.OCP.dataX = [0.00, 0.02, 0.04, 0.06, 0.08, 0.10, 0.12, 0.14, 0.16, 0.18, 0.20, 0.22, 0.24, 0.26, 0.28, 0.30, 0.32, 0.34, 0.36, 0.38, 0.40, 0.42, 0.44, 0.46, 0.48, 0.50, 0.52, 0.54, 0.56, 0.58, 0.60, 0.62, 0.64, 0.66, 0.68, 0.70, 0.72, 0.74, 0.76, 0.78, 0.80, 0.82, 0.84, 0.86, 0.88, 0.90, 0.92, 0.94, 0.96, 0.98, 1.00]; + p.OCP.dataY = [3.126000, 3.273000, 3.380000, 3.447000, 3.470000, 3.480000, 3.489000, 3.500000, 3.515000, 3.533000, 3.551000, 3.567000, 3.582000, 3.594000, 3.602000, 3.609000, 3.614000, 3.620000, 3.625000, 3.631000, 3.637000, 3.643000, 3.649000, 3.657000, 3.665000, 3.674000, 3.684000, 3.696000, 3.710000, 3.727000, 3.750000, 3.775000, 3.798000, 3.818000, 3.838000, 3.858000, 3.877000, 3.897000, 3.917000, 3.937000, 3.958000, 3.979000, 4.000000, 4.022000, 4.044000, 4.066000, 4.089000, 4.112000, 4.135000, 4.159000, 4.182000]; + + R0_tot = [0.0074814,0.00375135,0.00199926,0.000808379,0.000529892;0.0074814,0.00375135,0.00185474,0.000738337,0.00049771;0.0074814,0.00373265,0.00151744,0.000673703,0.000468303;0.0074814,0.00342697,0.00130825,0.000635903,0.000448193;0.0074814,0.00272078,0.00117717,0.000611909,0.000432971;0.0069826,0.00229452,0.00109773,0.00058627,0.000423262;0.00612671,0.00206227,0.00105217,0.000571475,0.000417812;0.00536849,0.00193337,0.00101405,0.000562212,0.000409251;0.00480278,0.00185103,0.000983687,0.000549631,0.000410476;0.00440167,0.00178588,0.000957044,0.000543895,0.000405368;0.00411575,0.00172477,0.000934531,0.000533558,0.0003994;0.00390187,0.00167949,0.000914869,0.000526768,0.000397983;0.00374658,0.00164546,0.000905014,0.000521591,0.000396355;0.0036334,0.00160835,0.000895256,0.000517151,0.000393944;0.00353281,0.00158063,0.0008793,0.000516464,0.000391954;0.00344791,0.0015608,0.000869252,0.00051227,0.000391406;0.00338977,0.00154377,0.000859461,0.000509065,0.000390664;0.00333872,0.00152388,0.000857998,0.000507685,0.00038969;0.00329657,0.00150926,0.000852066,0.000501813,0.000385995;0.00326141,0.00149673,0.000846844,0.000505136,0.000388294;0.00322799,0.00149188,0.00084232,0.00050115,0.000384211;0.00320237,0.00148102,0.000839376,0.00050073,0.000386191;0.00318137,0.00147045,0.000836217,0.000500865,0.000386833;0.00315527,0.00146685,0.000833235,0.000497888,0.000386251;0.00311928,0.00145992,0.000832092,0.000498747,0.000382866;0.00310424,0.00145576,0.000832276,0.000494884,0.000386295;0.00309834,0.00144722,0.000829437,0.000496968,0.000386992;0.00307867,0.00144351,0.000828243,0.000495848,0.000387399;0.00306673,0.00144798,0.000830089,0.000500481,0.000387638;0.00305254,0.00144507,0.000828833,0.000496636,0.000382113;0.00304161,0.00144131,0.000827309,0.000497478,0.000382439;0.0030343,0.00143125,0.00082455,0.000492916,0.000389673;0.00302615,0.00142702,0.000818782,0.000500463,0.00038546;0.00302297,0.0014205,0.000814619,0.000491503,0.0003858;0.00301535,0.00141463,0.000814864,0.000492482,0.000383651;0.00300786,0.00140957,0.000816839,0.000493931,0.000382345;0.00299892,0.00140604,0.000813864,0.000493379,0.000386187;0.00299335,0.00140043,0.000808475,0.000490454,0.000383351;0.0029867,0.00139212,0.000804998,0.00048877,0.000377422;0.00298272,0.00138903,0.000800962,0.000487252,0.000380367;0.00298442,0.00138775,0.000805805,0.000489012,0.000381427;0.00297911,0.00138254,0.000802961,0.000485793,0.000382455;0.0029803,0.00137618,0.000800822,0.000484924,0.000378882;0.00299019,0.0013734,0.000797575,0.000489591,0.000378052;0.00300121,0.00137477,0.000797567,0.00048441,0.000379647;0.00301775,0.00136948,0.000800597,0.000486859,0.000378768;0.00304181,0.00136736,0.000799303,0.000482635,0.000384018;0.00307807,0.00136092,0.000796617,0.000485164,0.000379421;0.00313868,0.00136725,0.000791133,0.000485946,0.000378392;0.00319877,0.00136795,0.000797194,0.000486639,0.000378749;0.00319877,0.00136795,0.000798674,0.000488616,0.00038044]; + R1_tot = [0.0109518,0.00514849,0.00241913,0.000853767,0.000473136;0.0109518,0.00514849,0.00220007,0.000747948,0.000445283;0.0109518,0.00501759,0.00156633,0.000667537,0.000418055;0.0109518,0.00370479,0.00129477,0.000664379,0.000438201;0.0109518,0.00216514,0.00106482,0.000640611,0.000507597;0.00728213,0.00153036,0.000911211,0.000570174,0.000443806;0.00408898,0.00127121,0.000813038,0.000536964,0.000411202;0.00271216,0.00113093,0.000750563,0.000509078,0.00039686;0.00214054,0.00103047,0.000707957,0.00049128,0.000367506;0.00183948,0.000958516,0.000681101,0.000487984,0.000367436;0.00161631,0.000900906,0.00065605,0.000471902,0.000371049;0.00143562,0.000870054,0.000636391,0.000458368,0.000371823;0.00131615,0.000855439,0.000624021,0.000446287,0.000361092;0.00125805,0.000850114,0.000654326,0.000460409,0.000374051;0.00124731,0.000866505,0.000676043,0.000473626,0.000375777;0.00125896,0.000869833,0.00068294,0.000483716,0.000367479;0.00126551,0.000878883,0.000685236,0.000493259,0.000394196;0.00126671,0.000886034,0.000680167,0.000500364,0.000410745;0.00126922,0.000890457,0.000679576,0.00050932,0.000430389;0.00126189,0.000899098,0.000680305,0.000511495,0.000425145;0.00124402,0.000907082,0.000684716,0.000517428,0.000411099;0.00121491,0.000919332,0.000686495,0.00051921,0.000407804;0.00118622,0.000932087,0.000694829,0.000514845,0.000406396;0.00115509,0.000939311,0.000704886,0.000514129,0.000393472;0.00114401,0.000955183,0.000713811,0.000514671,0.000395118;0.00112723,0.000972638,0.00072264,0.000515396,0.000390676;0.00112656,0.000993509,0.000733566,0.000513148,0.000391352;0.00115195,0.00101087,0.000750226,0.000523251,0.000395854;0.0011324,0.00100369,0.000766549,0.000531557,0.000404817;0.00109983,0.000968004,0.000760227,0.000547237,0.000416704;0.00106196,0.000866661,0.000703058,0.000514665,0.000397599;0.00103714,0.000801982,0.000640288,0.00051886,0.000418776;0.00103666,0.000804728,0.000635207,0.000576802,0.000470891;0.00106092,0.000845936,0.000689149,0.000600804,0.00046715;0.00114402,0.000924794,0.000766016,0.000589737,0.000446383;0.00127793,0.00100958,0.000815813,0.000568085,0.000426063;0.0014493,0.00106102,0.000822573,0.000544546,0.00040995;0.00160054,0.001059,0.000798524,0.000529874,0.000401881;0.00167306,0.00102742,0.000775283,0.000512792,0.000385689;0.0016553,0.000997881,0.000747467,0.000499305,0.000378046;0.0015888,0.000969107,0.0007177,0.000486984,0.000365287;0.00153194,0.000943878,0.000695967,0.000473997,0.000353817;0.00148159,0.000922603,0.000676232,0.000464284,0.000347846;0.0014293,0.00090407,0.000661561,0.000451668,0.00033672;0.00138895,0.000883109,0.000647176,0.000443648,0.00033387;0.00135019,0.000864676,0.000632936,0.000434683,0.000323263;0.00130261,0.000848781,0.000617294,0.00042993,0.000315138;0.00125322,0.000834973,0.00060042,0.000418737,0.000311007;0.00119786,0.000809893,0.000597469,0.000413424,0.000304475;0.00114129,0.000782067,0.000524998,0.000399362,0.000286744;0.00114129,0.000782067,0.000511038,0.000349277,0.000271429]; + tau1_tot=[4.9,11,12,17,15;4.9,11,12,18,17;4.9,13,14,17,18;4.9,13,20,16,19;4.9,12,17,16,18;3.7,13,16,14,15;3.2,14,17,15,14;3.2,16,18,15,15;4.3,16,17,15,15;5.9,16,17,15,13;8.8,17,17,15,15;14,18,17,16,15;18,20,18,16,14;21,23,22,18,16;23,24,24,19,17;23,26,25,20,16;23,26,24,23,21;22,26,23,23,23;22,24,22,24,22;21,23,21,22,20;19,23,20,21,19;19,23,20,20,17;18,23,19,19,16;17,23,19,18,15;16,23,19,18,15;15,24,19,17,14;15,25,20,17,14;17,26,20,18,15;17,27,22,18,14;20,27,24,20,16;25,29,25,26,21;23,27,22,22,23;21,22,20,27,24;21,21,21,24,21;21,23,25,22,18;22,26,26,20,18;24,28,23,21,18;27,27,21,20,18;31,25,21,21,17;31,24,20,21,18;29,23,19,21,17;26,23,20,20,18;24,23,20,21,19;22,23,20,20,19;21,23,20,22,19;19,23,21,23,17;19,24,21,21,18;18,24,20,20,18;17,24,21,23,19;17,23,15,22,19;17,23,14,15,17]; + R2_tot = [0.00465917,0.00296137,0.00671937,0.00670953,0.00638318;0.00465917,0.00296137,0.0049358,0.00369832,0.00345545;0.00465917,0.00360752,0.00337953,0.00232382,0.00217768;0.00465917,0.00385934,0.00299688,0.00165026,0.00141804;0.00465917,0.00279048,0.00167985,0.000868944,0.00062913;0.00344996,0.00198995,0.00127818,0.000602667,0.000490827;0.00243925,0.00173739,0.00100247,0.000502448,0.000400199;0.00195174,0.0016195,0.00117991,0.000653807,0.000522623;0.00196127,0.00143448,0.00090563,0.000593035,0.000515238;0.00210982,0.0013075,0.000808184,0.000489967,0.000453458;0.00266036,0.00121596,0.000742684,0.000416342,0.000378914;0.00344871,0.00125469,0.00074347,0.000401257,0.000350742;0.00389231,0.00126885,0.000696434,0.000379638,0.000337394;0.00392122,0.00188842,0.00105021,0.000512664,0.000344475;0.00377053,0.00200223,0.00129526,0.00051701,0.000287237;0.00349905,0.00192787,0.00134563,0.000629481,0.000378992;0.00314826,0.00181955,0.00130909,0.000850741,0.000526779;0.00285045,0.00175597,0.00125321,0.000856617,0.000718745;0.00255589,0.00161788,0.00119344,0.000887228,0.000640541;0.00227795,0.00154717,0.00113057,0.000748643,0.000532106;0.00204687,0.00150239,0.00108992,0.000741976,0.00046544;0.00179717,0.00147648,0.00104268,0.000650899,0.000408277;0.00166159,0.00146907,0.00104974,0.000648536,0.000425937;0.00154212,0.0014479,0.00100709,0.000584535,0.000365088;0.00142178,0.00147631,0.000971539,0.000574467,0.000359385;0.00138252,0.00149249,0.000966175,0.000567393,0.00034615;0.00135213,0.00150963,0.000946023,0.000509391,0.00032513;0.00161978,0.00157701,0.000973037,0.000512528,0.000321927;0.0018797,0.00176393,0.00107455,0.000552517,0.000360338;0.00233758,0.00204923,0.00128005,0.000575555,0.000345665;0.00303428,0.00251452,0.00155042,0.00101129,0.000740204;0.00397936,0.0023618,0.00194514,0.00167676,0.00124506;0.00494525,0.00268869,0.00239193,0.00146285,0.000987005;0.00575857,0.00294003,0.00244425,0.0011978,0.000867504;0.00596764,0.00278034,0.00213064,0.00106735,0.000733011;0.00569216,0.00241023,0.00177033,0.00095135,0.00066496;0.00494613,0.00199842,0.00153419,0.000885069,0.000648523;0.00403546,0.00176577,0.00135467,0.000898499,0.000675215;0.00326917,0.00161032,0.00130149,0.000889941,0.000633382;0.00274124,0.00153081,0.001227,0.000919164,0.000697074;0.00244927,0.00150333,0.0012413,0.000897664,0.000705495;0.00220336,0.00150596,0.00124922,0.000957534,0.000664272;0.00199516,0.00153147,0.00127479,0.000917949,0.0007466;0.00184088,0.0015754,0.00130798,0.000984517,0.000675919;0.00173782,0.00161791,0.00131662,0.000949267,0.000741;0.0017013,0.00158804,0.00133677,0.000999377,0.000704756;0.00162086,0.00170525,0.00138481,0.000946432,0.000736538;0.00152086,0.00166043,0.00128942,0.00097447,0.000683679;0.00150632,0.00176676,0.00137887,0.00096943,0.000718854;0.00128553,0.00167337,0.0016092,0.000944572,0.000726751;0.00128553,0.00167337,0.00170989,0.00124241,0.000828314]; + tau2_tot = [95,3.2e+02,6.7e+02,7.2e+02,7.2e+02;95,3.2e+02,6.7e+02,7.2e+02,7.2e+02;95,4.4e+02,7.2e+02,7.2e+02,7.2e+02;95,5.3e+02,7e+02,6.9e+02,7.2e+02;95,5.3e+02,6.4e+02,5.8e+02,6.9e+02;1.2e+02,5.4e+02,5.8e+02,4.5e+02,5.5e+02;1.8e+02,5.8e+02,5.8e+02,4.4e+02,5.4e+02;2.5e+02,6.2e+02,6.2e+02,5.6e+02,6e+02;3.7e+02,6e+02,5.7e+02,5.6e+02,6.2e+02;5.3e+02,5.6e+02,5.1e+02,5e+02,6.4e+02;6.7e+02,5.6e+02,4.8e+02,4.3e+02,6e+02;7.2e+02,5.9e+02,4.6e+02,3.8e+02,5.7e+02;7.2e+02,6.4e+02,4.9e+02,3.7e+02,5.2e+02;7.2e+02,6.9e+02,6.4e+02,4.6e+02,4.3e+02;7.2e+02,7.2e+02,6.9e+02,5e+02,3.9e+02;7.2e+02,7e+02,7e+02,6.4e+02,5.5e+02;7.2e+02,7e+02,7.1e+02,7.2e+02,7.2e+02;7.2e+02,7.1e+02,7e+02,6.9e+02,7.1e+02;7.2e+02,7e+02,7e+02,7e+02,6.8e+02;7.1e+02,7e+02,6.9e+02,6.5e+02,6e+02;6.8e+02,7.1e+02,6.8e+02,6.4e+02,5e+02;6.3e+02,7.1e+02,6.6e+02,6.1e+02,4.8e+02;6e+02,7e+02,6.6e+02,5.6e+02,5e+02;5.6e+02,6.9e+02,6.4e+02,5.3e+02,4.3e+02;5.3e+02,6.9e+02,6.1e+02,5.4e+02,4.2e+02;5.2e+02,6.8e+02,6e+02,5.2e+02,4.1e+02;5.7e+02,6.7e+02,5.9e+02,4.6e+02,3.9e+02;7e+02,6.9e+02,6e+02,4.9e+02,4.1e+02;7.2e+02,7.1e+02,6.3e+02,5e+02,4.2e+02;7.2e+02,7.2e+02,6.6e+02,5.9e+02,5.5e+02;7.2e+02,7.2e+02,7e+02,6.3e+02,6.2e+02;7.2e+02,7.2e+02,7.1e+02,7.2e+02,7.2e+02;7.2e+02,7.2e+02,7.1e+02,7.2e+02,7.2e+02;7.2e+02,7.2e+02,7.2e+02,7.2e+02,7e+02;7.2e+02,7.2e+02,7.2e+02,7.2e+02,6.5e+02;7.2e+02,7.2e+02,7.2e+02,6.9e+02,6.1e+02;7.2e+02,7.2e+02,7.1e+02,6.8e+02,6.3e+02;7.2e+02,7e+02,6.8e+02,6.9e+02,6.6e+02;7e+02,6.7e+02,6.8e+02,7e+02,6.4e+02;6.6e+02,6.5e+02,6.7e+02,7.2e+02,7.1e+02;6e+02,6.4e+02,6.7e+02,7.2e+02,6.9e+02;5.6e+02,6.5e+02,6.7e+02,7.2e+02,6.9e+02;5.4e+02,6.6e+02,6.8e+02,7.2e+02,7.2e+02;5.3e+02,6.7e+02,7e+02,7.2e+02,6.5e+02;5.4e+02,6.7e+02,7e+02,7.2e+02,7.1e+02;5.5e+02,6.5e+02,7.1e+02,7.2e+02,6.4e+02;5.3e+02,6.7e+02,7e+02,6.9e+02,6.8e+02;5.2e+02,6.4e+02,6.3e+02,6.8e+02,5.9e+02;4.9e+02,6.2e+02,6.5e+02,6.7e+02,5.8e+02;4.3e+02,5.6e+02,5.6e+02,6e+02,4.9e+02;4.3e+02,5.6e+02,5.6e+02,6e+02,4.9e+02]; + + R0 = R0_tot(:,3)'; + R1 = R1_tot(:,3)'; + tau1 = tau1_tot(:,3)'; + R2 = R2_tot(:,3)'; + tau2 = tau2_tot(:,3)'; + C1 = tau1 ./ R1; + C2 = tau2 ./ R2; + + + p.R0.functionFormat = 'tabulated'; + p.R0.argumentList = {'SOC'}; + p.R0.dataX = p.OCP.dataX; % same axis as SOC + p.R0.dataY = R0; + + p.R1.functionFormat = 'tabulated'; + p.R1.argumentList = {'SOC'}; + p.R1.dataX = p.OCP.dataX; + p.R1.dataY = R1; + + p.C1.functionFormat = 'tabulated'; + p.C1.argumentList = {'SOC'}; + p.C1.dataX = p.OCP.dataX; + p.C1.dataY = C1; + + p.R2.functionFormat = 'tabulated'; + p.R2.argumentList = {'SOC'}; + p.R2.dataX = p.OCP.dataX; + p.R2.dataY = R2; + + p.C2.functionFormat = 'tabulated'; + p.C2.argumentList = {'SOC'}; + p.C2.dataX = p.OCP.dataX; + p.C2.dataY = C2; + +end diff --git a/Battery/EquivalentCircuitModel/utils/plotExperience.m b/Battery/EquivalentCircuitModel/utils/plotExperience.m new file mode 100644 index 000000000..d18f714cc --- /dev/null +++ b/Battery/EquivalentCircuitModel/utils/plotExperience.m @@ -0,0 +1,76 @@ +function experience = plotexperience() + + + json_jp3 = 'C:\Users\Alexandre Fichter\Documents\stage_3A\contenu stage\data_August\jp3-params\jp3-opt-1d-full.json'; + jsonstruct = parseBattmoJson(json_jp3); + + % 1. Construction du modèle + [model, inputparams] = setupModelFromJson(jsonstruct); + + model.NegativeElectrode.Coating.ActiveMaterial.Interface.computeOCPFunc.argumentList = {'c'}; + model.NegativeElectrode.Coating.ActiveMaterial.Interface.computeOCPFunc.numberOfArguments = 1; + + model.PositiveElectrode.Coating.ActiveMaterial.Interface.computeOCPFunc.argumentList = {'c'}; + model.PositiveElectrode.Coating.ActiveMaterial.Interface.computeOCPFunc.numberOfArguments = 1; + + model.Electrolyte.computeDiffusionCoefficientFunc.argumentList = {'c', 'T'}; + model.Electrolyte.computeConductivityFunc.argumentList = {'c', 'T'}; + + % 3. Initialisation de l'état (maintenant le modèle sait quoi envoyer !) + state0 = setupInitialState(model); + + inputparams.Control.controlPolicy = 'CCDischarge'; + inputparams.Control.DRate = 2; + inputparams.Control.tmax = 30; + inputparams.Control.rampupTime = 0.1; + inputparams.Control.lowerCutoffVoltage = 2.5; + + schedule1 = setupSchedule(model, inputparams); + [~, states1] = simulateBattery(model, state0, schedule1); + state_after_pulse = states1{end}; + + inputparams.Control.DRate = 0; + inputparams.Control.tmax = 600; + inputparams.Control.rampupTime = 0.1; + + + schedule2 = setupSchedule(model, inputparams); + disp('Simulation du Repos en cours...'); + [~, states2] = simulateBattery(model, state_after_pulse, schedule2); + + all_states = [states1(:); states2(:)]; + N = length(all_states); + + time = zeros(N, 1); + voltage = zeros(N, 1); + current = zeros(N, 1); + + + for i = 1:N + time(i) = all_states{i}.time; + try + % Dans BattMo, les données du circuit extérieur sont dans "Control" + voltage(i) = all_states{i}.Control.E; + current(i) = all_states{i}.Control.I; + catch + % Alternative si les variables sont stockées à la racine + voltage(i) = all_states{i}.voltage; + current(i) = all_states{i}.current; + end + end + + figure; + subplot(2,1,1); + plot(time, voltage); + xlabel('Time (s)'); + ylabel('Voltage (V)'); + title('Voltage vs Time'); + + subplot(2,1,2); + plot(time, current); + xlabel('Time (s)'); + ylabel('Current (A)'); + title('Current vs Time'); + + experience = [current, voltage, time]; +end diff --git a/Documentation/exampleECM.nblink b/Documentation/exampleECM.nblink new file mode 100644 index 000000000..dafcad3e0 --- /dev/null +++ b/Documentation/exampleECM.nblink @@ -0,0 +1,3 @@ +{ + "path" : "../Battery/EquivalentCircuitModel/notebooks/exampleECM.ipynb" +} diff --git a/Examples/Advanced/Impedance/runImpedanceChen.m b/Examples/Advanced/Impedance/runImpedanceChen.m index 43966c1be..debee4b4b 100644 --- a/Examples/Advanced/Impedance/runImpedanceChen.m +++ b/Examples/Advanced/Impedance/runImpedanceChen.m @@ -1,6 +1,9 @@ % clear all % close all + + + % We define some shorthand names for simplicity. ne = 'NegativeElectrode'; pe = 'PositiveElectrode'; @@ -19,7 +22,7 @@ jsonstruct = mergeJsonStructs({jsonstruct_material, ... jsonstruct_geometry}); -includeDoubleLayer = false; +includeDoubleLayer = true; if includeDoubleLayer @@ -44,23 +47,32 @@ impsolv = ImpedanceSolver(inputparams, options, extrastructs); -%% + + + + + set(0, 'defaultlinelinewidth', 3); -omegas = linspace(-4, 2, 30); -omegas = 10.^omegas; -Z = impsolv.computeImpedance(omegas); + figure + + +omegas = logspace(-4, 2, 30); +Z = impsolv.computeImpedance(omegas); hold on plot(real(Z), -imag(Z), 'displayname', 'battmo'); +axis equal; + + docompare = true; if docompare p = fileparts(mfilename('fullpath')); - data = load(fullfile(p, 'utils', 'pybamm_chen_impedances.mat')); + data = load('C:\Users\Alexandre Fichter\Documents\stage_3A\contenu stage\matlab\BattMo\Examples\Advanced\Impedance\utils\pybamm_chen_impedances.mat'); Zpybamm = data.impedances; plot(real(Zpybamm), -imag(Zpybamm), 'displayname', 'pybamm'); end diff --git a/Tests/TestExamples/TestRunExamples.m b/Tests/TestExamples/TestRunExamples.m index 8cc100900..be4d19c01 100644 --- a/Tests/TestExamples/TestRunExamples.m +++ b/Tests/TestExamples/TestRunExamples.m @@ -15,6 +15,8 @@ 'runChen2020' , ... 'runImpedanceChen' , ... 'runImpedanceScript' , ... + 'exampleECM' , ... + 'runClassFitting' , ... 'runCR' , ... 'runGittTest' , ... 'runJellyRoll' , ... diff --git a/Utilities/FunctionInterface/NamedFunction.m b/Utilities/FunctionInterface/NamedFunction.m index 1918c35ec..2bc75df9f 100644 --- a/Utilities/FunctionInterface/NamedFunction.m +++ b/Utilities/FunctionInterface/NamedFunction.m @@ -20,7 +20,6 @@ fdnames = {'functionName'}; fn = dispatchParams(fn, jsonstruct, fdnames); - fn.functionHandler = str2func(fn.functionName); end