Implement GITT example

CHANGELOG.md

@@ -2,6 +2,7 @@
 
 ## Features
 
+- [#249](https://github.com/pybop-team/PyBOP/pull/249) - Add WeppnerHuggins model and GITT example.
 - [#304](https://github.com/pybop-team/PyBOP/pull/304) - Decreases the testing suite completion time.
 - [#301](https://github.com/pybop-team/PyBOP/pull/301) - Updates default echem solver to "fast with events" mode.
 - [#251](https://github.com/pybop-team/PyBOP/pull/251) - Increment PyBaMM > v23.5, remove redundant tests within integration tests, increment citation version, fix examples with incorrect model definitions.

examples/scripts/gitt.py

@@ -0,0 +1,70 @@
+import numpy as np
+
+import pybop
+
+# Define model
+parameter_set = pybop.ParameterSet.pybamm("Xu2019")
+model = pybop.lithium_ion.SPM(
+    parameter_set=parameter_set, options={"working electrode": "positive"}
+)
+
+# Generate data
+sigma = 0.005
+t_eval = np.arange(0, 150, 2)
+values = model.predict(t_eval=t_eval)
+corrupt_values = values["Voltage [V]"].data + np.random.normal(0, sigma, len(t_eval))
+
+# Form dataset
+dataset = pybop.Dataset(
+    {
+        "Time [s]": t_eval,
+        "Current function [A]": values["Current [A]"].data,
+        "Voltage [V]": corrupt_values,
+    }
+)
+
+# Define parameter set
+parameter_set.update(
+    {
+        "Reference OCP [V]": 4.1821,
+        "Derivative of the OCP wrt stoichiometry [V]": -1.38636,
+    },
+    check_already_exists=False,
+)
+
+# Define the cost to optimise
+model = pybop.lithium_ion.WeppnerHuggins(parameter_set=parameter_set)
+
+parameters = [
+    pybop.Parameter(
+        "Positive electrode diffusivity [m2.s-1]",
+        prior=pybop.Gaussian(5e-14, 1e-13),
+        bounds=[1e-16, 1e-11],
+        true_value=parameter_set["Positive electrode diffusivity [m2.s-1]"],
+    ),
+]
+
+problem = pybop.FittingProblem(
+    model,
+    parameters,
+    dataset,
+    signal=["Voltage [V]"],
+)
+
+cost = pybop.RootMeanSquaredError(problem)
+
+# Build the optimisation problem
+optim = pybop.Optimisation(cost=cost, optimiser=pybop.PSO, verbose=True)
+
+# Run the optimisation problem
+x, final_cost = optim.run()
+print("Estimated parameters:", x)
+
+# Plot the timeseries output
+pybop.quick_plot(problem, parameter_values=x, title="Optimised Comparison")
+
+# Plot convergence
+pybop.plot_convergence(optim)
+
+# Plot the parameter traces
+pybop.plot_parameters(optim)

pybop/__init__.py

@@ -46,8 +46,7 @@
 from ._utils import is_numeric
 
 #
-# Cost class
-# Problem class
+# Problem classes
 #
 from .problems.base_problem import BaseProblem
 from .problems.fitting_problem import FittingProblem

pybop/models/lithium_ion/__init__.py

@@ -2,4 +2,4 @@
 # Import lithium ion based models
 #
 from .base_echem import EChemBaseModel
-from .echem import SPM, SPMe, DFN, MPM, MSMR
+from .echem import SPM, SPMe, DFN, MPM, MSMR, WeppnerHuggins

pybop/models/lithium_ion/base_echem.py

@@ -73,16 +73,24 @@
         """
         parameter_set = parameter_set or self._parameter_set
 
-        electrode_params = [
-            (
-                "Negative electrode active material volume fraction",
-                "Negative electrode porosity",
-            ),
-            (
-                "Positive electrode active material volume fraction",
-                "Positive electrode porosity",
-            ),
-        ]
+        if self.pybamm_model.options["working electrode"] == "positive":
+            electrode_params = [
+                (
+                    "Positive electrode active material volume fraction",
+                    "Positive electrode porosity",
+                ),
+            ]
+        else:
+            electrode_params = [
+                (
+                    "Negative electrode active material volume fraction",
+                    "Negative electrode porosity",
+                ),
+                (
+                    "Positive electrode active material volume fraction",
+                    "Positive electrode porosity",
+                ),
+            ]
 
         related_parameters = {
             key: inputs.get(key) if inputs and key in inputs else parameter_set[key]

pybop/models/lithium_ion/echem.py

@@ -1,6 +1,7 @@
 import pybamm
 
 from .base_echem import EChemBaseModel
+from .weppner_huggins import BaseWeppnerHuggins
 
 
 class SPM(EChemBaseModel):
@@ -270,3 +271,50 @@
             spatial_methods=spatial_methods,
             solver=solver,
         )
+
+
+class WeppnerHuggins(EChemBaseModel):
+    """
+    Represents the Weppner & Huggins model to fit diffusion coefficients to GITT data.
+
+    Parameters
+    ----------
+    name: str, optional
+        A name for the model instance, defaults to "Weppner & Huggins model".
+    parameter_set: pybamm.ParameterValues or dict, optional
+        A dictionary or a ParameterValues object containing the parameters for the model. If None, the default parameters are used.
+    geometry: dict, optional
+        A dictionary defining the model's geometry. If None, the default geometry is used.
+    submesh_types: dict, optional
+        A dictionary defining the types of submeshes to use. If None, the default submesh types are used.
+    var_pts: dict, optional
+        A dictionary specifying the number of points for each variable for discretization. If None, the default variable points are used.
+    spatial_methods: dict, optional
+        A dictionary specifying the spatial methods for discretization. If None, the default spatial methods are used.
+    solver: pybamm.Solver, optional
+        The solver to use for simulating the model. If None, the default solver is used.
+    """
+
+    def __init__(
+        self,
+        name="Weppner & Huggins model",
+        parameter_set=None,
+        geometry=None,
+        submesh_types=None,
+        var_pts=None,
+        spatial_methods=None,
+        solver=None,
+    ):
+        self.pybamm_model = BaseWeppnerHuggins()
+        self._unprocessed_model = self.pybamm_model
+
+        super().__init__(
+            model=self.pybamm_model,
+            name=name,
+            parameter_set=parameter_set,
+            geometry=geometry,
+            submesh_types=submesh_types,
+            var_pts=var_pts,
+            spatial_methods=spatial_methods,
+            solver=solver,
+        )

pybop/models/lithium_ion/weppner_huggins.py

@@ -0,0 +1,100 @@
+#
+# Weppner Huggins Model
+#
+import numpy as np
+import pybamm
+
+
+class BaseWeppnerHuggins(pybamm.lithium_ion.BaseModel):
+    """WeppnerHuggins Model for GITT.
+
+    Parameters
+    ----------
+    name : str, optional
+        The name of the model.
+    """
+
+    def __init__(self, name="Weppner & Huggins model"):
+        super().__init__({}, name)
+
+        pybamm.citations.register("""
+            @article{Weppner1977,
+            title={{Determination of the kinetic parameters
+            of mixed-conducting electrodes and application to the system Li3Sb}},
+            author={Weppner, W and Huggins, R A},
+            journal={Journal of The Electrochemical Society},
+            volume={124},
+            number={10},
+            pages={1569},
+            year={1977},
+            publisher={IOP Publishing}
+            }
+        """)
+
+        # `self.param` is a class containing all the relevant parameters and functions for
+        # this model. These are purely symbolic at this stage, and will be set by the
+        # `ParameterValues` class when the model is processed.
+        self.options["working electrode"] = "positive"
+
+        t = pybamm.t
+        ######################
+        # Parameters
+        ######################
+
+        d_s = pybamm.Parameter("Positive electrode diffusivity [m2.s-1]")
+
+        c_s_max = pybamm.Parameter(
+            "Maximum concentration in positive electrode [mol.m-3]"
+        )
+
+        i_app = self.param.current_density_with_time
+
+        U = pybamm.Parameter("Reference OCP [V]")
+
+        U_prime = pybamm.Parameter("Derivative of the OCP wrt stoichiometry [V]")
+
+        epsilon = pybamm.Parameter("Positive electrode active material volume fraction")
+
+        r_particle = pybamm.Parameter("Positive particle radius [m]")
+
+        a = 3 * (epsilon / r_particle)
+
+        l_w = self.param.p.L
+
+        ######################
+        # Governing equations
+        ######################
+        u_surf = (
+            (2 / (np.pi**0.5))
+            * (i_app / ((d_s**0.5) * a * self.param.F * l_w))
+            * (t**0.5)
+        )
+        # Linearised voltage
+        V = U + (U_prime * u_surf) / c_s_max
+        ######################
+        # (Some) variables
+        ######################
+        self.variables = {
+            "Voltage [V]": V,
+            "Time [s]": t,
+        }
+
+    @property
+    def default_geometry(self):
+        return {}
+
+    @property
+    def default_submesh_types(self):
+        return {}
+
+    @property
+    def default_var_pts(self):
+        return {}
+
+    @property
+    def default_spatial_methods(self):
+        return {}
+
+    @property
+    def default_solver(self):
+        return pybamm.DummySolver()