Allow linked design parameters

CHANGELOG.md

@@ -2,6 +2,7 @@
 
 ## Features
 
+- [#532](https://github.com/pybop-team/PyBOP/issues/532) - Allows the user to update linked parameters during design optimisation.
 - [#529](https://github.com/pybop-team/PyBOP/issues/529) - Adds `GravimetricPowerDensity` and `VolumetricPowerDensity` costs, along with the mathjax extension for Sphinx.
 
 ## Optimisations

examples/scripts/linked_parameters.py

@@ -0,0 +1,69 @@
+import pybop
+
+# The aim of this script is to show how to systematically update
+# design parameters which depend on the optimisation parameters.
+
+# Define parameter set and model
+parameter_set = pybop.ParameterSet.pybamm("Chen2020", formation_concentrations=True)
+model = pybop.lithium_ion.SPMe(parameter_set=parameter_set)
+
+# Fitting parameters
+parameters = pybop.Parameters(
+    pybop.Parameter(
+        "Positive electrode thickness [m]",
+        prior=pybop.Gaussian(7.56e-05, 0.1e-05),
+        bounds=[65e-06, 10e-05],
+    ),
+    pybop.Parameter(
+        "Positive electrode active material volume fraction",
+        prior=pybop.Gaussian(0.6, 0.15),
+        bounds=[0.1, 0.9],
+    ),
+)
+
+
+# Define a function to update the linked parameters
+def update_porosity(parameter_set):
+    parameter_set["Positive electrode porosity"] = (
+        1 - parameter_set["Positive electrode active material volume fraction"]
+    )
+
+
+model.update_linked_parameters = update_porosity
+
+# Define test protocol
+experiment = pybop.Experiment(
+    [
+        "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]"]
+
+# Generate problem
+problem = pybop.DesignProblem(
+    model,
+    parameters,
+    experiment,
+    signal=signal,
+    initial_state={"Initial SoC": 1.0},
+    update_capacity=True,
+)
+
+# Define the cost
+cost = pybop.GravimetricEnergyDensity(problem)
+
+# 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 energy density: {cost(optim.x0):.2f} Wh.kg-1")
+print(f"Optimised gravimetric energy density: {cost(x):.2f} Wh.kg-1")
+
+# 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/models/base_model.py

@@ -417,6 +417,20 @@ def classify_parameters(
 
         return standard_parameters
 
+    def update_linked_parameters(self, parameter_set: None):
+        """
+        A placeholder function that may be overwritten by the user in order to
+        update model parameters as a function of the optimisation parameters.
+
+        ***Warning: Currently only for use in design problems.
+
+        Parameters
+        ----------
+        parameter_set : pybamm.ParameterValues
+            A PyBaMM parameter set containing standard lithium ion parameters.
+        """
+        pass
+
     def reinit(
         self, inputs: Inputs, t: float = 0.0, x: Optional[np.ndarray] = None
     ) -> TimeSeriesState:

pybop/problems/design_problem.py

@@ -101,9 +101,9 @@ def set_initial_state(self, initial_state: dict):
 
         self.initial_state = initial_state
 
-    def evaluate(self, inputs: Inputs):
+    def update_parameter_set(self, inputs: Inputs):
         """
-        Evaluate the model with the given parameters and return the signal.
+        Update any linked parameters in the active parameter set.
 
         Parameters
         ----------
@@ -112,20 +112,46 @@ def evaluate(self, inputs: Inputs):
 
         Returns
         -------
-        y : np.ndarray
-            The model output y(t) simulated with inputs.
+        parameter_set : pybamm.ParameterValues
+            An updated version of the parameter set.
         """
-        inputs = self.parameters.verify(inputs)
-
-        # Update the active parameter set
         parameter_set = self.model.parameter_set
+
+        # Reset the initial concentrations to the formation concentrations
         if isinstance(self._model, EChemBaseModel):
             set_formation_concentrations(parameter_set)
+
+        # Update the parameter set and any linked parameters
         parameter_set.update(inputs)
+        if self._model is not None:
+            self._model.update_linked_parameters(parameter_set)
+
+        # Update the nominal capacity in order to update the C-rate
         if self.update_capacity:
             approximate_capacity = self.model.approximate_capacity(parameter_set)
             parameter_set.update({"Nominal cell capacity [A.h]": approximate_capacity})
 
+        return parameter_set
+
+    def evaluate(self, inputs: Inputs):
+        """
+        Evaluate the model with the given parameters and return the signal.
+
+        Parameters
+        ----------
+        inputs : Inputs
+            Parameters for evaluation of the model.
+
+        Returns
+        -------
+        y : np.ndarray
+            The model output y(t) simulated with inputs.
+        """
+        inputs = self.parameters.verify(inputs)
+
+        # Update the active parameter set including any linked parameters
+        parameter_set = self.update_parameter_set(inputs)
+
         try:
             with warnings.catch_warnings():
                 for pattern in self.warning_patterns:

pybop/problems/fitting_problem.py

@@ -101,11 +101,7 @@ def set_initial_state(self, initial_state: Optional[dict] = None):
 
         self.initial_state = initial_state
 
-    def evaluate(
-        self,
-        inputs: Inputs,
-        update_capacity=False,
-    ) -> dict[str, np.ndarray[np.float64]]:
+    def evaluate(self, inputs: Inputs) -> dict[str, np.ndarray[np.float64]]:
         """
         Evaluate the model with the given parameters and return the signal.
 
@@ -121,7 +117,7 @@ def evaluate(
         """
         inputs = self.parameters.verify(inputs)
         if self.eis:
-            return self._evaluateEIS(inputs, update_capacity=update_capacity)
+            return self._evaluateEIS(inputs)
         else:
             try:
                 sol = self._model.simulate(
@@ -139,9 +135,7 @@ def evaluate(
                 for signal in (self.signal + self.additional_variables)
             }
 
-    def _evaluateEIS(
-        self, inputs: Inputs, update_capacity=False
-    ) -> dict[str, np.ndarray[np.float64]]:
+    def _evaluateEIS(self, inputs: Inputs) -> dict[str, np.ndarray[np.float64]]:
         """
         Evaluate the model with the given parameters and return the signal.
 

pybop/problems/multi_fitting_problem.py

@@ -76,7 +76,7 @@ def set_initial_state(self, initial_state: Optional[dict] = None):
         for problem in self.problems:
             problem.set_initial_state(initial_state)
 
-    def evaluate(self, inputs: Inputs, eis=False):
+    def evaluate(self, inputs: Inputs):
         """
         Evaluate the model with the given parameters and return the signal.