Add MultiFittingProblem, example and test

examples/scripts/multi_fitting.py

@@ -0,0 +1,79 @@
+import numpy as np
+
+import pybop
+
+# Parameter set and model definition
+parameter_set = pybop.ParameterSet.pybamm("Chen2020")
+model_1 = pybop.lithium_ion.SPM(parameter_set=parameter_set)
+model_2 = pybop.lithium_ion.SPM(parameter_set=parameter_set.copy())
+
+# Fitting parameters
+parameters = pybop.Parameters(
+    pybop.Parameter(
+        "Negative electrode active material volume fraction",
+        prior=pybop.Gaussian(0.68, 0.05),
+        true_value=parameter_set["Negative electrode active material volume fraction"],
+    ),
+    pybop.Parameter(
+        "Positive electrode active material volume fraction",
+        prior=pybop.Gaussian(0.58, 0.05),
+        true_value=parameter_set["Positive electrode active material volume fraction"],
+    ),
+)
+
+# Generate a dataset
+sigma = 0.001
+experiment_1 = pybop.Experiment([("Discharge at 0.5C for 2 minutes (4 second period)")])
+values_1 = model_1.predict(experiment=experiment_1)
+dataset_1 = pybop.Dataset(
+    {
+        "Time [s]": values_1["Time [s]"].data,
+        "Current function [A]": values_1["Current [A]"].data,
+        "Voltage [V]": values_1["Voltage [V]"].data
+        + np.random.normal(0, sigma, len(values_1["Voltage [V]"].data)),
+    }
+)
+
+# Generate a second dataset
+experiment_2 = pybop.Experiment([("Discharge at 1C for 2 minutes (4 second period)")])
+values_2 = model_2.predict(experiment=experiment_2)
+dataset_2 = pybop.Dataset(
+    {
+        "Time [s]": values_2["Time [s]"].data,
+        "Current function [A]": values_2["Current [A]"].data,
+        "Voltage [V]": values_2["Voltage [V]"].data
+        + np.random.normal(0, sigma, len(values_2["Voltage [V]"].data)),
+    }
+)
+
+# Generate a problem for each dataset and combine into one
+problem_1 = pybop.FittingProblem(model_1, parameters, dataset_1)
+problem_2 = pybop.FittingProblem(model_2, parameters, dataset_2)
+problem = pybop.MultiFittingProblem(problem_list=[problem_1, problem_2])
+
+# Generate the cost function and optimisation class
+cost = pybop.SumSquaredError(problem)
+optim = pybop.IRPropMin(
+    cost,
+    # sigma0=0.011,
+    verbose=True,
+    max_iterations=12,
+)
+
+# Run optimisation
+x, final_cost = optim.run()
+print("Estimated parameters:", x)
+
+# Plot the timeseries output
+pybop.quick_plot(problem_1, inputs=x, title="Optimised Comparison")
+pybop.quick_plot(problem_2, inputs=x, title="Optimised Comparison")
+
+# Plot convergence
+pybop.plot_convergence(optim)
+
+# Plot the parameter traces
+pybop.plot_parameters(optim)
+
+# Plot the cost landscape with optimisation path
+bounds = np.array([[0.5, 0.8], [0.4, 0.7]])
+pybop.plot2d(optim, bounds=bounds, steps=15)

pybop/__init__.py

@@ -75,7 +75,7 @@
 # Problem class
 #
 from .problems.base_problem import BaseProblem
-from .problems.fitting_problem import FittingProblem
+from .problems.fitting_problem import FittingProblem, MultiFittingProblem
 from .problems.design_problem import DesignProblem
 
 #

pybop/problems/base_problem.py

@@ -49,7 +49,8 @@ def __init__(
             )
 
         self.parameters = parameters
-        self._model = model
+        if model is not None:
+            self._model = model
         self.check_model = check_model
         if isinstance(signal, str):
             signal = [signal]

pybop/problems/fitting_problem.py

@@ -1,7 +1,9 @@
 import numpy as np
 
 from pybop import BaseProblem
+from pybop._dataset import Dataset
 from pybop.models.base_model import Inputs
+from pybop.parameters.parameter import Parameters
 
 
 class FittingProblem(BaseProblem):
@@ -134,3 +136,116 @@ def evaluateS1(self, inputs: Inputs):
         )
 
         return (y, np.asarray(dy))
+
+
+class MultiFittingProblem(BaseProblem):
+    """
+    Problem class for joining mulitple fitting problems.
+
+    Extends `FittingProblem` to multiple fitting problems.
+    """
+
+    def __init__(self, problem_list, weights=None):
+        self.problem_list = problem_list
+        self.weights = weights
+
+        # Compile the set of parameters, ignoring duplicates
+        combined_parameters = Parameters()
+        for problem in self.problem_list:
+            combined_parameters.join(problem.parameters)
+
+        # Combine the target datasets
+        combined_dataset = Dataset(
+            {"Time [s]": np.asarray([]), "Combined signal": np.asarray([])}
+        )
+        for problem in self.problem_list:
+            for signal in problem.signal:
+                combined_dataset["Time [s]"] = np.concatenate(
+                    (combined_dataset["Time [s]"], problem._time_data)
+                )
+                combined_dataset["Combined signal"] = np.concatenate(
+                    (combined_dataset["Combined signal"], problem._target[signal])
+                )
+
+        super().__init__(
+            parameters=combined_parameters,
+            model=None,
+            signal=["Combined signal"],
+        )
+        self._dataset = combined_dataset.data
+        self.parameters.initial_value()
+
+        # Unpack time and target data
+        self._time_data = self._dataset["Time [s]"]
+        self.n_time_data = len(self._time_data)
+        self.set_target(combined_dataset)
+
+    @property
+    def n_problems(self):
+        return len(self.problem_list)
+
+    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 given inputs.
+        """
+        inputs = self.parameters.verify(inputs)
+        self.parameters.update(values=list(inputs.values()))
+
+        y = {"Combined signal": np.asarray([])}
+        for problem in self.problem_list:
+            problem_inputs = problem.parameters.as_dict()
+            for signal in problem.signal:
+                yi = problem.evaluate(problem_inputs)
+                y["Combined signal"] = np.concatenate(
+                    (y["Combined signal"], yi[signal])
+                )
+
+        return y
+
+    def evaluateS1(self, inputs: Inputs):
+        """
+        Evaluate the model with the given parameters and return the signal and its derivatives.
+
+        Parameters
+        ----------
+        inputs : Inputs
+            Parameters for evaluation of the model.
+
+        Returns
+        -------
+        tuple
+            A tuple containing the simulation result y(t) and the sensitivities dy/dx(t) evaluated
+            with given inputs.
+        """
+
+        inputs = self.parameters.verify(inputs)
+        self.parameters.update(values=list(inputs.values()))
+
+        # y = np.empty((self._target_length))
+        # dy = np.empty((self._target_length, self.n_parameters))
+
+        y = {"Combined signal": np.asarray([])}
+        dy = None
+        for problem in self.problem_list:
+            problem_inputs = problem.parameters.as_dict()
+            for signal in problem.signal:
+                yi, dyi = problem.evaluateS1(problem_inputs)
+                y["Combined signal"] = np.concatenate(
+                    (y["Combined signal"], yi[signal])
+                )
+                if dy is None:
+                    dy = dyi
+                else:
+                    dy = np.concatenate((dy, dyi))
+
+        return (y, dy)

tests/integration/test_model_experiment_changes.py

@@ -107,3 +107,53 @@ def final_cost(self, solution, model, parameters, init_soc):
         optim = pybop.PSO(cost)
         x, final_cost = optim.run()
         return final_cost
+
+    @pytest.mark.integration
+    def test_multi_fitting_problem(self):
+        parameter_set = pybop.ParameterSet.pybamm("Chen2020")
+        parameters = pybop.Parameter(
+            "Negative electrode active material volume fraction",
+            prior=pybop.Gaussian(0.68, 0.05),
+            true_value=parameter_set[
+                "Negative electrode active material volume fraction"
+            ],
+        )
+
+        model_1 = pybop.lithium_ion.SPM(parameter_set=parameter_set)
+        experiment_1 = pybop.Experiment(
+            ["Discharge at 1C until 3 V (4 seconds period)"]
+        )
+        solution_1 = model_1.predict(experiment=experiment_1)
+        dataset_1 = pybop.Dataset(
+            {
+                "Time [s]": solution_1["Time [s]"].data,
+                "Current function [A]": solution_1["Current [A]"].data,
+                "Voltage [V]": solution_1["Voltage [V]"].data,
+            }
+        )
+
+        model_2 = pybop.lithium_ion.SPMe(parameter_set=parameter_set.copy())
+        experiment_2 = pybop.Experiment(
+            ["Discharge at 3C until 3 V (4 seconds period)"]
+        )
+        solution_2 = model_2.predict(experiment=experiment_2)
+        dataset_2 = pybop.Dataset(
+            {
+                "Time [s]": solution_2["Time [s]"].data,
+                "Current function [A]": solution_2["Current [A]"].data,
+                "Voltage [V]": solution_2["Voltage [V]"].data,
+            }
+        )
+
+        # Define a problem for each dataset and combine them into one
+        problem_1 = pybop.FittingProblem(model_1, parameters, dataset_1)
+        problem_2 = pybop.FittingProblem(model_2, parameters, dataset_2)
+        problem = pybop.MultiFittingProblem(problem_list=[problem_1, problem_2])
+        cost = pybop.RootMeanSquaredError(problem)
+
+        # Test with a gradient and non-gradient-based optimiser
+        for optimiser in [pybop.SNES, pybop.IRPropMin]:
+            optim = optimiser(cost)
+            x, final_cost = optim.run()
+            np.testing.assert_allclose(x, parameters.true_value, atol=2e-5)
+            np.testing.assert_allclose(final_cost, 0, atol=2e-5)