From 8f28cb601de1d674f328dd1e6f4a0decb3aec9ff Mon Sep 17 00:00:00 2001 From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com> Date: Fri, 30 Oct 2026 15:07:36 +0000 Subject: [PATCH 01/18] Add get time in hours --- pybop/processing/dataset.py | 2 ++ tests/unit/test_dataset.py | 3 +++ 2 files changed, 5 insertions(+) diff --git a/pybop/processing/dataset.py b/pybop/processing/dataset.py index f2b2b45cd..b37cf6324 100644 --- a/pybop/processing/dataset.py +++ b/pybop/processing/dataset.py @@ -62,6 +62,8 @@ def __setitem__(self, key, value): def __getitem__(self, key): """Return the data corresponding to a particular key.""" if key not in self.data.keys(): + if key == "Time [h]": + return self["Time [s]"] / 3600 raise ValueError(f"The key {key} does not exist in this dataset.") return self.data[key] diff --git a/tests/unit/test_dataset.py b/tests/unit/test_dataset.py index a64d69687..c0393caec 100644 --- a/tests/unit/test_dataset.py +++ b/tests/unit/test_dataset.py @@ -36,6 +36,9 @@ def test_dataset(self): assert dataset.data == data_dictionary assert np.all(dataset["Time [s]"] == solution["Time [s]"].data) + # Test time conversion + np.testing.assert_allclose(dataset["Time [h]"], solution["Time [h]"].data) + # Test exception for non-dictionary inputs with pytest.raises( TypeError, match="The input to pybop.Dataset must be a dictionary." From 747e2b87450c27488230fada308bfc5eeaafcc3d Mon Sep 17 00:00:00 2001 From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com> Date: Fri, 30 Oct 2026 16:56:42 +0000 Subject: [PATCH 02/18] Update shift ocv and add downsample --- pybop/applications/utils.py | 41 ++++++++++++++++++---------- tests/unit/test_application_utils.py | 24 ++++++++++++---- 2 files changed, 45 insertions(+), 20 deletions(-) diff --git a/pybop/applications/utils.py b/pybop/applications/utils.py index c4a9f681e..a0a2a54a4 100644 --- a/pybop/applications/utils.py +++ b/pybop/applications/utils.py @@ -161,7 +161,7 @@ def get_ocv_function(parameter_values, OCP_type: str = None): def shift_ocv_to(voltage_points, parameter_values, naive_soc, OCP_type: str = None): - """Returns a new lithium-ion parameter set with shifted OCP functions and SOC points.""" + """Returns a new parameter set with shifted OCPs, the SOC points and voltage offset.""" OCP_type = "OCP" if OCP_type is None else OCP_type param = parameter_values.copy() @@ -170,24 +170,37 @@ def shift_ocv_to(voltage_points, parameter_values, naive_soc, OCP_type: str = No def mean_absolute_error(values): shift, stretch, offset = values return np.mean( - np.abs(voltage_points - ocv_function(shift + stretch * naive_soc) - offset) + np.abs(voltage_points - ocv_function(shift + stretch * naive_soc) + offset) ) x0 = [0, 1, 0] - scipy_result = minimize(mean_absolute_error, x0=x0) + scipy_result = minimize(mean_absolute_error, x0=x0, method="trust-constr") SOC_vec = scipy_result.x[0] + scipy_result.x[1] * naive_soc + if not scipy_result.success: + print("WARNING:", scipy_result.message) # Shift each of the OCP functions by half the offset to remove series resistance offset = scipy_result.x[2] - param["Positive electrode OCP [V]"] = Interpolant( - param[f"Positive electrode {OCP_type} [V]"].x, - param[f"Positive electrode {OCP_type} [V]"].y + offset / 2, - name="Positive OCP", - ) - param["Negative electrode OCP [V]"] = Interpolant( - param[f"Negative electrode {OCP_type} [V]"].x, - param[f"Negative electrode {OCP_type} [V]"].y - offset / 2, - name="Negative OCP", - ) + for electrode, sign in zip( + ["Positive electrode", "Negative electrode"], [1, -1], strict=False + ): + param[f"{electrode} {OCP_type} [V]"] = Interpolant( + param[f"{electrode} {OCP_type} [V]"].x, + param[f"{electrode} {OCP_type} [V]"].y - sign * offset / 2, + name=f"{electrode} {OCP_type}", + ) + + return param, SOC_vec, offset + - return param, SOC_vec +def downsample(x, y, tol_x=300, tol_y=1e-3): + """Iterate over the data and discard points that are too close to the previous point.""" + new_x, new_y = [x[0]], [y[0]] + for xi, yi in zip(x[1:-1], y[1:-1], strict=False): + if np.abs(xi - new_x[-1]) < tol_x and np.abs(yi - new_y[-1]) < tol_y: + continue + new_x.append(xi) + new_y.append(yi) + new_x.append(x[-1]) + new_y.append(y[-1]) + return np.asarray(new_x), np.asarray(new_y) diff --git a/tests/unit/test_application_utils.py b/tests/unit/test_application_utils.py index 0f44c4ead..0dba5197c 100644 --- a/tests/unit/test_application_utils.py +++ b/tests/unit/test_application_utils.py @@ -13,6 +13,7 @@ import pybop from pybop.applications.utils import ( OpenCircuitVoltage, + downsample, filter_with_preceding_row, get_cells, get_ocp_functions, @@ -77,12 +78,15 @@ def test_get_ocv_function(self, parameter_values): def test_shift_ocv_to(self, parameter_values): ocv_function = get_ocv_function(parameter_values, OCP_type="OCP") - soc = np.linspace(0, 1, 5) - voltage_points = ocv_function(soc) + 0.1 - param, soc_points = shift_ocv_to(voltage_points, parameter_values, soc) - ocv_function = get_ocv_function(param) - assert np.all(ocv_function(soc) > 0.1) - np.testing.assert_allclose(soc, soc_points, atol=2.5e-4) + soc = np.linspace(0, 1, 10) + lower_voltages = ocv_function(soc) - 0.1 + shifted_param, soc_points, offset = shift_ocv_to( + lower_voltages, parameter_values, soc + ) + shifted_ocv = get_ocv_function(shifted_param) + np.testing.assert_allclose(offset, 0.1, atol=1e-6) + np.testing.assert_allclose(soc, soc_points, atol=1e-8) + np.testing.assert_allclose(shifted_ocv(soc), lower_voltages) def test_make_voltage_monotonic(self): # Test ascending @@ -99,6 +103,14 @@ def test_make_voltage_monotonic(self): data = make_voltage_monotonic(data) assert np.all(data.lf["Voltage [V]"].to_numpy() >= voltage) + def test_downsample(self): + tol_x, tol_y = 300, 1e-3 + x = np.array([0, 1, 1, 1, 3]) * tol_x + y = np.array([1, 1, 3, 3, 4]) * tol_y + x_new, y_new = downsample(x, y, tol_x=tol_x, tol_y=tol_y) + assert (x_new == np.array([0, 1, 1, 3]) * tol_x).all() + assert (y_new == np.array([1, 1, 3, 4]) * tol_y).all() + @pytest.mark.skipif( sys.version_info < (3, 11), reason="requires a python version >= 3.11" ) From 1f614ef1763d65c7e5966c15a639deb2967e3bb0 Mon Sep 17 00:00:00 2001 From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com> Date: Fri, 30 Oct 2026 15:11:40 +0000 Subject: [PATCH 03/18] Update 0_plot_data.py --- .../dfn_parameterisation/0_plot_data.py | 28 ++++++++++++------- 1 file changed, 18 insertions(+), 10 deletions(-) diff --git a/examples/scripts/dfn_parameterisation/0_plot_data.py b/examples/scripts/dfn_parameterisation/0_plot_data.py index db9369f89..8312b98a0 100644 --- a/examples/scripts/dfn_parameterisation/0_plot_data.py +++ b/examples/scripts/dfn_parameterisation/0_plot_data.py @@ -3,6 +3,7 @@ from pathlib import Path import matplotlib.pyplot as plt +import scienceplots # noqa: F401 from pybop.applications.utils import get_cells @@ -19,9 +20,15 @@ """ parser = argparse.ArgumentParser() +parser.add_argument("--draft", action="store_true", help="Draft mode") parser.add_argument("--show-plots", action="store_true", help="Show plots") args = parser.parse_args() +plt.style.use("science") + +if args.draft: + plt.rcParams.update({"text.usetex": False}) + for cell in get_cells(match=""): cell_label = cell.info["Cell label"] cell_type = cell.info["Cell type"] @@ -38,12 +45,10 @@ step_descriptions = procedure.readme_dict[experiment_name][ "Step Descriptions" ] - title = f"Cell: {cell.info['Name']}, Experiment: {experiment_name}\n" + str( - step_descriptions - ) + title = f"Cell: {cell.info['Name']}, Experiment: {experiment_name}" # Plot the experiment - fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 6)) + fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 6), dpi=300) ax_properties = dict(x="Time [s]", y="Voltage [V]", grid=True) # Plot the cycles over the full dataset @@ -80,13 +85,14 @@ label=name, ls="--", ) - ax2.legend() + ax2.legend(loc="center left", bbox_to_anchor=(1, 0.5)) - fig.suptitle(title, fontsize=10, wrap=True) - fig.savefig(cell_path / f"{cell_label}_{experiment_name}.svg") + fig.suptitle(title, fontsize=14, wrap=True) + fig.tight_layout() + fig.savefig(cell_path / f"{cell_label}_{experiment_name}.png") if "EIS" in experiment_name: - fig, ax = plt.subplots(figsize=(6, 6)) + fig, ax = plt.subplots(figsize=(6, 6), dpi=300) ax_properties = dict( ax=ax, x="Impedance (real) [Ohm]", @@ -108,8 +114,10 @@ ax.yaxis.set_inverted(True) ax.set_aspect("equal", adjustable="box") - fig.suptitle(title, fontsize=10, wrap=True) - fig.savefig(cell_path / f"{cell_label}_{experiment_name}_Nyquist.svg") + title += "\n" + str(step_descriptions) + fig.suptitle(title, fontsize=14, wrap=True) + fig.tight_layout() + fig.savefig(cell_path / f"{cell_label}_{experiment_name}_Nyquist.png") if args.show_plots: plt.show() From c9dcb9692349002004ad76ff55921d65ec95f560 Mon Sep 17 00:00:00 2001 From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com> Date: Fri, 30 Oct 2026 15:12:18 +0000 Subject: [PATCH 04/18] Update 1_create_parameters.py --- .../dfn_parameterisation/1_create_parameters.py | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/examples/scripts/dfn_parameterisation/1_create_parameters.py b/examples/scripts/dfn_parameterisation/1_create_parameters.py index bc206a9be..2704264ab 100644 --- a/examples/scripts/dfn_parameterisation/1_create_parameters.py +++ b/examples/scripts/dfn_parameterisation/1_create_parameters.py @@ -1,3 +1,4 @@ +import argparse import pickle from pathlib import Path @@ -16,6 +17,10 @@ dictionary for charge and discharge. """ +parser = argparse.ArgumentParser() +parser.add_argument("--draft", action="store_true", help="Draft mode") +args = parser.parse_args() + for cell in get_cells(): cell_label = cell.info["Cell label"] cell_type = cell.info["Cell type"] @@ -63,6 +68,16 @@ "Negative electrode relative thickness": 85.2 / three_layer_thickness, } ) + elif cell_type == "Enertech SPB655060 Synthetic": + three_layer_thickness = 68 + 25 + 76.5 + param.update( + { + "Lower voltage cut-off [V]": 3.0, + "Upper voltage cut-off [V]": 4.2, + "Positive electrode relative thickness": 68 / three_layer_thickness, + "Negative electrode relative thickness": 76.5 / three_layer_thickness, + } + ) else: raise ValueError(f"Unrecognised cell type: {cell_type}") From d47afef12026d94cce528bf48c14704c4d16da82 Mon Sep 17 00:00:00 2001 From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com> Date: Fri, 30 Oct 2026 15:13:52 +0000 Subject: [PATCH 05/18] Update 2_electrode_ocps.py --- .../dfn_parameterisation/2_electrode_ocps.py | 95 ++++++++++++------- 1 file changed, 61 insertions(+), 34 deletions(-) diff --git a/examples/scripts/dfn_parameterisation/2_electrode_ocps.py b/examples/scripts/dfn_parameterisation/2_electrode_ocps.py index 94cb2083d..39bbc0050 100644 --- a/examples/scripts/dfn_parameterisation/2_electrode_ocps.py +++ b/examples/scripts/dfn_parameterisation/2_electrode_ocps.py @@ -9,7 +9,11 @@ from pyprobe.analysis import smoothing import pybop -from pybop.applications.utils import get_cells, make_voltage_monotonic +from pybop.applications.utils import ( + filter_with_preceding_row, + get_cells, + make_voltage_monotonic, +) """ Our data is stored in this folder structure: @@ -21,21 +25,23 @@ - metadata.json The pseudo-OCP procedure for positive (/negative) half-cells is: -0. Rest -1. Charge (/discharge) at C/30 until voltage limit (direction corresponding to full-cell charging) -2. Rest for 2 hours -3. Discharge (/charge) at C/30 until voltage limit (direction corresponding to full-cell discharging) -4. Rest for 2 hours +0. Charge (/discharge) at C/30 until voltage limit (direction corresponding to full-cell charging) +1. Rest for 2 hours +2. Discharge (/charge) at C/30 until voltage limit (direction corresponding to full-cell discharging) +3. Rest for 2 hours For each full-cell, load the half-cell pOCPs and add as interpolants to the parameter set. """ parser = argparse.ArgumentParser() +parser.add_argument("--draft", action="store_true", help="Draft mode") 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 + +if args.draft: + plt.rcParams.update({"text.usetex": False}) # Options smoothing_options = { @@ -44,6 +50,8 @@ "monotonic": True, } account_for_overpotential = True +margin = 0.01 # to keep the stoichiometry within [margin, 1-margin] + for cell in get_cells(): cell_label = cell.info["Cell label"] @@ -67,14 +75,18 @@ half_cell = pyprobe.load_archive( str(cell.info["Archive"] / cell.info["Associated"][electrode]) ) - experiment = half_cell.procedure["pOCP"].experiment("Pseudo OCP") + procedure = half_cell.procedure["pOCP"] + experiment = procedure.experiment("Pseudo OCP") cycle = 1 # use the second cycle for fitting # Plot data for validation - validation_data = pybop.import_pyprobe_result(experiment) + validation_data = pybop.import_pyprobe_result( + filter_with_preceding_row(procedure, experiment="Pseudo OCP") + ) ax_v.plot( - (validation_data["Time [s]"] - validation_data["Time [s]"][0]) / 3600, + validation_data["Time [h]"] - validation_data["Time [h]"][0], validation_data["Voltage [V]"], + ls=":", label=f"{electrode}", ) @@ -82,11 +94,11 @@ if "Positive" in electrode: # Starts fully lithiated and then delithiates direction_order = ["charge", "discharge"] - x0_order = [1.0, 0.0] + x0_order = [1.0 - margin, margin] else: # starts fully delithiated and then lithiates direction_order = ["discharge", "charge"] - x0_order = [0.0, 1.0] + x0_order = [margin, 1.0 - margin] for j, (direction, x0) in enumerate( zip(direction_order, x0_order, strict=False) @@ -102,27 +114,41 @@ ) ) - # Ensure that the charge throughput starts at zero - Q_start = dataset["Discharge capacity [A.h]"][0] + # Ensure that the charge throughput starts from zero + Q_start = -filter_with_preceding_row( + procedure, experiment="Pseudo OCP", cycle=cycle, phase=direction + ).get("Capacity [Ah]")[0] charge_throughput = dataset["Discharge capacity [A.h]"] - Q_start Q_meas = np.abs(charge_throughput[-1]) - # Compute the Coulombic efficiency to take account of dynamic overshoot and make a closed loop - CE = 1.0 + # Compute the dynamic efficiency to take account of voltage overshoot + DE = 1.0 if account_for_overpotential: + rest_phase_voltage = ( + experiment.rest().step(2 * cycle + j).get("Voltage [V]") + ) inverse_ocv_interpolant = pybop.Interpolant( dataset["Voltage [V]"], charge_throughput / charge_throughput[-1], name="Inverse OCV", ) - V_relax = experiment.rest().step(2 * cycle + j).get("Voltage [V]")[-1] + V_relax = rest_phase_voltage[-1] + + # Get the local overpotential at the current switch-off and add it + # back onto the rest voltage to make it comparable with pseudo-OCV + V_relax += ( + (dataset["Voltage [V]"][-1] - rest_phase_voltage[1]) + / dataset["Current [A]"][-1] + * np.mean(dataset["Current [A]"][1:101]) + ) # Find the scaling needed to make the rest voltage correspond to a fully de/lithiated state - CE = np.minimum(inverse_ocv_interpolant(V_relax), 1.0) - Q_cell = CE * Q_meas + DE = np.minimum(inverse_ocv_interpolant(V_relax), 1.0) + Q_th = Q_meas / (1.0 - 2 * margin) + Q_scale = Q_th * DE # Define state of lithiation - dataset["State of lithiation"] = x0 + charge_throughput / Q_cell + dataset["State of lithiation"] = x0 + charge_throughput / Q_scale if direction == "charge": assert ( @@ -146,26 +172,27 @@ ) # Plot - sol = np.linspace(0, 1, 501) + sol = np.linspace(margin, 1.0 - margin, 501) colour = "tab:red" if direction == "charge" else "tab:blue" ax[i].scatter( dataset["State of lithiation"], dataset["Voltage [V]"], - marker="x", - s=2, + marker=".", + s=4, color=colour, label=direction, ) - ax[i].plot(sol, ocp_branch(sol), color=colour) + ax[i].plot(sol, ocp_branch(sol), linewidth=1, color=colour, alpha=0.8) # Plot validation, align the dis/charge branches at the start of the first dis/charge - offset = experiment.rest().step(j - 1).df["Capacity [Ah]"][-1] / Q_meas + offset = ( + experiment.charge() if direction == "charge" else experiment.discharge() + ).get("Capacity [Ah]")[0] + state = x0 + (validation_data["Discharge capacity [A.h]"] + offset) / Q_th ax_v.plot( - (validation_data["Time [s]"] - validation_data["Time [s]"][0]) / 3600, - ocp_branch( - x0 + offset + validation_data["Discharge capacity [A.h]"] / Q_meas - ), - ls="--", + validation_data["Time [h]"] - validation_data["Time [h]"][0], + ocp_branch(state), + linewidth=0.7, label=f"Fit to {direction}", ) @@ -175,8 +202,8 @@ param = pickle.load(file) param.update( { - f"{electrode} theoretical capacity (half-cell) [A.h]": Q_cell, - f"{electrode} Coulombic efficiency (half-cell)": CE, + f"{electrode} theoretical capacity (half-cell) [A.h]": Q_th, + f"{electrode} pOCP dynamic efficiency (half-cell)": DE, f"{electrode} pOCP (half-cell) [V]": ocp_branch, } ) @@ -186,9 +213,9 @@ ax[i].title.set_text(electrode) ax[i].legend() - fig.savefig(cell_path / f"{cell_label}_electrode_pOCPs.svg") + fig.savefig(cell_path / f"{cell_label}_electrode_pOCPs.png") ax_v.legend() - fig_v.savefig(cell_path / f"{cell_label}_pOCP_validation.svg") + fig_v.savefig(cell_path / f"{cell_label}_pOCP_validation.png") if args.show_plots: plt.show() From cc094a6e065791624fd463e7aefc04281dc86592 Mon Sep 17 00:00:00 2001 From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com> Date: Fri, 30 Oct 2026 15:16:55 +0000 Subject: [PATCH 06/18] Update 3_cell_balancing.py --- .../dfn_parameterisation/3_cell_balancing.py | 269 ++++++++++-------- 1 file changed, 147 insertions(+), 122 deletions(-) diff --git a/examples/scripts/dfn_parameterisation/3_cell_balancing.py b/examples/scripts/dfn_parameterisation/3_cell_balancing.py index acd05271c..e43b3aca1 100644 --- a/examples/scripts/dfn_parameterisation/3_cell_balancing.py +++ b/examples/scripts/dfn_parameterisation/3_cell_balancing.py @@ -8,7 +8,12 @@ from pyprobe.analysis import smoothing import pybop -from pybop.applications.utils import OpenCircuitVoltage, get_cells, get_ocp_functions +from pybop.applications.utils import ( + OpenCircuitVoltage, + filter_with_preceding_row, + get_cells, + get_ocp_functions, +) """ Our data is stored in this folder structure: @@ -20,22 +25,29 @@ - metadata.json Our pseudo-OCV procedure differs slightly between cell types: -0. Rest for 2 hours -1. Charge at C/30 until V_max -(2. Short hold at V_max for the LithiumWerks cells, not the Molicel) -2/3. Rest for 2 hours -3/4. Discharge at C/30 until V_max -4/5. Rest for 2 hours +0. Charge at C/30 until V_max +(1. Short hold at V_max for the LithiumWerks cells, not the Molicel) +1/2. Rest for 2 hours +2/3. Discharge at C/30 until V_max +3/4. Rest for 2 hours For each full-cell, perform electrode balancing to obtain the stoichiometry limits. """ parser = argparse.ArgumentParser() +parser.add_argument("--draft", action="store_true", help="Draft mode") 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 + +if args.draft: + plt.rcParams.update({"text.usetex": False}) + maxiter = 50 +else: + maxiter = 2500 + +np.random.seed(8) # users can remove this line # Options smoothing_options = { @@ -45,14 +57,14 @@ } OCP_type = "pOCP (half-cell)" account_for_overpotential = True -add_extra_resistance = True +apply_offset = True class OCVBalance(pybop.BaseSimulator): def __init__(self, parameter_values, dataset, direction: str, OCP_type: str): self._parameter_values = parameter_values.copy() super().__init__(parameters=self._parameter_values) - self.domain_data = dataset["Naive SoC"] + self.domain_data = dataset["SoC"] self.input_data = dataset["Current [A]"] self.direction = direction self.positive_ocp_function, self.negative_ocp_function = get_ocp_functions( @@ -63,25 +75,22 @@ def get_model_parameters(self, inputs): param = self._parameter_values.copy() param.update(inputs) Q = param["Cell capacity [A.h]"] - sto_n_0 = param["Minimum negative stoichiometry " + self.direction] - sto_n_100 = sto_n_0 + Q / param["Negative electrode theoretical capacity [A.h]"] - sto_p_100 = param["Minimum positive stoichiometry " + self.direction] - sto_p_0 = sto_p_100 + Q / param["Positive electrode theoretical capacity [A.h]"] - R0 = param["Series resistance [Ohm]"] - return Q, sto_n_0, sto_n_100, sto_p_100, sto_p_0, R0 + x_0 = param[f"Minimum negative stoichiometry ({self.direction})"] + x_100 = param["Maximum negative stoichiometry"] + y_100 = param["Minimum positive stoichiometry"] + y_0 = param[f"Maximum positive stoichiometry ({self.direction})"] + Rc = param[f"Cell resistance ({self.direction}) [Ohm]"] + return Q, x_0, x_100, y_100, y_0, Rc def solve_batch(self, inputs, calculate_sensitivities: bool = False): solutions = [] for x in inputs: - Q, sto_n_0, sto_n_100, sto_p_100, sto_p_0, R0 = self.get_model_parameters(x) + Q, x_0, x_100, y_100, y_0, Rc = self.get_model_parameters(x) soc = self.domain_data - positive_voltage = self.positive_ocp_function( - sto_p_0 + (sto_p_100 - sto_p_0) * soc - ) - negative_voltage = self.negative_ocp_function( - sto_n_0 + (sto_n_100 - sto_n_0) * soc - ) - voltage_offset = R0 * self.input_data + positive_voltage = self.positive_ocp_function(y_0 + (y_100 - y_0) * soc) + negative_voltage = self.negative_ocp_function(x_0 + (x_100 - x_0) * soc) + voltage_offset = Rc * self.input_data + voltage = positive_voltage - negative_voltage - voltage_offset sol = pybop.Solution() sol.set_solution_variable( @@ -90,9 +99,7 @@ def solve_batch(self, inputs, calculate_sensitivities: bool = False): sol.set_solution_variable( "Negative electrode voltage [V]", data=negative_voltage ) - sol.set_solution_variable( - "Voltage [V]", data=positive_voltage - negative_voltage - voltage_offset - ) + sol.set_solution_variable("Voltage [V]", data=voltage) solutions.append(sol) return solutions @@ -110,9 +117,8 @@ def solve_batch(self, inputs, calculate_sensitivities: bool = False): # Define parameters that are shared between charge and discharge branches Q_nom = cell.info["Nominal cell capacity [A.h]"] - Q_th_n = pybop.Parameter(initial_value=1.2 * Q_nom, bounds=[Q_nom, 1.5 * Q_nom]) - Q_th_p = pybop.Parameter(initial_value=1.2 * Q_nom, bounds=[Q_nom, 1.5 * Q_nom]) - R0 = pybop.Parameter(initial_value=0, bounds=[-1, 1]) + sto_n_100 = pybop.Parameter(initial_value=0.98, bounds=[0.7, 1]) + sto_p_100 = pybop.Parameter(initial_value=0.02, bounds=[0, 0.1]) # Prepare a plot fig, ax = plt.subplots(1, 1, figsize=(5, 3), dpi=300) @@ -125,79 +131,112 @@ def solve_batch(self, inputs, calculate_sensitivities: bool = False): ax_v.set_ylabel("Voltage [V]") ax_v.set_xlabel("Time [h]") + # Load the full-cell pOCV + procedure = cell.procedure["pOCV"] + experiment = procedure.experiment("Pseudo OCV") + # Plot data for validation - experiment = cell.procedure["pOCV"].experiment("Pseudo OCV") - validation_data = pybop.import_pyprobe_result(experiment) + validation_data = pybop.import_pyprobe_result( + filter_with_preceding_row(procedure, experiment="Pseudo OCV") + ) ax_v.plot( - validation_data["Time [s]"] / 3600, + validation_data["Time [h]"], validation_data["Voltage [V]"], + ls=":", + linewidth=1.4, + color="black", label="Full-cell data", ) + directions = ["charge", "discharge"] params = [] datasets = [] problems = [] - CEs = [] - for i, direction in enumerate(["charge", "discharge"]): + DEs = [] + currents = [] + for i, direction in enumerate(directions): # Load the parameters and experiment with open(cell_path / f"params_{direction}.pickle", "rb") as file: param = pickle.load(file) model_param = param.copy() + # Account for one unexpected point of discharge that appears in the C08 data + step = 1 if cell_label == "C08" and direction == "discharge" else 0 + # Load the relevant branch of the full-cell data dataset = pybop.import_pyprobe_result( smoothing.downsample( - experiment.charge().step(0) + experiment.charge().step(step) if direction == "charge" - else experiment.discharge().step(1 if cell_label == "C08" else 0), + else experiment.discharge().step(step), **smoothing_options, ) ) - # Define a naive SOC based on the measured charge throughput - charge_throughput = ( - dataset["Discharge capacity [A.h]"] - dataset["Discharge capacity [A.h]"][0] - ) + # Ensure that the charge throughput starts from zero + Q_start = -filter_with_preceding_row( + procedure, experiment="Pseudo OCV", phase=direction, step=step + ).get("Capacity [Ah]")[0] + charge_throughput = dataset["Discharge capacity [A.h]"] - Q_start Q_meas = np.abs(charge_throughput[-1]) - # Compute the Coulombic efficiency to take account of dynamic overshoot and make a closed loop - CE = 1.0 + # Compute the dynamic efficiency to take account of voltage overshoot + DE = 1.0 + current = np.mean(dataset["Current [A]"][1:101]) if account_for_overpotential: + rest_phase_voltage = experiment.rest().step(i).get("Voltage [V]") inverse_ocv_interpolant = pybop.Interpolant( dataset["Voltage [V]"], charge_throughput / charge_throughput[-1], name="Inverse pOCV", ) - V_relax = experiment.rest().step(i).get("Voltage [V]")[-1] + V_relax = rest_phase_voltage[-1] + + # Get the local overpotential at the current switch-off and add it + # back onto the rest voltage to make it comparable with pseudo-OCV + V_relax += ( + (dataset["Voltage [V]"][-1] - rest_phase_voltage[1]) + / dataset["Current [A]"][-1] + * current + ) + # Find the scaling needed to make the rest voltage correspond to full dis/charge - CE = np.minimum(inverse_ocv_interpolant(V_relax), 1.0) - Q_cell = CE * Q_meas + DE = np.minimum(inverse_ocv_interpolant(V_relax), 1.0) + Q_scale = DE * Q_meas - # Define state of charge - dataset["Naive SoC"] = (0.0 if i == 0 else 1.0) - charge_throughput / Q_cell + # Define the state of charge based on the measured charge throughput + dataset["SoC"] = ( + 0.0 if direction == "charge" else 1.0 + ) - charge_throughput / Q_scale + dataset["dSoC/dV [V-1]"] = np.gradient(dataset["SoC"], dataset["Voltage [V]"]) # Define optimisation parameters model_param.update( { - "Negative electrode theoretical capacity [A.h]": Q_th_n, - "Positive electrode theoretical capacity [A.h]": Q_th_p, - "Cell capacity [A.h]": Q_cell, - "Minimum negative stoichiometry " + direction: pybop.Parameter( - initial_value=0, bounds=[0, 0.4] + "Cell capacity [A.h]": Q_scale, + f"Minimum negative stoichiometry ({direction})": pybop.Parameter( + initial_value=0.02, bounds=[0, 0.1] + ), + "Maximum negative stoichiometry": sto_n_100, + "Minimum positive stoichiometry": sto_p_100, + f"Maximum positive stoichiometry ({direction})": pybop.Parameter( + initial_value=0.98, bounds=[0.7, 1] ), - "Minimum positive stoichiometry " + direction: pybop.Parameter( - initial_value=0, bounds=[0, 0.4] + f"Cell resistance ({direction}) [Ohm]": pybop.Parameter( + initial_value=0, bounds=[-1, 1] ), - "Series resistance [Ohm]": R0, } ) # Set up an optimisation problem + subset = dataset.get_subset( + np.where(np.logical_and(dataset["SoC"] >= 0, dataset["SoC"] <= 1)) + ) cost = pybop.RootMeanSquaredError( - dataset, target="Voltage [V]", weighting="domain" + subset, target="Voltage [V]", weighting="domain" ) simulator = OCVBalance( - model_param, dataset=dataset, direction=direction, OCP_type=OCP_type + model_param, dataset=subset, direction=direction, OCP_type=OCP_type ) problem = pybop.Problem( simulator=simulator, @@ -207,17 +246,18 @@ def solve_batch(self, inputs, calculate_sensitivities: bool = False): params.append(model_param) datasets.append(dataset) problems.append(problem) - CEs.append(CE) + DEs.append(DE) + currents.append(current) # Optimise the charge and discharge branches together ocv_problem = pybop.MetaProblem(*problems) - options = pybop.SciPyDifferentialEvolutionOptions(maxiter=100) + options = pybop.SciPyDifferentialEvolutionOptions(maxiter=maxiter) optim = pybop.SciPyDifferentialEvolution(ocv_problem, options=options) result = optim.run() - print(result, "\n") + print(result, end="\n\n") - for i, (direction, model_param, dataset, problem, CE) in enumerate( - zip(["charge", "discharge"], params, datasets, problems, CEs, strict=False) + for direction, model_param, dataset, problem, DE, current in zip( + directions, params, datasets, problems, DEs, currents, strict=False ): print("Direction:", direction) @@ -226,64 +266,43 @@ def solve_batch(self, inputs, calculate_sensitivities: bool = False): model_param.update(inputs) # Obtain the parameters - Q_cell, sto_n_0, sto_n_100, sto_p_100, sto_p_0, R0 = ( - problem.simulator.get_model_parameters(inputs) + Q_scale, x_0, x_100, y_100, y_0, Rc = problem.simulator.get_model_parameters( + inputs ) print( "Stoichiometry limits:", - [float(x) for x in [sto_n_0, sto_n_100, sto_p_100, sto_p_0]], + [float(np.round(x * 100, decimals=3)) for x in [x_0, x_100, y_100, y_0]], ) - print("Cell capacity [A.h]:", Q_cell) - print("Coulombic efficiency:", CE) - print("Series resistance [Ohm]:", R0) + print("Cell capacity [A.h]:", Q_scale) + print("pOCV dynamic efficiency:", DE) + print("Cell resistance [Ohm]:", Rc) rmse_mV = problem.evaluate(inputs).values.item() - print("RMSE [mV]:", rmse_mV, "\n") + print("RMSE [mV]:", rmse_mV, end="\n\n") # Optional step to assign and subtract the extra resistance from one or other electrode - positive_resistance, negative_resistance = 0, 0 - if add_extra_resistance: - mean_current = np.mean(problem.simulator.input_data) - if R0 > 0: - positive_resistance = -R0 * mean_current - elif R0 < 0: - negative_resistance = R0 * mean_current + positive_offset, negative_offset = 0, 0 + if apply_offset: + if (direction == "charge" and Rc > 0) or ( + direction == "discharge" and Rc < 0 + ): + positive_offset = Rc * current + elif (direction == "charge" and Rc < 0) or ( + direction == "discharge" and Rc > 0 + ): + negative_offset = -Rc * current # Re-create the pOCP and pOCV functions positive_ocp_function = pybop.Interpolant( model_param[f"Positive electrode {OCP_type} [V]"].x, - model_param[f"Positive electrode {OCP_type} [V]"].y + positive_resistance, - name=f"Positive {OCP_type}", + model_param[f"Positive electrode {OCP_type} [V]"].y - positive_offset, + name="Positive electrode pOCP", ) negative_ocp_function = pybop.Interpolant( model_param[f"Negative electrode {OCP_type} [V]"].x, - model_param[f"Negative electrode {OCP_type} [V]"].y + negative_resistance, - name=f"Negative {OCP_type}", + model_param[f"Negative electrode {OCP_type} [V]"].y - negative_offset, + name="Negative electrode pOCP", ) - # Find the stoichiometry limits and capacity corresponding to the voltage limits - # (rather than the voltage range of the measurement, although should be similar) - measured_ocv_function = OpenCircuitVoltage( - positive_ocp_function, - sto_p_0, - sto_p_100, - negative_ocp_function, - sto_n_0, - sto_n_100, - ) - - # Obtain naive SoC values corresponding to the voltages limits - inverse_measured_ocv = pybop.InverseOCV(measured_ocv_function) - inverse_measured_ocv.parameters["Root"] = pybop.Parameter(initial_value=0.5) - lower_soc = inverse_measured_ocv(param["Lower voltage cut-off [V]"]) - upper_soc = inverse_measured_ocv(param["Upper voltage cut-off [V]"]) - x_0 = sto_n_0 + (sto_n_100 - sto_n_0) * lower_soc - x_100 = sto_n_0 + (sto_n_100 - sto_n_0) * upper_soc - y_100 = sto_p_0 + (sto_p_100 - sto_p_0) * upper_soc - y_0 = sto_p_0 + (sto_p_100 - sto_p_0) * lower_soc - - # Convert from the naive SOC to SOC based on the voltage limits - dataset["SoC"] = (dataset["Naive SoC"] - lower_soc) / (upper_soc - lower_soc) - Q_soc = Q_cell * (upper_soc - lower_soc) ocv_function = OpenCircuitVoltage( positive_ocp_function, y_0, @@ -298,26 +317,25 @@ def solve_batch(self, inputs, calculate_sensitivities: bool = False): ax.plot( (negative_ocp_function.x - x_0) / (x_100 - x_0), negative_ocp_function.y, - label="Negative pOCP, " + direction, + label="Negative electrode pOCP, " + direction, ls="-.", color=colour, ) ax.plot( (positive_ocp_function.x - y_0) / (y_100 - y_0), positive_ocp_function.y, - label="Positive pOCP, " + direction, + label="Positive electrode pOCP, " + direction, ls="--", color=colour, ) # Plot the data and the result - colour = "tab:red" if direction == "charge" else "tab:blue" ax.scatter( dataset["SoC"], dataset["Voltage [V]"], label="pOCV data, " + direction, marker=".", - s=8, + s=4, facecolor="none", edgecolor=colour, ) @@ -326,21 +344,26 @@ def solve_batch(self, inputs, calculate_sensitivities: bool = False): soc, ocv_function(soc), label="pOCV fit, " + direction, + linewidth=0.7, color="black", + alpha=0.8, ) # Plot validation, align the dis/charge branches at the start of the first dis/charge - inverse_ocv_function = pybop.InverseOCV(ocv_function) - Q_meas = Q_soc / CE - offset = experiment.rest().step(i - 1).df["Capacity [Ah]"][-1] / Q_meas + offset = ( + experiment.charge() + if direction == "charge" + else experiment.discharge().step(1 if cell_label == "C08" else 0) + ).get("Capacity [Ah]")[0] + Q_meas = Q_scale / DE + state = (0.0 if direction == "charge" else 1.0) - ( + validation_data["Discharge capacity [A.h]"] + offset + ) / Q_meas ax_v.plot( - validation_data["Time [s]"] / 3600, - ocv_function( - dataset["SoC"][0] - - offset - - validation_data["Discharge capacity [A.h]"] / Q_meas - ), + validation_data["Time [h]"], + ocv_function(state), ls="--", + linewidth=1, color=colour, label=f"Fit to {direction}", ) @@ -354,11 +377,13 @@ def solve_batch(self, inputs, calculate_sensitivities: bool = False): "Maximum negative stoichiometry": x_100, "Minimum positive stoichiometry": y_100, "Maximum positive stoichiometry": y_0, - "Measured cell capacity [A.h]": Q_soc / CE, - "Coulombic efficiency": CE, + "Measured cell capacity [A.h]": Q_meas, + "pOCV constant current [A]": current, + "pOCV dynamic efficiency": DE, "Negative electrode pOCP [V]": negative_ocp_function, "Positive electrode pOCP [V]": positive_ocp_function, "pOCV RMSE [mV]": rmse_mV, + "Reference temperature [K]": 298.15, } ) with open(cell_path / f"params_{direction}.pickle", "wb") as file: @@ -366,9 +391,9 @@ def solve_batch(self, inputs, calculate_sensitivities: bool = False): ax.legend(fontsize="small", bbox_to_anchor=(0.4, 0.4), loc="center left") plt.tight_layout() - fig.savefig(cell_path / f"{cell_label}_pOCV_balance.svg") + fig.savefig(cell_path / f"{cell_label}_pOCV_balance.png") ax_v.legend() - fig_v.savefig(cell_path / f"{cell_label}_pOCV_validation.svg") + fig_v.savefig(cell_path / f"{cell_label}_pOCV_validation.png") if args.show_plots: plt.show() From 3ffe3c202a816c8ec21abbbee641db088f75b56f Mon Sep 17 00:00:00 2001 From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com> Date: Fri, 30 Oct 2026 15:19:04 +0000 Subject: [PATCH 07/18] Update 4_thermal_response.py --- .../4_thermal_response.py | 214 ++++++++++-------- 1 file changed, 121 insertions(+), 93 deletions(-) diff --git a/examples/scripts/dfn_parameterisation/4_thermal_response.py b/examples/scripts/dfn_parameterisation/4_thermal_response.py index e163b7795..2b20f5c2d 100644 --- a/examples/scripts/dfn_parameterisation/4_thermal_response.py +++ b/examples/scripts/dfn_parameterisation/4_thermal_response.py @@ -24,23 +24,28 @@ - {Procedure}.parquet - metadata.json -Our temperature experiment include nested cycles of the following steps: -0. Rest for 10 seconds -1. Charge at C/2 for 12 minutes -2. Rest for 4 hours -3. Charge at 1C for 50 seconds -4. Discharge at 1C for 50 seconds -5. Rest for 1 hour +Our temperature experiment consists of nested cycles of the following steps: +0. Charge at C/2 for 12 minutes +1. Rest for 4 hours +2. Charge at 1C for 50 seconds +3. Discharge at 1C for 50 seconds +4. Rest for 1 hour For each full-cell, estimate the thermal parameters by fitting the temperature experiment. """ parser = argparse.ArgumentParser() +parser.add_argument("--draft", action="store_true", help="Draft mode") 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 + +if args.draft: + plt.rcParams.update({"text.usetex": False}) + maxiter = 5 +else: + maxiter = 250 # Options OCP_type = "pOCP" @@ -63,16 +68,10 @@ procedure = cell.procedure["Temperature"] N_cycles = procedure.experiment("Temperature experiment").cycle_info[0][-1] voltage_points = np.zeros(N_cycles) - T_init = np.zeros(N_cycles) + charge_throughput = np.zeros(N_cycles) for i in range(N_cycles): - dataset = pybop.import_pyprobe_result( - filter_with_preceding_row( - procedure, - experiment="Temperature experiment", - cycle=i, - phase="charge", - step=1, - ), + d = pybop.import_pyprobe_result( + procedure.experiment("Temperature experiment").cycle(i).rest(0), variables=[ "Time [s]", "Current [A]", @@ -88,8 +87,8 @@ "Temperature [degC]", ], ) - voltage_points[i] = dataset["Voltage [V]"][0] - T_init[i] = dataset["Surface temperature [degC]"][0] + voltage_points[i] = d["Voltage [V]"][-1] + charge_throughput[i] = d["Discharge capacity [A.h]"][-1] # Define the thermal model model = pybop.lithium_ion.CellTemperature() @@ -102,30 +101,35 @@ "Positive electrode OCP entropic change [V.K-1]": ( pybamm.Parameter("OCV entropic change [V.K-1]") / 2 ), - "Reference temperature [K]": np.mean(T_init) + kelvin, } ) # Plot the given OCV - fig_v, ax_v = plt.subplots(1, 1, figsize=(5, 3)) + fig_v, ax_v = plt.subplots(1, 1, figsize=(5, 3), dpi=300) ax_v.set_xlabel("State of charge") ax_v.set_ylabel("Voltage [V]") soc = np.linspace(0, 1, 501) ocv_function = get_ocv_function(model_param, OCP_type) - ax_v.plot(soc, ocv_function(soc), label="Pseudo OCV") + ax_v.plot(soc, ocv_function(soc), ls="--", label="Pseudo OCV") # Shift the OCV to match the voltage points - naive_soc = 0.1 + 0.1 * np.arange(N_cycles) - model_param, SOC_vec = shift_ocv_to( + naive_soc = 0.1 + 0.8 * (charge_throughput[0] - charge_throughput) / ( + charge_throughput[0] - charge_throughput[-1] + ) + shift_param, SOC_vec, _ = shift_ocv_to( voltage_points, model_param, naive_soc, OCP_type ) # Validate the alignment by plotting ax_v.scatter(SOC_vec, voltage_points, marker="o", label="Temperature points") - ocv_function = get_ocv_function(model_param) + ocv_function = get_ocv_function(shift_param, OCP_type) ax_v.plot(soc, ocv_function(soc), label="OCV") ax_v.legend() - fig_v.savefig(cell_path / f"{cell_label}_thermal_points.svg") + fig_v.savefig(cell_path / f"{cell_label}_thermal_points.png") + + # Define the OCP functions + for electrode in ["Positive electrode", "Negative electrode"]: + model_param[f"{electrode} OCP [V]"] = shift_param[f"{electrode} {OCP_type} [V]"] # Define parameters that are fixed over state of charge cell_thermal_mass = pybop.Parameter(initial_value=50, bounds=[0, 200]) @@ -191,14 +195,9 @@ param_i["Initial SoC"] = SOC_vec[i] input_data = pybop.generate_consistent_current(dataset) input_data = pybop.downsample_constant_current(dataset) - param_i.update( - { - "Current function [A]": pybop.Interpolant( - input_data["Time [s]"], - input_data["Current [A]"], - name="Current function", - ) - } + print(f"Downsampled {len(dataset)} to {len(input_data)} points") + param_i["Current function [A]"] = pybop.Interpolant( + input_data["Time [s]"], input_data["Current [A]"], name="Current function" ) dataset.control_functions = ["Voltage function [V]"] @@ -218,55 +217,72 @@ # Provide cell-specific initial values if cell_type == "LithiumWerks M1B": - initial_values = [ - 3.43e01, - 5.29e01, - 3.96e-01, - 2.60e-01, - 2.95e-04, - 6.88e-02, - -2.11e-04, - 4.50e-04, - 9.49e-04, - 2.71e-04, - 3.29e-04, - -2.99e-04, - 1.24e-04, - 1.95e-04, + initial_values = [ # C08 + 2.51393602e01, + 3.73928787e01, + 3.75403053e-01, + 3.63267815e-01, + -6.70529360e-04, + 6.30304646e-02, + -5.71735448e-04, + -3.29303289e-04, + -3.47597781e-04, + -2.48047033e-04, + -2.24761317e-04, + -3.29335560e-04, + -1.70693066e-04, + -3.38005422e-04, ] elif cell_type == "Molicel P45B": - initial_values = [ - 1.09e00, - 1.76e02, - 3.42e-02, - 7.35e-01, - -1.27e-03, - 7.24e-02, - -9.58e-04, - -9.29e-04, - -7.29e-04, - -7.04e-04, - -6.06e-06, - -7.19e-04, - -7.34e-04, - -3.13e-04, + initial_values = [ # C07 + 1.39820755e00, + 1.60002272e02, + 3.99926243e-02, + 7.35816442e-01, + -1.27415147e-03, + 6.66433357e-02, + -9.73136847e-04, + -9.02538780e-04, + -7.60357406e-04, + -7.05476160e-04, + -8.03138690e-04, + -7.41881187e-04, + -7.39788016e-04, + -8.47107788e-04, ] elif cell_type == "LG M50 Synthetic": initial_values = [ - 5.30e00, - 4.17e01, - 2.04e-02, - 5.86e-02, - 3.11e-03, - 1.20e-01, - 4.93e-03, - 4.39e-03, - 4.47e-03, - 4.11e-03, - 4.27e-03, - 4.02e-03, - 4.17e-03, - 5.00e-03, + 5.07221910e00, + 4.07944338e01, + 1.56687823e-01, + 5.93269911e-02, + -4.04693011e-04, + 1.19999997e-01, + 5.19878867e-04, + 3.66850038e-04, + 3.99739715e-04, + 2.35419785e-04, + 3.39030068e-04, + 2.79801649e-04, + -1.79491896e-04, + -5.54414458e-06, + ] + elif cell_type == "Enertech SPB655060 Synthetic": + initial_values = [ + 1.80947723e-01, + 4.21659345e01, + 6.30496445e-01, + 2.17994595e-01, + -5.24874801e-04, + 6.00037766e-02, + -4.76336097e-04, + -3.72620475e-04, + -3.06295635e-04, + -2.55021359e-04, + -2.54506310e-04, + -2.09846799e-04, + -1.55028106e-04, + -9.53824722e-05, ] else: initial_values = None @@ -275,15 +291,14 @@ meta_problem = pybop.MetaProblem(*problems, weights=1e-4 * np.ones(len(problems))) if initial_values is not None: meta_problem.parameters.update(initial_values=initial_values) - options = pybop.SciPyMinimizeOptions( - method="trust-constr", maxiter=10 - ) # increase maxiter for accuracy + options = pybop.SciPyMinimizeOptions(method="trust-constr", maxiter=maxiter) optim = pybop.SciPyMinimize(meta_problem, options=options) result = optim.run() print(result) inputs = result.best_inputs init_soc = np.zeros(N_cycles) + init_temperature = np.zeros(N_cycles) entropic_change = np.zeros(N_cycles) rmse_K = np.zeros(N_cycles) for i, (dataset, problem) in enumerate(zip(datasets, problems, strict=False)): @@ -291,31 +306,38 @@ solution = problem.simulate(inputs) # Plot the simulation versus the data - fig, ax = plt.subplots(1, 1, figsize=(5, 3)) + fig, ax = plt.subplots(1, 1, figsize=(5, 3), dpi=300) ax.set_title(f"{cell_type}, Thermal cycle {i}") - t0 = dataset["Time [s]"][0] / 3600 + t0 = dataset["Time [h]"][0] ax.scatter( - dataset["Time [s]"] / 3600 - t0, + dataset["Time [h]"] - t0, dataset["Surface temperature [degC]"], marker=".", - s=4, + s=2, color="tab:red", label="Data", ) ax.plot( solution["Time [s]"].data / 3600 - t0, solution["Surface temperature [K]"].data - kelvin, + color="black", label="Simulation", ) ax.set_xlabel("Time [h]") ax.set_ylabel("Temperature [degC]") ax.legend() plt.grid(True) - fig.savefig(cell_path / f"{cell_label}_thermal_response_{i:02d}.svg") + fig.savefig(cell_path / f"{cell_label}_thermal_response_{i:02d}.png") init_soc[i] = problem.simulator.parameter_values["Initial SoC"] + init_temperature[i] = problem.simulator.parameter_values[ + "Initial temperature [K]" + ] entropic_change[i] = inputs[f"OCV entropic change coefficient {i} [V.K-1]"] - rmse_K[i] = problem.cost.evaluate(solution, inputs).values.item() + cost = pybop.RootMeanSquaredError( + dataset, target="Surface temperature [K]", weighting="domain" + ) + rmse_K[i] = cost.evaluate(solution, inputs).values.item() # Update and save the parameters inputs["OCV entropic change [V.K-1]"] = pybop.Interpolant( @@ -329,20 +351,26 @@ "OCV entropic change [V.K-1]", ]: param.update({key: inputs[key]}) - param["Mean thermal RMSE [K]"] = np.mean(rmse_K) + param.update( + { + "Thermal noise variance [K]": inputs["Sigma for output 1"], + "Thermal test temperature [K]": np.mean(init_temperature), + "Mean thermal RMSE [K]": np.mean(rmse_K), + } + ) with open(cell_path / f"params_{direction}.pickle", "wb") as file: pickle.dump(param, file) # Plot the entropic change versus stoichiometry - fig, ax = plt.subplots(1, 1, figsize=(5, 3)) + fig, ax = plt.subplots(1, 1, figsize=(5, 3), dpi=300) ax.set_title(cell_type) soc = np.linspace(0, 1, 101) - ax.plot(soc, inputs["OCV entropic change [V.K-1]"](soc)) - ax.scatter(init_soc, entropic_change, marker="x") + ax.plot(soc, inputs["OCV entropic change [V.K-1]"](soc) * 1e3) + ax.scatter(init_soc, entropic_change * 1e3, marker="x") # V to mV ax.set_xlabel("State of charge") - ax.set_ylabel("OCV entropic change [V.K-1]") + ax.set_ylabel("OCV entropic change [mV.K-1]") - fig.savefig(cell_path / f"{cell_label}_entropic_coefficient.svg") + fig.savefig(cell_path / f"{cell_label}_entropic_coefficient.png") if args.show_plots: plt.show() From 7337beceda1882b8ef438808091b242d1968beba Mon Sep 17 00:00:00 2001 From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com> Date: Fri, 30 Oct 2026 15:19:53 +0000 Subject: [PATCH 08/18] Update 5_fit_resistances.py --- .../dfn_parameterisation/5_fit_resistances.py | 76 ++++++++++++------- 1 file changed, 48 insertions(+), 28 deletions(-) diff --git a/examples/scripts/dfn_parameterisation/5_fit_resistances.py b/examples/scripts/dfn_parameterisation/5_fit_resistances.py index 33e7aa918..c8b464cbc 100644 --- a/examples/scripts/dfn_parameterisation/5_fit_resistances.py +++ b/examples/scripts/dfn_parameterisation/5_fit_resistances.py @@ -24,15 +24,26 @@ 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). +The full-cell procedure is 9 cycles of the following steps: +0. Charge at C/10 for 1 hour or until 4.2 V +1. Rest for 4 hours +2. Potentiostatic EIS with 0.0042 V amplitude + For each full-cell, estimate the charge transfer and series resistance by fitting the EIS data. """ parser = argparse.ArgumentParser() +parser.add_argument("--draft", action="store_true", help="Draft mode") 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 + +if args.draft: + plt.rcParams.update({"text.usetex": False}) + maxiter = 50 +else: + maxiter = 2500 # Options OCP_type = "pOCP" @@ -53,7 +64,7 @@ def get_model_parameters(self, inputs): a1 = param.evaluate(param["CPE exponent"]) # Get the equivalent capacitance for the CPE - Q1 = C1**a1 / (R1 ** (1 - a1)) + Q1 = C1**a1 / (R1 ** (1 - a1)) if R1 > 0 else np.inf return R0, R1, Q1, a1 def solve_batch(self, inputs, calculate_sensitivities: bool = False): @@ -104,29 +115,28 @@ def solve_batch(self, inputs, calculate_sensitivities: bool = False): N_SOC = len(voltage_points) # Plot the given OCV - fig_v, ax_v = plt.subplots(1, 1, figsize=(5, 3)) + fig_v, ax_v = plt.subplots(1, 1, figsize=(5, 3), dpi=300) 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") + ax_v.plot(soc, ocv_function(soc), ls="--", 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 + shift_param, SOC_vec, offset = shift_ocv_to( + voltage_points, param, naive_soc, OCP_type + ) # Validate the alignment by plotting ax_v.scatter(SOC_vec, voltage_points, marker="o", label="EIS points") - ocv_function = get_ocv_function(shift_param) + ocv_function = get_ocv_function(shift_param, OCP_type) ax_v.plot(soc, ocv_function(soc), label="OCV") ax_v.legend() - fig_v.savefig(cell_path / f"{cell_label}_EIS_points.svg") + fig_v.savefig(cell_path / f"{cell_label}_EIS_points.png") # Set up Nyquist plot - fig, ax = plt.subplots(1, 1, figsize=(5, 3)) + fig, ax = plt.subplots(1, 1, figsize=(5, 3), dpi=300) ax.set_title(f"{cell_type}, Nyquist plot") ax.set_xlabel(r"Re(Z) [m$\Omega$]") ax.set_ylabel(r"-Im(Z) [m$\Omega$]") @@ -143,7 +153,7 @@ def solve_batch(self, inputs, calculate_sensitivities: bool = False): series_resistance = pybop.Parameter( initial_value=initial_series_resistance, bounds=[0, 0.1] ) - dl_capacitance = pybop.Parameter(initial_value=0.1, bounds=[0, 1]) + dl_capacitance = pybop.Parameter(initial_value=0.1, bounds=[0, 10]) # Estimate charge transfer resistances over SOC for fixed R0 problems = [] @@ -176,23 +186,34 @@ def solve_batch(self, inputs, calculate_sensitivities: bool = False): color=colourmap(i / N_SOC), marker=".", s=4, + alpha=0.8, + label=f"_{label}", + ) + ax.scatter( + Z[[idx_fmax, idx_fmin - 1]].real * 1e3, + -Z[[idx_fmax, idx_fmin - 1]].imag * 1e3, + color=colourmap(i / N_SOC), + marker="*", + s=8, + alpha=0.8, + edgecolors="red", + linewidth=0.5, 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, + "Frequency [Hz]": frequencies[idx_fmax:idx_fmin], + "Impedance": Z[idx_fmax:idx_fmin], }, domain="Frequency [Hz]", ) # Define optimisation parameters - initial_cpe_resistance = Z_semicircle[-1].real - initial_series_resistance + initial_cpe_resistance = ( + dataset["Impedance"][-1].real - initial_series_resistance + ) model_param = pybamm.ParameterValues( { "Series resistance [Ohm]": series_resistance, @@ -219,8 +240,8 @@ def solve_batch(self, inputs, calculate_sensitivities: bool = False): problems.append(problem) meta_problem = pybop.MetaProblem(*problems) - options = pybop.SciPyDifferentialEvolutionOptions(maxiter=500) - optim = pybop.SciPyDifferentialEvolution(meta_problem, options=options) + options = pybop.SciPyMinimizeOptions(method="trust-constr", maxiter=maxiter) + optim = pybop.SciPyMinimize(meta_problem, options=options) result = optim.run() print(result) inputs = result.best_inputs @@ -238,8 +259,8 @@ def solve_batch(self, inputs, calculate_sensitivities: bool = False): solution["Impedance"].data.real * 1e3, -solution["Impedance"].data.imag * 1e3, color=colourmap(i / N_SOC), - ls="-", - linewidth=1.5, + alpha=0.8, + linewidth=1, label=label, ) @@ -252,20 +273,19 @@ def solve_batch(self, inputs, calculate_sensitivities: bool = False): ax.axis("equal") # Update and save parameters + Rct_interpolant = pybop.Interpolant(SOC_vec, Rct, name="Charge transfer resistance") param.update( { - "Series resistance [Ohm]": inputs["Series resistance [Ohm]"], - "Charge transfer resistance [Ohm]": pybop.Interpolant( - SOC_vec, Rct, name="Charge transfer resistance" - ), + "Series resistance (ECM) [Ohm]": inputs["Series resistance [Ohm]"], + "Charge transfer resistance [Ohm]": Rct_interpolant, "Double-layer capacitance [F]": inputs["CPE capacitance [F]"], - "Mean EIS RMSE [mOhm]": np.mean(rmse_mOhm), + "Mean high-freq. 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") + fig.savefig(cell_path / f"{cell_label}_EIS_resistances.png") if args.show_plots: plt.show() From 95feda8f057d645cac13fb172ac3e6d45873a364 Mon Sep 17 00:00:00 2001 From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com> Date: Fri, 30 Oct 2026 15:21:38 +0000 Subject: [PATCH 09/18] Update 5b_charge_transfer.py --- .../5b_charge_transfer.py | 83 +++++++++++++++---- 1 file changed, 69 insertions(+), 14 deletions(-) diff --git a/examples/scripts/dfn_parameterisation/5b_charge_transfer.py b/examples/scripts/dfn_parameterisation/5b_charge_transfer.py index 5a8510ab1..36a7ddf02 100644 --- a/examples/scripts/dfn_parameterisation/5b_charge_transfer.py +++ b/examples/scripts/dfn_parameterisation/5b_charge_transfer.py @@ -24,11 +24,17 @@ """ parser = argparse.ArgumentParser() +parser.add_argument("--draft", action="store_true", help="Draft mode") 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 + +if args.draft: + plt.rcParams.update({"text.usetex": False}) + maxiter = 50 +else: + maxiter = 2500 class ChargeTransferResistance(pybop.BaseSimulator): @@ -59,13 +65,13 @@ def solve_batch(self, inputs, calculate_sensitivities: bool = False): * (1 - sto) ** ((1 - alpha) * omega) * sto_e ** (1 - alpha) ) - Rct = Rct_typ / (2 * j0) + Rct = np.divide( + Rct_typ, 2 * j0, out=np.full_like(j0, np.inf), where=j0 != 0 + ) 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) - ) + sol.set_solution_variable("Charge transfer resistance [Ohm]", data=Rct) solutions.append(sol) return solutions @@ -132,32 +138,34 @@ def solve_batch(self, inputs, calculate_sensitivities: bool = False): cost=pybop.WeightedCost(cost, weights=[1e3]), # convert Ohm to mOhm ) - options = pybop.SciPyDifferentialEvolutionOptions(maxiter=250) - optim = pybop.SciPyDifferentialEvolution(problem, options=options) + options = pybop.SciPyMinimizeOptions(method="trust-constr", maxiter=maxiter) + optim = pybop.SciPyMinimize(problem, options=options) result = optim.run() + print(result) + inputs = result.best_inputs # 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) + alpha, omega, sto_e, Rct_typ = simulator.get_model_parameters(inputs) 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") + print(r"RMSE [m$\Omega$]:", result.best_cost, end="\n\n") # Use optimisation result to plot the fitted function - fig, ax = plt.subplots(1, 1, figsize=(5, 3)) + fig, ax = plt.subplots(1, 1, figsize=(5, 3), dpi=300) 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") + Rct_fit = simulator.solve(inputs)["Charge transfer resistance [Ohm]"].data sto_range = np.linspace(0.03, 0.97, 95) simulator.domain_data = sto_range - solution = simulator.solve(result.best_inputs) + solution = simulator.solve(inputs) ax.plot( sto_range, solution["Charge transfer resistance [Ohm]"].data * 1e3, label="BV fit", linewidth=1, + color="black", ) ax.set_title(f"{cell_type}, Butler-Volmer fit") ax.legend() @@ -179,6 +187,7 @@ def solve_batch(self, inputs, calculate_sensitivities: bool = False): zero = np.finfo(np.float64).eps # almost zero param.update( { + "Typical charge transfer resistance [Ohm]": Rct_typ, dominant_electrode + " electrode charge transfer time scale [s]": tau_ct, dominant_electrode + " electrode dimensionless exchange rate": ( BaseGroupedModel.get_multiphase_butler_volmer(dominant_electrode) @@ -193,12 +202,58 @@ def solve_batch(self, inputs, calculate_sensitivities: bool = False): BaseGroupedModel.symmetric_butler_volmer ), other_electrode + " electrode capacitance [F]": zero, + "Charge transfer resistance RMSE [mOhm]": result.best_cost, + "Charge transfer resistance max. error [mOhm]": np.max( + np.abs(Rct_fit - Rct_values) * 1e3 + ), } ) + R0 = param["Series resistance (ECM) [Ohm]"] + for electrode, sign in zip( + ["Positive electrode", "Negative electrode"], [1, -1], strict=False + ): + resistance = R0 / 2 + if electrode.split(" ")[0] == dominant_electrode: + dataset = {"Stoichiometry": param[f"{electrode} pOCP [V]"].x} + Rct_sim = ChargeTransferResistance(param, dataset, dominant_electrode) + Rct_sol = Rct_sim.solve( + {"Initial electrolyte mole fraction": 1.0, **inputs} + ) + resistance += Rct_sol["Charge transfer resistance [Ohm]"].data + param[f"{electrode} OCP [V]"] = pybop.Interpolant( + param[f"{electrode} pOCP [V]"].x, + param[f"{electrode} pOCP [V]"].y + + (sign * param["pOCV constant current [A]"] * resistance), + name=f"{electrode} OCP", + ) 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") + fig.savefig(cell_path / f"{cell_label}_charge_transfer.png") + + """ Without discharge EIS data, mirror the shift in the charge branch. """ + other_direction = "discharge" if direction == "charge" else "charge" + with open(cell_path / f"params_{other_direction}.pickle", "rb") as file: + param = pickle.load(file) + for electrode, sign in zip( + ["Positive electrode", "Negative electrode"], [1, -1], strict=False + ): + resistance = R0 / 2 + if electrode.split(" ")[0] == dominant_electrode: + dataset = {"Stoichiometry": param[f"{electrode} pOCP [V]"].x} + Rct_sim = ChargeTransferResistance(param, dataset, dominant_electrode) + Rct_sol = Rct_sim.solve( + {"Initial electrolyte mole fraction": 1.0, **inputs} + ) + resistance += Rct_sol["Charge transfer resistance [Ohm]"].data + param[f"{electrode} OCP [V]"] = pybop.Interpolant( + param[f"{electrode} pOCP [V]"].x, + param[f"{electrode} pOCP [V]"].y + + (sign * param["pOCV constant current [A]"] * resistance), + name=f"{electrode} OCP", + ) + with open(cell_path / f"params_{other_direction}.pickle", "wb") as file: + pickle.dump(param, file) if args.show_plots: plt.show() From fba33ff3123c0f9e1da04a788432e9f6d4d3335d Mon Sep 17 00:00:00 2001 From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com> Date: Fri, 30 Oct 2026 17:10:03 +0000 Subject: [PATCH 10/18] Test pipeline in draft mode --- tests/examples/test_examples.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/tests/examples/test_examples.py b/tests/examples/test_examples.py index 114adfb38..8fa11879f 100644 --- a/tests/examples/test_examples.py +++ b/tests/examples/test_examples.py @@ -52,7 +52,9 @@ def test_example_scripts(self, example: Path, tmp_path): ) def test_pipeline_scripts(self, list_of_pipeline_scripts: list[Path], tmp_path): # The pipeline scripts must be run sequentially + sys.argv.append("--draft") for example in sorted(list_of_pipeline_scripts): new_example_path = tmp_path / example.name shutil.copy(example, new_example_path) runpy.run_path(new_example_path) + sys.argv.pop(-1) From c0c1946a7241c6d34f8eb1de7c97022a08055e7d Mon Sep 17 00:00:00 2001 From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com> Date: Fri, 30 Oct 2026 18:13:12 +0000 Subject: [PATCH 11/18] Update CHANGELOG.md --- CHANGELOG.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index caf6a26cf..2a54a119c 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,7 @@ ## Features +- [#996](https://github.com/pybop-team/PyBOP/pull/996) - Updates the DFN parameterisation pipeline, adds draft mode for testing. - [#995](https://github.com/pybop-team/PyBOP/pull/995) - Updates the cell models, including adding the thermodynamic factor to the GroupedSPMe. - [#976](https://github.com/pybop-team/PyBOP/pull/976) - Updates experimental procedures and synthetic validation data. - [#991](https://github.com/pybop-team/PyBOP/pull/991) - Adds support for Python 3.15 (EP-BOLFI optimiser and PyProBE still restricted to Python 3.12 or below). From 82c2069dc66b792b55f531bda7324230bd4d857f Mon Sep 17 00:00:00 2001 From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com> Date: Fri, 30 Oct 2026 20:14:02 +0000 Subject: [PATCH 12/18] Update test tolerances --- tests/unit/test_application_utils.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/unit/test_application_utils.py b/tests/unit/test_application_utils.py index 0dba5197c..1f4e5845c 100644 --- a/tests/unit/test_application_utils.py +++ b/tests/unit/test_application_utils.py @@ -84,8 +84,8 @@ def test_shift_ocv_to(self, parameter_values): lower_voltages, parameter_values, soc ) shifted_ocv = get_ocv_function(shifted_param) - np.testing.assert_allclose(offset, 0.1, atol=1e-6) - np.testing.assert_allclose(soc, soc_points, atol=1e-8) + np.testing.assert_allclose(offset, 0.1, atol=5e-5) + np.testing.assert_allclose(soc, soc_points, atol=1e-6) np.testing.assert_allclose(shifted_ocv(soc), lower_voltages) def test_make_voltage_monotonic(self): From ab9ea4925bb1c2cceda67ac92a3aa2a867a13e95 Mon Sep 17 00:00:00 2001 From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com> Date: Fri, 30 Oct 2026 20:22:00 +0000 Subject: [PATCH 13/18] Update test tolerance --- tests/unit/test_application_utils.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/unit/test_application_utils.py b/tests/unit/test_application_utils.py index 1f4e5845c..14ce80acf 100644 --- a/tests/unit/test_application_utils.py +++ b/tests/unit/test_application_utils.py @@ -78,14 +78,14 @@ def test_get_ocv_function(self, parameter_values): def test_shift_ocv_to(self, parameter_values): ocv_function = get_ocv_function(parameter_values, OCP_type="OCP") - soc = np.linspace(0, 1, 10) + soc = np.linspace(0, 1, 11) lower_voltages = ocv_function(soc) - 0.1 shifted_param, soc_points, offset = shift_ocv_to( lower_voltages, parameter_values, soc ) shifted_ocv = get_ocv_function(shifted_param) np.testing.assert_allclose(offset, 0.1, atol=5e-5) - np.testing.assert_allclose(soc, soc_points, atol=1e-6) + np.testing.assert_allclose(soc, soc_points, atol=5e-5) np.testing.assert_allclose(shifted_ocv(soc), lower_voltages) def test_make_voltage_monotonic(self): From 8783e94d634a53fddee916ca5325c64173a64363 Mon Sep 17 00:00:00 2001 From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com> Date: Fri, 30 Oct 2026 21:08:11 +0000 Subject: [PATCH 14/18] Update test tolerances --- tests/unit/test_application_utils.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/unit/test_application_utils.py b/tests/unit/test_application_utils.py index 14ce80acf..6ab190533 100644 --- a/tests/unit/test_application_utils.py +++ b/tests/unit/test_application_utils.py @@ -84,9 +84,9 @@ def test_shift_ocv_to(self, parameter_values): lower_voltages, parameter_values, soc ) shifted_ocv = get_ocv_function(shifted_param) - np.testing.assert_allclose(offset, 0.1, atol=5e-5) + np.testing.assert_allclose(offset, 0.1, atol=5e-4) np.testing.assert_allclose(soc, soc_points, atol=5e-5) - np.testing.assert_allclose(shifted_ocv(soc), lower_voltages) + np.testing.assert_allclose(shifted_ocv(soc), lower_voltages, atol=1e-6) def test_make_voltage_monotonic(self): # Test ascending From d951f518fe86d476d08c7d7feb704354fe21b8cb Mon Sep 17 00:00:00 2001 From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com> Date: Fri, 30 Oct 2026 21:16:36 +0000 Subject: [PATCH 15/18] Update test tolerances --- tests/unit/test_application_utils.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/unit/test_application_utils.py b/tests/unit/test_application_utils.py index 6ab190533..c792d9807 100644 --- a/tests/unit/test_application_utils.py +++ b/tests/unit/test_application_utils.py @@ -84,8 +84,8 @@ def test_shift_ocv_to(self, parameter_values): lower_voltages, parameter_values, soc ) shifted_ocv = get_ocv_function(shifted_param) - np.testing.assert_allclose(offset, 0.1, atol=5e-4) - np.testing.assert_allclose(soc, soc_points, atol=5e-5) + np.testing.assert_allclose(offset, 0.1, atol=1e-3) + np.testing.assert_allclose(soc, soc_points, atol=1e-3) np.testing.assert_allclose(shifted_ocv(soc), lower_voltages, atol=1e-6) def test_make_voltage_monotonic(self): From 9472647b9b190fe72d3fafa39413e77b7e70d239 Mon Sep 17 00:00:00 2001 From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com> Date: Fri, 30 Oct 2026 21:49:04 +0000 Subject: [PATCH 16/18] Update test tolerance --- tests/unit/test_application_utils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/unit/test_application_utils.py b/tests/unit/test_application_utils.py index c792d9807..2c5e1cef5 100644 --- a/tests/unit/test_application_utils.py +++ b/tests/unit/test_application_utils.py @@ -86,7 +86,7 @@ def test_shift_ocv_to(self, parameter_values): shifted_ocv = get_ocv_function(shifted_param) np.testing.assert_allclose(offset, 0.1, atol=1e-3) np.testing.assert_allclose(soc, soc_points, atol=1e-3) - np.testing.assert_allclose(shifted_ocv(soc), lower_voltages, atol=1e-6) + np.testing.assert_allclose(shifted_ocv(soc), lower_voltages, atol=5e-4) def test_make_voltage_monotonic(self): # Test ascending From ed9774d0c1b1b3a61d06dabfcf2d989661676c5a Mon Sep 17 00:00:00 2001 From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com> Date: Sun, 1 Nov 2026 16:18:10 +0000 Subject: [PATCH 17/18] Add capacity column to EIS data --- .../Full cell/EIS/EIS charge.parquet | Bin 124623 -> 124551 bytes pybop/pybamm/synthetic_utils.py | 18 +++++++++++------- 2 files changed, 11 insertions(+), 7 deletions(-) diff --git a/examples/data/LG M50 Synthetic/Full cell/EIS/EIS charge.parquet b/examples/data/LG M50 Synthetic/Full cell/EIS/EIS charge.parquet index d345765d5ce48ede5f3510b0f89d793090d7c6cc..979c9dfe2503fa6a77a8572ed4f4145b8495c99c 100644 GIT binary patch delta 12330 zcmeI&Vc_b95 zH@N(9bWO)){_R9%Z2TIp$cC{dBX1r50jAnE3^%>6QoD)#-4BpVdgoqzi{`S8zVww7 zIfD$MVI&N9=JTY@TLUk7>Um$$g&=C2_aQ^EePoU))hFBc)OA~_3vDC~w8m^We{JJR zb4=5_`kut6eZs`j4WXGoe0+Z?Eve2SvYvkEu)@QEbLUa>q@v81ujT%i3d%az0B|$O zL_lL}UI)sM$1f&+)#Tp!m3xr$C5>xGays+?G3BDkyzFT0kFDBLAQRf3Qg;$glHh9| ze3k-`xPSVSTzjlaFfy!___6&fAL(ej$q4`V=K^WGT!L&hq^~|dCm;Zxnep{Ex$bW* zn$_Nxq{7MZp<I7LZoRlBiGi?)>%9x=BpK4$ANA20myvkZBoGoq-O*jL6 zLo_|_?lum-_2y4S32SZK-l+-CeZ)tivWq@sUL$V!Kr=wFcZTq|QLO^ZDt-=b>?o_KNy^Fu@t zavJIaTo!USj^-_kf5Y(|KIEzSEpf3rK{;b%#W(c<%j*D2kHNAB1pl;6XbOxWXIfZ28eesMeDnF_r+ z^tZ|cshlwF>DrFtXTD=c_tFkRRsfOuDw0Z>X4RfwqZ@QmXldf44`~jGa#lt-5qS0s zy(L)(ASXGez!M#YahkUz`gFelBQ>(3Ce7X%Vc=7>%qZ6Ctl{zZvj$=3YzB|zxUFZ_Vk)Hq!qcsPBlY3VT`>_f=J%cTcERO) z$zl}u1s-}GZw%}*lX9JPrK`7afI-if2j0YQG(D&gdjotwN>}DpNwyI_=!J+E%X%0R z=%sYBHu%{RBRPx28U=6FnjNn((LcAU`~1So1O5mi@_f+=dieRZN4_E~@BHpfK4K2BpmAwAcIRo2=b(xGwnNC>&OcgS(}^+IAE$4|8@_ zHK@GHkr}<@JWosf@4eLg93!#9TV18ePn}IYmhLs_GpKxEcx;FX%SIu}35xBo0diVz zOkZfwxv2%wQ;`(#x_+5MtOWAjU4C=1%PQ6s2vfDG82gv=$Vn@OchnymmJ6!`k=443BoQoCD zodk*0D2sp8yiCm9o+@#6@Coj+4xMu9QHqcaDch`)LOYv{KKJ2G{sIx1OVSAL%%qH| z&eQk?#>0;O-ucP}Uvjw>o;hh9kh$Sv%V;kCR`y(DJ}ldGM1z9_F>DiIal3AS*@uj! zB+hYdg)d4Xy4p>C(n7{+hpWuf!=RU)h2H^i>5hv_aL;juk3~IPG}BshrdbSkgWez3n8hdxJyb^h zOwEa~e&MurUno0lAwab085!-PENgzZ@Xb-y(*qy-q z(mE__QLD?mb5`#C2XS14{-0GrN6|DfK7>ZUENH@A z|4bS_vxufTQM16AQon1q{-b`l*qde109qDN$GV;!6}@y1104zBch}3ljBkx3e&iVY zy=>1f5VC&`?$)4nAQfl`&t!?Ki-UPD5z?cpm?i0h5n2ZCo#PxC$L~CBsn*Lyh;ND1 zr|vLJLmd0zEm_(U-<+YL87+Hy4UwBb`}%(ph$pQ4lf-+8E*|P__VLV1)ygMZc;T3%H_z z0!9|Nqb|y=_u1LU^#jW$iR6{Vg=R2tqi3nXtUyuOT#A#>t>Epip|TecM-D0c++ODLBl_WHR*uU)7+9j1Koz zAA2QI#;VjoOrL>?jc0O0-Mz+OrSkaMxib;Br0v>&^<9WBCWP?&suJC6VN@9&PwPHy zqhj?9RhttHv1j=QX0tKWNryos(;E)VN zp-kT``KP#+6qz~IhL#LN1f$W2MKA3XA#n+@mY-bhX@D?T%SMsZyGB_FkW2f>Qm@7?MqusDY8p#(LgRvD81c}dVWd08mt#R;WYU5?pI{jqNegp!jtJ{FLv(GOHzF69P z@TboqJCD%t6QU;{%V7Ea%IcIeoQ8X~SdHe6=5j^Tc`!JZ2GgKlyb5C>|FO3P7n$yj zhobZ<%YF#_KHc&Fv2<3Xpi<158UdIy`*bl0Uen zGb z3UV7~ImE;8?awFhjdiItl^f4LziDx;_ZLZwyi-Ph$2jw0xX)LM^$FEi(e^qNC}Z-z z1CaqQ5sILdu_rwU{d;y5lN^4{^>zgvR)qojVx0J*iRnCva6it_hjsH3cfabOdT zToi)@)YK$KwmKpBd$2w$Kx4f;V|hZ8$BsL=7UjR(qYXn)9Bmcg$kX0gGQPyczisj4 zclMLv3G>v-4@UBI4u2mLX%m~ z!7}9qb=Oyfq_9vt#&*{AeCsaso2)o@oWM4yvt#w~u%plyb&>DTy{!W1eg%;Of79m_ zB^vKY&!$b{ksPh60MDaJh{Ze7X@XSJe0>>ecoLaLf#idJ96&@l@$)3WQg~{}R#Ns{ zul*;;yW-JZTOzlX(%`12w|DRNtcX-emWzQ6(&4oT=TY|z&r~WsM(x8@#s$n```pIM zTOWamsbTH$go#-@cog*E7-9`eI-F41{Sk)tBh%c2xPn@qy~EhVTT`!-LEcVszZXUn?}(|d#MJexu^ z#ys9s{FQSZ`M9KSCq3i~t6JNK%>cgx1GmjCt_*g}CmGgVgKZlQcg*fE?u&O^F<6gP zY&x9onaNwtf%uiAACr!BeanZF9;&JuMYWv?AiheWW1l-n_3w}lkZ?a_Tn{hlhO;`w<(95BH)B$BVO^C^ekKy_il3@6XM&Tev{Gr z$}MGp>t%EEIw$q)ZQ85jArA0EYzWu(&&-=YKsQ=K86;?;QjG#Eec_f zFKAAY$M~v#x~O`0QS&5ucS?sYah9{iUmbD$vJ0pXXrdF-!J}e&Auzitw=LE2vK(^ zSoix_R}&W{Jz-)ZsM3@&a&c}SNN$-DV`3r#kav*3Zl?xqwWCc;LT12Jt_MH4!A`9i zMXspzxfY(Lv(|K& zM3dTW5ji?!eoW+ziCKc4S|psA1$Ua5H1XXxQ{3^GsJ+AOO0rye0-3`rp;Ke5MLs6z7=`>C)HH*I2>iSC;yE`LkK*#V*0C5s$ z^+24u!G@GttsR>zlRl-0BdlF((`Avv5*GY{HzVkdq5Q9z&PTE$Y}_B3F8cBa+BYzF z)Ah=2lNPMsQ}_gpoQXb8m|e*_^@v8jIG7u#ad{5-faiMQ5kO{#B(hvs4^#_2Og*ax z@Y3Frkwp_fWKu)`h@%$RTUfq*6xqfoGzy55RW&BSP!gbv5av9>;<9FYu2?OE(g;`{6z*ZPca=qL zS1?a3rMC1ITr}ySY{!s2ELYnv{Q6QClj98JqW=Oq@-2e7Xu_Aa7(bt1E!xEfLTJ-E zd@hPp#swpFZvHN7GtyYHgdz5J8_?_aHH;!~q1*^_ap4WH&{}D3!fQvWzP@rIkbXt@ zI0lUk_~xQRTJ7>Gf1187_f(W}-WAAPy0AWbL{wDuTa=Q2UbFNZ^WuTVgU0s@J6!iB z_dCu#kmd;D4A$y}9r?9l==vVE^ubRAIOKl>E zU5v8{ST(u#YY)Ieisqe54RHnacYbQ zfzyNtaNzJ~C=+;qPVt9Cex~3)Hvze6M#0seKfB7A`fe^vly=O%6N>C@Vi;Y{Z8v)x zh}U@Bxno1+P9nDG_IuF|9wU$FzOz|&x}4;QwN3Zrbm^f6tIvYD4ifOn5LaT#dv9hc(BT2i$fgh9u*-tj`q(lLgE$f2npv^So&jlbYl0+Cdk2^gC zkwP*H1F35uBEVP>?%>#UrW3i&U&VlZP+(q7fC}&oihGn#D1}6+T%S)d1deNiHnGU; zB_J*3wVB@t`i-F92>Okn-w67Rpx+4kjiBEM`i-F92>Okn-w67Rpx+4kjiBEM`i-F9 z2>Okn-w67Rp#Q%Tbmmzg`~M1hZx!*ip#Sw}IMQ(z*nljY1(F~|F~AWZ0124`s+e9g zI~)K2P}9LT(+FUzm7ssqwXnms{}Zlx9X1L1x5=A_z*?044cEX9Tm4Vi#1g_LJO8uk z`qr?6@|ub{peqm@zDok5f)8&Zx#ocn!P(b|1|I(d=`jyf1lRwa27e&i=YfyFMlfH|YIm5&A1J9a8QK zkn`Vh1Rb>Kx~9AS;rj1uG;Wa4A^(q;NDc4&dui|WQfUBuy@?h+DZ>tir(tRYmw_p` OU=l(AKuhk4;Qs(ouwc0W delta 12426 zcmeI2S5#A7`|U|+fzUz;y>|gAQWX)US5Z1BT{?m^>7ghfO+X}wNC)Ym2k9LIks4a) zy%P{3RSti@_x;YfIhW^LakIymd#*Xwv&O#J`!{#Q0O8;OVfmkLG)XnMG~*QjL7MRg zG&nQ_Gz3hNiA3-x^8hFr4quJ%rG>gPsWqTzXdIqd3J`&<5F?5Eh#yOVy=SC#2r8m< zD1Tx{St@Qzc|M;;NEirXTiIznbIAp(>YXZ`JZ|vWzH!ho5d^svJnji7vQ;yjPd!c@__<7z0(~ zV;ut}BR|4;lNJZ;ehHw8D zH8KHG`y-LFZ_c|hjsO@w`$o?WyTnp+{Ip%0Fe`@1;jvS}X2+E86fPfbr@bJ`{1Dc1 z8zS+XF~^5)rZra*x!#dt`(*=ONS-s5O{wdCE)r4XE?uCH#LfYQYPX-(~s`)$$!}Sn&6DTG7G{8 zw%mQEnx;H>$PqdPS`D>Ja7wby+6aQu-?J{KB5FB&9%4)4)-qy%!*fS!-*B8i-i91# z{#My#NX=oXxN$K}ge=qNW3}}dy=}Bgf8D(o-?9!i?jS@5Dekg-d76nW{c*-GJOmws zLG2%`J}ed8pqPKmz8T_T*7M%o!HU#^`jXF_j!o$U&IM8BJymV|!HfdBL`y2&BVIDS zY)_EfrSBxR`epj@dM_+X;4cQq_D@!V}c$r2R{FBbstACG*mvo6Zs zQKO@MRyMM5U(F1XJ3+yfo~sE-O2+H{z@?*;=ijre=I*Y$+|cycN{)$BPS^KN06seNF((9rKU%hVrOH%2a8d! zZORFp`l(b9$7oNvrvB0B4N6A2j>!H)U2e=l?98#bB* zU;5yek^}N@UDlPHknE}}K2Yu;{Q{$*W zVFg&wPZYfbhfM#F%RrxsPEBD?36c4%XZeG6F;-ey!Cp&!v>zp(Sn}G6RZ_DbDYK4O zEDr!lzLCE-X6W2s&f%FHapOqr>we9EYuTLXcdO zA&YL)z=^&~DV9lhj-Uyc@X1TRy^UR!@Sb)(f%!{SGvPV>H@k;6nN-~NX3u{v>h$pd z9N?d+mw8a~$93R1?;J?o9@t{oxrZ1u+dIJUq7tuU4n1ma!A~i^W zm%lG)V5_!P#mCv!XExb638V*#HtS#BhVqowOgx=R)#uMTHl}ZIn9aKP6zUy<1E{ zS)`*TI7`SQLXnFeX^vW#0P8&}|DEanD*A&yN?+K{OGER_Q$HoyEXDdm`qpP&oA9pK zICWgmOR#Z7(5(U+SneC{qWu$fJQf1|y?j6%9c=>cZf?_aE(On{`Z)Sj`*Bo$dWXVU z-jQqFGR8eRMBo*cQ94&z#9i3)$s8*yG*;kx zT)(0-2$1ENV1X4D0`V#ExtT{+0*D7sx8kh#W%9&cxN@mcxADE6*06{>3r->3i?(bL zmvo!;q_j;h12+}%0^?`^pHuExk=6Io$yG5!o}WCm!=EOnFPRMPT$>~X&(o$*Al?V> zR99ees-h?P4^z&KvS`bO7EKu|M<0UTe%!?NVjpTf!h<H&id z?hEz49_1rY%gCIvKCz;!&UoS@$~Ev43GYHF=CkW=kniCP&cVv>mkiQI05uF>6UKdi zDnB;-&FxRJ<_D}FG=XkCj=%P0&5_qVO}lgIB-qe!A!CWkb8Z%yK#5&Kj`zb%N|$~k z(R%DWe!a(84~iJH=cfqyeo8YhBxLXFKo1l}oF-qfJ;N@2$#3}anE)v}xkT2no>z>F zMw{x1Nz?HhPU$+*rBSjP4-d+wv90+^H`#RUF|rgm@v_HvWKO1m%;ii0*z-k5Tw}g1 zF06P0V+HmtZsnG1eHS8OIo74!7I9~cg${ke1pADW z{B~&JqmZK!ThQjDy4UkOIl;XrEYF;F{R(8Bk_@;f*-HOrSi61*FHGYgo(mJYOj%~} zZi7&?Iy8Rpg%8<`($qw^AnEs=a5W9G_u+P{Efcul;q54_-*$c>!(ek;5XY#|J z!th5LY-%Zt14UZEn+1%-bJGnoYlXdS$!EgcO4T7Y>Q{+j(bgZ2PlQ1%-*w1RF8Xni zR%Gl**f)h@g7Kphct?6I%*x@6mRjP`10L6_Y~Khk_kbmHd41H&5dVvQG9XOaZ}94n zfv>PVdDGn8p@GmC38N1p&h92Y-AW|*{-MqpA;%t?(~aFj-ST}Ksd%@sUWv3;%%-y# zCFp$jhg@_0m*B{XlGWlDM#93a-p$i%v2hb#M<{@L+K7q1d1)9i?9;EDW0Y*Ht&zc$ z;kiy``9@f8F)T`TwzwX_;`XdHFJiwq%soa#pEO7!L}5Xmr5nK#CBoCyak0h0#?V*D zhma)xnH7>6R1;#T-#rRJ$Zau;q?MV+=5i1@6TUO#g)Fh`ruf#3Fy%~-br!=rD+9DV z%m&^poK))_qBK)sB@*iMh?mn2uZ(XkptPbwVspPCywwvhA!COqZ9gq~)p@gB8@znr z=1?)Y>H_H~1a=&vEIKBDGGsMk&Ks#Xb;R4`7!v3-NfAN75@t)yWz=IgGktl{wFATdy zcgbQ&U@MZTl4;FnYIvd)Kl?aI8iHSm%20`a)^yUJkbtD`ZIaw6MJjzs7^0+ut*KwfpOP}M;^93jI|Tp5*hWEu6Y;J~Qx zPpaZ?;ZB^>8ub|->etu_p%!D8Ld76b&NM(i)$a+MwWJJQADEU$V-$4b+$rrvRn*Kr z!dE?9vMU^R>G4!Kzx8B7p07ngrMCt8h&=Iv@2b94J(yK;&77*j)1*1+V!Gf(F>G+< zVw>~&JI*S}oVW0IKvX8IsXF`<*!a|o3a*()F@yNB#3WU)ba5un{&VvA76S4to3C;T zv;B31B}pg8oVQnB)HDP;p~ENt#i7$&<-u8sPWSVlI49vwVe82;BT0=9$z3s4#D6PguJp#jABqoB(b+;xAq9l>< zU*@L1SY;*0Khs4?`en+i69yHTJk9g&J2iDECE-#Q*|fu@|5T+Y1UrQMDC6 zKAtete2RD>-qLK6iTCNn@`NF2@KEZ|wOez?g!)=qN&*SeF3yM_nTpu$Oc>fmXi2P@ zo2%JR;ai+;3y42bq5F+MjHtp?NK_r-mH!;myrrRMSTM+A)4M$Wm@4%CDML$h5hfxg zUlYg?&8es#$9lb011$Dn$)GJR2eSks=WS`W`lVHiv#=Gb*@{3+3*Mo1F^V|}R$F*v z$Y+?V{cdGs2eU08g_0??q3q!BQTZI!N(yacXfh00IvEjrcI(etugd}LvO?-zT(rr0 zfsKZ>wGZAOBby+)tb++d*M4gEGM!}Gx2Ad=y5%o;eNc5n#<01A9G691#Vry%Nn>Rk zGhmK=tLbX5o9Cma`ph8i<#cTxE!Jvn*H&?0$>UDmyW+ABm+AJ02*`*Z8F#+PvwgHn zrtwZk_8r3`F52(%h4fPZawZW80as;6IdJjGnHA5tpIjf%@hf8&8nMvhdps2= zNdL0QIpco2SUBe0qTOS&?~nb42lu%CEK{!8MYtPV1uy%8(73hzQOKW3t?e4uP zh3*kI|JA`YuVEZ5>b-cEIaU5w(rS9fwiKPI$GbO?1ZnxHrp_#*zfuJg3~KV zhlIq8zIrf7?*sYq=|jf zi8oaoNcyR2OCVdQ7^~3%5<;TFp+y{cuMNaE=QD-OA|Uw`aO=W*>^y;V`lObc{(U{_ zH(FOkQ7sS{&=7@b3Mt!08{S{%uA4_9mMkJ2C1CU2MGrUoDVNd4q3y;Z)y)GqCF)&a zXjiErV{u`p6YnsyXH_&}^W#fx=&sX6v?+LjZtTR>wGER;`2~g!8B7D}=D6NtJ3x!hcGTM4cWIc}$(tO00PQ>LP} z+EuOTR?+fB09r^Qqm_CIy;oR9EAo2C!{Bfoy;n5jxmWcul1*Nz988vD{S32~6%K0Yo|Tkoj4 zY?*8*$%&HatrXF452LRT`$?!|?BvMWD;B$P6X%6Dsd-(**0S`T-AGAb6=0}?m=H$t zP#*0IMhgD~OcKN7O^8A!3fMuwJ9b>p_4syU0K*ZUlon~oKUY_-DzUmR_BVb$U8#Fu zibosZq4lwyCC6Km{r#D-e+1*FlSHQ?>UWY|&UHI<@MDq)oX|Jw7IE6T-!^{*AW0_QFE7j6Qt1G0085xYPM-LC%C#Da>%4Ha4TO~tY8HR0~zNKKh z6}!c~XQ6GXT)--*@r2Ti^OG{EyhJ{C{;`%LSDonFDNQ3!V996r@EyP@UfK|mEd1&YfF8Ln z4u}H+6UX3M3jjuBq79-( z4Rybv?l;u^hPvNS_Z#YdL)~wv`wex!q3$=-{f4^VQ1=_^enZ`FsQV3dzoG6o)cuCK z-%$7e_tf3_|0i{a&rJb%{)^b-YDxc#*x`KBfCjkIG=Lm_Fbx<6UjK*Qfd;5qT%&j> z4h{}I6SU+rA&z;c0gM6^Ua+G~-17I2&_BG_A`2h`M}Ynw5XWHH4=b@BTM* zy*o{HCDM5Y;06FCVaRCcpp%>Mm05rysN`BPKvTBhB69#GQ1@Roz6B4Q1E_)){;H&H z_`n=Mjp6KCF+s<+$!O$pc#2X9nf}`=zz?PXRR2Bd|IsGUk|!VhPbg{UZzz5c{(2rD z@K3wINCB<<>yG`$jU0s!&I4}!@6Lf7`t#anfF|w!O>e-t7XUo}R1i$$8-KmA|9Ay? zaGwQ$5MTp--AFJ|LPHP#2KW93({_oO;DV;t8~>-MV1lImwa5Nrht`lX!C6ZI{Qou= zjF8%E|9a~Ezw!1BGN%9Q!=$Af`MY-P?^ pl.LazyFrame: """Create a PyProBE-compatible dataframe for one EIS sweep.""" row_count = len(frequencies) + initial_capacity, initial_voltage = initial_point return pl.LazyFrame( { "Time [s]": [None] * row_count, "Step": [step_number] * row_count, "Event": [event_number] * row_count, - "Initial voltage [V]": [initial_voltage] * row_count, + "Initial voltage [V]": [float(initial_voltage)] * row_count, "Frequency [Hz]": frequencies, "Impedance (real) [Ohm]": np.real(impedance), "Impedance (imag) [Ohm]": np.imag(impedance), "Current [A]": [None] * row_count, "Voltage [V]": [None] * row_count, - "Capacity [Ah]": [None] * row_count, + "Capacity [Ah]": [float(initial_capacity)] * row_count, } ) @@ -351,17 +352,20 @@ def _solve_eis_experiment( eis_model, f_eval=frequencies, parameter_values=parameter_values ) impedance = np.asarray(eis_simulator.solve()["Impedance"].data) - initial_voltage = ( - None + initial_point = ( + (None, None) if latest_solution is None - else float(latest_solution["Terminal voltage [V]"].data[-1]) + else ( + latest_solution["Discharge capacity [A.h]"].data[-1] * -1, + latest_solution["Terminal voltage [V]"].data[-1], + ) ) sweep_frame = _build_eis_sweep_dataframe( step_number=step, event_number=step_offset + i, frequencies=frequencies, impedance=impedance, - initial_voltage=initial_voltage, + initial_point=initial_point, ) sweep_frames.append(sweep_frame) From 9861dab701e9c4279c7f69b8d511c2c160ea402f Mon Sep 17 00:00:00 2001 From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com> Date: Sun, 1 Nov 2026 21:51:57 +0000 Subject: [PATCH 18/18] Update 4_thermal_response.py --- .../4_thermal_response.py | 68 ++++++++++--------- 1 file changed, 36 insertions(+), 32 deletions(-) diff --git a/examples/scripts/dfn_parameterisation/4_thermal_response.py b/examples/scripts/dfn_parameterisation/4_thermal_response.py index 2b20f5c2d..89c6bb896 100644 --- a/examples/scripts/dfn_parameterisation/4_thermal_response.py +++ b/examples/scripts/dfn_parameterisation/4_thermal_response.py @@ -136,7 +136,7 @@ surface_thermal_mass = pybop.Parameter(initial_value=50, bounds=[0, 200]) cell_heat_transfer = pybop.Parameter(initial_value=0.25, bounds=[0, 1]) surface_heat_transfer = pybop.Parameter(initial_value=0.25, bounds=[0, 1]) - noise_variance = pybop.Parameter(initial_value=9e-2, bounds=[6e-2, 12e-2]) + noise_variance = pybop.Parameter(initial_value=9e-2, bounds=[3e-2, 12e-2]) datasets = [] problems = [] @@ -191,7 +191,7 @@ }, ) - # Set the initial SoC + # Set the initial SoC and the protocol param_i["Initial SoC"] = SOC_vec[i] input_data = pybop.generate_consistent_current(dataset) input_data = pybop.downsample_constant_current(dataset) @@ -199,11 +199,15 @@ param_i["Current function [A]"] = pybop.Interpolant( input_data["Time [s]"], input_data["Current [A]"], name="Current function" ) - dataset.control_functions = ["Voltage function [V]"] + param_i["Voltage function [V]"] = pybop.Interpolant( + dataset["Time [s]"], dataset["Voltage [V]"], name="Voltage function" + ) # Set up an optimisation problem simulator = pybop.pybamm.Simulator( - model, parameter_values=param_i, protocol=dataset + model, + parameter_values=param_i, + protocol=(input_data["Time [s]"], dataset["Time [s]"]), ) cost = pybop.GaussianLogLikelihood( dataset, @@ -252,37 +256,37 @@ ] elif cell_type == "LG M50 Synthetic": initial_values = [ - 5.07221910e00, - 4.07944338e01, - 1.56687823e-01, - 5.93269911e-02, - -4.04693011e-04, - 1.19999997e-01, - 5.19878867e-04, - 3.66850038e-04, - 3.99739715e-04, - 2.35419785e-04, - 3.39030068e-04, - 2.79801649e-04, - -1.79491896e-04, - -5.54414458e-06, + 5.07188089e00, + 4.07949072e01, + 1.56556488e-01, + 5.93242898e-02, + -4.04751633e-04, + 1.19999087e-01, + 5.19896327e-04, + 3.67652039e-04, + 4.00373608e-04, + 2.36077551e-04, + 3.39727570e-04, + 2.80494186e-04, + -1.80232589e-04, + -6.01453524e-06, ] elif cell_type == "Enertech SPB655060 Synthetic": initial_values = [ - 1.80947723e-01, - 4.21659345e01, - 6.30496445e-01, - 2.17994595e-01, - -5.24874801e-04, - 6.00037766e-02, - -4.76336097e-04, - -3.72620475e-04, - -3.06295635e-04, - -2.55021359e-04, - -2.54506310e-04, - -2.09846799e-04, - -1.55028106e-04, - -9.53824722e-05, + 4.02487923e-02, + 4.22616038e01, + 6.70578752e-01, + 2.18039161e-01, + -5.28150497e-04, + 3.00461298e-02, + -4.73122270e-04, + -3.73105240e-04, + -3.08568057e-04, + -2.51359780e-04, + -2.56819393e-04, + -2.21431954e-04, + -1.61524511e-04, + -8.42979972e-05, ] else: initial_values = None