Merge pull request #167 from pybop-team/162-update-design-optimisatio…

…n-functionality

Design Optimisation Refactor and Additions

examples/scripts/spme_max_energy.py

@@ -1,6 +1,4 @@
 import pybop
-import numpy as np
-import warnings
 
 ## NOTE: This is a brittle example, the classes and methods below will be
 ## integrated into pybop in a future release.
@@ -31,90 +29,6 @@
         bounds=[2e-06, 9e-06],
     ),
 ]
-# Define stoichiometries
-sto = model._electrode_soh.get_min_max_stoichiometries(parameter_set)
-mean_sto_neg = np.mean(sto[0:2])
-mean_sto_pos = np.mean(sto[2:4])
-
-
-# Define functions
-def nominal_capacity(
-    x, model
-):  # > 50% of the simulation time is spent in this function (~0.7 sec per iteration vs ~1.1 for the forward simulation)
-    """
-    Update the nominal capacity based on the theoretical energy density and an
-    average voltage.
-    """
-    inputs = {
-        key: x[i] for i, key in enumerate([param.name for param in model.parameters])
-    }
-    model._parameter_set.update(inputs)
-
-    theoretical_energy = model._electrode_soh.calculate_theoretical_energy(  # All of the computational time is in this line (~0.7)
-        model._parameter_set
-    )
-    average_voltage = model._parameter_set["Positive electrode OCP [V]"](
-        mean_sto_pos
-    ) - model._parameter_set["Negative electrode OCP [V]"](mean_sto_neg)
-
-    theoretical_capacity = theoretical_energy / average_voltage
-    model._parameter_set.update({"Nominal cell capacity [A.h]": theoretical_capacity})
-
-
-def cell_mass(parameter_set):  # This is very low compute time
-    """
-    Compute the total cell mass [kg] for the current parameter set.
-    """
-
-    # Approximations due to SPM(e) parameter set limitations
-    electrolyte_density = parameter_set["Separator density [kg.m-3]"]
-
-    # Electrode mass densities [kg.m-3]
-    positive_mass_density = (
-        parameter_set["Positive electrode active material volume fraction"]
-        * parameter_set["Positive electrode density [kg.m-3]"]
-    )
-    +(parameter_set["Positive electrode porosity"] * electrolyte_density)
-    negative_mass_density = (
-        parameter_set["Negative electrode active material volume fraction"]
-        * parameter_set["Negative electrode density [kg.m-3]"]
-    )
-    +(parameter_set["Negative electrode porosity"] * electrolyte_density)
-
-    # Area densities [kg.m-2]
-    positive_area_density = (
-        parameter_set["Positive electrode thickness [m]"] * positive_mass_density
-    )
-    negative_area_density = (
-        parameter_set["Negative electrode thickness [m]"] * negative_mass_density
-    )
-    separator_area_density = (
-        parameter_set["Separator thickness [m]"]
-        * parameter_set["Separator porosity"]
-        * parameter_set["Separator density [kg.m-3]"]
-    )
-    positive_current_collector_area_density = (
-        parameter_set["Positive current collector thickness [m]"]
-        * parameter_set["Positive current collector density [kg.m-3]"]
-    )
-    negative_current_collector_area_density = (
-        parameter_set["Negative current collector thickness [m]"]
-        * parameter_set["Negative current collector density [kg.m-3]"]
-    )
-
-    # Cross-sectional area [m2]
-    cross_sectional_area = (
-        parameter_set["Electrode height [m]"] * parameter_set["Electrode width [m]"]
-    )
-
-    return cross_sectional_area * (
-        positive_area_density
-        + separator_area_density
-        + negative_area_density
-        + positive_current_collector_area_density
-        + negative_current_collector_area_density
-    )
-
 
 # Define test protocol
 experiment = pybop.Experiment(
@@ -128,58 +42,12 @@ def cell_mass(parameter_set):  # This is very low compute time
     model, parameters, experiment, signal=signal, init_soc=init_soc
 )
 
-# Update the C-rate and the example dataset
-nominal_capacity(problem.x0, model)
-sol = problem.evaluate(problem.x0)
-problem._time_data = sol[:, -1]
-problem._target = sol[:, 0:-1]
-dt = sol[1, -1] - sol[0, -1]
-
-
-# Define cost function as a subclass
-class GravimetricEnergyDensity(pybop.BaseCost):
-    """
-    Defines the (negative*) gravimetric energy density corresponding to a
-    normalised 1C discharge from upper to lower voltage limits.
-    *The energy density is maximised by minimising the negative energy density.
-    """
-
-    def __init__(self, problem):
-        super().__init__(problem)
-
-    def _evaluate(self, x, grad=None):
-        """
-        Compute the cost
-        """
-        with warnings.catch_warnings(record=True) as w:
-            # Update the C-rate and run the simulation
-            nominal_capacity(x, self.problem._model)
-            sol = self.problem.evaluate(x)
-
-            if any(w) and issubclass(w[-1].category, UserWarning):
-                # Catch infeasible parameter combinations e.g. E_Li > Q_p
-                print(f"ignoring this sample due to: {w[-1].message}")
-                return np.inf
-
-            else:
-                voltage = sol[:, 0]
-                current = sol[:, 1]
-                gravimetric_energy_density = -np.trapz(
-                    voltage * current, dx=dt
-                ) / (  # trapz over-estimates compares to pybamm (~0.5%)
-                    3600 * cell_mass(self.problem._model._parameter_set)
-                )
-            # Return the negative energy density, as the optimiser minimises
-            # this function, to carry out maximisation of the energy density
-            return gravimetric_energy_density
-
-
 # Generate cost function and optimisation class
-cost = GravimetricEnergyDensity(problem)
+cost = pybop.GravimetricEnergyDensity(problem)
 optim = pybop.Optimisation(
     cost, optimiser=pybop.PSO, verbose=True, allow_infeasible_solutions=False
 )
-optim.set_max_iterations(5)
+optim.set_max_iterations(15)
 
 # Run optimisation
 x, final_cost = optim.run()
@@ -188,7 +56,7 @@ def _evaluate(self, x, grad=None):
 print(f"Optimised gravimetric energy density: {-final_cost:.2f} Wh.kg-1")
 
 # Plot the timeseries output
-nominal_capacity(x, cost.problem._model)
+# model.approximate_capacity(x)
 pybop.quick_plot(x, cost, title="Optimised Comparison")
 
 # Plot the cost landscape with optimisation path

pybop/__init__.py

@@ -26,7 +26,13 @@
 #
 # Cost function class
 #
-from ._costs import BaseCost, RootMeanSquaredError, SumSquaredError, ObserverCost
+from ._costs import (
+    BaseCost,
+    RootMeanSquaredError,
+    SumSquaredError,
+    ObserverCost,
+    GravimetricEnergyDensity,
+)
 
 #
 # Dataset class

pybop/_costs.py

@@ -1,4 +1,5 @@
 import numpy as np
+import warnings
 
 from pybop.observers.observer import Observer
 
@@ -288,6 +289,91 @@ def set_fail_gradient(self, de):
         self._de = de
 
 
+class GravimetricEnergyDensity(BaseCost):
+    """
+    Represents the gravimetric energy density of a battery cell, calculated based
+    on a normalized discharge from upper to lower voltage limits. The goal is to
+    maximize the energy density, which is achieved by minimizing the negative energy
+    density reported by this class.
+
+    Attributes:
+        problem (object): The associated problem containing model and evaluation methods.
+        parameter_set (object): The set of parameters from the problem's model.
+        dt (float): The time step size used in the simulation.
+    """
+
+    def __init__(self, problem, update_capacity=False):
+        """
+        Initializes the gravimetric energy density calculator with a problem.
+
+        Args:
+            problem (object): The problem instance containing the model and data.
+        """
+        super().__init__(problem)
+        self.problem = problem
+        if update_capacity is True:
+            nominal_capacity_warning = (
+                "The nominal capacity is approximated for each iteration."
+            )
+        else:
+            nominal_capacity_warning = (
+                "The nominal capacity is fixed at the initial model value."
+            )
+        warnings.warn(nominal_capacity_warning, UserWarning)
+        self.update_capacity = update_capacity
+        self.parameter_set = problem._model._parameter_set
+        self.update_simulation_data(problem.x0)
+
+    def update_simulation_data(self, initial_conditions):
+        """
+        Updates the simulation data based on the initial conditions.
+
+        Args:
+            initial_conditions (array): The initial conditions for the simulation.
+        """
+        if self.update_capacity:
+            self.problem._model.approximate_capacity(self.problem.x0)
+        solution = self.problem.evaluate(initial_conditions)
+        self.problem._time_data = solution[:, -1]
+        self.problem._target = solution[:, 0:-1]
+        self.dt = solution[1, -1] - solution[0, -1]
+
+    def _evaluate(self, x, grad=None):
+        """
+        Computes the cost function for the energy density.
+
+        Args:
+            x (array): The parameter set for which to compute the cost.
+            grad (array, optional): Gradient information, not used in this method.
+
+        Returns:
+            float: The negative gravimetric energy density or infinity in case of infeasible parameters.
+        """
+        try:
+            with warnings.catch_warnings():
+                # Convert UserWarning to an exception
+                warnings.filterwarnings("error", category=UserWarning)
+
+                if self.update_capacity:
+                    self.problem._model.approximate_capacity(x)
+                solution = self.problem.evaluate(x)
+
+                voltage, current = solution[:, 0], solution[:, 1]
+                negative_energy_density = -np.trapz(voltage * current, dx=self.dt) / (
+                    3600 * self.problem._model.cell_mass(self.parameter_set)
+                )
+
+                return negative_energy_density
+
+        except UserWarning as e:
+            print(f"Ignoring this sample due to: {e}")
+            return np.inf
+
+        except Exception as e:
+            print(f"An error occurred during the evaluation: {e}")
+            return np.inf
+
+
 class ObserverCost(BaseCost):
     """
     Observer cost function.

pybop/models/base_model.py

@@ -441,6 +441,44 @@ def _check_params(self, inputs=None, allow_infeasible_solutions=True):
         """
         return True
 
+    def cell_mass(self, parameter_set=None):
+        """
+        Calculate the cell mass in kilograms.
+
+        This method must be implemented by subclasses.
+
+        Parameters
+        ----------
+        parameter_set : dict, optional
+            A dictionary containing the parameter values necessary for the mass
+            calculations.
+
+        Raises
+        ------
+        NotImplementedError
+            If the method has not been implemented by the subclass.
+        """
+        raise NotImplementedError
+
+    def approximate_capacity(self, x):
+        """
+        Calculate a new estimate for the nominal capacity based on the theoretical energy density
+        and an average voltage.
+
+        This method must be implemented by subclasses.
+
+        Parameters
+        ----------
+        x : array-like
+            An array of values representing the model inputs.
+
+        Raises
+        ------
+        NotImplementedError
+            If the method has not been implemented by the subclass.
+        """
+        raise NotImplementedError
+
     @property
     def built_model(self):
         return self._built_model

pybop/models/empirical/ecm.py

@@ -1,8 +1,8 @@
 import pybamm
-from ..base_model import BaseModel
+from .ecm_base import ECircuitModel
 
 
-class Thevenin(BaseModel):
+class Thevenin(ECircuitModel):
     """
     The Thevenin class represents an equivalent circuit model based on the Thevenin model in PyBaMM.
 

pybop/models/empirical/ecm_base.py

@@ -0,0 +1,29 @@
+from ..base_model import BaseModel
+
+
+class ECircuitModel(BaseModel):
+    """
+    Overwrites and extends `BaseModel` class for circuit-based PyBaMM models.
+    """
+
+    def __init__(self):
+        super().__init__()
+
+    def _check_params(self, inputs=None, allow_infeasible_solutions=True):
+        """
+        Check the compatibility of the model parameters.
+
+        Parameters
+        ----------
+        inputs : dict
+            The input parameters for the simulation.
+        allow_infeasible_solutions : bool, optional
+            If True, infeasible parameter values will be allowed in the optimisation (default: True).
+
+        Returns
+        -------
+        bool
+            A boolean which signifies whether the parameters are compatible.
+
+        """
+        return True