Adds ExponentialDecay Model (#583)

* feat: adds ExponentialDecay model, updates problem class for non-battery models, updates tests.

* refactor: update model variables, rename model, add standalone model test

* tests: updates for decay model

examples/scripts/comparison_examples/exp_UKF.py

@@ -1,11 +1,11 @@
 import numpy as np
+import pybamm
 
 import pybop
-from examples.standalone.model import ExponentialDecay
 
 # Parameter set and model definition
-parameter_set = {"k": 0.1, "y0": 1.0}
-model = ExponentialDecay(parameter_set=parameter_set, n_states=1)
+parameter_set = pybamm.ParameterValues({"k": 0.1, "y0": 1.0})
+model = pybop.ExponentialDecayModel(parameter_set=parameter_set, n_states=1)
 
 # Fitting parameters
 parameters = pybop.Parameters(
@@ -23,50 +23,30 @@
     ),
 )
 
+
+def noise(sigma):
+    return np.random.normal(0, sigma, len(values["y_0"].data))
+
+
 # Make a prediction with measurement noise
 sigma = 1e-2
 t_eval = np.linspace(0, 20, 10)
-true_inputs = parameters.as_dict("true")
 values = model.predict(t_eval=t_eval)
-values = values["2y"].data
-corrupt_values = values + np.random.normal(0, sigma, len(t_eval))
-
-# Verification step: compute the analytical solution for 2y
-expected_values = (
-    2 * parameters["y0"].true_value * np.exp(-parameters["k"].true_value * t_eval)
-)
-
-# Verification step: make another prediction using the Observer class
-model.build(parameters=parameters)
-simulator = pybop.Observer(parameters, model, signal=["2y"])
-simulator.domain_data = t_eval
-measurements = simulator.evaluate(true_inputs)
-
-# Verification step: Compare by plot
-go = pybop.plot.PlotlyManager().go
-line1 = go.Scatter(x=t_eval, y=corrupt_values, name="Corrupt values", mode="markers")
-line2 = go.Scatter(
-    x=t_eval, y=expected_values, name="Expected trajectory", mode="lines"
-)
-line3 = go.Scatter(
-    x=t_eval, y=measurements["2y"], name="Observed values", mode="markers"
-)
-fig = go.Figure(data=[line1, line2, line3])
 
 # Form dataset
 dataset = pybop.Dataset(
     {
         "Time [s]": t_eval,
         "Current function [A]": 0 * t_eval,  # placeholder
-        "2y": corrupt_values,
+        "y_0": values["y_0"].data + noise(sigma),
     }
 )
 
 # Build the model to get the number of states
 model.build(dataset=dataset.data, parameters=parameters)
 
 # Define the UKF observer
-signal = ["2y"]
+signal = ["y_0"]
 n_states = model.n_states
 n_signals = len(signal)
 covariance = np.diag([sigma**2] * n_states)
@@ -82,15 +62,6 @@
     signal=signal,
 )
 
-# Verification step: Find the maximum likelihood estimate given the true parameters
-estimation = observer.evaluate(true_inputs)
-
-# Verification step: Add the estimate to the plot
-line4 = go.Scatter(
-    x=t_eval, y=estimation["2y"], name="Estimated trajectory", mode="lines"
-)
-fig.add_trace(line4)
-fig.show()
 
 # Generate problem, cost function, and optimisation class
 cost = pybop.ObserverCost(observer)

examples/scripts/comparison_examples/exponential_decay.py

@@ -0,0 +1,71 @@
+import numpy as np
+import pybamm
+
+import pybop
+
+# Define model and use high-performant solver for sensitivities
+solver = pybamm.IDAKLUSolver
+parameter_set = pybamm.ParameterValues({"k": 1, "y0": 0.5})
+model = pybop.ExponentialDecayModel(
+    parameter_set=parameter_set, solver=solver, n_states=2
+)
+
+# Fitting parameters
+parameters = pybop.Parameters(
+    pybop.Parameter(
+        "k",
+        prior=pybop.Gaussian(0.5, 0.05),
+    ),
+    pybop.Parameter(
+        "y0",
+        prior=pybop.Gaussian(0.2, 0.05),
+    ),
+)
+
+# Generate data
+sigma = 0.003
+t_eval = np.linspace(0, 10, 100)
+values = model.predict(t_eval=t_eval)
+
+
+def noise(sigma):
+    return np.random.normal(0, sigma, len(values["y_0"].data))
+
+
+# Form dataset
+dataset = pybop.Dataset(
+    {
+        "Time [s]": t_eval,
+        "Current function [A]": 0 * t_eval,
+        "y_0": values["y_0"].data + noise(sigma),
+        "y_1": values["y_1"].data + noise(sigma),
+    }
+)
+
+signal = ["y_0", "y_1"]
+# Generate problem, cost function, and optimisation class
+problem = pybop.FittingProblem(model, parameters, dataset, signal=signal)
+cost = pybop.Minkowski(problem, p=2)
+optim = pybop.AdamW(
+    cost,
+    verbose=True,
+    allow_infeasible_solutions=True,
+    sigma0=0.02,
+    max_iterations=100,
+    max_unchanged_iterations=20,
+)
+
+# Run optimisation
+results = optim.run()
+
+# Plot the timeseries output
+pybop.plot.quick(problem, problem_inputs=results.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
+pybop.plot.surface(optim)

pybop/__init__.py

@@ -80,6 +80,7 @@
 from .models.base_model import BaseModel
 from .models import lithium_ion
 from .models import empirical
+from .models._exponential_decay import ExponentialDecayModel
 from .models.base_model import TimeSeriesState
 from .models.base_model import Inputs
 

pybop/models/_exponential_decay.py

@@ -0,0 +1,82 @@
+import pybamm
+
+from pybop.models.base_model import BaseModel
+
+
+class ExponentialDecayModel(BaseModel):
+    """
+    Exponential decay model defined by the equation:
+
+        dy/dt = -k * y,  y(0) = y0
+
+    Note: The output variables are named "y_{i}" for each state.
+        For example, the first state is "y_0", the second is "y_1", etc.
+    Attributes:
+        n_states (int): Number of states in the system (default is 1).
+        pybamm_model (pybamm.BaseModel): PyBaMM model representation.
+        default_parameter_values (pybamm.ParameterValues): Default parameter values
+            for the model, with "k" (decay rate) and "y0" (initial condition).
+
+    Parameters:
+        name (str): Name of the model (default: "Experimental Decay Model").
+        parameter_set (pybamm.ParameterValues): Parameter values for the model.
+        n_states (int): Number of states in the system. Must be >= 1.
+    """
+
+    def __init__(
+        self,
+        name: str = "Experimental Decay Model",
+        parameter_set: pybamm.ParameterValues = None,
+        n_states: int = 1,
+        solver=None,
+    ):
+        if n_states < 1:
+            raise ValueError("The number of states (n_states) must be at least 1.")
+
+        super().__init__(name=name, parameter_set=parameter_set)
+
+        self.n_states = n_states
+        if solver is None:
+            self._solver = pybamm.CasadiSolver
+            self._solver.mode = "fast with events"
+            self._solver.max_step_decrease_count = 1
+        else:
+            self._solver = solver
+
+        # Initialise the PyBaMM model, variables, parameters
+        self.pybamm_model = pybamm.BaseModel()
+        ys = [pybamm.Variable(f"y_{i}") for i in range(n_states)]
+        k = pybamm.Parameter("k")
+        y0 = pybamm.Parameter("y0")
+
+        # Set up the right-hand side and initial conditions
+        self.pybamm_model.rhs = {y: -k * y for y in ys}
+        self.pybamm_model.initial_conditions = {y: y0 for y in ys}
+
+        # Define model outputs and set default values
+        self.pybamm_model.variables = {f"y_{en}": i for en, i in enumerate(ys)} | {
+            "Time [s]": pybamm.t
+        }
+        self.default_parameter_values = parameter_set or pybamm.ParameterValues(
+            {"k": 0.1, "y0": 1}
+        )
+
+        # Store model attributes to be used by the solver
+        self._unprocessed_model = self.pybamm_model
+        self._parameter_set = self.default_parameter_values
+        self._unprocessed_parameter_set = self._parameter_set
+
+        # Geometry and solver setup
+        self._geometry = self.pybamm_model.default_geometry
+        self._submesh_types = self.pybamm_model.default_submesh_types
+        self._var_pts = self.pybamm_model.default_var_pts
+        self._spatial_methods = self.pybamm_model.default_spatial_methods
+        self._solver = pybamm.CasadiSolver(mode="fast")
+
+        # Additional attributes for solver and discretisation
+        self._model_with_set_params = None
+        self._built_model = None
+        self._built_initial_soc = None
+        self._mesh = None
+        self._disc = None
+        self.geometric_parameters = {}

pybop/problems/base_problem.py

@@ -79,8 +79,6 @@ def __init__(
 
         # Add domain and remove duplicates
         self.additional_variables = additional_variables or []
-        self.additional_variables.extend([self.domain, "Current [A]"])
-        self.additional_variables = list(set(self.additional_variables))
 
         # If model.solver is IDAKLU, set output vars for improved performance
         self.output_vars = tuple(self.signal + self.additional_variables)

pybop/problems/design_problem.py

@@ -56,6 +56,10 @@ def __init__(
             "Non-physical point encountered",
         ]
 
+        # Add "Current" to the variable list
+        self.additional_variables.extend([self.domain, "Current [A]"])
+        self.additional_variables = list(set(self.additional_variables))
+
         # Set whether to update the nominal capacity along with the design parameters
         if update_capacity is True:
             nominal_capacity_warning = (

pybop/problems/fitting_problem.py

@@ -59,7 +59,7 @@ def __init__(
         self._n_parameters = len(self.parameters)
 
         # Check that the dataset contains necessary variables
-        dataset.check(domain=self.domain, signal=[*self.signal, "Current function [A]"])
+        dataset.check(domain=self.domain, signal=self.signal)
 
         # Unpack domain and target data
         self._domain_data = self._dataset[self.domain]

tests/integration/test_observer_parameterisation.py

@@ -1,8 +1,8 @@
 import numpy as np
+import pybamm
 import pytest
 
 import pybop
-from examples.standalone.model import ExponentialDecay
 
 
 class TestObservers:
@@ -20,11 +20,15 @@ def setup(self):
 
     @pytest.fixture
     def parameter_set(self):
-        return {"k": self.ground_truth[0], "y0": self.ground_truth[1]}
+        return pybamm.ParameterValues(
+            {"k": self.ground_truth[0], "y0": self.ground_truth[1]}
+        )
 
     @pytest.fixture
     def model(self, parameter_set):
-        return ExponentialDecay(parameter_set=parameter_set, n_states=1)
+        return pybop.ExponentialDecayModel(
+            parameter_set=parameter_set, solver=pybamm.IDAKLUSolver, n_states=1
+        )
 
     @pytest.fixture
     def parameters(self, parameter_set):
@@ -51,20 +55,20 @@ def test_observer_exponential_decay(self, parameters, model):
         # Make a prediction with measurement noise
         sigma = 1e-2
         t_eval = np.linspace(0, 20, 10)
-        values = model.predict(t_eval=t_eval)["2y"].data
+        values = model.predict(t_eval=t_eval)["y_0"].data
         corrupt_values = values + self.noise(sigma, len(t_eval))
 
         # Form dataset
         dataset = pybop.Dataset(
             {
                 "Time [s]": t_eval,
                 "Current function [A]": 0 * t_eval,  # placeholder
-                "2y": corrupt_values,
+                "y_0": corrupt_values,
             }
         )
 
         # Define the UKF observer
-        signal = ["2y"]
+        signal = ["y_0"]
         n_states = model.n_states
         n_signals = len(signal)
         covariance = np.diag([sigma**2] * n_states)