diff --git a/CHANGELOG.md b/CHANGELOG.md index 8261b0641..72df32aa5 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. - [#998](https://github.com/pybop-team/PyBOP/pull/998) - Add `remove_voltage_events` option, `default_variables` property and the option to pass both `t_eval` and `t_interp` to the PyBaMM simulators. - [#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. 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 d345765d5..979c9dfe2 100644 Binary files a/examples/data/LG M50 Synthetic/Full cell/EIS/EIS charge.parquet and b/examples/data/LG M50 Synthetic/Full cell/EIS/EIS charge.parquet differ 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() 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}") 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() 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() diff --git a/examples/scripts/dfn_parameterisation/4_thermal_response.py b/examples/scripts/dfn_parameterisation/4_thermal_response.py index e163b7795..89c6bb896 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,37 +101,42 @@ "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]) 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 = [] @@ -187,24 +191,23 @@ }, ) - # 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) - 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" + ) + param_i["Voltage function [V]"] = pybop.Interpolant( + dataset["Time [s]"], dataset["Voltage [V]"], name="Voltage function" ) - dataset.control_functions = ["Voltage function [V]"] # 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, @@ -218,55 +221,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.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 = [ + 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 @@ -275,15 +295,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 +310,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 +355,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() 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() 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() 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/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/pybop/pybamm/synthetic_utils.py b/pybop/pybamm/synthetic_utils.py index 170eea8d1..ce5c2479b 100644 --- a/pybop/pybamm/synthetic_utils.py +++ b/pybop/pybamm/synthetic_utils.py @@ -285,22 +285,23 @@ def _build_eis_sweep_dataframe( event_number: int, frequencies: np.ndarray, impedance: np.ndarray, - initial_voltage: float | None, + initial_point: tuple, ) -> 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) 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) diff --git a/tests/unit/test_application_utils.py b/tests/unit/test_application_utils.py index 0f44c4ead..2c5e1cef5 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, 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=1e-3) + np.testing.assert_allclose(soc, soc_points, atol=1e-3) + np.testing.assert_allclose(shifted_ocv(soc), lower_voltages, atol=5e-4) 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" ) 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."