From dfd96d98fd986e1bd8b4ec7237fa66e8891ada3d Mon Sep 17 00:00:00 2001 From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com> Date: Thu, 6 Aug 2026 17:44:47 +0100 Subject: [PATCH 1/4] Remove "negative" --- examples/data/LG M50 Synthetic/Full cell/EIS/metadata.json | 2 +- pybop/pybamm/synthetic_utils.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/data/LG M50 Synthetic/Full cell/EIS/metadata.json b/examples/data/LG M50 Synthetic/Full cell/EIS/metadata.json index 8487da385..b4e751942 100644 --- a/examples/data/LG M50 Synthetic/Full cell/EIS/metadata.json +++ b/examples/data/LG M50 Synthetic/Full cell/EIS/metadata.json @@ -28,7 +28,7 @@ "Initial voltage": "The terminal voltage prior to the EIS measurement.", "Frequency": "The input perturbation frequency.", "Impedance (real)": "The real component of the complex impedance.", - "Impedance (imag)": "The negative imaginary component of the complex impedance.", + "Impedance (imag)": "The imaginary component of the complex impedance.", "Procedure Time": "Time elapsed since beginning of procedure.", "Procedure Capacity": "The net charge passed since beginning of procedure." }, diff --git a/pybop/pybamm/synthetic_utils.py b/pybop/pybamm/synthetic_utils.py index 4b49f2c67..f1e54983d 100644 --- a/pybop/pybamm/synthetic_utils.py +++ b/pybop/pybamm/synthetic_utils.py @@ -51,7 +51,7 @@ "Initial voltage": "The terminal voltage prior to the EIS measurement.", "Frequency": "The input perturbation frequency.", "Impedance (real)": "The real component of the complex impedance.", - "Impedance (imag)": "The negative imaginary component of the complex impedance.", + "Impedance (imag)": "The imaginary component of the complex impedance.", # "Temperature": "The temperature of the cell.", } From 9d5d8d8eea971c6db15224dc75ceb05ffefeb0b2 Mon Sep 17 00:00:00 2001 From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com> Date: Thu, 6 Aug 2026 17:44:50 +0100 Subject: [PATCH 2/4] Create 5_fit_resistances.py --- .../dfn_parameterisation/5_fit_resistances.py | 272 ++++++++++++++++++ 1 file changed, 272 insertions(+) create mode 100644 examples/scripts/dfn_parameterisation/5_fit_resistances.py diff --git a/examples/scripts/dfn_parameterisation/5_fit_resistances.py b/examples/scripts/dfn_parameterisation/5_fit_resistances.py new file mode 100644 index 000000000..33e7aa918 --- /dev/null +++ b/examples/scripts/dfn_parameterisation/5_fit_resistances.py @@ -0,0 +1,272 @@ +import argparse +import pickle +from pathlib import Path + +import matplotlib.pyplot as plt +import numpy as np +import pybamm +import pyprobe +import scienceplots # noqa: F401 +from matplotlib import colormaps + +import pybop +from pybop.applications.utils import get_cells, get_ocv_function, shift_ocv_to + +""" +Our data is stored in this folder structure: +examples/data + - {Battery type} + - {Cell type} + - {Cell_number} + - {Procedure}.parquet + - metadata.json + +The EIS procedure differs for the full-cells (evenly spread SOC points) and half-cells +(not evenly spaced but can be found from the measured charge throughput). + +For each full-cell, estimate the charge transfer and series resistance by fitting the EIS data. +""" + +parser = argparse.ArgumentParser() +parser.add_argument("--show-plots", action="store_true", help="Show plots") +args = parser.parse_args() + +plt.style.use("science") +plt.rcParams.update({"text.usetex": False}) # users can remove this line + +# Options +OCP_type = "pOCP" + + +class R_CPE(pybop.BaseSimulator): + def __init__(self, parameter_values, dataset): + self._parameter_values = parameter_values.copy() + super().__init__(parameters=self._parameter_values) + self.domain_data = dataset["Frequency [Hz]"] + + def get_model_parameters(self, inputs): + param = self._parameter_values.copy() + param.update(inputs) + R0 = param.evaluate(param["Series resistance [Ohm]"]) + R1 = param.evaluate(param["CPE resistance [Ohm]"]) # soc-dependent + C1 = param.evaluate(param["CPE capacitance [F]"]) + a1 = param.evaluate(param["CPE exponent"]) + + # Get the equivalent capacitance for the CPE + Q1 = C1**a1 / (R1 ** (1 - a1)) + return R0, R1, Q1, a1 + + def solve_batch(self, inputs, calculate_sensitivities: bool = False): + solutions = [] + for x in inputs: + R0, R1, Q1, a1 = self.get_model_parameters(x) + f_eval = self.domain_data + Z = R0 + R1 / ( + 1 + R1 * Q1 * np.asarray(1j * 2 * np.pi * f_eval, np.complex128) ** a1 + ) + + sol = pybop.Solution() + sol.set_solution_variable( + "Impedance", data=np.asarray(Z, dtype=np.complex128) + ) + solutions.append(sol) + return solutions + + +for cell in get_cells(): + cell_label = cell.info["Cell label"] + cell_type = cell.info["Cell type"] + cell_path = Path(__file__).parent / "results" / cell_type / cell_label + + """ Choose the charge branch for now. Later decide how to deal with hysteresis """ + direction = "charge" + + # Load parameters + with open(cell_path / f"params_{direction}.pickle", "rb") as file: + param = pickle.load(file) + + # Get data + eis_cell = pyprobe.load_archive( + str(cell.info["Archive"] / (cell.info["Associated"]["Full-cell EIS"])) + ) + experiment = eis_cell.procedure["EIS charge"].experiment("EIS charge") + + # Fit the DC voltages to the OCV function to obtain SOC, assume evenly spaced + # but allow a small voltage offset between initial (rest) voltage and OCV + voltage_points = ( + experiment.lf.select("Initial voltage [V]") + .drop_nulls() + .unique(maintain_order=True) + .collect() + .to_numpy() + .flatten() + ) + N_SOC = len(voltage_points) + + # Plot the given OCV + fig_v, ax_v = plt.subplots(1, 1, figsize=(5, 3)) + ax_v.set_xlabel("State of charge") + ax_v.set_ylabel("Voltage [V]") + soc = np.linspace(0, 1, 501) + ocv_function = get_ocv_function(param, OCP_type) + ax_v.plot(soc, ocv_function(soc), label="Pseudo OCV") + + # Shift the OCV to match the voltage points + naive_soc = 0.1 + 0.1 * np.arange(N_SOC) + shift_param, SOC_vec = shift_ocv_to(voltage_points, param, naive_soc, OCP_type) + + """ Overwrite the OCV functions in the saved *charge* params. """ + param = shift_param + + # Validate the alignment by plotting + ax_v.scatter(SOC_vec, voltage_points, marker="o", label="EIS points") + ocv_function = get_ocv_function(shift_param) + ax_v.plot(soc, ocv_function(soc), label="OCV") + ax_v.legend() + fig_v.savefig(cell_path / f"{cell_label}_EIS_points.svg") + + # Set up Nyquist plot + fig, ax = plt.subplots(1, 1, figsize=(5, 3)) + ax.set_title(f"{cell_type}, Nyquist plot") + ax.set_xlabel(r"Re(Z) [m$\Omega$]") + ax.set_ylabel(r"-Im(Z) [m$\Omega$]") + colourmap = colormaps["viridis"].resampled(N_SOC) + + # Define shared parameters + initial_series_resistance = ( + experiment.lf.select("Impedance (real) [Ohm]") + .drop_nulls() + .first() + .collect() + .item() + ) + series_resistance = pybop.Parameter( + initial_value=initial_series_resistance, bounds=[0, 0.1] + ) + dl_capacitance = pybop.Parameter(initial_value=0.1, bounds=[0, 1]) + + # Estimate charge transfer resistances over SOC for fixed R0 + problems = [] + for i in range(N_SOC): + eis_data = experiment.cycle(i).step(2) + frequencies = eis_data.get("Frequency [Hz]") + Zr = eis_data.get("Impedance (real) [Ohm]") + Zj = eis_data.get("Impedance (imag) [Ohm]") + + Z = Zr + 1j * Zj # impedance data [Ohm] + Z = np.asarray(Z, dtype=np.complex128) + + # Select the frequency range of semi-circle by first finding the peak of the semicircle + # and then finding the mid-frequency where imaginary part is smallest + peak_idx = np.where(np.diff(Z.imag) > 0)[0][0] + idx_fmin = peak_idx + np.argmin(np.abs(Z.imag[peak_idx:])) + + # Remove some low-frequency points that may contain diffusion phenomena and + # some high frequencies that may be affected by inductance + idx_fmin -= 3 + idx_to_the_left = np.where(Z.imag[:peak_idx] > Z.imag[idx_fmin])[0] + idx_fmax = idx_to_the_left[-3] if len(idx_to_the_left) > 2 else 0 + + # Add to plot, include a few points outside of the fitting region + label = rf"{(100 * SOC_vec[i]):.0f} \%" + Z_plot = Z[idx_fmax : idx_fmin + 5] + ax.scatter( + Z_plot.real * 1e3, + -Z_plot.imag * 1e3, + color=colourmap(i / N_SOC), + marker=".", + s=4, + label=f"_{label}", + ) + + # Select range for fitting + f_semicircle = frequencies[idx_fmax:idx_fmin] + Z_semicircle = Z[idx_fmax:idx_fmin] + + dataset = pybop.Dataset( + { + "Frequency [Hz]": f_semicircle, + "Impedance": Z_semicircle, + }, + domain="Frequency [Hz]", + ) + + # Define optimisation parameters + initial_cpe_resistance = Z_semicircle[-1].real - initial_series_resistance + model_param = pybamm.ParameterValues( + { + "Series resistance [Ohm]": series_resistance, + "CPE capacitance [F]": dl_capacitance, + f"CPE resistance {i} [Ohm]": pybop.Parameter( + initial_value=initial_cpe_resistance, + bounds=[0, 4 * initial_cpe_resistance], + ), + f"CPE exponent {i}": pybop.Parameter( + initial_value=0.9, bounds=[0.4, 1] + ), + "CPE resistance [Ohm]": pybamm.Parameter(f"CPE resistance {i} [Ohm]"), + "CPE exponent": pybamm.Parameter(f"CPE exponent {i}"), + } + ) + + # Set up an optimisation problem + simulator = R_CPE(parameter_values=model_param, dataset=dataset) + cost = pybop.RootMeanSquaredError(dataset, target="Impedance") + problem = pybop.Problem( + simulator=simulator, + cost=pybop.WeightedCost(cost, weights=[1e3]), # convert Ohm to mOhm + ) + problems.append(problem) + + meta_problem = pybop.MetaProblem(*problems) + options = pybop.SciPyDifferentialEvolutionOptions(maxiter=500) + optim = pybop.SciPyDifferentialEvolution(meta_problem, options=options) + result = optim.run() + print(result) + inputs = result.best_inputs + + Rct = np.zeros(N_SOC) + alpha = np.zeros(N_SOC) + rmse_mOhm = np.zeros(N_SOC) + for i, problem in enumerate(problems): + # Simulate the best fit parameters + solution = problem.simulate(inputs=inputs) + + # Add simluation to plot + label = rf"{(100 * SOC_vec[i]):.0f} \%" + ax.plot( + solution["Impedance"].data.real * 1e3, + -solution["Impedance"].data.imag * 1e3, + color=colourmap(i / N_SOC), + ls="-", + linewidth=1.5, + label=label, + ) + + Rct[i] = inputs[f"CPE resistance {i} [Ohm]"] + alpha[i] = inputs[f"CPE exponent {i}"] + rmse_mOhm[i] = problem.cost.evaluate(solution, inputs).values.item() + + ax.legend(fontsize="small", loc="center left", bbox_to_anchor=(1, 0.5)) + ax.grid(True) + ax.axis("equal") + + # Update and save parameters + param.update( + { + "Series resistance [Ohm]": inputs["Series resistance [Ohm]"], + "Charge transfer resistance [Ohm]": pybop.Interpolant( + SOC_vec, Rct, name="Charge transfer resistance" + ), + "Double-layer capacitance [F]": inputs["CPE capacitance [F]"], + "Mean EIS RMSE [mOhm]": np.mean(rmse_mOhm), + } + ) + with open(cell_path / f"params_{direction}.pickle", "wb") as file: + pickle.dump(param, file) + + fig.savefig(cell_path / f"{cell_label}_EIS_resistances.svg") + +if args.show_plots: + plt.show() +plt.close() From b69bebd1efe10e8601444b61506624e564bec223 Mon Sep 17 00:00:00 2001 From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com> Date: Mon, 17 Aug 2026 22:49:57 +0100 Subject: [PATCH 3/4] Create 5b_charge_transfer.py --- .../5b_charge_transfer.py | 205 ++++++++++++++++++ 1 file changed, 205 insertions(+) create mode 100644 examples/scripts/dfn_parameterisation/5b_charge_transfer.py diff --git a/examples/scripts/dfn_parameterisation/5b_charge_transfer.py b/examples/scripts/dfn_parameterisation/5b_charge_transfer.py new file mode 100644 index 000000000..5a8510ab1 --- /dev/null +++ b/examples/scripts/dfn_parameterisation/5b_charge_transfer.py @@ -0,0 +1,205 @@ +import argparse +import pickle +from pathlib import Path + +import matplotlib.pyplot as plt +import numpy as np +import pybamm +import scienceplots # noqa: F401 + +import pybop +from pybop.applications.utils import get_cells +from pybop.models.lithium_ion.base_model import BaseGroupedModel + +""" +Our data is stored in this folder structure: +examples/data + - {Battery type} + - {Cell type} + - {Cell_number} + - {Procedure}.parquet + - metadata.json + +For each full-cell, fit the Butler-Volmer relation to the charge transfer resistance estimates. +""" + +parser = argparse.ArgumentParser() +parser.add_argument("--show-plots", action="store_true", help="Show plots") +args = parser.parse_args() + +plt.style.use("science") +plt.rcParams.update({"text.usetex": False}) # users can remove this line + + +class ChargeTransferResistance(pybop.BaseSimulator): + def __init__(self, parameter_values, dataset, electrode: str): + self._parameter_values = parameter_values.copy() + super().__init__(parameters=self._parameter_values) + self.domain_data = dataset["Stoichiometry"] + self.electrode = electrode + + def get_model_parameters(self, inputs): + param = self._parameter_values.copy() + param.update(inputs) + Domain = self.electrode.capitalize() + alpha = param[f"{Domain} electrode charge transfer coefficient"] + omega = param[f"{Domain} electrode charge transfer ideality factor"] + sto_e = param["Initial electrolyte mole fraction"] + Rct_typ = param["Typical charge transfer resistance [Ohm]"] + return alpha, omega, sto_e, Rct_typ + + def solve_batch(self, inputs, calculate_sensitivities: bool = False): + solutions = [] + for x in inputs: + alpha, omega, sto_e, Rct_typ = self.get_model_parameters(x) + sto = np.clip(self.domain_data, 0, 1) + + j0 = ( + sto ** (alpha * omega) + * (1 - sto) ** ((1 - alpha) * omega) + * sto_e ** (1 - alpha) + ) + Rct = Rct_typ / (2 * j0) + + sol = pybop.Solution() + sol.set_solution_variable("Dimensionless exchange rate", data=j0) + sol.set_solution_variable( + "Charge transfer resistance [Ohm]", data=np.nan_to_num(Rct, nan=1e6) + ) + solutions.append(sol) + return solutions + + +for cell in get_cells(): + cell_label = cell.info["Cell label"] + cell_type = cell.info["Cell type"] + cell_path = Path(__file__).parent / "results" / cell_type / cell_label + + """ Use the charge values for both charge and discharge parameter sets. """ + direction = "charge" + + """ Select the dominant electrode to attribute the charge transfer resistance to. """ + dominant_electrode = "Negative" + + # Load parameters + with open(cell_path / f"params_{direction}.pickle", "rb") as file: + param = pickle.load(file) + + # Fit the Butler-Volmer relation + soc_values = param["Charge transfer resistance [Ohm]"].x + Rct_values = param["Charge transfer resistance [Ohm]"].y + + if dominant_electrode == "Negative": + x_0 = param["Minimum negative stoichiometry"] + x_100 = param["Maximum negative stoichiometry"] + sto_values = x_0 + (x_100 - x_0) * soc_values + else: + y_100 = param["Minimum positive stoichiometry"] + y_0 = param["Maximum positive stoichiometry"] + sto_values = y_0 + (y_100 - y_0) * soc_values + + dataset = pybop.Dataset( + { + "Stoichiometry": sto_values, + "Charge transfer resistance [Ohm]": Rct_values, + }, + domain="Stoichiometry", + ) + + # Define parameter values + model_param = { + f"{dominant_electrode} electrode charge transfer coefficient": pybop.Parameter( + initial_value=0.75, bounds=[0, 1] + ), + f"{dominant_electrode} electrode charge transfer ideality factor": pybop.Parameter( + initial_value=0.75, bounds=[0, 1] + ), + "Initial electrolyte mole fraction": 1.0, + "Typical charge transfer resistance [Ohm]": pybop.Parameter( + initial_value=np.mean(Rct_values), bounds=[0, 1] + ), + } + + # Fit the reaction rate to the exchange current density + simulator = ChargeTransferResistance( + model_param, dataset=dataset, electrode=dominant_electrode + ) + cost = pybop.RootMeanSquaredError( + dataset, target="Charge transfer resistance [Ohm]" + ) + problem = pybop.Problem( + simulator=simulator, + cost=pybop.WeightedCost(cost, weights=[1e3]), # convert Ohm to mOhm + ) + + options = pybop.SciPyDifferentialEvolutionOptions(maxiter=250) + optim = pybop.SciPyDifferentialEvolution(problem, options=options) + result = optim.run() + + # Round and display results + alpha, omega, sto_e, Rct_typ = simulator.get_model_parameters(result.best_inputs) + alpha = np.round(alpha, decimals=2) + omega = np.round(omega, decimals=2) + print("Typical charge transfer resistance [mOhm]:", Rct_typ) + print("Charge transfer coefficient:", alpha) + print("Charge transfer ideality factor:", omega) + print(r"RMSE [m$\Omega$]:", result.best_cost, "\n") + + # Use optimisation result to plot the fitted function + fig, ax = plt.subplots(1, 1, figsize=(5, 3)) + ax.set_xlabel(f"{dominant_electrode} electrode stoichiometry") + ax.set_ylabel(r"Charge transfer resistance [m$\Omega$]") + ax.scatter(sto_values, Rct_values * 1e3, label="ECM fits") + sto_range = np.linspace(0.03, 0.97, 95) + simulator.domain_data = sto_range + solution = simulator.solve(result.best_inputs) + ax.plot( + sto_range, + solution["Charge transfer resistance [Ohm]"].data * 1e3, + label="BV fit", + linewidth=1, + ) + ax.set_title(f"{cell_type}, Butler-Volmer fit") + ax.legend() + + # Convert charge transfer resistance into a charge transfer timescale + Q_th = (param["Measured cell capacity [A.h]"] * 3600) / ( + param["Maximum " + dominant_electrode.lower() + " stoichiometry"] + - param["Minimum " + dominant_electrode.lower() + " stoichiometry"] + ) + RT_F = ( + pybamm.constants.R.value + * param["Ambient temperature [K]"] + / pybamm.constants.F.value + ) + tau_ct = (3 * Q_th * Rct_typ) / (2 * RT_F) + + # Update and save parameters + other_electrode = "Positive" if dominant_electrode == "Negative" else "Negative" + zero = np.finfo(np.float64).eps # almost zero + param.update( + { + dominant_electrode + " electrode charge transfer time scale [s]": tau_ct, + dominant_electrode + " electrode dimensionless exchange rate": ( + BaseGroupedModel.get_multiphase_butler_volmer(dominant_electrode) + ), + dominant_electrode + " electrode charge transfer coefficient": alpha, + dominant_electrode + " electrode charge transfer ideality factor": omega, + dominant_electrode + " electrode capacitance [F]": param[ + "Double-layer capacitance [F]" + ], + other_electrode + " electrode charge transfer time scale [s]": 10, # fast + other_electrode + " electrode dimensionless exchange rate": ( + BaseGroupedModel.symmetric_butler_volmer + ), + other_electrode + " electrode capacitance [F]": zero, + } + ) + with open(cell_path / f"params_{direction}.pickle", "wb") as file: + pickle.dump(param, file) + + fig.savefig(cell_path / f"{cell_label}_charge_transfer.svg") + +if args.show_plots: + plt.show() +plt.close() From edd2622fb5f4af4108d7c958bda49bc9c792dbcd Mon Sep 17 00:00:00 2001 From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com> Date: Mon, 17 Aug 2026 22:59:05 +0100 Subject: [PATCH 4/4] Update CHANGELOG.md --- CHANGELOG.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index b75e88531..641fdcd6c 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,7 @@ ## Features +- [#973](https://github.com/pybop-team/PyBOP/pull/973) - Adds example scripts for EIS parameterisation. - [#974](https://github.com/pybop-team/PyBOP/pull/974) - Adds voltage components to each grouped model as well as asymmetric and multiphase Butler-Volmer kinetics. - [#969](https://github.com/pybop-team/PyBOP/pull/969) - Updates synthetic data and adds example script for thermal parameterisation. - [#965](https://github.com/pybop-team/PyBOP/pull/965) - Adds synthetic data and example scripts for OCV parameterisation.