Add calc nom. voltage from half-cell ocps, np.trapz dx variable, prin…

…t userwarning, reduce iterations for tests

examples/scripts/spme_max_energy.py

@@ -2,6 +2,9 @@
 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.
+
 # A design optimisation example loosely based on work by L.D. Couto
 # available at https://doi.org/10.1016/j.energy.2022.125966.
 
@@ -20,18 +23,24 @@
     pybop.Parameter(
         "Positive electrode thickness [m]",
         prior=pybop.Gaussian(7.56e-05, 0.05e-05),
-        bounds=[7e-05, 10e-05],
+        bounds=[65e-06, 10e-05],
     ),
     pybop.Parameter(
         "Positive particle radius [m]",
         prior=pybop.Gaussian(5.22e-06, 0.05e-06),
-        bounds=[3e-06, 9e-06],
+        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):
+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.
@@ -41,18 +50,18 @@ def nominal_capacity(x, model):
     }
     model._parameter_set.update(inputs)
 
-    theoretical_energy = model._electrode_soh.calculate_theoretical_energy(
+    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["Upper voltage cut-off [V]"]
-        + model._parameter_set["Lower voltage cut-off [V]"]
-    ) / 2
+    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):
+def cell_mass(parameter_set):  # This is very low compute time
     """
     Compute the total cell mass [kg] for the current parameter set.
     """
@@ -124,6 +133,7 @@ def cell_mass(parameter_set):
 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
@@ -148,23 +158,26 @@ def _evaluate(self, x, grad=None):
 
             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_Ah = np.trapz(voltage * current) / (
+                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_Ah
+            # 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)
 optim = pybop.Optimisation(cost, optimiser=pybop.PSO, verbose=True)
-optim.set_max_iterations(15)
+optim.set_max_iterations(5)
 
 # Run optimisation
 x, final_cost = optim.run()