Update calculation of cyclable lithium capacity (#349)

* Update calculation of cyclable lithium capacity

* Add test

* Increase coverage

* Add formation_concentrations flag

CHANGELOG.md

@@ -25,6 +25,7 @@
 
 ## Bug Fixes
 
+- [#339](https://github.com/pybop-team/PyBOP/issues/339) - Updates the calculation of the cyclable lithium capacity in the spme_max_energy example.
 - [#387](https://github.com/pybop-team/PyBOP/issues/387) - Adds keys to ParameterSet and updates ECM OCV check.
 - [#380](https://github.com/pybop-team/PyBOP/pull/380) - Restore self._boundaries construction for `pybop.PSO`
 - [#372](https://github.com/pybop-team/PyBOP/pull/372) - Converts `np.array` to `np.asarray` for Numpy v2.0 support.

examples/scripts/spme_max_energy.py

@@ -16,7 +16,7 @@
 # print(f"Optimised volumetric energy density: {final_cost:.2f} Wh.m-3")
 
 # Define parameter set and model
-parameter_set = pybop.ParameterSet.pybamm("Chen2020")
+parameter_set = pybop.ParameterSet.pybamm("Chen2020", formation_concentrations=True)
 model = pybop.lithium_ion.SPMe(parameter_set=parameter_set)
 
 # Fitting parameters

pybop/models/lithium_ion/base_echem.py

@@ -282,16 +282,6 @@ def approximate_capacity(self, x):
         None
             The nominal cell capacity is updated directly in the model's parameter set.
         """
-        # Extract stoichiometries and compute mean values
-        (
-            min_sto_neg,
-            max_sto_neg,
-            min_sto_pos,
-            max_sto_pos,
-        ) = self._electrode_soh.get_min_max_stoichiometries(self._parameter_set)
-        mean_sto_neg = (min_sto_neg + max_sto_neg) / 2
-        mean_sto_pos = (min_sto_pos + max_sto_pos) / 2
-
         inputs = {
             key: x[i] for i, key in enumerate([param.name for param in self.parameters])
         }
@@ -302,6 +292,16 @@ def approximate_capacity(self, x):
             self._parameter_set
         )
 
+        # Extract stoichiometries and compute mean values
+        (
+            min_sto_neg,
+            max_sto_neg,
+            min_sto_pos,
+            max_sto_pos,
+        ) = self._electrode_soh.get_min_max_stoichiometries(self._parameter_set)
+        mean_sto_neg = (min_sto_neg + max_sto_neg) / 2
+        mean_sto_pos = (min_sto_pos + max_sto_pos) / 2
+
         # Calculate average voltage
         positive_electrode_ocp = self._parameter_set["Positive electrode OCP [V]"]
         negative_electrode_ocp = self._parameter_set["Negative electrode OCP [V]"]

pybop/parameters/parameter_set.py

@@ -2,7 +2,7 @@
 import types
 from typing import List
 
-from pybamm import ParameterValues, parameter_sets
+from pybamm import LithiumIonParameters, ParameterValues, parameter_sets
 
 
 class ParameterSet:
@@ -174,14 +174,16 @@ def is_json_serializable(self, value):
             return False
 
     @classmethod
-    def pybamm(cls, name):
+    def pybamm(cls, name, formation_concentrations=False):
         """
         Retrieves a PyBaMM parameter set by name.
 
         Parameters
         ----------
         name : str
             The name of the PyBaMM parameter set to retrieve.
+        set_formation_concentrations : bool, optional
+            If True, re-calculates the initial concentrations of lithium in the active material (default: False).
 
         Returns
         -------
@@ -194,4 +196,44 @@ def pybamm(cls, name):
         if name not in list(parameter_sets):
             raise ValueError(msg)
 
-        return ParameterValues(name).copy()
+        parameter_set = ParameterValues(name).copy()
+
+        if formation_concentrations:
+            set_formation_concentrations(parameter_set)
+
+        return parameter_set
+
+
+def set_formation_concentrations(parameter_set):
+    """
+    Compute the concentration of lithium in the positive electrode assuming that
+    all lithium in the active material originated from the positive electrode.
+
+    Parameters
+    ----------
+    parameter_set : pybamm.ParameterValues
+        A PyBaMM parameter set containing standard lithium ion parameters.
+    """
+    # Obtain the total amount of lithium in the active material
+    Q_Li_particles_init = parameter_set.evaluate(
+        LithiumIonParameters().Q_Li_particles_init
+    )
+
+    # Convert this total amount to a concentration in the positive electrode
+    c_init = (
+        Q_Li_particles_init
+        * 3600
+        / (
+            parameter_set["Positive electrode active material volume fraction"]
+            * parameter_set["Positive electrode thickness [m]"]
+            * parameter_set["Electrode height [m]"]
+            * parameter_set["Electrode width [m]"]
+            * parameter_set["Faraday constant [C.mol-1]"]
+        )
+    )
+
+    # Update the initial lithium concentrations
+    parameter_set.update({"Initial concentration in negative electrode [mol.m-3]": 0})
+    parameter_set.update(
+        {"Initial concentration in positive electrode [mol.m-3]": c_init}
+    )

tests/unit/test_optimisation.py

@@ -403,6 +403,7 @@ def test_halting(self, cost):
         optim = pybop.Optimisation(cost=cost)
 
         # Trigger threshold
+        optim.set_threshold(None)
         optim.set_threshold(np.inf)
         optim.run()
         optim.set_max_unchanged_iterations()

tests/unit/test_parameter_sets.py

@@ -123,3 +123,16 @@ def test_bpx_parameter_sets(self):
             match="Parameter set already constructed, or path to bpx file not provided.",
         ):
             bpx_parameters.import_from_bpx()
+
+    @pytest.mark.unit
+    def test_set_formation_concentrations(self):
+        parameter_set = pybop.ParameterSet.pybamm(
+            "Chen2020", formation_concentrations=True
+        )
+
+        assert (
+            parameter_set["Initial concentration in negative electrode [mol.m-3]"] == 0
+        )
+        assert (
+            parameter_set["Initial concentration in positive electrode [mol.m-3]"] > 0
+        )