Merge pull request #3456 from abillscmu/issue-3224-initial_soc

make initial soc work with half cell models

CHANGELOG.md

@@ -37,6 +37,7 @@
 - Error generated when invalid parameter values are passed ([#3132](https://github.com/pybamm-team/PyBaMM/pull/3132))
 - Parameters in `Prada2013` have been updated to better match those given in the paper, which is a 2.3 Ah cell, instead of the mix-and-match with the 1.1 Ah cell from Lain2019 ([#3096](https://github.com/pybamm-team/PyBaMM/pull/3096))
 - The `OneDimensionalX` thermal model has been updated to account for edge/tab cooling and account for the current collector volumetric heat capacity. It now gives the correct behaviour compared with a lumped model with the correct total heat transfer coefficient and surface area for cooling. ([#3042](https://github.com/pybamm-team/PyBaMM/pull/3042))
+- Fixed a bug where supplying an initial soc did not work with half cell models ([#3456](https://github.com/pybamm-team/PyBaMM/pull/3456))
 
 ## Optimizations
 

pybamm/models/full_battery_models/base_battery_model.py

@@ -551,7 +551,7 @@ def __init__(self, extra_options):
             )
         if options["working electrode"] == "negative":
             raise pybamm.OptionError(
-                "The 'negative' working elecrtrode option has been removed because "
+                "The 'negative' working electrode option has been removed because "
                 "the voltage - and therefore the energy stored - would be negative."
                 "Use the 'positive' working electrode option instead and set whatever "
                 "would normally be the negative electrode as the positive electrode."

pybamm/models/full_battery_models/lithium_ion/__init__.py

@@ -9,7 +9,10 @@
     get_initial_ocps,
     get_min_max_ocps,
 )
-from .electrode_soh_half_cell import ElectrodeSOHHalfCell
+from .electrode_soh_half_cell import (
+    ElectrodeSOHHalfCell,
+    get_initial_stoichiometry_half_cell
+)
 from .spm import SPM
 from .spme import SPMe
 from .dfn import DFN

pybamm/models/full_battery_models/lithium_ion/electrode_soh_half_cell.py

@@ -56,3 +56,94 @@ def __init__(self, name="Electrode-specific SOH model"):
     def default_solver(self):
         # Use AlgebraicSolver as CasadiAlgebraicSolver gives unnecessary warnings
         return pybamm.AlgebraicSolver()
+
+
+def get_initial_stoichiometry_half_cell(
+    initial_value,
+    parameter_values,
+    param=None,
+    known_value="cyclable lithium capacity",
+    options=None,
+):
+    """
+    Calculate initial stoichiometry to start off the simulation at a particular
+    state of charge, given voltage limits, open-circuit potential, etc defined by
+    parameter_values
+
+    Parameters
+    ----------
+    initial_value : float
+        Target initial value.
+        If integer, interpreted as SOC, 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 : pybamm.ParameterValues
+        The parameter values to use in the calculation
+
+    Returns
+    -------
+    x
+        The initial stoichiometry that give the desired initial state of charge
+    """
+    param = pybamm.LithiumIonParameters(options)
+    x_0, x_100 = get_min_max_stoichiometries(parameter_values)
+
+    if isinstance(initial_value, str) and initial_value.endswith("V"):
+        V_init = float(initial_value[:-1])
+        V_min = parameter_values.evaluate(param.voltage_low_cut)
+        V_max = parameter_values.evaluate(param.voltage_high_cut)
+
+        if not V_min < V_init < V_max:
+            raise ValueError(
+                f"Initial voltage {V_init}V is outside the voltage limits "
+                f"({V_min}, {V_max})"
+            )
+
+        # Solve simple model for initial soc based on target voltage
+        soc_model = pybamm.BaseModel()
+        soc = pybamm.Variable("soc")
+        Up = param.p.prim.U
+        T_ref = parameter_values["Reference temperature [K]"]
+        x = x_0 + soc * (x_100 - x_0)
+
+        soc_model.algebraic[soc] = Up(x, T_ref) - V_init
+        # initial guess for soc linearly interpolates between 0 and 1
+        # based on V linearly interpolating between V_max and V_min
+        soc_model.initial_conditions[soc] = (V_init - V_min) / (V_max - V_min)
+        soc_model.variables["soc"] = soc
+        parameter_values.process_model(soc_model)
+        initial_soc = pybamm.AlgebraicSolver().solve(soc_model, [0])["soc"].data[0]
+    elif isinstance(initial_value, (int, float)):
+        initial_soc = initial_value
+        if not 0 <= initial_soc <= 1:
+            raise ValueError("Initial SOC should be between 0 and 1")
+
+    else:
+        raise ValueError(
+            "Initial value must be a float between 0 and 1, "
+            "or a string ending in 'V'"
+        )
+
+    x = x_0 + initial_soc * (x_100 - x_0)
+
+    return x
+
+
+def get_min_max_stoichiometries(
+    parameter_values, options={"working electrode": "positive"}
+):
+    """
+    Get the minimum and maximum stoichiometries from the parameter values
+
+    Parameters
+    ----------
+    parameter_values : pybamm.ParameterValues
+        The parameter values to use in the calculation
+    """
+    esoh_model = pybamm.lithium_ion.ElectrodeSOHHalfCell(options)
+    param = pybamm.LithiumIonParameters(options)
+    esoh_sim = pybamm.Simulation(esoh_model, parameter_values=parameter_values)
+    Q_w = parameter_values.evaluate(param.p.Q_init)
+    esoh_sol = esoh_sim.solve([0], inputs={"Q_w": Q_w})
+    x_0, x_100 = esoh_sol["x_0"].data[0], esoh_sol["x_100"].data[0]
+    return x_0, x_100

pybamm/parameters/parameter_values.py

@@ -271,6 +271,39 @@ def update(self, values, check_conflict=False, check_already_exists=True, path="
         # reset processed symbols
         self._processed_symbols = {}
 
+    def set_initial_stoichiometry_half_cell(
+        self,
+        initial_value,
+        param=None,
+        known_value="cyclable lithium capacity",
+        inplace=True,
+        options=None,
+    ):
+        """
+        Set the initial stoichiometry of the working electrode, based on the initial
+        SOC or voltage
+        """
+        param = param or pybamm.LithiumIonParameters(options)
+        x = pybamm.lithium_ion.get_initial_stoichiometry_half_cell(
+            initial_value, self, param=param, known_value=known_value, options=options
+        )
+        if inplace:
+            parameter_values = self
+        else:
+            parameter_values = self.copy()
+
+        c_max = self.evaluate(param.p.prim.c_max)
+
+        parameter_values.update(
+            {
+                "Initial concentration in {} electrode [mol.m-3]".format(
+                    options["working electrode"]
+                ): x
+                * c_max
+            }
+        )
+        return parameter_values
+
     def set_initial_stoichiometries(
         self,
         initial_value,

pybamm/simulation.py

@@ -290,9 +290,10 @@ def update_new_model_events(self, new_model, op):
                 # figure out whether the voltage event is greater than the starting
                 # voltage (charge) or less (discharge) and set the sign of the
                 # event accordingly
-                if (isinstance(op.value, pybamm.Interpolant) or
-                    isinstance(op.value, pybamm.Multiplication)):
-                    inpt = {"start time":0}
+                if isinstance(op.value, pybamm.Interpolant) or isinstance(
+                    op.value, pybamm.Multiplication
+                ):
+                    inpt = {"start time": 0}
                     init_curr = op.value.evaluate(t=0, inputs=inpt).flatten()[0]
                     sign = np.sign(init_curr)
                 else:
@@ -373,8 +374,16 @@ def set_initial_soc(self, initial_soc):
         options = self.model.options
         param = self._model.param
         if options["open-circuit potential"] == "MSMR":
-            self._parameter_values = self._unprocessed_parameter_values.set_initial_ocps(  # noqa: E501
-                initial_soc, param=param, inplace=False, options=options
+            self._parameter_values = (
+                self._unprocessed_parameter_values.set_initial_ocps(  # noqa: E501
+                    initial_soc, param=param, inplace=False, options=options
+                )
+            )
+        elif options["working electrode"] == "positive":
+            self._parameter_values = (
+                self._unprocessed_parameter_values.set_initial_stoichiometry_half_cell(
+                    initial_soc, param=param, inplace=False, options=options
+                )
             )
         else:
             self._parameter_values = (

tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_electrode_soh.py

@@ -346,6 +346,9 @@ def test_initial_soc_cell_capacity(self):
 
     def test_error(self):
         parameter_values = pybamm.ParameterValues("Chen2020")
+        parameter_values_half_cell = pybamm.lithium_ion.DFN(
+            {"working electrode": "positive"}
+        ).default_parameter_values
 
         with self.assertRaisesRegex(
             ValueError, "Initial SOC should be between 0 and 1"
@@ -358,6 +361,23 @@ def test_error(self):
         with self.assertRaisesRegex(ValueError, "must be a float"):
             pybamm.lithium_ion.get_initial_stoichiometries("5 A", parameter_values)
 
+        with self.assertRaisesRegex(ValueError, "outside the voltage limits"):
+            pybamm.lithium_ion.get_initial_stoichiometry_half_cell(
+                "1 V", parameter_values_half_cell
+            )
+
+        with self.assertRaisesRegex(ValueError, "must be a float"):
+            pybamm.lithium_ion.get_initial_stoichiometry_half_cell(
+                "5 A", parameter_values_half_cell
+            )
+
+        with self.assertRaisesRegex(
+            ValueError, "Initial SOC should be between 0 and 1"
+        ):
+            pybamm.lithium_ion.get_initial_stoichiometry_half_cell(
+                2, parameter_values_half_cell
+            )
+
 
 class TestGetInitialOCP(TestCase):
     def test_get_initial_ocp(self):

tests/unit/test_parameters/test_parameter_values.py

@@ -15,6 +15,7 @@
     lico2_ocp_Dualfoil1998,
     lico2_diffusivity_Dualfoil1998,
 )
+from pybamm.expression_tree.exceptions import OptionError
 import casadi
 
 
@@ -119,6 +120,56 @@ def test_set_initial_stoichiometries(self):
         y_100 = param_100["Initial concentration in positive electrode [mol.m-3]"]
         self.assertAlmostEqual(y, y_0 - 0.4 * (y_0 - y_100))
 
+    def test_set_initial_stoichiometry_half_cell(self):
+        param = pybamm.lithium_ion.DFN(
+            {"working electrode": "positive"}
+        ).default_parameter_values
+        param = param.set_initial_stoichiometry_half_cell(
+            0.4, inplace=False, options={"working electrode": "positive"}
+        )
+        param_0 = param.set_initial_stoichiometry_half_cell(
+            0, inplace=False, options={"working electrode": "positive"}
+        )
+        param_100 = param.set_initial_stoichiometry_half_cell(
+            1, inplace=False, options={"working electrode": "positive"}
+        )
+
+        y = param["Initial concentration in positive electrode [mol.m-3]"]
+        y_0 = param_0["Initial concentration in positive electrode [mol.m-3]"]
+        y_100 = param_100["Initial concentration in positive electrode [mol.m-3]"]
+        self.assertAlmostEqual(y, y_0 - 0.4 * (y_0 - y_100))
+
+        # inplace for 100% coverage
+        param_t = pybamm.lithium_ion.DFN(
+            {"working electrode": "positive"}
+        ).default_parameter_values
+        param_t.set_initial_stoichiometry_half_cell(
+            0.4, inplace=True, options={"working electrode": "positive"}
+        )
+        y = param_t["Initial concentration in positive electrode [mol.m-3]"]
+        param_0 = pybamm.lithium_ion.DFN(
+            {"working electrode": "positive"}
+        ).default_parameter_values
+        param_0.set_initial_stoichiometry_half_cell(
+            0, inplace=True, options={"working electrode": "positive"}
+        )
+        y_0 = param_0["Initial concentration in positive electrode [mol.m-3]"]
+        param_100 = pybamm.lithium_ion.DFN(
+            {"working electrode": "positive"}
+        ).default_parameter_values
+        param_100.set_initial_stoichiometry_half_cell(
+            1, inplace=True, options={"working electrode": "positive"}
+        )
+        y_100 = param_100["Initial concentration in positive electrode [mol.m-3]"]
+        self.assertAlmostEqual(y, y_0 - 0.4 * (y_0 - y_100))
+
+        # test error
+        param = pybamm.ParameterValues("Chen2020")
+        with self.assertRaisesRegex(OptionError, "working electrode"):
+            param.set_initial_stoichiometry_half_cell(
+                0.1, options={"working electrode": "negative"}
+            )
+
     def test_set_initial_ocps(self):
         options = {
             "open-circuit potential": "MSMR",

tests/unit/test_simulation.py

@@ -204,6 +204,33 @@ def test_solve_with_initial_soc(self):
         sim.build(initial_soc=0.5)
         self.assertEqual(sim._built_initial_soc, 0.5)
 
+        # Test whether initial_soc works with half cell (solve)
+        options = {"working electrode": "positive"}
+        model = pybamm.lithium_ion.DFN(options)
+        sim = pybamm.Simulation(model)
+        sim.solve([0,1], initial_soc = 0.9)
+        self.assertEqual(sim._built_initial_soc, 0.9)
+
+        # Test whether initial_soc works with half cell (build)
+        options = {"working electrode": "positive"}
+        model = pybamm.lithium_ion.DFN(options)
+        sim = pybamm.Simulation(model)
+        sim.build(initial_soc = 0.9)
+        self.assertEqual(sim._built_initial_soc, 0.9)
+
+        # Test whether initial_soc works with half cell when it is a voltage
+        model = pybamm.lithium_ion.SPM({"working electrode": "positive"})
+        parameter_values = model.default_parameter_values
+        ucv = parameter_values["Open-circuit voltage at 100% SOC [V]"]
+        parameter_values["Open-circuit voltage at 100% SOC [V]"] = ucv + 1e-12
+        parameter_values["Upper voltage cut-off [V]"] = ucv + 1e-12
+        options = {"working electrode": "positive"}
+        parameter_values["Current function [A]"] = 0.0
+        sim = pybamm.Simulation(model, parameter_values=parameter_values)
+        sol = sim.solve([0,1], initial_soc = "{} V".format(ucv))
+        voltage = sol["Terminal voltage [V]"].entries
+        self.assertAlmostEqual(voltage[0], ucv, places=5)
+
         # test with MSMR
         model = pybamm.lithium_ion.MSMR({"number of MSMR reactions": ("6", "4")})
         param = pybamm.ParameterValues("MSMR_Example")