Merge branch 'develop' into 441-functional-parameters

.github/workflows/lychee_links.yaml

@@ -44,6 +44,7 @@ jobs:
             --exclude "https://a.tile.openstreetmap.org/*"
             --exclude "https://openstreetmap.org/*|https://www.openstreetmap.org/*"
             --exclude "https://cdn.plot.ly/*"
+            --exclude "http://www.w3.org/*|https://www.w3.org/*"
             --exclude "https://doi.org/*"
             --exclude-path ./CHANGELOG.md
             --exclude-path asv.conf.json

.pre-commit-config.yaml

@@ -4,7 +4,7 @@ ci:
 
 repos:
   - repo: https://github.com/astral-sh/ruff-pre-commit
-    rev: "v0.5.7"
+    rev: "v0.6.1"
     hooks:
       - id: ruff
         args: [--fix, --show-fixes]

CHANGELOG.md

@@ -3,6 +3,7 @@
 ## Features
 
 - [#441](https://github.com/pybop-team/PyBOP/issues/441) - Adds an example for estimating constants within a `pybamm.FunctionalParameter`.
+- [#405](https://github.com/pybop-team/PyBOP/pull/405) - Adds frequency-domain based EIS prediction methods via `model.simulateEIS` and updates to `problem.evaluate` with examples and tests.
 - [#460](https://github.com/pybop-team/PyBOP/pull/460) - Notebook example files added for ECM and folder structure updated.
 - [#450](https://github.com/pybop-team/PyBOP/pull/450) - Adds support for IDAKLU with output variables, and corresponding examples, tests.
 - [#364](https://github.com/pybop-team/PyBOP/pull/364) - Adds the MultiFittingProblem class and the multi_fitting example script.
@@ -23,6 +24,7 @@
 
 ## Breaking Changes
 
+- [#436](https://github.com/pybop-team/PyBOP/pull/436) - **API Change:** The functionality from `BaseCost.evaluate/S1` & `BaseCost._evaluate/S1` is represented in `BaseCost.__call__` & `BaseCost.compute`. `BaseCost.compute` directly acts on the predictions, while `BaseCost.__call__` calls `BaseProblem.evaluate/S1` before `BaseCost.compute`. `compute` has optional args for gradient cost calculations.
 - [#424](https://github.com/pybop-team/PyBOP/issues/424) - Replaces the `init_soc` input to `FittingProblem` with the option to pass an initial OCV value, updates `BaseModel` and fixes `multi_model_identification.ipynb` and `spm_electrode_design.ipynb`.
 
 # [v24.6.1](https://github.com/pybop-team/PyBOP/tree/v24.6.1) - 2024-07-31

examples/scripts/eis_fitting.py

@@ -0,0 +1,82 @@
+import numpy as np
+
+import pybop
+
+# Define model
+parameter_set = pybop.ParameterSet.pybamm("Chen2020")
+parameter_set["Contact resistance [Ohm]"] = 0.0
+initial_state = {"Initial SoC": 0.5}
+n_frequency = 20
+sigma0 = 1e-4
+f_eval = np.logspace(-4, 5, n_frequency)
+model = pybop.lithium_ion.SPM(
+    parameter_set=parameter_set,
+    eis=True,
+    options={"surface form": "differential", "contact resistance": "true"},
+)
+
+# Create synthetic data for parameter inference
+sim = model.simulateEIS(
+    inputs={
+        "Negative electrode active material volume fraction": 0.531,
+        "Positive electrode active material volume fraction": 0.732,
+    },
+    f_eval=f_eval,
+    initial_state=initial_state,
+)
+
+# Fitting parameters
+parameters = pybop.Parameters(
+    pybop.Parameter(
+        "Negative electrode active material volume fraction",
+        prior=pybop.Uniform(0.4, 0.75),
+        bounds=[0.375, 0.75],
+    ),
+    pybop.Parameter(
+        "Positive electrode active material volume fraction",
+        prior=pybop.Uniform(0.4, 0.75),
+        bounds=[0.375, 0.75],
+    ),
+)
+
+
+def noise(sigma, values):
+    # Generate real part noise
+    real_noise = np.random.normal(0, sigma, values)
+
+    # Generate imaginary part noise
+    imag_noise = np.random.normal(0, sigma, values)
+
+    # Combine them into a complex noise
+    return real_noise + 1j * imag_noise
+
+
+# Form dataset
+dataset = pybop.Dataset(
+    {
+        "Frequency [Hz]": f_eval,
+        "Current function [A]": np.ones(n_frequency) * 0.0,
+        "Impedance": sim["Impedance"] + noise(sigma0, len(sim["Impedance"])),
+    }
+)
+
+signal = ["Impedance"]
+# Generate problem, cost function, and optimisation class
+problem = pybop.FittingProblem(model, parameters, dataset, signal=signal)
+cost = pybop.GaussianLogLikelihoodKnownSigma(problem, sigma0=sigma0)
+optim = pybop.CMAES(cost, max_iterations=100, sigma0=0.25, max_unchanged_iterations=30)
+
+x, final_cost = optim.run()
+print("Estimated parameters:", x)
+
+# Plot the nyquist
+pybop.nyquist(problem, problem_inputs=x, title="Optimised Comparison")
+
+# Plot convergence
+pybop.plot_convergence(optim)
+
+# Plot the parameter traces
+pybop.plot_parameters(optim)
+
+# Plot 2d landscape
+pybop.plot2d(optim, steps=10)

examples/scripts/exp_UKF.py

@@ -39,7 +39,7 @@
 # Verification step: make another prediction using the Observer class
 model.build(parameters=parameters)
 simulator = pybop.Observer(parameters, model, signal=["2y"])
-simulator.time_data = t_eval
+simulator.domain_data = t_eval
 measurements = simulator.evaluate(true_inputs)
 
 # Verification step: Compare by plotting

examples/scripts/spm_IRPropMin.py

@@ -35,7 +35,7 @@
 
 # Generate problem, cost function, and optimisation class
 problem = pybop.FittingProblem(model, parameters, dataset)
-cost = pybop.SumSquaredError(problem)
+cost = pybop.Minkowski(problem, p=2)
 optim = pybop.IRPropMin(cost, max_iterations=100)
 
 x, final_cost = optim.run()

examples/scripts/spm_MLE.py

@@ -2,8 +2,14 @@
 
 import pybop
 
-# Define model
+# Define model and set initial parameter values
 parameter_set = pybop.ParameterSet.pybamm("Chen2020")
+parameter_set.update(
+    {
+        "Negative electrode active material volume fraction": 0.63,
+        "Positive electrode active material volume fraction": 0.51,
+    }
+)
 model = pybop.lithium_ion.SPM(parameter_set=parameter_set)
 
 # Fitting parameters
@@ -19,31 +25,39 @@
     ),
 )
 
-# Set initial parameter values
-parameter_set.update(
-    {
-        "Negative electrode active material volume fraction": 0.63,
-        "Positive electrode active material volume fraction": 0.51,
-    }
-)
 # Generate data
-sigma = 0.005
-t_eval = np.arange(0, 900, 3)
-values = model.predict(t_eval=t_eval)
-corrupt_values = values["Voltage [V]"].data + np.random.normal(0, sigma, len(t_eval))
+sigma = 0.002
+experiment = pybop.Experiment(
+    [
+        (
+            "Discharge at 0.5C for 3 minutes (3 second period)",
+            "Charge at 0.5C for 3 minutes (3 second period)",
+        ),
+    ]
+)
+values = model.predict(initial_state={"Initial SoC": 0.5}, experiment=experiment)
+
+
+def noise(sigma):
+    return np.random.normal(0, sigma, len(values["Voltage [V]"].data))
+
 
 # Form dataset
 dataset = pybop.Dataset(
     {
-        "Time [s]": t_eval,
+        "Time [s]": values["Time [s]"].data,
         "Current function [A]": values["Current [A]"].data,
-        "Voltage [V]": corrupt_values,
+        "Voltage [V]": values["Voltage [V]"].data + noise(sigma),
+        "Bulk open-circuit voltage [V]": values["Bulk open-circuit voltage [V]"].data
+        + noise(sigma),
     }
 )
 
+
+signal = ["Voltage [V]", "Bulk open-circuit voltage [V]"]
 # Generate problem, cost function, and optimisation class
-problem = pybop.FittingProblem(model, parameters, dataset)
-likelihood = pybop.GaussianLogLikelihoodKnownSigma(problem, sigma0=sigma)
+problem = pybop.FittingProblem(model, parameters, dataset, signal=signal)
+likelihood = pybop.GaussianLogLikelihood(problem, sigma0=sigma * 4)
 optim = pybop.IRPropMin(
     likelihood,
     max_unchanged_iterations=20,

examples/standalone/cost.py

@@ -1,3 +1,5 @@
+import numpy as np
+
 import pybop
 
 
@@ -43,22 +45,18 @@ def __init__(self, problem=None):
         )
         self.x0 = self.parameters.initial_value()
 
-    def compute(self, inputs):
+    def compute(
+        self, y: dict = None, dy: np.ndarray = None, calculate_grad: bool = False
+    ):
         """
         Compute the cost for a given parameter value.
 
         The cost function is defined as cost(x) = x^2 + 42, where x is the
         parameter value.
 
-        Parameters
-        ----------
-        inputs : Dict
-            The parameters for which to evaluate the cost.
-
         Returns
         -------
         float
             The calculated cost value for the given parameter.
         """
-
-        return inputs["x"] ** 2 + 42
+        return self.parameters["x"].value ** 2 + 42

examples/standalone/problem.py

@@ -16,7 +16,7 @@ def __init__(
         check_model=True,
         signal=None,
         additional_variables=None,
-        init_soc=None,
+        initial_state=None,
     ):
         super().__init__(parameters, model, check_model, signal, additional_variables)
         self._dataset = dataset.data
@@ -26,15 +26,15 @@ def __init__(
             if name not in self._dataset:
                 raise ValueError(f"expected {name} in list of dataset")
 
-        self._time_data = self._dataset["Time [s]"]
-        self.n_time_data = len(self._time_data)
-        if np.any(self._time_data < 0):
+        self._domain_data = self._dataset[self.domain]
+        self.n_data = len(self._domain_data)
+        if np.any(self._domain_data < 0):
             raise ValueError("Times can not be negative.")
-        if np.any(self._time_data[:-1] >= self._time_data[1:]):
+        if np.any(self._domain_data[:-1] >= self._domain_data[1:]):
             raise ValueError("Times must be increasing.")
 
         for signal in self.signal:
-            if len(self._dataset[signal]) != self.n_time_data:
+            if len(self._dataset[signal]) != self.n_data:
                 raise ValueError(
                     f"Time data and {signal} data must be the same length."
                 )
@@ -56,7 +56,7 @@ def evaluate(self, inputs, **kwargs):
         """
 
         return {
-            signal: inputs["Gradient"] * self._time_data + inputs["Intercept"]
+            signal: inputs["Gradient"] * self._domain_data + inputs["Intercept"]
             for signal in self.signal
         }
 
@@ -78,7 +78,7 @@ def evaluateS1(self, inputs):
 
         y = self.evaluate(inputs)
 
-        dy = np.zeros((self.n_time_data, self.n_outputs, self.n_parameters))
-        dy[:, 0, 0] = self._time_data
+        dy = np.zeros((self.n_data, self.n_outputs, self.n_parameters))
+        dy[:, 0, 0] = self._domain_data
 
         return (y, dy)

pybop/__init__.py

@@ -43,7 +43,7 @@
 #
 # Utilities
 #
-from ._utils import is_numeric
+from ._utils import is_numeric, SymbolReplacer
 
 #
 # Experiment class
@@ -157,6 +157,7 @@
 from .plotting.plot_convergence import plot_convergence
 from .plotting.plot_parameters import plot_parameters
 from .plotting.plot_problem import quick_plot
+from .plotting.nyquist import nyquist
 
 #
 # Remove any imported modules, so we don't expose them as part of pybop