Merge pull request pybop-team#405 from pybop-team/269-add-eis-paramet…

…er-identification-methods

Adds EIS prediction methods

.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

CHANGELOG.md

@@ -2,6 +2,7 @@
 
 ## Features
 
+- [#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.

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/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

pybop/_dataset.py

@@ -1,3 +1,5 @@
+from typing import Union
+
 import numpy as np
 from pybamm import solvers
 
@@ -76,41 +78,68 @@ def __getitem__(self, key):
 
         return self.data[key]
 
-    def check(self, signal=None):
+    def check(self, domain: str = None, signal: Union[str, list[str]] = None) -> bool:
         """
         Check the consistency of a PyBOP Dataset against the expected format.
 
+        Parameters
+        ----------
+        domain : str, optional
+            The domain of the dataset. Defaults to "Time [s]".
+        signal : str or List[str], optional
+            The signal(s) to check. Defaults to ["Voltage [V]"].
+
         Returns
         -------
         bool
-            If True, the dataset has the expected attributes.
+            True if the dataset has the expected attributes.
 
         Raises
         ------
         ValueError
             If the time series and the data series are not consistent.
         """
-        if signal is None:
-            signal = ["Voltage [V]"]
-        if isinstance(signal, str):
-            signal = [signal]
-
-        # Check that the dataset contains time and chosen signal
-        for name in ["Time [s]", *signal]:
-            if name not in self.names:
-                raise ValueError(f"expected {name} in list of dataset")
-
-        # Check for increasing times
-        time_data = self.data["Time [s]"]
+        self.domain = domain or "Time [s]"
+        signals = [signal] if isinstance(signal, str) else (signal or ["Voltage [V]"])
+
+        # Check that the dataset contains domain and chosen signals
+        missing_attributes = set([self.domain, *signals]) - set(self.names)
+        if missing_attributes:
+            raise ValueError(
+                f"Expected {', '.join(missing_attributes)} in list of dataset"
+            )
+
+        domain_data = self.data[self.domain]
+
+        # Check domain-specific constraints
+        if self.domain == "Time [s]":
+            self._check_time_constraints(domain_data)
+        elif self.domain == "Frequency [Hz]":
+            self._check_frequency_constraints(domain_data)
+
+        # Check for consistent data length
+        self._check_data_consistency(domain_data, signals)
+
+        return True
+
+    @staticmethod
+    def _check_time_constraints(time_data: np.ndarray) -> None:
         if np.any(time_data < 0):
-            raise ValueError("Times can not be negative.")
+            raise ValueError("Times cannot be negative.")
         if np.any(time_data[:-1] >= time_data[1:]):
             raise ValueError("Times must be increasing.")
 
-        # Check for consistent data
-        n_time_data = len(time_data)
-        for s in signal:
-            if len(self.data[s]) != n_time_data:
-                raise ValueError(f"Time data and {s} data must be the same length.")
-
-        return True
+    @staticmethod
+    def _check_frequency_constraints(freq_data: np.ndarray) -> None:
+        if np.any(freq_data < 0):
+            raise ValueError("Frequencies cannot be negative.")
+
+    def _check_data_consistency(
+        self, domain_data: np.ndarray, signals: list[str]
+    ) -> None:
+        n_domain_data = len(domain_data)
+        for s in signals:
+            if len(self.data[s]) != n_domain_data:
+                raise ValueError(
+                    f"{self.domain} data and {s} data must be the same length."
+                )