Add PowerDensity costs (#530)

* Add PowerDensity costs

* Add integrals to docstrings

* Use mathjax for sphinx

* Insert line breaks

* Update CHANGELOG.md

---------

Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>

CHANGELOG.md

@@ -2,6 +2,8 @@
 
 ## Features
 
+- [#529](https://github.com/pybop-team/PyBOP/issues/529) - Adds `GravimetricPowerDensity` and `VolumetricPowerDensity` costs, along with the mathjax extension for Sphinx.
+
 ## Optimisations
 
 - [#512](https://github.com/pybop-team/PyBOP/pull/513) - Refactors `LogPosterior` with attributes pointing to composed likelihood object.
@@ -16,7 +18,6 @@
 
 # [v24.9.1](https://github.com/pybop-team/PyBOP/tree/v24.9.0) - 2024-09-16
 
-
 ## Features
 
 ## Bug Fixes

docs/conf.py

@@ -20,12 +20,13 @@
     "sphinx.ext.napoleon",
     "sphinx.ext.autodoc",
     "sphinx.ext.autosummary",
+    "sphinx.ext.mathjax",
     "sphinx.ext.todo",
     "sphinx.ext.viewcode",
     "sphinx_design",
     "sphinx_copybutton",
     "autoapi.extension",
-    # custom extentions
+    # custom extensions
     "_extension.gallery_directive",
     # For extension examples and demos
     "myst_parser",

examples/scripts/spme_max_energy.py → examples/scripts/maximising_energy.py

@@ -3,8 +3,8 @@
 # A design optimisation example loosely based on work by L.D. Couto
 # available at https://doi.org/10.1016/j.energy.2022.125966.
 
-# The target is to maximise the gravimetric energy density over a
-# range of possible design parameter values, including for example:
+# The target is to maximise the energy density over a range of
+# possible design parameter values, including for example:
 # cross-sectional area = height x width (only need change one)
 # electrode widths, particle radii, volume fractions and
 # separator width.
@@ -29,7 +29,10 @@
 
 # Define test protocol
 experiment = pybop.Experiment(
-    ["Discharge at 1C until 2.5 V (5 seconds period)"],
+    [
+        "Discharge at 1C until 2.5 V (5 seconds period)",
+        "Hold at 2.5 V for 30 minutes or until 10 mA (5 seconds period)",
+    ],
 )
 signal = ["Voltage [V]", "Current [A]"]
 

examples/scripts/maximising_power.py

@@ -0,0 +1,62 @@
+import pybop
+
+# Define parameter set and model
+parameter_set = pybop.ParameterSet.pybamm("Chen2020", formation_concentrations=True)
+model = pybop.lithium_ion.SPMe(parameter_set=parameter_set)
+
+# Define useful quantities
+nominal_capacity = parameter_set["Nominal cell capacity [A.h]"]
+target_c_rate = 2
+discharge_rate = target_c_rate * nominal_capacity
+
+# Fitting parameters
+parameters = pybop.Parameters(
+    pybop.Parameter(
+        "Positive electrode thickness [m]",
+        prior=pybop.Gaussian(7.56e-05, 0.5e-05),
+        bounds=[65e-06, 10e-05],
+    ),
+    pybop.Parameter(
+        "Nominal cell capacity [A.h]",  # controls the C-rate in the experiment
+        prior=pybop.Gaussian(discharge_rate, 0.2),
+        bounds=[0.8 * discharge_rate, 1.2 * discharge_rate],
+    ),
+)
+
+# Define test protocol
+experiment = pybop.Experiment(
+    ["Discharge at 1C for 30 minutes or until 2.5 V (5 seconds period)"],
+)
+signal = ["Voltage [V]", "Current [A]"]
+
+# Generate problem
+problem = pybop.DesignProblem(
+    model,
+    parameters,
+    experiment,
+    signal=signal,
+    initial_state={"Initial SoC": 1.0},
+)
+
+# Generate multiple cost functions and combine them
+cost1 = pybop.GravimetricPowerDensity(problem, target_time=3600 / target_c_rate)
+cost2 = pybop.VolumetricPowerDensity(problem, target_time=3600 / target_c_rate)
+cost = pybop.WeightedCost(cost1, cost2, weights=[1, 1e-3])
+
+# Run optimisation
+optim = pybop.XNES(
+    cost, verbose=True, allow_infeasible_solutions=False, max_iterations=10
+)
+x, final_cost = optim.run()
+print("Estimated parameters:", x)
+print(f"Initial gravimetric power density: {cost1(optim.x0):.2f} W.kg-1")
+print(f"Optimised gravimetric power density: {cost1(x):.2f} W.kg-1")
+print(f"Initial volumetric power density: {cost2(optim.x0):.2f} W.m-3")
+print(f"Optimised volumetric power density: {cost2(x):.2f} W.m-3")
+print(f"Optimised discharge rate: {x[-1]:.2f} A = {x[-1]/nominal_capacity:.2f} C")
+
+# Plot the timeseries output
+pybop.quick_plot(problem, problem_inputs=x, title="Optimised Comparison")
+
+# Plot the cost landscape with optimisation path
+pybop.plot2d(optim, steps=5)

pybop/__init__.py

@@ -105,6 +105,8 @@
     DesignCost,
     GravimetricEnergyDensity,
     VolumetricEnergyDensity,
+    GravimetricPowerDensity,
+    VolumetricPowerDensity,
 )
 from .costs._likelihoods import (
     BaseLikelihood,

pybop/costs/design_costs.py

@@ -30,11 +30,20 @@ def __init__(self, problem):
 
 class GravimetricEnergyDensity(DesignCost):
     """
-    Represents the gravimetric energy density of a battery cell, calculated based
-    on a normalised discharge from upper to lower voltage limits. The goal is to
-    maximise the energy density, which is achieved by setting minimising = False
+    Calculates the gravimetric energy density (specific energy) of a battery cell,
+    when applied to a normalised discharge from upper to lower voltage limits. The
+    goal of maximising the energy density is achieved by setting minimising = False
     in the optimiser settings.
 
+    The gravimetric energy density [Wh.kg-1] is calculated as
+
+    .. math::
+        \\frac{1}{3600 m} \\int_{t=0}^{t=T} I(t) V(t) \\mathrm{d}t
+
+    where m is the cell mass, t is the time, T is the total time, I is the current
+    and V is the voltage. The factor of 1/3600 is included to convert from seconds
+    to hours.
+
     Inherits all parameters and attributes from ``DesignCost``.
     """
 
@@ -79,10 +88,19 @@ def compute(
 
 class VolumetricEnergyDensity(DesignCost):
     """
-    Represents the volumetric energy density of a battery cell, calculated based
-    on a normalised discharge from upper to lower voltage limits. The goal is to
-    maximise the energy density, which is achieved by setting minimising = False
-    in the optimiser settings.
+    Calculates the (volumetric) energy density of a battery cell, when applied to a
+    normalised discharge from upper to lower voltage limits. The goal of maximising
+    the energy density is achieved by setting minimising = False in the optimiser
+    settings.
+
+    The volumetric energy density [Wh.m-3] is calculated as
+
+    .. math::
+        \\frac{1}{3600 v} \\int_{t=0}^{t=T} I(t) V(t) \\mathrm{d}t
+
+    where v is the cell volume, t is the time, T is the total time, I is the current
+    and V is the voltage. The factor of 1/3600 is included to convert from seconds
+    to hours.
 
     Inherits all parameters and attributes from ``DesignCost``.
     """
@@ -124,3 +142,130 @@ def compute(
         )
 
         return energy_density
+
+
+class GravimetricPowerDensity(DesignCost):
+    """
+    Calculates the gravimetric power density (specific power) of a battery cell,
+    when applied to a discharge from upper to lower voltage limits. The goal of
+    maximising the power density is achieved by setting minimising = False in the
+    optimiser settings.
+
+    The time-averaged gravimetric power density [W.kg-1] is calculated as
+
+    .. math::
+        \\frac{1}{3600 m T} \\int_{t=0}^{t=T} I(t) V(t) \\mathrm{d}t
+
+    where m is the cell mass, t is the time, T is the total time, I is the current
+    and V is the voltage. The factor of 1/3600 is included to convert from seconds
+    to hours.
+
+    Inherits all parameters and attributes from ``DesignCost``.
+
+    Additional parameters
+    ---------------------
+    target_time : int
+        The length of time (seconds) over which the power should be sustained.
+    """
+
+    def __init__(self, problem, target_time: int = 3600):
+        super().__init__(problem)
+        self.target_time = target_time
+
+    def compute(
+        self,
+        y: dict,
+        dy: np.ndarray = None,
+        calculate_grad: bool = False,
+    ) -> float:
+        """
+        Computes the cost function for the given predictions.
+
+        Parameters
+        ----------
+        y : dict
+            The dictionary of predictions with keys designating the signals for fitting.
+        dy : np.ndarray, optional
+            The corresponding gradient with respect to the parameters for each signal.
+            Note: not used in design optimisation classes.
+        calculate_grad : bool, optional
+            A bool condition designating whether to calculate the gradient.
+
+        Returns
+        -------
+        float
+            The gravimetric power density or -infinity in case of infeasible parameters.
+        """
+        if not any(np.isfinite(y[signal][0]) for signal in self.signal):
+            return -np.inf
+
+        voltage, current = y["Voltage [V]"], y["Current [A]"]
+        dt = y["Time [s]"][1] - y["Time [s]"][0]
+        time_averaged_power_density = np.trapz(voltage * current, dx=dt) / (
+            self.target_time * 3600 * self.problem.model.cell_mass()
+        )
+
+        return time_averaged_power_density
+
+
+class VolumetricPowerDensity(DesignCost):
+    """
+    Calculates the (volumetric) power density of a battery cell, when applied to a
+    discharge from upper to lower voltage limits. The goal of maximising the power
+    density is achieved by setting minimising = False in the optimiser settings.
+
+    The time-averaged volumetric power density [W.m-3] is calculated as
+
+    .. math::
+        \\frac{1}{3600 v T} \\int_{t=0}^{t=T} I(t) V(t) \\mathrm{d}t
+
+    where v is the cell volume, t is the time, T is the total time, I is the current
+    and V is the voltage. The factor of 1/3600 is included to convert from seconds
+    to hours.
+
+    Inherits all parameters and attributes from ``DesignCost``.
+
+    Additional parameters
+    ---------------------
+    target_time : int
+        The length of time (seconds) over which the power should be sustained.
+    """
+
+    def __init__(self, problem, target_time: int = 3600):
+        super().__init__(problem)
+        self.target_time = target_time
+
+    def compute(
+        self,
+        y: dict,
+        dy: np.ndarray = None,
+        calculate_grad: bool = False,
+    ) -> float:
+        """
+        Computes the cost function for the given predictions.
+
+        Parameters
+        ----------
+        y : dict
+            The dictionary of predictions with keys designating the signals for fitting.
+        dy : np.ndarray, optional
+            The corresponding gradient with respect to the parameters for each signal.
+            Note: not used in design optimisation classes.
+        calculate_grad : bool, optional
+            A bool condition designating whether to calculate the gradient.
+
+        Returns
+        -------
+        float
+            The volumetric power density or -infinity in case of infeasible parameters.
+        """
+        if not any(np.isfinite(y[signal][0]) for signal in self.signal):
+            return -np.inf
+
+        voltage, current = y["Voltage [V]"], y["Current [A]"]
+        dt = y["Time [s]"][1] - y["Time [s]"][0]
+        time_averaged_power_density = np.trapz(voltage * current, dx=dt) / (
+            self.target_time * 3600 * self.problem.model.cell_volume()
+        )
+
+        return time_averaged_power_density

tests/unit/test_cost.py

@@ -219,6 +219,8 @@ def design_problem(self, parameters, experiment, signal):
             pybop.DesignCost,
             pybop.GravimetricEnergyDensity,
             pybop.VolumetricEnergyDensity,
+            pybop.GravimetricPowerDensity,
+            pybop.VolumetricPowerDensity,
         ],
     )
     @pytest.mark.unit