diff --git a/.all-contributorsrc b/.all-contributorsrc index f732a780..4d250cdf 100644 --- a/.all-contributorsrc +++ b/.all-contributorsrc @@ -141,6 +141,15 @@ "test", "code" ] + }, + { + "login": "noelhallemans", + "name": "Noël Hallemans", + "avatar_url": "https://avatars.githubusercontent.com/u/90609505?v=4", + "profile": "https://eng.ox.ac.uk/people/noel-hallemans/", + "contributions": [ + "example" + ] } ], "contributorsPerLine": 7, diff --git a/README.md b/README.md index d7428803..a935f350 100644 --- a/README.md +++ b/README.md @@ -141,6 +141,7 @@ Thanks goes to these wonderful people ([emoji key](https://allcontributors.org/d MarkBlyth
MarkBlyth
💻 f-g-r-i-m-m
f-g-r-i-m-m
💡 Dibyendu-IITKGP
Dibyendu-IITKGP
💡 ⚠️ 💻 + Noël Hallemans
Noël Hallemans
💡 diff --git a/papers/Hallemans et al/Data/Z_SPMegrouped_SOC_chen2020.mat b/papers/Hallemans et al/Data/Z_SPMegrouped_SOC_chen2020.mat new file mode 100644 index 00000000..27aec468 Binary files /dev/null and b/papers/Hallemans et al/Data/Z_SPMegrouped_SOC_chen2020.mat differ diff --git a/papers/Hallemans et al/Data/timeDomainSimulation_SPMegrouped.mat b/papers/Hallemans et al/Data/timeDomainSimulation_SPMegrouped.mat new file mode 100644 index 00000000..2cc08d88 Binary files /dev/null and b/papers/Hallemans et al/Data/timeDomainSimulation_SPMegrouped.mat differ diff --git a/papers/Hallemans et al/Fig11_tauSensitivity.py b/papers/Hallemans et al/Fig11_tauSensitivity.py new file mode 100644 index 00000000..f95874d7 --- /dev/null +++ b/papers/Hallemans et al/Fig11_tauSensitivity.py @@ -0,0 +1,65 @@ +import matplotlib.pyplot as plt +import numpy as np +from scipy.io import savemat + +import pybop +from pybop.models.lithium_ion.basic_SPMe import convert_physical_to_grouped_parameters + +SOC = 0.2 + +factor = 2 +Nparams = 11 +Nfreq = 60 +fmin = 2e-4 +fmax = 1e3 + +# Get grouped parameters +R0 = 0.01 + +parameter_set = pybop.ParameterSet.pybamm("Chen2020") +parameter_set["Electrolyte diffusivity [m2.s-1]"] = 1.769e-10 +parameter_set["Electrolyte conductivity [S.m-1]"] = 1e16 +parameter_set["Negative electrode conductivity [S.m-1]"] = 1e16 +parameter_set["Positive electrode conductivity [S.m-1]"] = 1e16 + +grouped_parameters = convert_physical_to_grouped_parameters(parameter_set) +grouped_parameters["Series resistance [Ohm]"] = R0 +model_options = {"surface form": "differential", "contact resistance": "true"} +var_pts = {"x_n": 100, "x_s": 20, "x_p": 100, "r_n": 100, "r_p": 100} + +## Change parameters +parameter_name = "Negative particle diffusion time scale [s]" +param0 = grouped_parameters[parameter_name] +params = np.logspace(np.log10(param0 / factor), np.log10(param0 * factor), Nparams) + +# Simulate impedance at these parameter values +frequencies = np.logspace(np.log10(fmin), np.log10(fmax), Nfreq) + +impedances = 1j * np.zeros((Nfreq, Nparams)) +for ii, param in enumerate(params): + grouped_parameters[parameter_name] = param + model = pybop.lithium_ion.GroupedSPMe( + parameter_set=grouped_parameters, + eis=True, + options=model_options, + var_pts=var_pts, + ) + model.build( + initial_state={"Initial SoC": SOC}, + ) + simulation = model.simulateEIS(inputs=None, f_eval=frequencies) + impedances[:, ii] = simulation["Impedance"] + +fig, ax = plt.subplots() +for ii in range(Nparams): + ax.plot( + np.real(impedances[:, ii]), + -np.imag(impedances[:, ii]), + ) +ax.set(xlabel=r"$Z_r(\omega)$ [$\Omega$]", ylabel=r"$-Z_j(\omega)$ [$\Omega$]") +ax.grid() +ax.set_aspect("equal", "box") +plt.show() + +mdic = {"Z": impedances, "f": frequencies, "name": parameter_name} +savemat("Data/Z_SPMegrouped_taudn_20.mat", mdic) diff --git a/papers/Hallemans et al/Fig12_EIS_Fitting_Simulation.py b/papers/Hallemans et al/Fig12_EIS_Fitting_Simulation.py new file mode 100644 index 00000000..a4c38fbe --- /dev/null +++ b/papers/Hallemans et al/Fig12_EIS_Fitting_Simulation.py @@ -0,0 +1,303 @@ +import matplotlib.pyplot as plt +import numpy as np +import scipy +from scipy.io import savemat + +import pybop +from pybop.models.lithium_ion.basic_SPMe import convert_physical_to_grouped_parameters + +# To duplicate paper results, modify the below: +n_runs = 1 # 10 +max_iterations = 10 # 1000 + +## Define true parameters +R0 = 0.01 + +parameter_set = pybop.ParameterSet.pybamm("Chen2020") +parameter_set["Electrolyte diffusivity [m2.s-1]"] = 1.769e-10 +parameter_set["Electrolyte conductivity [S.m-1]"] = 1e16 +parameter_set["Negative electrode conductivity [S.m-1]"] = 1e16 +parameter_set["Positive electrode conductivity [S.m-1]"] = 1e16 + +grouped_parameters = convert_physical_to_grouped_parameters(parameter_set) +grouped_parameters["Series resistance [Ohm]"] = R0 +model_options = {"surface form": "differential", "contact resistance": "true"} +var_pts = {"x_n": 100, "x_s": 20, "x_p": 100, "r_n": 100, "r_p": 100} + +## Construct model +model = pybop.lithium_ion.GroupedSPMe( + parameter_set=grouped_parameters, eis=True, var_pts=var_pts, options=model_options +) + +## Parameter bounds for optimisation +tau_d_bounds = [5e2, 1e4] +tau_e_bounds = [2e2, 1e3] +zeta_bounds = [0.5, 1.5] +Qe_bounds = [5e2, 1e3] +tau_ct_bounds = [1e3, 5e4] +C_bounds = [1e-5, 1] +c0p_bounds = [0.8, 0.9] +c0n_bounds = [1e-5, 0.1] +c100p_bounds = [0.2, 0.3] +c100n_bounds = [0.85, 0.95] +t_plus_bounds = [0.2, 0.5] +R0_bounds = [1e-5, 0.05] + + +# Create the parameters object +parameters = pybop.Parameters( + pybop.Parameter( + "Series resistance [Ohm]", + bounds=R0_bounds, + prior=pybop.Uniform(*R0_bounds), + initial_value=np.mean(R0_bounds), + true_value=grouped_parameters["Series resistance [Ohm]"], + transformation=pybop.LogTransformation(), + ), + pybop.Parameter( + "Positive particle diffusion time scale [s]", + bounds=tau_d_bounds, + initial_value=np.mean(tau_d_bounds), + prior=pybop.Uniform(*tau_d_bounds), + true_value=grouped_parameters["Positive particle diffusion time scale [s]"], + transformation=pybop.LogTransformation(), + ), + pybop.Parameter( + "Negative particle diffusion time scale [s]", + bounds=tau_d_bounds, + initial_value=np.mean(tau_d_bounds), + prior=pybop.Uniform(*tau_d_bounds), + true_value=grouped_parameters["Negative particle diffusion time scale [s]"], + transformation=pybop.LogTransformation(), + ), + pybop.Parameter( + "Cation transference number", + bounds=t_plus_bounds, + initial_value=np.mean(t_plus_bounds), + prior=pybop.Uniform(*t_plus_bounds), + true_value=grouped_parameters["Cation transference number"], + transformation=pybop.LogTransformation(), + ), + pybop.Parameter( + "Positive electrode electrolyte diffusion time scale [s]", + bounds=tau_e_bounds, + initial_value=np.mean(tau_e_bounds), + prior=pybop.Uniform(*tau_e_bounds), + true_value=grouped_parameters[ + "Positive electrode electrolyte diffusion time scale [s]" + ], + transformation=pybop.LogTransformation(), + ), + pybop.Parameter( + "Negative electrode electrolyte diffusion time scale [s]", + bounds=tau_e_bounds, + initial_value=np.mean(tau_e_bounds), + prior=pybop.Uniform(*tau_e_bounds), + true_value=grouped_parameters[ + "Negative electrode electrolyte diffusion time scale [s]" + ], + transformation=pybop.LogTransformation(), + ), + pybop.Parameter( + "Separator electrolyte diffusion time scale [s]", + bounds=tau_e_bounds, + initial_value=np.mean(tau_e_bounds), + prior=pybop.Uniform(*tau_e_bounds), + true_value=grouped_parameters["Separator electrolyte diffusion time scale [s]"], + transformation=pybop.LogTransformation(), + ), + pybop.Parameter( + "Positive electrode charge transfer time scale [s]", + bounds=tau_ct_bounds, + initial_value=np.mean(tau_ct_bounds), + prior=pybop.Uniform(*tau_ct_bounds), + true_value=grouped_parameters[ + "Positive electrode charge transfer time scale [s]" + ], + transformation=pybop.LogTransformation(), + ), + pybop.Parameter( + "Negative electrode charge transfer time scale [s]", + bounds=tau_ct_bounds, + initial_value=np.mean(tau_ct_bounds), + prior=pybop.Uniform(*tau_ct_bounds), + true_value=grouped_parameters[ + "Negative electrode charge transfer time scale [s]" + ], + transformation=pybop.LogTransformation(), + ), + pybop.Parameter( + "Positive electrode capacitance [F]", + bounds=C_bounds, + initial_value=np.mean(C_bounds), + prior=pybop.Uniform(*C_bounds), + true_value=grouped_parameters["Positive electrode capacitance [F]"], + transformation=pybop.LogTransformation(), + ), + pybop.Parameter( + "Negative electrode capacitance [F]", + bounds=C_bounds, + initial_value=np.mean(C_bounds), + prior=pybop.Uniform(*C_bounds), + true_value=grouped_parameters["Negative electrode capacitance [F]"], + transformation=pybop.LogTransformation(), + ), + pybop.Parameter( + "Positive electrode relative porosity", + bounds=zeta_bounds, + initial_value=np.mean(zeta_bounds), + prior=pybop.Uniform(*zeta_bounds), + true_value=grouped_parameters["Positive electrode relative porosity"], + transformation=pybop.LogTransformation(), + ), + pybop.Parameter( + "Negative electrode relative porosity", + bounds=zeta_bounds, + initial_value=np.mean(zeta_bounds), + prior=pybop.Uniform(*zeta_bounds), + true_value=grouped_parameters["Negative electrode relative porosity"], + transformation=pybop.LogTransformation(), + ), + pybop.Parameter( + "Reference electrolyte capacity [A.s]", + bounds=Qe_bounds, + initial_value=np.mean(Qe_bounds), + prior=pybop.Uniform(*Qe_bounds), + true_value=grouped_parameters["Reference electrolyte capacity [A.s]"], + transformation=pybop.LogTransformation(), + ), + pybop.Parameter( + "Minimum positive stoichiometry", + bounds=c100p_bounds, + initial_value=np.mean(c100p_bounds), + prior=pybop.Uniform(*c100p_bounds), + true_value=grouped_parameters["Minimum positive stoichiometry"], + transformation=pybop.LogTransformation(), + ), + pybop.Parameter( + "Maximum positive stoichiometry", + bounds=c0p_bounds, + initial_value=np.mean(c0p_bounds), + prior=pybop.Uniform(*c0p_bounds), + true_value=grouped_parameters["Maximum positive stoichiometry"], + transformation=pybop.LogTransformation(), + ), + pybop.Parameter( + "Minimum negative stoichiometry", + bounds=c0n_bounds, + initial_value=np.mean(c0n_bounds), + prior=pybop.Uniform(*c0n_bounds), + true_value=grouped_parameters["Minimum negative stoichiometry"], + transformation=pybop.LogTransformation(), + ), + pybop.Parameter( + "Maximum negative stoichiometry", + bounds=c100n_bounds, + initial_value=np.mean(c100n_bounds), + prior=pybop.Uniform(*c100n_bounds), + true_value=grouped_parameters["Maximum negative stoichiometry"], + transformation=pybop.LogTransformation(), + ), +) + + +## Read simulated impedance data +EIS_data = scipy.io.loadmat("Data/Z_SPMegrouped_SOC_chen2020.mat") + +impedances = EIS_data.get("Z") +frequencies = EIS_data.get("f") +frequencies = frequencies.flatten() +SOCs = EIS_data.get("SOC") +SOCs = SOCs.flatten() + +# Form dataset +signal = ["Impedance"] +datasets = [None] * len(SOCs) +models = [None] * len(SOCs) +problems = [None] * len(SOCs) + +for ii in range(len(SOCs)): + datasets[ii] = pybop.Dataset( + { + "Frequency [Hz]": frequencies, + "Current function [A]": np.ones(len(frequencies)) * 0.0, + "Impedance": impedances[:, ii], + } + ) + models[ii] = model.new_copy() + models[ii].build(initial_state={"Initial SoC": SOCs[ii]}) + problems[ii] = pybop.FittingProblem( + models[ii], + parameters, + datasets[ii], + signal=signal, + ) + +problem = pybop.MultiFittingProblem(*problems) +cost = pybop.SumSquaredError(problem) +optim = pybop.PSO( + cost, + parallel=True, + multistart=n_runs, + max_iterations=max_iterations, + max_unchanged_iterations=max_iterations, + # polish=False, # For SciPyDifferential + # popsize=5, # For SciPyDifferential +) + +results = optim.run() + +# Print optimised parameters +print("True grouped parameters", parameters.true_value()) +grouped_parameters.update(parameters.as_dict(results.best_x)) + +# Plot convergence +pybop.plot.convergence(optim) + +# Plot the parameter traces +pybop.plot.parameters(optim) + +## Save estimated impedance data +model_hat = pybop.lithium_ion.GroupedSPMe( + parameter_set=grouped_parameters, + eis=True, + var_pts=var_pts, + options=model_options, +) + +Nfreq = len(frequencies) +NSOC = len(SOCs) + +# Compute impedance for estimated parameters +impedances_hat = 1j * np.zeros((Nfreq, NSOC)) +for ii in range(len(SOCs)): + model_hat.build(initial_state={"Initial SoC": SOCs[ii]}) + simulation = model_hat.simulateEIS(inputs=None, f_eval=frequencies) + impedances_hat[:, ii] = simulation["Impedance"] + +colors = plt.get_cmap("tab10") +fig, ax = plt.subplots() +for ii in range(len(SOCs)): + ax.plot(np.real(impedances[:, ii]), -np.imag(impedances[:, ii]), color=colors(ii)) + ax.plot( + np.real(impedances_hat[:, ii]), + -np.imag(impedances_hat[:, ii]), + color=colors(ii), + ) +ax.set(xlabel=r"$Z_r(\omega)$ [$\Omega$]", ylabel=r"$-Z_j(\omega)$ [$\Omega$]") +ax.grid() +ax.set_aspect("equal", "box") +plt.show() + +mdic = { + "Z": impedances, + "Zhat": impedances_hat, + "f": frequencies, + "SOC": SOCs, + "final_cost": results.final_cost, + "theta": parameters.true_value(), + "thetahat": results.x, + "thetahatbest": results.best_x, + "computationTime": results.time, +} +savemat("Data/Zhat_SOC_SPMe_Simulation.mat", mdic) diff --git a/papers/Hallemans et al/Fig2_timeDomainSimulation.py b/papers/Hallemans et al/Fig2_timeDomainSimulation.py new file mode 100644 index 00000000..fec25409 --- /dev/null +++ b/papers/Hallemans et al/Fig2_timeDomainSimulation.py @@ -0,0 +1,70 @@ +import matplotlib.pyplot as plt +import numpy as np +from scipy.io import savemat + +import pybop +from pybop.models.lithium_ion.basic_SPMe import convert_physical_to_grouped_parameters + +## Grouped parameter set +R0 = 0.01 + +parameter_set = pybop.ParameterSet.pybamm("Chen2020") +parameter_set["Electrolyte diffusivity [m2.s-1]"] = 1.769e-10 +parameter_set["Electrolyte conductivity [S.m-1]"] = 1e16 +parameter_set["Negative electrode conductivity [S.m-1]"] = 1e16 +parameter_set["Positive electrode conductivity [S.m-1]"] = 1e16 + +grouped_parameters = convert_physical_to_grouped_parameters(parameter_set) +grouped_parameters["Series resistance [Ohm]"] = R0 +model_options = {"surface form": "differential", "contact resistance": "true"} +var_pts = {"x_n": 100, "x_s": 20, "x_p": 100, "r_n": 100, "r_p": 100} + +## Create model +model = pybop.lithium_ion.GroupedSPMe( + parameter_set=grouped_parameters, eis=True, var_pts=var_pts, options=model_options +) + +## Test model in the time domain +SOC0 = 0.9 +model.build(initial_state={"Initial SoC": SOC0}) + +Ts = 10 # Sampling period +T = 2 * 60 * 60 # 2h +N = int(T / Ts) +time = np.linspace(0, T - Ts, N) + +i_relax7 = np.zeros([7 * int(60 / Ts)]) # 7 min +i_relax20 = np.zeros([20 * int(60 / Ts)]) # 20 min +i_discharge = 5 * np.ones([53 * int(60 / Ts)]) # 53 min +i_charge = -5 * np.ones([20 * int(60 / Ts)]) # 20 min +current = np.concatenate((i_relax7, i_discharge, i_relax20, i_charge, i_relax20)) + +experiment = pybop.Dataset( + { + "Time [s]": time, + "Current function [A]": current, + } +) +model.set_current_function(dataset=experiment) +simulation = model.predict(t_eval=time) + +# Plot traces +fig, ax = plt.subplots() +ax.plot(simulation["Time [s]"].data, simulation["Current [A]"].data) +ax.set(xlabel="time [s]", ylabel="Current [A]") +ax.grid() +plt.show() + +fig, ax = plt.subplots() +ax.plot(simulation["Time [s]"].data, simulation["Voltage [V]"].data) +ax.set(xlabel="time [s]", ylabel="Voltage [V]") +ax.grid() +plt.show() + +## Save data +t = simulation["Time [s]"].data +i = simulation["Current [A]"].data +v = simulation["Voltage [V]"].data + +mdic = {"t": t, "i": i, "v": v, "SOC0": SOC0} +savemat("Data/timeDomainSimulation_SPMegrouped.mat", mdic) diff --git a/papers/Hallemans et al/Fig4_comparisonBruteForceFreqDomain.py b/papers/Hallemans et al/Fig4_comparisonBruteForceFreqDomain.py new file mode 100644 index 00000000..75317265 --- /dev/null +++ b/papers/Hallemans et al/Fig4_comparisonBruteForceFreqDomain.py @@ -0,0 +1,127 @@ +import time as timer + +import numpy as np +import pybamm + +import pybop + +## Fixed parameters +n_runs = 1 + +SOC = 0.5 +Nfreq = 60 +fmin = 2e-4 +fmax = 1e3 +frequencies = np.logspace(np.log10(fmin), np.log10(fmax), Nfreq) + +parameter_set = pybop.ParameterSet.pybamm("Chen2020") +parameter_set["Contact resistance [Ohm]"] = 0.01 +model_options = {"surface form": "differential", "contact resistance": "true"} + + +## Time domain simulation +# Time domain +I_app = 100e-3 +number_of_periods = 10 +samples_per_period = 56 + + +def current_function(t): + return I_app * np.sin(2 * np.pi * pybamm.Parameter("Frequency [Hz]") * t) + + +parameter_set.update( + { + "Current function [A]": current_function, + "Frequency [Hz]": 10, + }, + check_already_exists=False, +) + +model = pybop.lithium_ion.DFN( + parameter_set=parameter_set, + options=model_options, +) +var_pts = model.pybamm_model.default_var_pts +var_pts["x_n"] = 100 +var_pts["x_p"] = 100 +var_pts["r_n"] = 100 +var_pts["r_p"] = 100 + +model_time = pybop.lithium_ion.DFN( + parameter_set=parameter_set, + solver=pybamm.ScipySolver(atol=1e-9), + options=model_options, + var_pts=var_pts, +) +model_freq = pybop.lithium_ion.DFN( + parameter_set=parameter_set, + options=model_options, + var_pts=var_pts, + eis=True, +) +parameters = pybop.Parameters(pybop.Parameter("Frequency [Hz]", initial_value=10)) +model_time.build(parameters=parameters, initial_state={"Initial SoC": SOC}) +model_freq.build(initial_state={"Initial SoC": SOC}) + +## Time domain simulation +# time_elapsed = np.zeros([Nruns, 1]) +# for ii in range(Nruns): +# start_time = timer.time() + +# impedance_time = [] +# for frequency in frequencies: +# # Solve +# period = 1 / frequency +# dt = period / samples_per_period +# t_eval = np.array(range(0, samples_per_period * number_of_periods)) * dt +# sol = model_time.simulate( +# inputs={"Frequency [Hz]": frequency}, +# t_eval=t_eval, +# ) +# # Extract final P periods of the solution +# P = 5 +# time = sol["Time [s]"].entries[-P * samples_per_period :] +# current = sol["Current [A]"].entries[-P * samples_per_period :] +# voltage = sol["Voltage [V]"].entries[-P * samples_per_period :] + +# # FFT +# current_fft = np.fft.fft(current) +# voltage_fft = np.fft.fft(voltage) + +# # Get index of first harmonic +# impedance = -voltage_fft[P] / current_fft[P] +# impedance_time.append(impedance) + +# end_time = timer.time() +# time_elapsed[ii] = end_time - start_time +# print("Time domain method: ", np.mean(time_elapsed), "s") + +# plt.plot(np.real(impedance_time), -np.imag(impedance_time)) +# plt.show() + +## Frequency domain simulation +time_elapsed = np.zeros([n_runs, 1]) +for ii in range(n_runs): + start_time = timer.time() + simulationFD = model_freq.simulateEIS( + inputs=None, + f_eval=frequencies, + ) + impedance_freq = simulationFD["Impedance"] + + end_time = timer.time() + time_elapsed[ii] = end_time - start_time +print("Frequency domain method: ", np.mean(time_elapsed), "s") + +# plt.plot(np.real(impedance_freq), -np.imag(impedance_freq)) +# plt.show() + + +# Write to matfile +# mdic = { +# "Ztime": impedance_time, +# "Zfreq": impedance_freq, +# "f": frequencies, +# } +# savemat("Data comparison time freq/Z_DFN_timeFreq.mat", mdic) diff --git a/papers/Hallemans et al/Fig5_impedance_models.py b/papers/Hallemans et al/Fig5_impedance_models.py new file mode 100644 index 00000000..d329a9bd --- /dev/null +++ b/papers/Hallemans et al/Fig5_impedance_models.py @@ -0,0 +1,62 @@ +import matplotlib.pyplot as plt +import numpy as np +from scipy.io import savemat + +import pybop + +Nfreq = 60 +SOC = 0.5 +fmin = 2e-4 +fmax = 1e3 + +frequencies = np.logspace(np.log10(fmin), np.log10(fmax), Nfreq) +impedances = 1j * np.zeros((Nfreq, 3)) + +# Define model +parameter_set = pybop.ParameterSet.pybamm("Chen2020") +parameter_set["Contact resistance [Ohm]"] = 0.01 + +model_options = {"surface form": "differential", "contact resistance": "true"} +var_pts = {"x_n": 100, "x_s": 20, "x_p": 100, "r_n": 100, "r_p": 100} + + +## SPM +model = pybop.lithium_ion.SPM( + parameter_set=parameter_set, options=model_options, eis=True, var_pts=var_pts +) +model.build(initial_state={"Initial SoC": SOC}) + +simulation = model.simulateEIS(inputs=None, f_eval=frequencies) +impedances[:, 0] = simulation["Impedance"] + +## SPMe +model = pybop.lithium_ion.SPMe( + parameter_set=parameter_set, options=model_options, eis=True, var_pts=var_pts +) +model.build(initial_state={"Initial SoC": SOC}) +simulation = model.simulateEIS(inputs=None, f_eval=frequencies) +impedances[:, 1] = simulation["Impedance"] + +## DFN +model = pybop.lithium_ion.DFN( + parameter_set=parameter_set, options=model_options, eis=True, var_pts=var_pts +) +model.build(initial_state={"Initial SoC": SOC}) +simulation = model.simulateEIS(inputs=None, f_eval=frequencies) +impedances[:, 2] = simulation["Impedance"] + +## Plot +fig, ax = plt.subplots() +for ii in range(3): + ax.plot( + np.real(impedances[:, ii]), + -np.imag(impedances[:, ii]), + ) +ax.set(xlabel=r"$Z_r(\omega)$ [$\Omega$]", ylabel=r"$-Z_j(\omega)$ [$\Omega$]") +ax.grid() +ax.set_aspect("equal", "box") +plt.show() + +## Save +mdic = {"Z": impedances, "f": frequencies} +savemat("Data/Z_SPM_SPMe_DFN_Pybop_chen2020.mat", mdic) diff --git a/papers/Hallemans et al/Fig7_impedance_SOC_SPMegrouped.py b/papers/Hallemans et al/Fig7_impedance_SOC_SPMegrouped.py new file mode 100644 index 00000000..3785f1fc --- /dev/null +++ b/papers/Hallemans et al/Fig7_impedance_SOC_SPMegrouped.py @@ -0,0 +1,53 @@ +import matplotlib.pyplot as plt +import numpy as np + +import pybop +from pybop.models.lithium_ion.basic_SPMe import convert_physical_to_grouped_parameters + +## Group parameter set +R0 = 0.01 + +parameter_set = pybop.ParameterSet.pybamm("Chen2020") +parameter_set["Electrolyte diffusivity [m2.s-1]"] = 1.769e-10 +parameter_set["Electrolyte conductivity [S.m-1]"] = 1e16 +parameter_set["Negative electrode conductivity [S.m-1]"] = 1e16 +parameter_set["Positive electrode conductivity [S.m-1]"] = 1e16 + +grouped_parameters = convert_physical_to_grouped_parameters(parameter_set) +grouped_parameters["Series resistance [Ohm]"] = R0 + +## Create model +model_options = {"surface form": "differential", "contact resistance": "true"} +var_pts = {"x_n": 100, "x_s": 20, "x_p": 100, "r_n": 100, "r_p": 100} +model = pybop.lithium_ion.GroupedSPMe( + parameter_set=grouped_parameters, eis=True, var_pts=var_pts, options=model_options +) + +## Simulate impedance +Nfreq = 60 +fmin = 2e-4 +fmax = 1e3 +NSOC = 9 +frequencies = np.logspace(np.log10(fmin), np.log10(fmax), Nfreq) +SOCs = np.linspace(0.1, 0.9, NSOC) + +impedances = 1j * np.zeros((Nfreq, NSOC)) +for ii, SOC in enumerate(SOCs): + model.build(initial_state={"Initial SoC": SOC}) + simulation = model.simulateEIS(inputs=None, f_eval=frequencies) + impedances[:, ii] = simulation["Impedance"] + +fig, ax = plt.subplots() +for ii in range(len(SOCs)): + ax.plot( + np.real(impedances[:, ii]), + -np.imag(impedances[:, ii]), + ) +ax.set(xlabel=r"$Z_r(\omega)$ [$\Omega$]", ylabel=r"$-Z_j(\omega)$ [$\Omega$]") +ax.grid() +ax.set_aspect("equal", "box") +plt.show() + +## Save data +# mdic = {"Z": impedances, "f": frequencies, "SOC": SOCs} +# savemat("Data/Z_SPMegrouped_SOC_chen2020.mat", mdic) diff --git a/papers/Hallemans et al/Fig8_groupedParamSensitivity.py b/papers/Hallemans et al/Fig8_groupedParamSensitivity.py new file mode 100644 index 00000000..1e8f36e9 --- /dev/null +++ b/papers/Hallemans et al/Fig8_groupedParamSensitivity.py @@ -0,0 +1,81 @@ +import matplotlib.pyplot as plt +import numpy as np +from scipy.io import savemat + +import pybop +from pybop.models.lithium_ion.basic_SPMe import convert_physical_to_grouped_parameters + +# +factor = 2 +Nparams = 11 +SOC = 0.5 +Nfreq = 60 +fmin = 2e-4 +fmax = 1e3 + +# Get grouped parameters +R0 = 0.01 + +parameter_set = pybop.ParameterSet.pybamm("Chen2020") +parameter_set["Electrolyte diffusivity [m2.s-1]"] = 1.769e-10 +parameter_set["Electrolyte conductivity [S.m-1]"] = 1e16 +parameter_set["Negative electrode conductivity [S.m-1]"] = 1e16 +parameter_set["Positive electrode conductivity [S.m-1]"] = 1e16 + +grouped_parameters = convert_physical_to_grouped_parameters(parameter_set) +grouped_parameters["Series resistance [Ohm]"] = R0 +model_options = {"surface form": "differential", "contact resistance": "true"} + +var_pts = {"x_n": 100, "x_s": 20, "x_p": 100, "r_n": 100, "r_p": 100} + +## Change parameters +parameter_name = "Negative electrode relative porosity" + +# "Positive particle diffusion time scale [s]" +# "Positive electrode electrolyte diffusion time scale [s]" +# "Separator electrolyte diffusion time scale [s]" +# "Positive electrode charge transfer time scale [s]" +# "Series resistance [Ohm]" +# "Positive electrode relative porosity" +# "Cation transference number" +# "Reference electrolyte capacity [A.s]" +# "Positive electrode capacitance [F]" +# "Positive theoretical electrode capacity [As]" +# "Positive electrode relative thickness" +# "Measured cell capacity [A.s]" + +param0 = grouped_parameters[parameter_name] + +params = np.logspace(np.log10(param0 / factor), np.log10(param0 * factor), Nparams) + +# Simulate impedance at these parameter values +frequencies = np.logspace(np.log10(fmin), np.log10(fmax), Nfreq) + +impedances = 1j * np.zeros((Nfreq, Nparams)) +for ii, param in enumerate(params): + grouped_parameters[parameter_name] = param + model = pybop.lithium_ion.GroupedSPMe( + parameter_set=grouped_parameters, + eis=True, + options=model_options, + var_pts=var_pts, + ) + model.build( + initial_state={"Initial SoC": SOC}, + ) + simulation = model.simulateEIS(inputs=None, f_eval=frequencies) + impedances[:, ii] = simulation["Impedance"] + +fig, ax = plt.subplots() +for ii in range(Nparams): + ax.plot( + np.real(impedances[:, ii]), + -np.imag(impedances[:, ii]), + ) +ax.set(xlabel=r"$Z_r(\omega)$ [$\Omega$]", ylabel=r"$-Z_j(\omega)$ [$\Omega$]") +ax.grid() +ax.set_aspect("equal", "box") +plt.show() + +mdic = {"Z": impedances, "f": frequencies, "name": parameter_name} +savemat("Data/Sensitivity/Z_SPMegrouped_zetan.mat", mdic) diff --git a/papers/Hallemans et al/README.md b/papers/Hallemans et al/README.md new file mode 100644 index 00000000..0af8c239 --- /dev/null +++ b/papers/Hallemans et al/README.md @@ -0,0 +1,5 @@ +# Physics-based battery model parametrisation from impedance data + +This directory contains the code used for [Hallemans et al.](http://arxiv.org/abs/2412.10896) article. + +Note: These scripts were ran using `PyBOP` v24.12, and are not maintained for future PyBOP releases. diff --git a/papers/Hallemans et al/timeDomain_Fitting_Simulation.py b/papers/Hallemans et al/timeDomain_Fitting_Simulation.py new file mode 100644 index 00000000..f479a771 --- /dev/null +++ b/papers/Hallemans et al/timeDomain_Fitting_Simulation.py @@ -0,0 +1,270 @@ +import numpy as np +import scipy +from scipy.io import savemat + +import pybop +from pybop.models.lithium_ion.basic_SPMe import convert_physical_to_grouped_parameters + +# To duplicate paper results, modify the below: +n_runs = 1 # 10 +max_iterations = 10 # 1000 + +## Grouped parameter set +R0 = 0.01 + +parameter_set = pybop.ParameterSet.pybamm("Chen2020") +parameter_set["Electrolyte diffusivity [m2.s-1]"] = 1.769e-10 +parameter_set["Electrolyte conductivity [S.m-1]"] = 1e16 +parameter_set["Negative electrode conductivity [S.m-1]"] = 1e16 +parameter_set["Positive electrode conductivity [S.m-1]"] = 1e16 + +grouped_parameters = convert_physical_to_grouped_parameters(parameter_set) +grouped_parameters["Series resistance [Ohm]"] = R0 +model_options = {"surface form": "differential", "contact resistance": "true"} +var_pts = {"x_n": 100, "x_s": 20, "x_p": 100, "r_n": 100, "r_p": 100} + +## Construct model +model = pybop.lithium_ion.GroupedSPMe( + parameter_set=grouped_parameters, var_pts=var_pts, options=model_options +) + +## Parameter bounds for optimisation +R0_bounds = [1e-5, 0.05] +tau_d_bounds = [5e2, 1e4] +t_plus_bounds = [0.2, 0.5] +tau_e_bounds = [2e2, 1e3] +tau_ct_bounds = [1e3, 5e4] +C_bounds = [1e-5, 1] +zeta_bounds = [0.5, 1.5] +Qe_bounds = [5e2, 1e3] +# c0p_bounds = [0.8, 0.9] +# c0n_bounds = [1e-5, 0.1] +# c100p_bounds = [0.2, 0.3] +# c100n_bounds = [0.85, 0.95] + + +# Create the parameters object +parameters = pybop.Parameters( + pybop.Parameter( + "Series resistance [Ohm]", + bounds=R0_bounds, + prior=pybop.Uniform(*R0_bounds), + initial_value=np.mean(R0_bounds), + true_value=grouped_parameters["Series resistance [Ohm]"], + transformation=pybop.LogTransformation(), + ), + pybop.Parameter( + "Positive particle diffusion time scale [s]", + bounds=tau_d_bounds, + initial_value=np.mean(tau_d_bounds), + prior=pybop.Uniform(*tau_d_bounds), + true_value=grouped_parameters["Positive particle diffusion time scale [s]"], + transformation=pybop.LogTransformation(), + ), + pybop.Parameter( + "Negative particle diffusion time scale [s]", + bounds=tau_d_bounds, + initial_value=np.mean(tau_d_bounds), + prior=pybop.Uniform(*tau_d_bounds), + true_value=grouped_parameters["Negative particle diffusion time scale [s]"], + transformation=pybop.LogTransformation(), + ), + pybop.Parameter( + "Cation transference number", + bounds=t_plus_bounds, + initial_value=np.mean(t_plus_bounds), + prior=pybop.Uniform(*t_plus_bounds), + true_value=grouped_parameters["Cation transference number"], + transformation=pybop.LogTransformation(), + ), + pybop.Parameter( + "Positive electrode electrolyte diffusion time scale [s]", + bounds=tau_e_bounds, + initial_value=np.mean(tau_e_bounds), + prior=pybop.Uniform(*tau_e_bounds), + true_value=grouped_parameters[ + "Positive electrode electrolyte diffusion time scale [s]" + ], + transformation=pybop.LogTransformation(), + ), + pybop.Parameter( + "Negative electrode electrolyte diffusion time scale [s]", + bounds=tau_e_bounds, + initial_value=np.mean(tau_e_bounds), + prior=pybop.Uniform(*tau_e_bounds), + true_value=grouped_parameters[ + "Negative electrode electrolyte diffusion time scale [s]" + ], + transformation=pybop.LogTransformation(), + ), + pybop.Parameter( + "Separator electrolyte diffusion time scale [s]", + bounds=tau_e_bounds, + initial_value=np.mean(tau_e_bounds), + prior=pybop.Uniform(*tau_e_bounds), + true_value=grouped_parameters["Separator electrolyte diffusion time scale [s]"], + transformation=pybop.LogTransformation(), + ), + pybop.Parameter( + "Positive electrode charge transfer time scale [s]", + bounds=tau_ct_bounds, + initial_value=np.mean(tau_ct_bounds), + prior=pybop.Uniform(*tau_ct_bounds), + true_value=grouped_parameters[ + "Positive electrode charge transfer time scale [s]" + ], + transformation=pybop.LogTransformation(), + ), + pybop.Parameter( + "Negative electrode charge transfer time scale [s]", + bounds=tau_ct_bounds, + initial_value=np.mean(tau_ct_bounds), + prior=pybop.Uniform(*tau_ct_bounds), + true_value=grouped_parameters[ + "Negative electrode charge transfer time scale [s]" + ], + transformation=pybop.LogTransformation(), + ), + pybop.Parameter( + "Positive electrode capacitance [F]", + bounds=C_bounds, + initial_value=np.mean(C_bounds), + prior=pybop.Uniform(*C_bounds), + true_value=grouped_parameters["Positive electrode capacitance [F]"], + transformation=pybop.LogTransformation(), + ), + pybop.Parameter( + "Negative electrode capacitance [F]", + bounds=C_bounds, + initial_value=np.mean(C_bounds), + prior=pybop.Uniform(*C_bounds), + true_value=grouped_parameters["Negative electrode capacitance [F]"], + transformation=pybop.LogTransformation(), + ), + pybop.Parameter( + "Positive electrode relative porosity", + bounds=zeta_bounds, + initial_value=np.mean(zeta_bounds), + prior=pybop.Uniform(*zeta_bounds), + true_value=grouped_parameters["Positive electrode relative porosity"], + transformation=pybop.LogTransformation(), + ), + pybop.Parameter( + "Negative electrode relative porosity", + bounds=zeta_bounds, + initial_value=np.mean(zeta_bounds), + prior=pybop.Uniform(*zeta_bounds), + true_value=grouped_parameters["Negative electrode relative porosity"], + transformation=pybop.LogTransformation(), + ), + pybop.Parameter( + "Reference electrolyte capacity [A.s]", + bounds=Qe_bounds, + initial_value=np.mean(Qe_bounds), + prior=pybop.Uniform(*Qe_bounds), + true_value=grouped_parameters["Reference electrolyte capacity [A.s]"], + transformation=pybop.LogTransformation(), + ), + # pybop.Parameter( + # "Minimum positive stoichiometry", + # bounds=c100p_bounds, + # initial_value=np.mean(c100p_bounds), + # prior=pybop.Uniform(*c100p_bounds), + # true_value=grouped_parameters["Minimum positive stoichiometry"], + # transformation=pybop.LogTransformation(), + # ), + # pybop.Parameter( + # "Maximum positive stoichiometry", + # bounds=c0p_bounds, + # initial_value=np.mean(c0p_bounds), + # prior=pybop.Uniform(*c0p_bounds), + # true_value=grouped_parameters["Maximum positive stoichiometry"], + # transformation=pybop.LogTransformation(), + # ), + # pybop.Parameter( + # "Minimum negative stoichiometry", + # bounds=c0n_bounds, + # initial_value=np.mean(c0n_bounds), + # prior=pybop.Uniform(*c0n_bounds), + # true_value=grouped_parameters["Minimum negative stoichiometry"], + # transformation=pybop.LogTransformation(), + # ), + # pybop.Parameter( + # "Maximum negative stoichiometry", + # bounds=c100n_bounds, + # initial_value=np.mean(c100n_bounds), + # prior=pybop.Uniform(*c100n_bounds), + # true_value=grouped_parameters["Maximum negative stoichiometry"], + # transformation=pybop.LogTransformation(), + # ), +) + +## Read simulated time domain data +timeDomainData = scipy.io.loadmat("Data/timeDomainSimulation_SPMegrouped.mat") +SOC0 = timeDomainData.get("SOC0") +t = timeDomainData.get("t").flatten() +i = timeDomainData.get("i").flatten() +v = timeDomainData.get("v").flatten() + +SOC0 = SOC0.flatten() +SOC0 = SOC0[0] + +model.build(initial_state={"Initial SoC": SOC0}) +dataset = pybop.Dataset( + { + "Time [s]": t, + "Current function [A]": i, + "Voltage [V]": v, + } +) +problem = pybop.FittingProblem(model, parameters, dataset) +cost = pybop.SumSquaredError(problem) +optim = pybop.PSO( + cost, + parallel=True, + multistart=n_runs, + max_iterations=max_iterations, + max_unchanged_iterations=max_iterations, + # polish=False, # For SciPyDifferential + # popsize=5, # For SciPyDifferential +) + +results = optim.run() + +# Print optimised parameters +print("True grouped parameters", parameters.true_value()) +grouped_parameters.update(parameters.as_dict(results.best_x)) + +# Plot traces +pybop.plot.quick(problem, results.best_x) + +# Plot convergence +pybop.plot.convergence(optim) + +# Plot the parameter traces +pybop.plot.parameters(optim) + +## Save estimated voltage +model_hat = pybop.lithium_ion.GroupedSPMe( + parameter_set=grouped_parameters, var_pts=var_pts, options=model_options +) + +# Grouped SPMe +model_hat.build(initial_state={"Initial SoC": SOC0}) +model_hat.set_current_function(dataset=dataset) +simulation_hat = model_hat.predict(t_eval=dataset["Time [s]"]) + +## Save data +mdic = { + "t": simulation_hat["Time [s]"].data, + "i": simulation_hat["Current [A]"].data, + "v": v, + "vhat": simulation_hat["Voltage [V]"].data, + "SOC0": SOC0, + "final_cost": results.final_cost, + "theta": parameters.true_value(), + "thetahat": results.x, + "thetahatbest": results.best_x, + "computationTime": results.time, +} +savemat("Data/Estimate_timeDomainSimulation.mat", mdic) diff --git a/tests/integration/test_spm_parameterisations.py b/tests/integration/test_spm_parameterisations.py index 861c97fc..fa7729ad 100644 --- a/tests/integration/test_spm_parameterisations.py +++ b/tests/integration/test_spm_parameterisations.py @@ -230,11 +230,10 @@ def test_multiple_signals(self, multi_optimiser, spm_two_signal_cost): else 0.02, } - if ( - isinstance(spm_two_signal_cost, pybop.BaseJaxCost) - and multi_optimiser is pybop.SciPyDifferentialEvolution - ): + if multi_optimiser is pybop.SciPyDifferentialEvolution: common_args["bounds"] = [[0.375, 0.775], [0.375, 0.775]] + if isinstance(spm_two_signal_cost, pybop.GaussianLogLikelihood): + common_args["bounds"].extend([[0.0, 0.1], [0.0, 0.1]]) # Test each optimiser optim = multi_optimiser(**common_args)