Merge pull request #2210 from pybamm-team/faster-electrode-soh

Faster electrode soh

CHANGELOG.md

@@ -4,6 +4,10 @@
 
 -   For experiments, the simulation now automatically checks and skips steps that cannot be performed (e.g. "Charge at 1C until 4.2V" from 100% SOC) ([#2212](https://github.com/pybamm-team/PyBaMM/pull/2212))
 
+## Optimizations
+
+-   Sped up calculations of Electrode SOH variables for summary variables ([#2210](https://github.com/pybamm-team/PyBaMM/pull/2210))
+
 ## Breaking changes
 
 -   Events must now be defined in such a way that they are positive at the initial conditions (events will be triggered when they become negative, instead of when they change sign in either direction) ([#2212](https://github.com/pybamm-team/PyBaMM/pull/2212))

pybamm/models/full_battery_models/lithium_ion/electrode_soh.py

@@ -53,10 +53,12 @@ def __init__(self, name="ElectrodeSOH model", param=None):
         C = Cn * (x_100 - x_0)
         y_0 = y_100 + C / Cp
 
-        self.algebraic = {
-            x_100: Up(y_100, T_ref) - Un(x_100, T_ref) - V_max,
-            x_0: Up(y_0, T_ref) - Un(x_0, T_ref) - V_min,
-        }
+        Un_0 = Un(x_0, T_ref)
+        Up_0 = Up(y_0, T_ref)
+        Un_100 = Un(x_100, T_ref)
+        Up_100 = Up(y_100, T_ref)
+
+        self.algebraic = {x_100: Up_100 - Un_100 - V_max, x_0: Up_0 - Un_0 - V_min}
 
         self.initial_conditions = {x_0: pybamm.Scalar(0.1), x_100: pybamm.Scalar(0.9)}
 
@@ -67,12 +69,12 @@ def __init__(self, name="ElectrodeSOH model", param=None):
             "y_100": y_100,
             "x_0": x_0,
             "y_0": y_0,
-            "Un(x_100)": Un(x_100, T_ref),
-            "Up(y_100)": Up(y_100, T_ref),
-            "Un(x_0)": Un(x_0, T_ref),
-            "Up(y_0)": Up(y_0, T_ref),
-            "Up(y_0) - Un(x_0)": Up(y_0, T_ref) - Un(x_0, T_ref),
-            "Up(y_100) - Un(x_100)": Up(y_100, T_ref) - Un(x_100, T_ref),
+            "Un(x_100)": Un_100,
+            "Up(y_100)": Up_100,
+            "Un(x_0)": Un_0,
+            "Up(y_0)": Up_0,
+            "Up(y_0) - Un(x_0)": Up_0 - Un_0,
+            "Up(y_100) - Un(x_100)": Up_100 - Un_100,
             "n_Li_100": 3600 / param.F * (y_100 * Cp + x_100 * Cn),
             "n_Li_0": 3600 / param.F * (y_0 * Cp + x_0 * Cn),
             "n_Li": n_Li,
@@ -117,12 +119,12 @@ def __init__(self, name="ElectrodeSOHx100 model", param=None):
         Cp = pybamm.InputParameter("C_p")
 
         x_100 = pybamm.Variable("x_100")
-
         y_100 = (n_Li * param.F / 3600 - x_100 * Cn) / Cp
 
-        self.algebraic = {
-            x_100: Up(y_100, T_ref) - Un(x_100, T_ref) - V_max,
-        }
+        Un_100 = Un(x_100, T_ref)
+        Up_100 = Up(y_100, T_ref)
+
+        self.algebraic = {x_100: Up_100 - Un_100 - V_max}
 
         self.initial_conditions = {x_100: pybamm.Scalar(0.9)}
 
@@ -170,7 +172,12 @@ def __init__(self, name="ElectrodeSOHx0 model", param=None):
         C = Cn * (x_100 - x_0)
         y_0 = y_100 + C / Cp
 
-        self.algebraic = {x_0: Up(y_0, T_ref) - Un(x_0, T_ref) - V_min}
+        Un_0 = Un(x_0, T_ref)
+        Up_0 = Up(y_0, T_ref)
+        Un_100 = Un(x_100, T_ref)
+        Up_100 = Up(y_100, T_ref)
+
+        self.algebraic = {x_0: Up_0 - Un_0 - V_min}
 
         self.initial_conditions = {x_0: pybamm.Scalar(0.1)}
 
@@ -179,12 +186,12 @@ def __init__(self, name="ElectrodeSOHx0 model", param=None):
             "Capacity [A.h]": C,
             "x_0": x_0,
             "y_0": y_0,
-            "Un(x_100)": Un(x_100, T_ref),
-            "Up(y_100)": Up(y_100, T_ref),
-            "Un(x_0)": Un(x_0, T_ref),
-            "Up(y_0)": Up(y_0, T_ref),
-            "Up(y_0) - Un(x_0)": Up(y_0, T_ref) - Un(x_0, T_ref),
-            "Up(y_100) - Un(x_100)": Up(y_100, T_ref) - Un(x_100, T_ref),
+            "Un(x_100)": Un_100,
+            "Up(y_100)": Up_100,
+            "Un(x_0)": Un_0,
+            "Up(y_0)": Up_0,
+            "Up(y_0) - Un(x_0)": Up_0 - Un_0,
+            "Up(y_100) - Un(x_100)": Up_100 - Un_100,
             "n_Li_100": 3600 / param.F * (y_100 * Cp + x_100 * Cn),
             "n_Li_0": 3600 / param.F * (y_0 * Cp + x_0 * Cn),
             "n_Li": n_Li,
@@ -206,7 +213,6 @@ class ElectrodeSOHSolver:
     def __init__(self, parameter_values, param=None):
         self.parameter_values = parameter_values
         self.param = param or pybamm.LithiumIonParameters()
-        self.sims = self.create_electrode_soh_sims(parameter_values, self.param)
 
         # Check whether each electrode OCP is a function (False) or data (True)
         OCPp_data = isinstance(parameter_values["Positive electrode OCP [V]"], tuple)
@@ -240,27 +246,48 @@ def __init__(self, parameter_values, param=None):
             - self.param.n.prim.U_dimensional(x, T)
         )
 
-    def create_electrode_soh_sims(self, parameter_values, param):
-        full_model = ElectrodeSOH(param=param)
-        full_sim = pybamm.Simulation(full_model, parameter_values=parameter_values)
-        x100_model = ElectrodeSOHx100(param=param)
-        x100_sim = pybamm.Simulation(x100_model, parameter_values=parameter_values)
-        x0_model = ElectrodeSOHx0(param=param)
-        x0_sim = pybamm.Simulation(x0_model, parameter_values=parameter_values)
-        return {"combined": full_sim, "split": [x100_sim, x0_sim]}
+    def _get_electrode_soh_sims_full(self):
+        try:
+            return self._full_sim
+        except AttributeError:
+            full_model = ElectrodeSOH(param=self.param)
+            self._full_sim = pybamm.Simulation(
+                full_model, parameter_values=self.parameter_values
+            )
+            return self._full_sim
+
+    def _get_electrode_soh_sims_split(self):
+        try:
+            return self._split_sims
+        except AttributeError:
+            x100_model = ElectrodeSOHx100(param=self.param)
+            x100_sim = pybamm.Simulation(
+                x100_model, parameter_values=self.parameter_values
+            )
+            x0_model = ElectrodeSOHx0(param=self.param)
+            x0_sim = pybamm.Simulation(x0_model, parameter_values=self.parameter_values)
+            self._split_sims = [x100_sim, x0_sim]
+            return self._split_sims
 
     def solve(self, inputs):
-        ics = self.set_up_solve(inputs)
+        ics = self._set_up_solve(inputs)
         try:
-            sol = self.solve_full(inputs, ics)
+            sol = self._solve_full(inputs, ics)
         except pybamm.SolverError:
             # just in case solving one by one works better
-            sol = self.solve_split(inputs, ics)
+            try:
+                sol = self._solve_split(inputs, ics)
+            except pybamm.SolverError as original_error:
+                # check if the error is due to the simulation not being feasible
+                self._check_esoh_feasible(inputs)
+                # if that didn't raise an error, raise the original error instead
+                raise original_error  # pragma: no cover (don't know how to get here)
+
         return sol
 
-    def set_up_solve(self, inputs):
-        sim = self.sims["combined"]
-        x0_min, x100_max, _, _ = self.check_esoh_feasible(inputs)
+    def _set_up_solve(self, inputs):
+        sim = self._get_electrode_soh_sims_full()
+        x0_min, x100_max, _, _ = self._get_lims(inputs)
 
         x100_init = x100_max
         x0_init = x0_min
@@ -274,15 +301,15 @@ def set_up_solve(self, inputs):
                 x0_init = x0_init_sol
         return {"x_100": np.array(x100_init), "x_0": np.array(x0_init)}
 
-    def solve_full(self, inputs, ics):
-        sim = self.sims["combined"]
+    def _solve_full(self, inputs, ics):
+        sim = self._get_electrode_soh_sims_full()
         sim.build()
         sim.built_model.set_initial_conditions_from(ics)
         sol = sim.solve([0], inputs=inputs)
         return sol
 
-    def solve_split(self, inputs, ics):
-        x100_sim, x0_sim = self.sims["split"]
+    def _solve_split(self, inputs, ics):
+        x100_sim, x0_sim = self._get_electrode_soh_sims_split()
         x100_sim.build()
         x100_sim.built_model.set_initial_conditions_from(ics)
         x100_sol = x100_sim.solve([0], inputs=inputs)
@@ -295,9 +322,10 @@ def solve_split(self, inputs, ics):
 
         return x0_sol
 
-    def check_esoh_feasible(self, inputs):
-        Vmax = inputs["V_max"]
-        Vmin = inputs["V_min"]
+    def _get_lims(self, inputs):
+        """
+        Get stoichiometry limits based on n_Li, C_n, and C_p
+        """
         Cp = inputs["C_p"]
         Cn = inputs["C_n"]
         n_Li = inputs["n_Li"]
@@ -316,6 +344,15 @@ def check_esoh_feasible(self, inputs):
         x0_min = max(x0_min_from_y0_max, x0_min)
         y100_min = max(y100_min_from_x100_max, y100_min)
         y0_max = min(y0_max_from_x0_min, y0_max)
+        return (x0_min, x100_max, y100_min, y0_max)
+
+    def _check_esoh_feasible(self, inputs):
+        """
+        Check that the electrode SOH calculation is feasible, based on voltage limits
+        """
+        x0_min, x100_max, y100_min, y0_max = self._get_lims(inputs)
+        Vmax = inputs["V_max"]
+        Vmin = inputs["V_min"]
 
         # Check stoich limits are between 0 and 1
         for x in ["x0_min", "x100_max", "y100_min", "y0_max"]:
@@ -348,8 +385,6 @@ def check_esoh_feasible(self, inputs):
                 )
             )
 
-        return (x0_min, x100_max, y100_min, y0_max)
-
 
 def get_initial_stoichiometries(initial_soc, parameter_values):
     """

pybamm/simulation.py

@@ -741,7 +741,9 @@ def solve(
                     cycle_sum_vars,
                     cycle_first_state,
                 ) = pybamm.make_cycle_solution(
-                    starting_solution.steps, esoh_solver, True
+                    starting_solution.steps,
+                    esoh_solver=esoh_solver,
+                    save_this_cycle=True
                 )
                 starting_solution_cycles = [cycle_solution]
                 starting_solution_summary_variables = [cycle_sum_vars]
@@ -896,9 +898,7 @@ def solve(
                                 "due to exceeded bounds at initial conditions."
                             )
                     cycle_sol = pybamm.make_cycle_solution(
-                        steps,
-                        esoh_solver,
-                        save_this_cycle=save_this_cycle,
+                        steps, esoh_solver=esoh_solver, save_this_cycle=save_this_cycle
                     )
                     cycle_solution, cycle_sum_vars, cycle_first_state = cycle_sol
                     all_cycle_solutions.append(cycle_solution)

pybamm/solvers/base_solver.py

@@ -796,8 +796,15 @@ def solve(
             ics_set_up = self.models_set_up[model]["initial conditions"]
             # Check that initial conditions have not been updated
             if ics_set_up != model.concatenated_initial_conditions:
-                # If the new initial conditions are different, set up again
-                self.set_up(model, ext_and_inputs_list[0], t_eval, ics_only=True)
+                if self.algebraic_solver is True:
+                    # For an algebraic solver, we don't need to set up the initial
+                    # conditions function and we can just evaluate
+                    # model.concatenated_initial_conditions
+                    model.y0 = model.concatenated_initial_conditions.evaluate()
+                else:
+                    # If the new initial conditions are different
+                    # and cannot be evaluated directly, set up again
+                    self.set_up(model, ext_and_inputs_list[0], t_eval, ics_only=True)
                 self.models_set_up[model][
                     "initial conditions"
                 ] = model.concatenated_initial_conditions

pybamm/solvers/solution.py

@@ -827,7 +827,7 @@ def make_cycle_solution(step_solutions, esoh_solver=None, save_this_cycle=True):
 
     cycle_solution.steps = step_solutions
 
-    cycle_summary_variables = get_cycle_summary_variables(cycle_solution, esoh_solver)
+    cycle_summary_variables = _get_cycle_summary_variables(cycle_solution, esoh_solver)
 
     cycle_first_state = cycle_solution.first_state
 
@@ -839,7 +839,7 @@ def make_cycle_solution(step_solutions, esoh_solver=None, save_this_cycle=True):
     return cycle_solution, cycle_summary_variables, cycle_first_state
 
 
-def get_cycle_summary_variables(cycle_solution, esoh_solver):
+def _get_cycle_summary_variables(cycle_solution, esoh_solver):
     model = cycle_solution.all_models[0]
     cycle_summary_variables = pybamm.FuzzyDict({})
 

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

@@ -30,11 +30,30 @@ def test_known_solution(self):
         self.assertAlmostEqual(sol["n_Li_0"].data[0], n_Li, places=5)
 
         # Solve with split esoh and check outputs
-        ics = esoh_solver.set_up_solve(inputs)
-        sol_split = esoh_solver.solve_split(inputs, ics)
+        ics = esoh_solver._set_up_solve(inputs)
+        sol_split = esoh_solver._solve_split(inputs, ics)
         for key in sol.all_models[0].variables:
             self.assertAlmostEqual(sol[key].data[0], sol_split[key].data[0], places=5)
 
+    def test_error(self):
+
+        param = pybamm.LithiumIonParameters()
+        parameter_values = pybamm.ParameterValues("Mohtat2020")
+
+        esoh_solver = pybamm.lithium_ion.ElectrodeSOHSolver(parameter_values, param)
+
+        Vmin = 3
+        Vmax = 4.2
+        Cn = parameter_values.evaluate(param.n.cap_init)
+        Cp = parameter_values.evaluate(param.p.cap_init)
+        n_Li = parameter_values.evaluate(param.n_Li_particles_init) * 10
+
+        inputs = {"V_max": Vmax, "V_min": Vmin, "n_Li": n_Li, "C_n": Cn, "C_p": Cp}
+
+        # Solve the model and check outputs
+        with self.assertRaisesRegex(ValueError, "should be between 0 and 1"):
+            esoh_solver.solve(inputs)
+
 
 class TestElectrodeSOHHalfCell(unittest.TestCase):
     def test_known_solution(self):
@@ -96,10 +115,10 @@ def test_error(self):
 
         inputs = {"V_min": 0, "V_max": 6, "C_n": C_n, "C_p": C_p, "n_Li": n_Li}
         with self.assertRaisesRegex(ValueError, "lower bound of the voltage"):
-            esoh_solver.check_esoh_feasible(inputs)
+            esoh_solver._check_esoh_feasible(inputs)
         inputs = {"V_min": 3, "V_max": 6, "C_n": C_n, "C_p": C_p, "n_Li": n_Li}
         with self.assertRaisesRegex(ValueError, "upper bound of the voltage"):
-            esoh_solver.check_esoh_feasible(inputs)
+            esoh_solver._check_esoh_feasible(inputs)
 
 
 if __name__ == "__main__":