Update pipeline scripts

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).

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()

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}")
 

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,26 +75,30 @@
         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}",
         )
 
         # Define initial state of lithiation and direction
         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()