Update grouped models

CHANGELOG.md

@@ -2,6 +2,7 @@
 
 ## Features
 
+- [#974](https://github.com/pybop-team/PyBOP/pull/974) - Adds voltage components to each grouped model as well as asymmetric and multiphase Butler-Volmer kinetics.
 - [#969](https://github.com/pybop-team/PyBOP/pull/969) - Updates synthetic data and adds example script for thermal parameterisation.
 - [#965](https://github.com/pybop-team/PyBOP/pull/965) - Adds synthetic data and example scripts for OCV parameterisation.
 - [#963](https://github.com/pybop-team/PyBOP/pull/963) - Adds an example for generating synthetic data from a specification and exporting it to a PyProBE-compatible parquet file.

examples/notebooks/battery_parameterisation/lgm50_pulse_validation.ipynb

@@ -223,7 +223,7 @@
    "id": "12",
    "metadata": {},
    "source": [
-    "We can construct the model and parameter values for a two-RC circuit model with initial values as listed below. Note, the \"Initial SOC\" is shifted slightly to better match the zero degree data."
+    "We can construct the model and parameter values for a two-RC circuit model with initial values as listed below. Note, the \"Initial SoC\" is shifted slightly to better match the zero degree data."
    ]
   },
   {

examples/scripts/battery_parameterisation/gitt_fitting.py

@@ -47,7 +47,7 @@
     )
 
 # Group the parameters
-grouped_parameter_values = pybop.lithium_ion.SPDiffusion.create_grouped_parameters(
+grouped_parameter_values = pybop.li_half_cell.SPDiffusion.create_grouped_parameters(
     parameter_values
 )
 
@@ -60,11 +60,13 @@
 gitt_parameter_data = gitt_fit()
 
 # Plot the functional parameters
-pybop.plot.dataset(gitt_parameter_data, signal=["Particle diffusion time scale [s]"])
+pybop.plot.dataset(
+    gitt_parameter_data, signal=["Positive particle diffusion time scale [s]"]
+)
 pybop.plot.dataset(gitt_parameter_data, signal=["Series resistance [Ohm]"])
 
 # Run the identified model
-identified_model = pybop.lithium_ion.SPDiffusion(build=True)
+identified_model = pybop.li_half_cell.SPDiffusion(build=True)
 grouped_parameter_values.update(gitt_fit.best_inputs)
 grouped_parameter_values["Current function [A]"] = pybamm.Interpolant(
     dataset["Time [s]"], dataset["Current [A]"], pybamm.t
@@ -76,7 +78,7 @@
 # Return to the original model and update the diffusivity value
 diffusivity = np.mean(
     parameter_values["Positive particle radius [m]"] ** 2
-    / gitt_parameter_data["Particle diffusion time scale [s]"]
+    / gitt_parameter_data["Positive particle diffusion time scale [s]"]
 )
 parameter_values.update({"Positive particle diffusivity [m2.s-1]": diffusivity})
 

examples/scripts/battery_parameterisation/gitt_pulse.py

@@ -30,7 +30,7 @@
 )
 
 # Group the parameters
-grouped_parameter_values = pybop.lithium_ion.SPDiffusion.create_grouped_parameters(
+grouped_parameter_values = pybop.li_half_cell.SPDiffusion.create_grouped_parameters(
     parameter_values
 )
 

examples/scripts/comparison_examples/gitt_models.py

@@ -2,6 +2,7 @@
 import pybamm
 
 import pybop
+from pybop.models.li_half_cell import SPDiffusion, WeppnerHuggins
 
 # Define model and parameter values
 model_options = {"working electrode": "positive"}
@@ -30,13 +31,13 @@
     }
 )
 
-for model in [pybop.lithium_ion.WeppnerHuggins(), pybop.lithium_ion.SPDiffusion()]:
+for model in [WeppnerHuggins(), SPDiffusion()]:
     # GITT target parameter
     diffusion_parameter = pybop.Parameter(pybop.Gaussian(5000, 1000))
-    if isinstance(model, pybop.lithium_ion.WeppnerHuggins):
+    if isinstance(model, WeppnerHuggins):
         # Group parameter values
-        grouped_parameter_values = (
-            pybop.lithium_ion.WeppnerHuggins.create_grouped_parameters(parameter_values)
+        grouped_parameter_values = WeppnerHuggins.create_grouped_parameters(
+            parameter_values
         )
 
         # We can fit only the duration of the pulse
@@ -49,7 +50,7 @@
                 dataset["Discharge capacity [A.h]"][-1]
                 - dataset["Discharge capacity [A.h]"][0]
             )
-            * (grouped_parameter_values["Theoretical electrode capacity [A.s]"] / 3600)
+            * grouped_parameter_values["Theoretical electrode capacity [A.h]"]
         )
         grouped_parameter_values.update(
             {
@@ -61,7 +62,7 @@
         # Fitting parameters
         grouped_parameter_values.update(
             {
-                "Particle diffusion time scale [s]": diffusion_parameter,
+                "Positive particle diffusion time scale [s]": diffusion_parameter,
                 "Reference voltage [V]": pybop.Parameter(
                     initial_value=grouped_parameter_values["Reference voltage [V]"],
                 ),
@@ -70,14 +71,14 @@
 
     else:
         # Group parameter values
-        grouped_parameter_values = (
-            pybop.lithium_ion.SPDiffusion.create_grouped_parameters(parameter_values)
+        grouped_parameter_values = SPDiffusion.create_grouped_parameters(
+            parameter_values
         )
 
         # Fitting parameters
         grouped_parameter_values.update(
             {
-                "Particle diffusion time scale [s]": diffusion_parameter,
+                "Positive particle diffusion time scale [s]": diffusion_parameter,
                 "Series resistance [Ohm]": pybop.Parameter(
                     initial_value=grouped_parameter_values["Series resistance [Ohm]"],
                 ),
@@ -87,7 +88,7 @@
     # Build the problem
     gitt_dataset = (
         dataset.get_subset(pulse_index)
-        if isinstance(model, pybop.lithium_ion.WeppnerHuggins)
+        if isinstance(model, WeppnerHuggins)
         else dataset
     )
     simulator = pybop.pybamm.Simulator(
@@ -97,13 +98,15 @@
     problem = pybop.Problem(simulator, cost)
 
     # Build the optimisation problem
-    optim = pybop.SciPyMinimize(problem)
+    options = pybop.SciPyMinimizeOptions(method="Nelder-Mead")
+    optim = pybop.SciPyMinimize(problem, options)
 
     # Run the optimisation problem
     result = optim.run()
     print(result)
     print(
-        "Diffusion time [s]:", result.best_inputs["Particle diffusion time scale [s]"]
+        "Diffusion time [s]:",
+        result.best_inputs["Positive particle diffusion time scale [s]"],
     )
 
     # Plot the timeseries output

examples/scripts/dfn_parameterisation/3_cell_balancing.py

@@ -354,7 +354,7 @@ 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.s]": Q_soc / CE * 3600,
+                "Measured cell capacity [A.h]": Q_soc / CE,
                 "Coulombic efficiency": CE,
                 "Negative electrode pOCP [V]": negative_ocp_function,
                 "Positive electrode pOCP [V]": positive_ocp_function,

pybop/__init__.py

@@ -56,6 +56,7 @@
 # Model classes
 #
 from .models import lithium_ion
+from .models import li_half_cell
 from .models._exponential_decay import ExponentialDecayModel
 from .models.lithium_ion.utils import Interpolant, InverseOCV
 

pybop/applications/gitt_methods.py

@@ -36,15 +36,17 @@ def __init__(
     ):
         self.parameter_values = parameter_values.copy()
         self.parameters = {
-            "Particle diffusion time scale [s]": pybop.Parameter(bounds=[0, np.inf]),
+            "Positive particle diffusion time scale [s]": pybop.Parameter(
+                bounds=[0, np.inf]
+            ),
             "Series resistance [Ohm]": pybop.Parameter(bounds=[0, np.inf]),
         }
         self.cost = cost or pybop.RootMeanSquaredError
         self.optimiser = optimiser or pybop.SciPyMinimize
         self.optimiser_options = optimiser_options or self.optimiser.default_options()
 
         # Create model
-        self.model = pybop.lithium_ion.SPDiffusion()
+        self.model = pybop.li_half_cell.SPDiffusion()
         self.problem = None
 
     def __call__(
@@ -110,7 +112,9 @@ def __init__(
         self.optimiser_options = optimiser_options or self.optimiser.default_options()
 
         # Set up OCV root-finding function
-        self.inverse_ocp = pybop.InverseOCV(parameter_values["Electrode OCP [V]"])
+        self.inverse_ocp = pybop.InverseOCV(
+            parameter_values["Positive electrode OCP [V]"]
+        )
 
         # Initialise single pulse fitter
         self.pulse_fit = GITTPulseFit(
@@ -148,7 +152,9 @@ def __call__(self) -> pybop.Dataset:
                 # Log the result
                 self.pulses.append(pulse_result)
                 diffusion_time.append(
-                    pulse_result.best_inputs["Particle diffusion time scale [s]"]
+                    pulse_result.best_inputs[
+                        "Positive particle diffusion time scale [s]"
+                    ]
                 )
                 series_resistance.append(
                     pulse_result.best_inputs["Series resistance [Ohm]"]
@@ -169,14 +175,16 @@ def __call__(self) -> pybop.Dataset:
         self.parameter_data = pybop.Dataset(
             {
                 "Stoichiometry": np.asarray(stoichiometry),
-                "Particle diffusion time scale [s]": np.asarray(diffusion_time),
+                "Positive particle diffusion time scale [s]": np.asarray(
+                    diffusion_time
+                ),
                 "Series resistance [Ohm]": np.asarray(series_resistance),
                 cost_name: np.asarray(best_cost),
             }
             if len(stoichiometry) > 1 and stoichiometry[-1] > stoichiometry[0]
             else {
                 "Stoichiometry": np.flipud(np.asarray(stoichiometry)),
-                "Particle diffusion time scale [s]": np.flipud(
+                "Positive particle diffusion time scale [s]": np.flipud(
                     np.asarray(diffusion_time)
                 ),
                 "Series resistance [Ohm]": np.flipud(np.asarray(series_resistance)),
@@ -187,8 +195,8 @@ def __call__(self) -> pybop.Dataset:
 
         # Compute mean values
         self.best_inputs = {
-            "Particle diffusion time scale [s]": np.mean(
-                self.parameter_data["Particle diffusion time scale [s]"]
+            "Positive particle diffusion time scale [s]": np.mean(
+                self.parameter_data["Positive particle diffusion time scale [s]"]
             ),
             "Series resistance [Ohm]": np.mean(
                 self.parameter_data["Series resistance [Ohm]"]

pybop/models/alternative_functions.py

@@ -0,0 +1,39 @@
+import pybamm
+
+""" Alternative functions written as classes to allow pickling. """
+
+
+class FunctionalDiffusionTime:
+    def __init__(self, r2_scale, D, c_scale):
+        self.r2_scale = r2_scale
+        self.D = D
+        self.c_scale = c_scale
+
+    def __call__(self, sto, T):
+        r2_scale, D, c_scale = self.r2_scale, self.D, self.c_scale
+        return r2_scale / D(sto * c_scale, T)
+
+
+class AsymmetricButlerVolmer:
+    def __init__(self, alpha):
+        self.alpha = alpha  # cathodic transfer coefficient
+
+    def __call__(self, sto_surf, sto_e, eta_RT_F):
+        alpha = self.alpha
+        j0 = sto_surf**alpha * (sto_e * (1 - sto_surf)) ** (1 - alpha)
+        return j0 * (pybamm.exp((1 - alpha) * eta_RT_F) - pybamm.exp(-alpha * eta_RT_F))
+
+
+class MultiphaseButlerVolmer:
+    def __init__(self, alpha, omega):
+        self.alpha = alpha
+        self.omega = omega
+
+    def __call__(self, sto_surf, sto_e, eta_RT_F):
+        alpha, omega = self.alpha, self.omega
+        j0 = (
+            sto_surf ** (alpha * omega)
+            * (1 - sto_surf) ** ((1 - alpha) * omega)
+            * sto_e ** (1 - alpha)
+        )
+        return j0 * (pybamm.exp((1 - alpha) * eta_RT_F) - pybamm.exp(-alpha * eta_RT_F))

pybop/models/li_half_cell/__init__.py

@@ -0,0 +1,5 @@
+#
+# Import lithium half cell models
+#
+from .sp_diffusion import SPDiffusion
+from .weppner_huggins import WeppnerHuggins

pybop/models/li_half_cell/base_model.py

@@ -0,0 +1,112 @@
+import pybamm
+from pybamm import FunctionParameter, Parameter
+from pybamm import lithium_ion as pybamm_lithium_ion
+
+from pybop.models.lithium_ion.utils import InverseOCV
+
+
+class BaseHalfCellModel(pybamm_lithium_ion.BaseModel):
+    """
+    A base model for PyBOP's grouped-parameter lithium-ion half-cell models.
+
+    Parameters
+    ----------
+    name : str, optional
+        The name of the model.
+    **model_kwargs : optional
+        Valid PyBaMM model option keys and their values, for example:
+        options : dict, optional
+            A dictionary of options to customise the behaviour of the PyBaMM model.
+        build : bool, optional
+            If True, the model is built upon creation (default: False).
+    """
+
+    def __init__(self, name="Base Model", **model_kwargs):
+        super().__init__(name=name, **model_kwargs)
+
+    def build_model(self):
+        """
+        Build model variables and equations
+        Credit: PyBaMM
+        """
+        self._build_model()
+
+        self._built = True
+        pybamm.logger.info(f"Finish building {self.name}")
+
+    @staticmethod
+    def set_initial_state(
+        initial_value,
+        parameter_values,
+        direction=None,
+        param=None,
+        inplace=True,
+        options=None,
+        inputs=None,
+        tol=1e-6,
+    ):
+        """
+        Set the value of the initial state of charge.
+
+        Parameters
+        ----------
+        initial_value : float
+            Target initial value.
+            If float, interpreted as stoichiometry, must be between 0 and 1.
+            If string e.g. "4 V", interpreted as voltage, must be between V_min and V_max.
+        parameter_values : :class:`pybamm.ParameterValues`
+            Parameters and their corresponding values.
+        param : :class:`pybamm.LithiumIonParameters`, optional
+            The symbolic parameter set to use for the simulation.
+            If not provided, the default parameter set will be used.
+        inplace: bool, optional
+            If True, replace the parameters values in place. Otherwise, return a new set of
+            parameter values. Default is True.
+        options : dict-like, optional
+            A dictionary of options to be passed to the model, see
+            :class:`pybamm.BatteryModelOptions`.
+        inputs : dict, optional
+            A dictionary of input parameters to pass to the model when solving.
+        tol : float, optional
+            The tolerance for the solver used to compute the initial stoichiometries.
+            A lower value results in higher precision but may increase computation time.
+            Default is 1e-6.
+        """
+        parameter_values = parameter_values if inplace else parameter_values.copy()
+
+        if isinstance(initial_value, str) and initial_value.endswith("V"):
+            V_init = float(initial_value[:-1])
+            V_min = parameter_values.evaluate(
+                Parameter("Lower voltage cut-off [V]"), inputs=inputs
+            )
+            V_max = parameter_values.evaluate(
+                Parameter("Upper voltage cut-off [V]"), inputs=inputs
+            )
+
+            if not V_min - tol <= V_init <= V_max + tol:
+                raise ValueError(
+                    f"Initial voltage {V_init}V is outside the voltage limits ({V_min}, {V_max})."
+                )
+
+            def ocv_function(sto_p):
+                U_p = FunctionParameter(
+                    "Positive electrode OCP [V]",
+                    {"Positive particle stoichiometry": sto_p},
+                )
+                return parameter_values.evaluate(U_p, inputs=inputs).squeeze()
+
+            inverse_ocv = InverseOCV(ocv_function)
+            sto_p = inverse_ocv(V_init)
+
+        elif isinstance(initial_value, int | float):
+            sto_p = initial_value
+
+        else:
+            raise ValueError("Initial value must be a float or a string ending in 'V'.")
+
+        if not 0 <= sto_p <= 1:
+            raise ValueError("Initial stoichiometry should be between 0 and 1.")
+
+        parameter_values["Initial stoichiometry"] = sto_p
+
+        return parameter_values