#1511 add function for initial soc

pybamm/models/full_battery_models/lithium_ion/__init__.py

@@ -2,7 +2,7 @@
 # Root of the lithium-ion models module.
 #
 from .base_lithium_ion_model import BaseModel
-from .electrode_soh import ElectrodeSOH
+from .electrode_soh import ElectrodeSOH, get_initial_stoichiometries
 from .spm import SPM
 from .spme import SPMe
 from .dfn import DFN

pybamm/models/full_battery_models/lithium_ion/electrode_soh.py

@@ -93,3 +93,57 @@ 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_stoichiometries(initial_soc, parameter_values):
+    """
+    Calculate initial stoichiometries to start off the simulation at a particular
+    state of charge, given voltage limits, open-circuit potentials, etc defined by
+    parameter_values
+
+    Parameters
+    ----------
+    initial_soc : float
+        Target initial SOC. Must be between 0 and 1.
+    parameter_values : :class:`pybamm.ParameterValues`
+        The parameter values class that will be used for the simulation. Required for
+        calculating appropriate initial stoichiometries.
+
+    Returns
+    -------
+    x, y
+        The initial stoichiometries that give the desired initial state of charge
+    """
+    if initial_soc < 0 or initial_soc > 1:
+        raise ValueError("Initial SOC should be between 0 and 1")
+
+    model = pybamm.lithium_ion.ElectrodeSOH()
+
+    param = pybamm.LithiumIonParameters()
+    sim = pybamm.Simulation(model, parameter_values=parameter_values)
+
+    V_min = parameter_values.evaluate(param.voltage_low_cut_dimensional)
+    V_max = parameter_values.evaluate(param.voltage_high_cut_dimensional)
+    C_n = parameter_values.evaluate(param.C_n_init)
+    C_p = parameter_values.evaluate(param.C_p_init)
+    n_Li = parameter_values.evaluate(param.n_Li_particles_init)
+
+    # Solve the model and check outputs
+    sol = sim.solve(
+        [0],
+        inputs={
+            "V_min": V_min,
+            "V_max": V_max,
+            "C_n": C_n,
+            "C_p": C_p,
+            "n_Li": n_Li,
+        },
+    )
+
+    x_0 = sol["x_0"].data[0]
+    y_0 = sol["y_0"].data[0]
+    C = sol["C"].data[0]
+    x = x_0 + initial_soc * C / C_n
+    y = y_0 - initial_soc * C / C_p
+
+    return x, y

pybamm/simulation.py

@@ -138,6 +138,7 @@ def __init__(
         # Initialize empty built states
         self._model_with_set_params = None
         self._built_model = None
+        self.op_conds_to_built_models = None
         self._mesh = None
         self._disc = None
         self._solution = None
@@ -329,7 +330,6 @@ def set_up_model_for_experiment_old(self, model):
             op_cond[:2]: (new_model, self.parameter_values)
             for op_cond in set(self.experiment.operating_conditions)
         }
-        self.op_conds_to_built_models = None
 
     def set_up_model_for_experiment_new(self, model):
         """
@@ -340,7 +340,6 @@ def set_up_model_for_experiment_new(self, model):
         reduces simulation time since the model formulation is efficient.
         """
         self.op_conds_to_model_and_param = {}
-        self.op_conds_to_built_models = None
         for op_cond, op_inputs in zip(
             self.experiment.operating_conditions, self._experiment_inputs
         ):
@@ -545,6 +544,7 @@ def solve(
         save_at_cycles=None,
         calc_esoh=True,
         starting_solution=None,
+        initial_soc=None,
         **kwargs,
     ):
         """
@@ -585,6 +585,10 @@ def solve(
         starting_solution : :class:`pybamm.Solution`
             The solution to start stepping from. If None (default), then self._solution
             is used. Must be None if not using an experiment.
+        initial_soc : float, optional
+            Initial State of Charge (SOC) for the simulation. Must be between 0 and 1.
+            If given, overwrites the initial concentrations provided in the parameter
+            set.
         **kwargs
             Additional key-word arguments passed to `solver.solve`.
             See :meth:`pybamm.BaseSolver.solve`.
@@ -593,6 +597,39 @@ def solve(
         if solver is None:
             solver = self.solver
 
+        if initial_soc is not None:
+            if (
+                hasattr(self, "_built_initial_soc")
+                and self._built_initial_soc != initial_soc
+            ):
+                # reset
+                self._model_with_set_params = None
+                self._built_model = None
+                self.op_conds_to_built_models = None
+
+            c_n_init = self.parameter_values[
+                "Initial concentration in negative electrode [mol.m-3]"
+            ]
+            c_p_init = self.parameter_values[
+                "Initial concentration in positive electrode [mol.m-3]"
+            ]
+            param = pybamm.LithiumIonParameters()
+            c_n_max = self.parameter_values.evaluate(param.c_n_max)
+            c_p_max = self.parameter_values.evaluate(param.c_p_max)
+            x, y = pybamm.lithium_ion.get_initial_stoichiometries(
+                initial_soc, self.parameter_values
+            )
+            self.parameter_values.update(
+                {
+                    "Initial concentration in negative electrode [mol.m-3]": x
+                    * c_n_max,
+                    "Initial concentration in positive electrode [mol.m-3]": y
+                    * c_p_max,
+                }
+            )
+            # Save solved initial SOC in case we need to re-build the model
+            self._built_initial_soc = initial_soc
+
         if self.operating_mode in ["without experiment", "drive cycle"]:
             self.build(check_model=check_model)
             if save_at_cycles is not None:
@@ -838,6 +875,15 @@ def solve(
                 "Finish experiment simulation, took {}".format(timer.time())
             )
 
+        # reset parameter values
+        if initial_soc is not None:
+            self.parameter_values.update(
+                {
+                    "Initial concentration in negative electrode [mol.m-3]": c_n_init,
+                    "Initial concentration in positive electrode [mol.m-3]": c_p_init,
+                }
+            )
+
         return self.solution
 
     def step(

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

@@ -11,7 +11,7 @@ def test_known_solution(self):
 
         param = pybamm.LithiumIonParameters()
         parameter_values = pybamm.ParameterValues(
-            chemistry=pybamm.parameter_sets.Marquis2019
+            chemistry=pybamm.parameter_sets.Mohtat2020
         )
         sim = pybamm.Simulation(model, parameter_values=parameter_values)
 
@@ -38,6 +38,42 @@ def test_known_solution(self):
         self.assertAlmostEqual(sol["n_Li_0"].data[0], n_Li, places=5)
 
 
+class TestSetInitialSOC(unittest.TestCase):
+    def test_known_solutions(self):
+        model = pybamm.lithium_ion.ElectrodeSOH()
+
+        param = pybamm.LithiumIonParameters()
+        parameter_values = pybamm.ParameterValues(
+            chemistry=pybamm.parameter_sets.Mohtat2020
+        )
+        sim = pybamm.Simulation(model, parameter_values=parameter_values)
+
+        V_min = parameter_values.evaluate(param.voltage_low_cut_dimensional)
+        V_max = parameter_values.evaluate(param.voltage_high_cut_dimensional)
+        C_n = parameter_values.evaluate(param.C_n_init)
+        C_p = parameter_values.evaluate(param.C_p_init)
+        n_Li = parameter_values.evaluate(param.n_Li_particles_init)
+
+        # Solve the model and check outputs
+        esoh_sol = sim.solve(
+            [0],
+            inputs={
+                "V_min": V_min,
+                "V_max": V_max,
+                "C_n": C_n,
+                "C_p": C_p,
+                "n_Li": n_Li,
+            },
+        )
+
+        x, y = pybamm.lithium_ion.get_initial_stoichiometries(1, parameter_values)
+        self.assertAlmostEqual(x, esoh_sol["x_100"].data[0])
+        self.assertAlmostEqual(y, esoh_sol["y_100"].data[0])
+        x, y = pybamm.lithium_ion.get_initial_stoichiometries(0, parameter_values)
+        self.assertAlmostEqual(x, esoh_sol["x_0"].data[0])
+        self.assertAlmostEqual(y, esoh_sol["y_0"].data[0])
+
+
 if __name__ == "__main__":
     print("Add -v for more debug output")
     import sys

tests/unit/test_simulation.py

@@ -215,6 +215,15 @@ def test_step(self):
         self.assertEqual(sim.solution.t[1], dt / tau)
         self.assertEqual(sim.solution.t[2], 2 * dt / tau)
 
+    def test_solve_with_initial_soc(self):
+        model = pybamm.lithium_ion.SPM()
+        param = model.default_parameter_values
+        sim = pybamm.Simulation(model, parameter_values=param)
+        sim.solve(t_eval=[0, 600], initial_soc=1)
+        self.assertEqual(sim._built_initial_soc, 1)
+        sim.solve(t_eval=[0, 600], initial_soc=0.5)
+        self.assertEqual(sim._built_initial_soc, 0.5)
+
     def test_solve_with_inputs(self):
         model = pybamm.lithium_ion.SPM()
         param = model.default_parameter_values