Update mass and capacity calculations for half-cell models (#543)

* Update cell_mass parameters

* Update theoretical capacity

* Update base_echem.py

* Update linked_parameters.py

* Update CHANGELOG.md

* Update test_weighted_cost.py

* Create test_half_cell_model.py

* More descriptive variable name

* Update test_half_cell_model.py

* Add comment on balancing

---------

Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>

CHANGELOG.md

@@ -2,6 +2,7 @@
 
 ## Features
 
+- [#452](https://github.com/pybop-team/PyBOP/issues/452) - Extends `cell_mass` and `approximate_capacity` for half-cell models.
 - [#544](https://github.com/pybop-team/PyBOP/issues/544) - Allows iterative plotting using `StandardPlot`.
 - [#541](https://github.com/pybop-team/PyBOP/pull/541) - Adds `ScaledLogLikelihood` and `BaseMetaLikelihood` classes.
 - [#409](https://github.com/pybop-team/PyBOP/pull/409) - Adds plotting and convergence methods for Monte Carlo sampling. Includes open-access Tesla 4680 dataset for Bayesian inference example. Fixes transformations for sampling.

examples/notebooks/energy_based_electrode_design.ipynb

@@ -76,6 +76,7 @@
    "outputs": [],
    "source": [
     "import numpy as np\n",
+    "from pybamm import Parameter\n",
     "\n",
     "import pybop\n",
     "\n",
@@ -125,6 +126,24 @@
    "outputs": [],
    "source": [
     "parameter_set = pybop.ParameterSet.pybamm(\"Chen2020\")\n",
+    "parameter_set.update(\n",
+    "    {\n",
+    "        \"Electrolyte density [kg.m-3]\": Parameter(\"Separator density [kg.m-3]\"),\n",
+    "        \"Negative electrode active material density [kg.m-3]\": Parameter(\n",
+    "            \"Negative electrode density [kg.m-3]\"\n",
+    "        ),\n",
+    "        \"Negative electrode carbon-binder density [kg.m-3]\": Parameter(\n",
+    "            \"Negative electrode density [kg.m-3]\"\n",
+    "        ),\n",
+    "        \"Positive electrode active material density [kg.m-3]\": Parameter(\n",
+    "            \"Positive electrode density [kg.m-3]\"\n",
+    "        ),\n",
+    "        \"Positive electrode carbon-binder density [kg.m-3]\": Parameter(\n",
+    "            \"Positive electrode density [kg.m-3]\"\n",
+    "        ),\n",
+    "    },\n",
+    "    check_already_exists=False,\n",
+    ")\n",
     "model = pybop.lithium_ion.SPMe(parameter_set=parameter_set)"
    ]
   },

examples/scripts/linked_parameters.py

@@ -1,18 +1,40 @@
-import pybamm
+from pybamm import Parameter
 
 import pybop
 
 # The aim of this script is to show how to systematically update
 # design parameters which depend on the optimisation parameters.
 
-# Define parameter set and model
+# Define parameter set and additional parameters needed for the cost function
 parameter_set = pybop.ParameterSet.pybamm("Chen2020", formation_concentrations=True)
+parameter_set.update(
+    {
+        "Electrolyte density [kg.m-3]": Parameter("Separator density [kg.m-3]"),
+        "Negative electrode active material density [kg.m-3]": Parameter(
+            "Negative electrode density [kg.m-3]"
+        ),
+        "Negative electrode carbon-binder density [kg.m-3]": Parameter(
+            "Negative electrode density [kg.m-3]"
+        ),
+        "Positive electrode active material density [kg.m-3]": Parameter(
+            "Positive electrode density [kg.m-3]"
+        ),
+        "Positive electrode carbon-binder density [kg.m-3]": Parameter(
+            "Positive electrode density [kg.m-3]"
+        ),
+    },
+    check_already_exists=False,
+)
+
+# Link parameters
 parameter_set.update(
     {
         "Positive electrode porosity": 1
-        - pybamm.Parameter("Positive electrode active material volume fraction")
+        - Parameter("Positive electrode active material volume fraction")
     }
 )
+
+# Define model
 model = pybop.lithium_ion.SPMe(parameter_set=parameter_set)
 
 # Fitting parameters

examples/scripts/maximising_energy.py

@@ -1,3 +1,5 @@
+from pybamm import Parameter
+
 import pybop
 
 # A design optimisation example loosely based on work by L.D. Couto
@@ -9,8 +11,28 @@
 # electrode widths, particle radii, volume fractions and
 # separator width.
 
-# Define parameter set and model
+# Define parameter set and additional parameters needed for the cost function
 parameter_set = pybop.ParameterSet.pybamm("Chen2020", formation_concentrations=True)
+parameter_set.update(
+    {
+        "Electrolyte density [kg.m-3]": Parameter("Separator density [kg.m-3]"),
+        "Negative electrode active material density [kg.m-3]": Parameter(
+            "Negative electrode density [kg.m-3]"
+        ),
+        "Negative electrode carbon-binder density [kg.m-3]": Parameter(
+            "Negative electrode density [kg.m-3]"
+        ),
+        "Positive electrode active material density [kg.m-3]": Parameter(
+            "Positive electrode density [kg.m-3]"
+        ),
+        "Positive electrode carbon-binder density [kg.m-3]": Parameter(
+            "Positive electrode density [kg.m-3]"
+        ),
+    },
+    check_already_exists=False,
+)
+
+# Define model
 model = pybop.lithium_ion.SPMe(parameter_set=parameter_set)
 
 # Fitting parameters

examples/scripts/maximising_power.py

@@ -1,7 +1,29 @@
+from pybamm import Parameter
+
 import pybop
 
-# Define parameter set and model
+# Define parameter set and additional parameters needed for the cost function
 parameter_set = pybop.ParameterSet.pybamm("Chen2020", formation_concentrations=True)
+parameter_set.update(
+    {
+        "Electrolyte density [kg.m-3]": Parameter("Separator density [kg.m-3]"),
+        "Negative electrode active material density [kg.m-3]": Parameter(
+            "Negative electrode density [kg.m-3]"
+        ),
+        "Negative electrode carbon-binder density [kg.m-3]": Parameter(
+            "Negative electrode density [kg.m-3]"
+        ),
+        "Positive electrode active material density [kg.m-3]": Parameter(
+            "Positive electrode density [kg.m-3]"
+        ),
+        "Positive electrode carbon-binder density [kg.m-3]": Parameter(
+            "Positive electrode density [kg.m-3]"
+        ),
+    },
+    check_already_exists=False,
+)
+
+# Define model
 model = pybop.lithium_ion.SPMe(parameter_set=parameter_set)
 
 # Define useful quantities

pybop/models/lithium_ion/base_echem.py

@@ -2,6 +2,7 @@
 import warnings
 from typing import Optional
 
+from pybamm import LithiumIonParameters
 from pybamm import lithium_ion as pybamm_lithium_ion
 
 from pybop.models.base_model import BaseModel, Inputs
@@ -86,6 +87,7 @@ def __init__(
         self._disc = None
 
         self._electrode_soh = pybamm_lithium_ion.electrode_soh
+        self._electrode_soh_half_cell = pybamm_lithium_ion.electrode_soh_half_cell
         self.geometric_parameters = self.set_geometric_parameters()
 
     def _check_params(
@@ -213,30 +215,36 @@ def cell_mass(self, parameter_set: Optional[ParameterSet] = None):
         parameter_set = parameter_set or self._parameter_set
 
         def mass_density(
-            active_material_vol_frac, density, porosity, electrolyte_density
+            active_material_vol_frac,
+            density,
+            porosity,
+            electrolyte_density,
+            carbon_binder_domain_density,
         ):
-            return (active_material_vol_frac * density) + (
-                porosity * electrolyte_density
+            return (
+                (active_material_vol_frac * density)
+                + (porosity * electrolyte_density)
+                + (1.0 - active_material_vol_frac - porosity)
+                * carbon_binder_domain_density
             )
 
         def area_density(thickness, mass_density):
             return thickness * mass_density
 
-        # Approximations due to SPM(e) parameter set limitations
-        electrolyte_density = parameter_set["Separator density [kg.m-3]"]
-
         # Calculate mass densities
         positive_mass_density = mass_density(
             parameter_set["Positive electrode active material volume fraction"],
-            parameter_set["Positive electrode density [kg.m-3]"],
+            parameter_set["Positive electrode active material density [kg.m-3]"],
             parameter_set["Positive electrode porosity"],
-            electrolyte_density,
+            parameter_set["Electrolyte density [kg.m-3]"],
+            parameter_set["Positive electrode carbon-binder density [kg.m-3]"],
         )
         negative_mass_density = mass_density(
             parameter_set["Negative electrode active material volume fraction"],
-            parameter_set["Negative electrode density [kg.m-3]"],
+            parameter_set["Negative electrode active material density [kg.m-3]"],
             parameter_set["Negative electrode porosity"],
-            electrolyte_density,
+            parameter_set["Electrolyte density [kg.m-3]"],
+            parameter_set["Negative electrode carbon-binder density [kg.m-3]"],
         )
 
         # Calculate area densities
@@ -248,7 +256,7 @@ def area_density(thickness, mass_density):
         )
         separator_area_density = area_density(
             parameter_set["Separator thickness [m]"],
-            parameter_set["Separator porosity"] * electrolyte_density,
+            parameter_set["Separator density [kg.m-3]"],
         )
         positive_cc_area_density = area_density(
             parameter_set["Positive current collector thickness [m]"],
@@ -279,8 +287,8 @@ def area_density(thickness, mass_density):
     def approximate_capacity(self, parameter_set: Optional[ParameterSet] = None):
         """
         Calculate an estimate for the nominal cell capacity. The estimate is computed
-        by dividing the theoretical energy (in watt-hours) by the average open circuit
-        potential (voltage) of the cell.
+        by estimating the capacity of the positive electrode that lies between the
+        stoichiometric limits corresponding to the upper and lower voltage limits.
 
         Parameters
         ----------
@@ -290,39 +298,31 @@ def approximate_capacity(self, parameter_set: Optional[ParameterSet] = None):
         Returns
         -------
         float
-            The estimate of the nominal cell capacity.
+            The estimate of the nominal cell capacity [A.h].
         """
         parameter_set = parameter_set or self._parameter_set
 
-        # Calculate theoretical energy density
-        theoretical_energy = self._electrode_soh.calculate_theoretical_energy(
-            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(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 = parameter_set["Positive electrode OCP [V]"]
-        negative_electrode_ocp = parameter_set["Negative electrode OCP [V]"]
-        try:
-            average_voltage = positive_electrode_ocp(
-                mean_sto_pos
-            ) - negative_electrode_ocp(mean_sto_neg)
-        except Exception as e:
-            raise ValueError(f"Error in average voltage calculation: {e}") from e
-
-        # Calculate the capacity estimate
-        approximate_capacity = theoretical_energy / average_voltage
-
-        return ParameterSet.evaluate_symbol(approximate_capacity, parameter_set)
+        # Calculate the theoretical capacity in the limit of low current
+        if self.pybamm_model.options["working electrode"] == "positive":
+            (
+                max_sto_p,
+                min_sto_p,
+            ) = self._electrode_soh_half_cell.get_min_max_stoichiometries(parameter_set)
+        else:
+            (
+                min_sto_n,
+                max_sto_n,
+                min_sto_p,
+                max_sto_p,
+            ) = self._electrode_soh.get_min_max_stoichiometries(parameter_set)
+            # Note that the stoichiometric limits correspond to 0 and 100% SOC.
+            # Stoichiometric balancing is performed within get_min_max_stoichiometries
+            # such that the capacity accessible between the limits should be the same
+            # for both electrodes, so we consider just the positive electrode below.
+
+        Q_p = LithiumIonParameters().p.prim.Q_init
+        theoretical_capacity = Q_p * (max_sto_p - min_sto_p)
+        return ParameterSet.evaluate_symbol(theoretical_capacity, parameter_set)
 
     def set_geometric_parameters(self):
         """