Merge pull request #2095 from pybamm-team/bound-ocv

bound sto for better eSOH results

CHANGELOG.md

@@ -2,6 +2,12 @@
 ## Features
 
 -   Added `__eq__` and `__hash__` methods for `Symbol` objects, using `.id` ([#1978](https://github.com/pybamm-team/PyBaMM/pull/1978))
+
+
+## Optimizations
+
+-   Stoichiometry inputs to OCP functions are now bounded between 1e-10 and 1-1e-10, with singularities at 0 and 1 so that OCP goes to +- infinity ([#2095](https://github.com/pybamm-team/PyBaMM/pull/2095))
+
 ## Breaking changes
 
 -   Changed some dictionary keys to `Symbol` instead of `Symbol.id` (internal change only, should not affect external facing functions) ([#1978](https://github.com/pybamm-team/PyBaMM/pull/1978))

pybamm/expression_tree/binary_operators.py

@@ -623,6 +623,10 @@ def _binary_new_copy(self, left, right):
         """See :meth:`pybamm.BinaryOperator._binary_new_copy()`."""
         return pybamm.minimum(left, right)
 
+    def _sympy_operator(self, left, right):
+        """Override :meth:`pybamm.BinaryOperator._sympy_operator`"""
+        return sympy.Min(left, right)
+
 
 class Maximum(BinaryOperator):
     """Returns the greater of two objects."""
@@ -655,6 +659,10 @@ def _binary_new_copy(self, left, right):
         """See :meth:`pybamm.BinaryOperator._binary_new_copy()`."""
         return pybamm.maximum(left, right)
 
+    def _sympy_operator(self, left, right):
+        """Override :meth:`pybamm.BinaryOperator._sympy_operator`"""
+        return sympy.Max(left, right)
+
 
 def simplify_elementwise_binary_broadcasts(left, right):
     left, right = preprocess_binary(left, right)

pybamm/models/full_battery_models/lithium_ion/electrode_soh_half_cell.py

@@ -58,7 +58,7 @@ def __init__(self, working_electrode, name="Electrode-specific SOH model"):
         x_100_init = 0.85
         # Make sure x_0 = x_100 - C/C_w > 0
         C_init = param.Q
-        C_init = pybamm.minimum(Cw * x_100_init - 0.1, C_init)
+        C_init = pybamm.minimum(Cw * x_100_init, C_init)
         self.initial_conditions = {x_100: x_100_init, C: C_init}
 
         self.variables = {

pybamm/parameters/lithium_ion_parameters.py

@@ -676,6 +676,11 @@ def j0_dimensional(self, c_e, c_s_surf, T):
 
     def U_dimensional(self, sto, T):
         """Dimensional open-circuit potential [V]"""
+        # bound stoichiometry between tol and 1-tol. Adding 1/sto + 1/(sto-1) later
+        # will ensure that ocp goes to +- infinity if sto goes into that region
+        # anyway
+        tol = 1e-10
+        sto = pybamm.maximum(pybamm.minimum(sto, 1 - tol), tol)
         inputs = {f"{self.domain} particle stoichiometry": sto}
         u_ref = pybamm.FunctionParameter(f"{self.domain} electrode OCP [V]", inputs)
         # add a term to ensure that the OCP goes to infinity at 0 and -infinity at 1