Update set_initial_state

CHANGELOG.md

@@ -2,6 +2,7 @@
 
 ## Features
 
+- [#889](https://github.com/pybop-team/PyBOP/pull/889) - Adds methods for setting the initial state from a voltage to the grouped models.
 - [#869](https://github.com/pybop-team/PyBOP/issues/869) - Adds methods for pre-processing current data for linear interpolation.
 - [#868](https://github.com/pybop-team/PyBOP/pull/868) - Adds support for Python3.13 (NumPy restricted to <2.4, EP-BOLFI optimiser and PyProBE do not support Python 3.13).
 - [#871](https://github.com/pybop-team/PyBOP/pull/871) - Adds a lumped thermal model called `CellTemperature`.

pybop/__init__.py

@@ -58,6 +58,7 @@
 #
 from .models import lithium_ion
 from .models._exponential_decay import ExponentialDecayModel
+from .models.lithium_ion.utils import Interpolant, InverseOCV
 
 #
 # PyBaMM utility classes
@@ -176,7 +177,7 @@
 #
 # Applications
 #
-from .applications.base_method import BaseApplication, Interpolant, InverseOCV
+from .applications.base_method import BaseApplication
 from .applications.ocp_methods import OCPMerge, OCPAverage, OCPCapacityToStoichiometry
 from .applications.gitt_methods import GITTPulseFit, GITTFit
 

pybop/_result.py

@@ -1,6 +1,8 @@
 import numpy as np
 
-from pybop import Logger, Problem, plot
+from pybop import plot
+from pybop._logging import Logger
+from pybop.problems.problem import Problem
 
 
 class Result:

pybop/applications/base_method.py

@@ -1,11 +1,6 @@
 import warnings
-from collections.abc import Callable
 
 import numpy as np
-from pybamm import Interpolant as PybammInterpolant
-from scipy import interpolate
-
-import pybop
 
 
 class BaseApplication:
@@ -27,141 +22,3 @@ def check_monotonicity(self, voltage: np.ndarray) -> None:
 
         if not (is_increasing or is_decreasing):
             warnings.warn("OCV is not strictly monotonic.", stacklevel=2)
-
-
-class Interpolant:
-    """
-    A class that returns a pybamm.Interpolant to pybamm models and otherwise
-    a numeric interpolant.
-
-    Parameters
-    ----------
-    x : array_like
-        Input coordinates.
-    y : array_like
-        Output values corresponding to x.
-    name : str, optional
-        Name for the interpolant when used in PyBaMM.
-    bounds_error : bool, optional
-        If True, raise error when interpolating outside bounds.
-    fill_value : str or float, optional
-        Value to use for out-of-bounds interpolation.
-    axis : int, optional
-        Axis along which to interpolate.
-    """
-
-    def __init__(
-        self,
-        x: np.ndarray,
-        y: np.ndarray,
-        name: str | None = None,
-        bounds_error: bool = False,
-        fill_value: str | float = "extrapolate",
-        axis: int = 0,
-    ):
-        self.x = np.asarray(x)
-        self.y = np.asarray(y)
-        self.name = name
-        self._interp_func = self._create_interpolant(bounds_error, fill_value, axis)
-
-    def _create_interpolant(
-        self, bounds_error: bool, fill_value: str | float, axis: int
-    ):
-        """Create the scipy interpolation function."""
-        return interpolate.interp1d(
-            self.x,
-            self.y,
-            bounds_error=bounds_error,
-            fill_value=fill_value,
-            axis=axis,
-        )
-
-    def __call__(self, x: float | np.ndarray):
-        """
-        Evaluate the interpolant at given points.
-
-        Parameters
-        ----------
-        x : float or array_like
-            Points at which to evaluate the interpolant.
-
-        Returns
-        -------
-        float, array_like, or PybammInterpolant
-            Interpolated values or PyBaMM interpolant object.
-        """
-        try:
-            # Try numeric evaluation first
-            return self._interp_func(x)
-        except Exception:
-            # Fall back to PyBaMM interpolant for symbolic evaluation
-            return PybammInterpolant(self.x, self.y, x, name=self.name)
-
-
-class InverseOCV:
-    """
-    A class to find the stoichiometry corresponding to a given open-circuit
-    voltage.
-
-    Parameters
-    ----------
-    ocv_function : Callable
-        The open-circuit voltage as a function of stoichiometry.
-    optimiser : pybop.BaseOptimiser, optional
-        The optimisation algorithm to use (default: pybop.SciPyMinimize).
-    optimiser_options : pybop.OptimiserOptions, optional
-        Options for the optimiser.
-    """
-
-    def __init__(
-        self,
-        ocv_function: Callable,
-        optimiser: pybop.BaseOptimiser | None = None,
-        optimiser_options: pybop.OptimiserOptions | None = None,
-    ):
-        self.optimiser = optimiser or pybop.SciPyMinimize
-        self.optimiser_options = optimiser_options or self.optimiser.default_options()
-
-        parameters = pybop.Parameters(
-            {"Root": pybop.Parameter(initial_value=0.5, bounds=[0, 1])}
-        )
-
-        # Set up a root-finding cost function
-        class OCVRoot(pybop.BaseSimulator):
-            def __init__(self, ocv_value: float):
-                super().__init__(parameters=parameters)
-                self.ocv_value = ocv_value
-
-            def solve_batch(self, inputs, calculate_sensitivities: bool = False):
-                solutions = []
-                for x in inputs:
-                    diff = np.abs(ocv_function(x["Root"]) - self.ocv_value)
-                    sol = pybop.Solution()
-                    sol.set_solution_variable("Difference", data=np.asarray([diff]))
-                    solutions.append(sol)
-                return solutions
-
-        self.ocv_root = OCVRoot
-
-        # Minimise to find the stoichiometry
-        self.cost = pybop.DesignCost(target="Difference")
-        self.cost.minimising = True
-
-    def __call__(self, ocv_value: float) -> float:
-        """
-        Estimate and return the stoichiometry.
-
-        Parameters
-        ----------
-        ocv_value : float
-            The open-circuit voltage value [V] for which to estimate the stoichiometry.
-
-        Returns
-        -------
-        float
-            The stoichiometry corresponding to the open-circuit voltage value.
-        """
-        problem = pybop.Problem(self.ocv_root(ocv_value), self.cost)
-        optim = self.optimiser(problem, options=self.optimiser_options)
-        result = optim.run()
-        return result.best_inputs["Root"]

pybop/models/lithium_ion/__init__.py

@@ -1,5 +1,5 @@
 #
-# Import lithium ion based models
+# Import lithium ion models
 #
 from .sp_diffusion import SPDiffusion
 from .grouped_spm import GroupedSPM

pybop/models/lithium_ion/base_model.py

@@ -0,0 +1,124 @@
+import pybamm
+from pybamm import FunctionParameter, Parameter
+from pybamm import lithium_ion as pybamm_lithium_ion
+
+from pybop.models.lithium_ion.utils import InverseOCV
+
+
+class BaseGroupedModel(pybamm_lithium_ion.BaseModel):
+    """
+    A base model for PyBOP's grouped-parameter lithium-ion battery models.
+
+    Parameters
+    ----------
+    name : str, optional
+        The name of the model.
+    **model_kwargs : optional
+        Valid PyBaMM model option keys and their values, for example:
+        options : dict, optional
+            A dictionary of options to customise the behaviour of the PyBaMM model.
+        build : bool, optional
+            If True, the model is built upon creation (default: False).
+    """
+
+    def __init__(self, name="Base Model", **model_kwargs):
+        super().__init__(name=name, **model_kwargs)
+
+    def build_model(self):
+        """
+        Build model variables and equations
+        Credit: PyBaMM
+        """
+        self._build_model()
+
+        self._built = True
+        pybamm.logger.info(f"Finish building {self.name}")
+
+    @staticmethod
+    def set_initial_state(
+        initial_value,
+        parameter_values,
+        direction=None,
+        param=None,
+        inplace=True,
+        options=None,
+        inputs=None,
+        tol=1e-6,
+    ):
+        """
+        Set the value of the initial state of charge.
+
+        Parameters
+        ----------
+        initial_value : float
+            Target initial value.
+            If float, 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 : :class:`pybamm.ParameterValues`
+            Parameters and their corresponding values.
+        param : :class:`pybamm.LithiumIonParameters`, optional
+            The symbolic parameter set to use for the simulation.
+            If not provided, the default parameter set will be used.
+        inplace: bool, optional
+            If True, replace the parameters values in place. Otherwise, return a new set of
+            parameter values. Default is True.
+        options : dict-like, optional
+            A dictionary of options to be passed to the model, see
+            :class:`pybamm.BatteryModelOptions`.
+        inputs : dict, optional
+            A dictionary of input parameters to pass to the model when solving.
+        tol : float, optional
+            The tolerance for the solver used to compute the initial stoichiometries.
+            A lower value results in higher precision but may increase computation time.
+            Default is 1e-6.
+        """
+        parameter_values = parameter_values if inplace else parameter_values.copy()
+
+        if isinstance(initial_value, str) and initial_value.endswith("V"):
+            V_init = float(initial_value[:-1])
+            V_min = parameter_values.evaluate(
+                pybamm.Parameter("Lower voltage cut-off [V]"), inputs=inputs
+            )
+            V_max = parameter_values.evaluate(
+                pybamm.Parameter("Upper voltage cut-off [V]"), inputs=inputs
+            )
+
+            if not V_min - tol <= V_init <= V_max + tol:
+                raise ValueError(
+                    f"Initial voltage {V_init}V is outside the voltage limits ({V_min}, {V_max})."
+                )
+
+            x_0 = Parameter("Minimum negative stoichiometry")
+            x_100 = Parameter("Maximum negative stoichiometry")
+            y_100 = Parameter("Minimum positive stoichiometry")
+            y_0 = Parameter("Maximum positive stoichiometry")
+
+            def ocv_function(soc):
+                sto_p = y_0 - soc * (y_0 - y_100)
+                sto_n = x_0 + soc * (x_100 - x_0)
+                U_p = FunctionParameter(
+                    "Positive electrode OCP [V]",
+                    {"Positive particle stoichiometry": sto_p},
+                )
+                U_n = FunctionParameter(
+                    "Negative electrode OCP [V]",
+                    {"Negative particle stoichiometry": sto_n},
+                )
+
+                return parameter_values.evaluate(U_p - U_n, inputs=inputs)
+
+            inverse_ocv = InverseOCV(ocv_function)
+            soc = inverse_ocv(V_init)
+
+        elif isinstance(initial_value, int | float):
+            soc = initial_value
+
+        else:
+            raise ValueError("Initial value must be a float or a string ending in 'V'.")
+
+        if not 0 <= soc <= 1:
+            raise ValueError("Initial SOC should be between 0 and 1.")
+
+        parameter_values["Initial SoC"] = soc
+
+        return parameter_values