Add fitting functionality for geometric parameters

CHANGELOG.md

@@ -2,6 +2,7 @@
 
 ## Features
 
+- [#18](https://github.com/pybop-team/PyBOP/pull/18) - Adds geometric parameter fitting capability, via `model.rebuild()` with `model.rebuild_parameters`.
 - [#203](https://github.com/pybop-team/PyBOP/pull/203) - Adds support for modern Python packaging via a `pyproject.toml` file and configures the `pytest` test runner and `ruff` linter to use their configurations stored as declarative metadata.
 - [#123](https://github.com/pybop-team/PyBOP/issues/123) - Configures scheduled tests to run against the last three PyPI releases of PyBaMM via dynamic GitHub Actions matrix generation.
 - [#187](https://github.com/pybop-team/PyBOP/issues/187) - Adds M1 Github runner to `test_on_push` workflow, updt. self-hosted supported python versions in scheduled tests.

examples/scripts/spm_CMAES.py

@@ -8,14 +8,16 @@
 # Fitting parameters
 parameters = [
     pybop.Parameter(
-        "Negative electrode active material volume fraction",
-        prior=pybop.Gaussian(0.6, 0.05),
-        bounds=[0.5, 0.8],
+        "Negative particle radius [m]",
+        prior=pybop.Gaussian(6e-06, 0.1e-6),
+        bounds=[1e-6, 9e-6],
+        true_value=parameter_set["Negative particle radius [m]"],
     ),
     pybop.Parameter(
-        "Positive electrode active material volume fraction",
-        prior=pybop.Gaussian(0.48, 0.05),
-        bounds=[0.4, 0.7],
+        "Positive particle radius [m]",
+        prior=pybop.Gaussian(4.5e-06, 0.1e-6),
+        bounds=[1e-6, 9e-6],
+        true_value=parameter_set["Positive particle radius [m]"],
     ),
 ]
 
@@ -42,6 +44,13 @@
 
 # Run the optimisation
 x, final_cost = optim.run()
+print(
+    "True parameters:",
+    [
+        parameters[0].true_value,
+        parameters[1].true_value,
+    ],
+)
 print("Estimated parameters:", x)
 
 # Plot the timeseries output
@@ -57,5 +66,4 @@
 pybop.plot_cost2d(cost, steps=15)
 
 # Plot the cost landscape with optimisation path and updated bounds
-bounds = np.array([[0.55, 0.75], [0.48, 0.68]])
-pybop.plot_cost2d(cost, optim=optim, bounds=bounds, steps=15)
+pybop.plot_cost2d(cost, optim=optim, steps=15)

examples/standalone/model.py

@@ -56,3 +56,4 @@ def __init__(
         self._built_initial_soc = None
         self._mesh = None
         self._disc = None
+        self.rebuild_parameters = {}

pybop/_problem.py

@@ -61,7 +61,7 @@ def __init__(
 
         # Add the initial values to the parameter definitions
         for i, param in enumerate(self.parameters):
-            param.update(value=self.x0[i])
+            param.update(initial_value=self.x0[i])
 
     def evaluate(self, x):
         """
@@ -118,6 +118,10 @@ def target(self):
         """
         return self._target
 
+    @property
+    def model(self):
+        return self._model
+
 
 class FittingProblem(BaseProblem):
     """
@@ -149,6 +153,7 @@ def __init__(
     ):
         super().__init__(parameters, model, check_model, signal, init_soc, x0)
         self._dataset = dataset.data
+        self.x = self.x0
 
         # Check that the dataset contains time and current
         for name in ["Time [s]", "Current function [A]"] + self.signal:
@@ -199,6 +204,13 @@ def evaluate(self, x):
         y : np.ndarray
             The model output y(t) simulated with inputs x.
         """
+        if (x != self.x).any() and self._model.matched_parameters:
+            for i, param in enumerate(self.parameters):
+                param.update(value=x[i])
+
+            self._model.rebuild(parameters=self.parameters)
+            self.x = x
+
         y = np.asarray(self._model.simulate(inputs=x, t_eval=self._time_data))
 
         return y
@@ -218,6 +230,10 @@ def evaluateS1(self, x):
             A tuple containing the simulation result y(t) and the sensitivities dy/dx(t) evaluated
             with given inputs x.
         """
+        if self._model.matched_parameters:
+            raise RuntimeError(
+                "Gradient not available when using geometric parameters."
+            )
 
         y, dy = self._model.simulateS1(
             inputs=x,

pybop/costs/design_costs.py

@@ -41,7 +41,7 @@ def __init__(self, problem, update_capacity=False):
             )
         warnings.warn(nominal_capacity_warning, UserWarning)
         self.update_capacity = update_capacity
-        self.parameter_set = problem._model._parameter_set
+        self.parameter_set = problem.model.parameter_set
         self.update_simulation_data(problem.x0)
 
     def update_simulation_data(self, initial_conditions):
@@ -54,7 +54,7 @@ def update_simulation_data(self, initial_conditions):
             The initial conditions for the simulation.
         """
         if self.update_capacity:
-            self.problem._model.approximate_capacity(self.problem.x0)
+            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]
@@ -116,12 +116,12 @@ def _evaluate(self, x, grad=None):
                 warnings.filterwarnings("error", category=UserWarning)
 
                 if self.update_capacity:
-                    self.problem._model.approximate_capacity(x)
+                    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)
+                    3600 * self.problem.model.cell_mass(self.parameter_set)
                 )
 
                 return negative_energy_density
@@ -166,12 +166,12 @@ def _evaluate(self, x, grad=None):
                 warnings.filterwarnings("error", category=UserWarning)
 
                 if self.update_capacity:
-                    self.problem._model.approximate_capacity(x)
+                    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_volume(self.parameter_set)
+                    3600 * self.problem.model.cell_volume(self.parameter_set)
                 )
 
                 return negative_energy_density

pybop/costs/fitting_costs.py

@@ -37,7 +37,6 @@ def _evaluate(self, x, grad=None):
             The root mean square error.
 
         """
-
         prediction = self.problem.evaluate(x)
 
         if len(prediction) < len(self._target):

pybop/models/base_model.py

@@ -2,6 +2,7 @@
 from dataclasses import dataclass
 from typing import Any, Dict, Optional
 import pybamm
+import copy
 import numpy as np
 import casadi
 
@@ -57,6 +58,9 @@ def __init__(self, name="Base Model"):
         self.parameters = None
         self.dataset = None
         self.signal = None
+        self.matched_parameters = {}
+        self.non_matched_parameters = {}
+        self.fit_keys = []
         self.param_check_counter = 0
         self.allow_infeasible_solutions = True
 
@@ -87,10 +91,11 @@ def build(
         """
         self.dataset = dataset
         self.parameters = parameters
-        if self.parameters is None:
-            self.fit_keys = []
-        else:
+        if self.parameters is not None:
+            self.set_parameter_classification(self.parameters)
             self.fit_keys = [param.name for param in self.parameters]
+        else:
+            self.fit_keys = []
 
         if init_soc is not None:
             self.set_init_soc(init_soc)
@@ -140,24 +145,24 @@ def set_init_soc(self, init_soc):
         # Save solved initial SOC in case we need to rebuild the model
         self._built_initial_soc = init_soc
 
-    def set_params(self):
+    def set_params(self, rebuild=False):
         """
         Assign the parameters to the model.
 
         This method processes the model with the given parameters, sets up
         the geometry, and updates the model instance.
         """
-        if self.model_with_set_params:
+        if self.model_with_set_params and not rebuild:
             return
 
         # Mark any simulation inputs in the parameter set
-        if self.parameters is not None:
+        if self.non_matched_parameters:
             for i in self.fit_keys:
                 self._parameter_set[i] = "[input]"
 
-        if self.dataset is not None and self.parameters is not None:
+        if self.dataset is not None and self.non_matched_parameters:
             if "Current function [A]" not in self.fit_keys:
-                self.parameter_set["Current function [A]"] = pybamm.Interpolant(
+                self._parameter_set["Current function [A]"] = pybamm.Interpolant(
                     self.dataset["Time [s]"],
                     self.dataset["Current function [A]"],
                     pybamm.t,
@@ -172,6 +177,93 @@ def set_params(self):
             self._parameter_set.process_geometry(self.geometry)
         self.pybamm_model = self._model_with_set_params
 
+    def rebuild(
+        self,
+        dataset=None,
+        parameters=None,
+        parameter_set=None,
+        check_model=True,
+        init_soc=None,
+    ):
+        """
+        Rebuild the PyBaMM model for a given parameter set.
+
+        This method requires the self.build() method to be called first, and
+        then rebuilds the model for a given parameter set. Specifically,
+        this method applies the given parameters, sets up the mesh and discretization if needed, and prepares the model
+        for simulations.
+
+        Parameters
+        ----------
+        dataset : pybamm.Dataset, optional
+            The dataset to be used in the model construction.
+        parameters : dict, optional
+            A dictionary containing parameter values to apply to the model.
+        parameter_set : pybop.parameter_set, optional
+            A PyBOP parameter set object or a dictionary containing the parameter values
+        check_model : bool, optional
+            If True, the model will be checked for correctness after construction.
+        init_soc : float, optional
+            The initial state of charge to be used in simulations.
+        """
+        self.dataset = dataset
+        self.parameters = parameters
+        if parameters is not None:
+            self.set_parameter_classification(parameters)
+
+        if init_soc is not None:
+            self.set_init_soc(init_soc)
+
+        if self._built_model is None:
+            raise ValueError("Model must be built before calling rebuild")
+
+        self.set_params(rebuild=True)
+        self._mesh = pybamm.Mesh(self.geometry, self.submesh_types, self.var_pts)
+        self._disc = pybamm.Discretisation(self.mesh, self.spatial_methods)
+        self._built_model = self._disc.process_model(
+            self._model_with_set_params, inplace=False, check_model=check_model
+        )
+
+        # Clear solver and setup model
+        self._solver._model_set_up = {}
+
+    def set_parameter_classification(self, parameters):
+        """
+        Set the parameter classification for the model.
+
+        Parameters
+        ----------
+        parameters : Pybop.ParameterSet
+
+        Returns
+        -------
+        None
+            The method updates attributes on self.
+
+        """
+        processed_parameters = {param.name: param.value for param in parameters}
+        matched_parameters = {
+            param: processed_parameters[param]
+            for param in processed_parameters
+            if param in self.rebuild_parameters
+        }
+        non_matched_parameters = {
+            param: processed_parameters[param]
+            for param in processed_parameters
+            if param not in self.rebuild_parameters
+        }
+
+        self.matched_parameters.update(matched_parameters)
+        self.non_matched_parameters.update(non_matched_parameters)
+
+        if self.matched_parameters:
+            self._parameter_set.update(self.matched_parameters)
+            self._unprocessed_parameter_set = self._parameter_set
+            self.geometry = self.pybamm_model.default_geometry
+
+        if self.non_matched_parameters:
+            self.fit_keys = list(self.non_matched_parameters.keys())
+
     def reinit(
         self, inputs: Inputs, t: float = 0.0, x: Optional[np.ndarray] = None
     ) -> TimeSeriesState:
@@ -243,25 +335,29 @@ def simulate(self, inputs, t_eval) -> np.ndarray[np.float64]:
         ValueError
             If the model has not been built before simulation.
         """
-
         if self._built_model is None:
             raise ValueError("Model must be built before calling simulate")
         else:
-            if not isinstance(inputs, dict):
-                inputs = {key: inputs[i] for i, key in enumerate(self.fit_keys)}
+            if not self.fit_keys and self.matched_parameters:
+                sol = self.solver.solve(self.built_model, t_eval=t_eval)
 
-            if self.check_params(
-                inputs=inputs,
-                allow_infeasible_solutions=self.allow_infeasible_solutions,
-            ):
-                sol = self._solver.solve(self.built_model, inputs=inputs, t_eval=t_eval)
+            else:
+                if not isinstance(inputs, dict):
+                    inputs = {key: inputs[i] for i, key in enumerate(self.fit_keys)}
 
-                predictions = [sol[signal].data for signal in self.signal]
+                if self.check_params(
+                    inputs=inputs,
+                    allow_infeasible_solutions=self.allow_infeasible_solutions,
+                ):
+                    sol = self.solver.solve(
+                        self.built_model, inputs=inputs, t_eval=t_eval
+                    )
+                else:
+                    return [np.inf]
 
-                return np.vstack(predictions).T
+            predictions = [sol[signal].data for signal in self.signal]
 
-            else:
-                return [np.inf]
+            return np.vstack(predictions).T
 
     def simulateS1(self, inputs, t_eval):
         """
@@ -449,6 +545,17 @@ def _check_params(self, inputs=None, allow_infeasible_solutions=True):
         """
         return True
 
+    def copy(self):
+        """
+        Return a copy of the model.
+
+        Returns
+        -------
+        BaseModel
+            A copy of the model.
+        """
+        return copy.copy(self)
+
     def cell_mass(self, parameter_set=None):
         """
         Calculate the cell mass in kilograms.

pybop/models/empirical/ecm.py

@@ -75,6 +75,7 @@ def __init__(
         self._built_initial_soc = None
         self._mesh = None
         self._disc = None
+        self.rebuild_parameters = {}
 
     def _check_params(self, inputs=None, allow_infeasible_solutions=True):
         """