diff --git a/.github/workflows/lychee_links.yaml b/.github/workflows/lychee_links.yaml
index fa82a4c78..1e1ed3df3 100644
--- a/.github/workflows/lychee_links.yaml
+++ b/.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
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 3f4d462e5..e09ef0e94 100644
--- a/CHANGELOG.md
+++ b/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.
diff --git a/examples/scripts/eis_fitting.py b/examples/scripts/eis_fitting.py
new file mode 100644
index 000000000..f86e7707a
--- /dev/null
+++ b/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)
diff --git a/examples/scripts/exp_UKF.py b/examples/scripts/exp_UKF.py
index f305443de..3845a9679 100644
--- a/examples/scripts/exp_UKF.py
+++ b/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
diff --git a/examples/scripts/spm_IRPropMin.py b/examples/scripts/spm_IRPropMin.py
index 1969f6f9d..2ef494c1c 100644
--- a/examples/scripts/spm_IRPropMin.py
+++ b/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()
diff --git a/examples/standalone/problem.py b/examples/standalone/problem.py
index 87fa99395..52e3a4643 100644
--- a/examples/standalone/problem.py
+++ b/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)
diff --git a/pybop/__init__.py b/pybop/__init__.py
index 49c7eb83e..6f8b1b828 100644
--- a/pybop/__init__.py
+++ b/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
diff --git a/pybop/_dataset.py b/pybop/_dataset.py
index 66bcb1f10..120fcb61d 100644
--- a/pybop/_dataset.py
+++ b/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."
+ )
diff --git a/pybop/_utils.py b/pybop/_utils.py
index 6fbfeaab5..426b0d277 100644
--- a/pybop/_utils.py
+++ b/pybop/_utils.py
@@ -1,4 +1,7 @@
+from typing import Optional
+
import numpy as np
+import pybamm
def is_numeric(x):
@@ -6,3 +9,168 @@ def is_numeric(x):
Check if a variable is numeric.
"""
return isinstance(x, (int, float, np.number))
+
+
+class SymbolReplacer:
+ """
+ Helper class to replace all instances of one or more symbols in an expression tree
+ with another symbol, as defined by the dictionary `symbol_replacement_map`
+ Originally developed by pybamm: https://github.com/pybamm-team/pybamm
+
+ Parameters
+ ----------
+ symbol_replacement_map : dict {:class:`pybamm.Symbol` -> :class:`pybamm.Symbol`}
+ Map of which symbols should be replaced by which.
+ processed_symbols: dict {:class:`pybamm.Symbol` -> :class:`pybamm.Symbol`}, optional
+ cached replaced symbols
+ process_initial_conditions: bool, optional
+ Whether to process initial conditions, default is True
+ """
+
+ def __init__(
+ self,
+ symbol_replacement_map: dict[pybamm.Symbol, pybamm.Symbol],
+ processed_symbols: Optional[dict[pybamm.Symbol, pybamm.Symbol]] = None,
+ process_initial_conditions: bool = True,
+ ):
+ self._symbol_replacement_map = symbol_replacement_map
+ self._processed_symbols = processed_symbols or {}
+ self._process_initial_conditions = process_initial_conditions
+
+ def process_model(self, unprocessed_model, inplace=True):
+ """
+ Replace all instances of a symbol in a PyBaMM model class.
+
+ Parameters
+ ----------
+ unprocessed_model : :class:`pybamm.BaseModel`
+ Model class to assign parameter values to
+ inplace: bool, optional
+ If True, replace the parameters in the model in place. Otherwise, return a
+ new model with parameter values set. Default is True.
+ """
+
+ model = unprocessed_model if inplace else unprocessed_model.new_copy()
+
+ for variable, equation in unprocessed_model.rhs.items():
+ pybamm.logger.verbose(f"Replacing symbols in {variable!r} (rhs)")
+ model.rhs[self.process_symbol(variable)] = self.process_symbol(equation)
+
+ for variable, equation in unprocessed_model.algebraic.items():
+ pybamm.logger.verbose(f"Replacing symbols in {variable!r} (algebraic)")
+ model.algebraic[self.process_symbol(variable)] = self.process_symbol(
+ equation
+ )
+
+ for variable, equation in unprocessed_model.initial_conditions.items():
+ pybamm.logger.verbose(
+ f"Replacing symbols in {variable!r} (initial conditions)"
+ )
+ if self._process_initial_conditions:
+ model.initial_conditions[self.process_symbol(variable)] = (
+ self.process_symbol(equation)
+ )
+ else:
+ model.initial_conditions[self.process_symbol(variable)] = equation
+
+ model.boundary_conditions = self.process_boundary_conditions(unprocessed_model)
+
+ for variable, equation in unprocessed_model.variables.items():
+ pybamm.logger.verbose(f"Replacing symbols in {variable!r} (variables)")
+ model.variables[variable] = self.process_symbol(equation)
+
+ model.events = self._process_events(unprocessed_model.events)
+ pybamm.logger.info(f"Finish replacing symbols in {model.name}")
+
+ return model
+
+ def _process_events(self, events: list) -> list:
+ new_events = []
+ for event in events:
+ pybamm.logger.verbose(f"Replacing symbols in event '{event.name}'")
+ new_events.append(
+ pybamm.Event(
+ event.name, self.process_symbol(event.expression), event.event_type
+ )
+ )
+ return new_events
+
+ def process_boundary_conditions(self, model):
+ """
+ Process boundary conditions for a PybaMM model class
+ Boundary conditions are dictionaries {"left": left bc, "right": right bc}
+ in general, but may be imposed on the tabs (or *not* on the tab) for a
+ small number of variables, e.g. {"negative tab": neg. tab bc,
+ "positive tab": pos. tab bc "no tab": no tab bc}.
+ """
+ boundary_conditions = {}
+ sides = ["left", "right", "negative tab", "positive tab", "no tab"]
+ for variable, bcs in model.boundary_conditions.items():
+ processed_variable = self.process_symbol(variable)
+ boundary_conditions[processed_variable] = {}
+
+ for side in sides:
+ try:
+ bc, typ = bcs[side]
+ pybamm.logger.verbose(
+ f"Replacing symbols in {variable!r} ({side} bc)"
+ )
+ processed_bc = (self.process_symbol(bc), typ)
+ boundary_conditions[processed_variable][side] = processed_bc
+ except KeyError as err:
+ # Don't raise if side is not in the boundary conditions
+ if err.args[0] in side:
+ pass
+ # Raise otherwise
+ else: # pragma: no cover
+ raise KeyError(err) from err
+
+ return boundary_conditions
+
+ def process_symbol(self, symbol):
+ """
+ This function recurses down the tree, replacing any symbols in
+ self._symbol_replacement_map.keys() with their corresponding value
+
+ Parameters
+ ----------
+ symbol : :class:`pybamm.Symbol`
+ The symbol to replace
+
+ Returns
+ -------
+ :class:`pybamm.Symbol`
+ Symbol with all replacements performed
+ """
+ if symbol in self._processed_symbols:
+ return self._processed_symbols[symbol]
+
+ processed_symbol = self._process_symbol(symbol)
+ self._processed_symbols[symbol] = processed_symbol
+ return processed_symbol
+
+ def _process_symbol(self, symbol: pybamm.Symbol) -> pybamm.Symbol:
+ if symbol in self._symbol_replacement_map:
+ return self._symbol_replacement_map[symbol]
+
+ if isinstance(symbol, pybamm.BinaryOperator):
+ # process children
+ new_left = self.process_symbol(symbol.left)
+ new_right = self.process_symbol(symbol.right)
+ return symbol._binary_new_copy(new_left, new_right) # noqa: SLF001
+
+ if isinstance(symbol, pybamm.UnaryOperator):
+ new_child = self.process_symbol(symbol.child)
+ return symbol._unary_new_copy(new_child) # noqa: SLF001
+
+ if isinstance(symbol, pybamm.Function):
+ new_children = [self.process_symbol(child) for child in symbol.children]
+ # Return a new copy with the replaced symbols
+ return symbol._function_new_copy(new_children) # noqa: SLF001
+
+ if isinstance(symbol, pybamm.Concatenation):
+ new_children = [self.process_symbol(child) for child in symbol.children]
+ return symbol._concatenation_new_copy(new_children) # noqa: SLF001
+
+ # Return leaf
+ return symbol
diff --git a/pybop/costs/_likelihoods.py b/pybop/costs/_likelihoods.py
index fa3eac27c..d651d89eb 100644
--- a/pybop/costs/_likelihoods.py
+++ b/pybop/costs/_likelihoods.py
@@ -15,7 +15,7 @@ class BaseLikelihood(BaseCost):
def __init__(self, problem: BaseProblem):
super().__init__(problem)
- self.n_time_data = problem.n_time_data
+ self.n_data = problem.n_data
class GaussianLogLikelihoodKnownSigma(BaseLikelihood):
@@ -36,7 +36,7 @@ def __init__(self, problem: BaseProblem, sigma0: Union[list[float], float]):
super().__init__(problem)
sigma0 = self.check_sigma0(sigma0)
self.sigma2 = sigma0**2.0
- self._offset = -0.5 * self.n_time_data * np.log(2 * np.pi * self.sigma2)
+ self._offset = -0.5 * self.n_data * np.log(2 * np.pi * self.sigma2)
self._multip = -1 / (2.0 * self.sigma2)
def compute(
@@ -59,7 +59,7 @@ def compute(
# Calculate residuals and error
r = np.asarray([self._target[signal] - y[signal] for signal in self.signal])
- e = np.sum(self._offset + self._multip * np.sum(r**2.0))
+ e = np.sum(self._offset + self._multip * np.sum(np.real(r * np.conj(r))))
if calculate_grad:
dl = np.sum((np.sum((r * dy.T), axis=2) / self.sigma2), axis=1)
@@ -107,7 +107,7 @@ def __init__(
):
super().__init__(problem)
self._dsigma_scale = dsigma_scale
- self._logpi = -0.5 * self.n_time_data * np.log(2 * np.pi)
+ self._logpi = -0.5 * self.n_data * np.log(2 * np.pi)
self.sigma = Parameters()
self._add_sigma_parameters(sigma0)
@@ -187,14 +187,14 @@ def compute(
r = np.asarray([self._target[signal] - y[signal] for signal in self.signal])
e = np.sum(
self._logpi
- - self.n_time_data * np.log(sigma)
- - np.sum(r**2.0, axis=1) / (2.0 * sigma**2.0)
+ - self.n_data * np.log(sigma)
+ - np.sum(np.real(r * np.conj(r)), axis=1) / (2.0 * sigma**2.0)
)
if calculate_grad:
dl = np.sum((np.sum((r * dy.T), axis=2) / (sigma**2.0)), axis=1)
dsigma = (
- -self.n_time_data / sigma + np.sum(r**2.0, axis=1) / (sigma**3.0)
+ -self.n_data / sigma + np.sum(r**2.0, axis=1) / (sigma**3.0)
) / self._dsigma_scale
dl = np.concatenate((dl.flatten(), dsigma))
return e, dl
diff --git a/pybop/costs/design_costs.py b/pybop/costs/design_costs.py
index e0dca790a..319202bec 100644
--- a/pybop/costs/design_costs.py
+++ b/pybop/costs/design_costs.py
@@ -61,7 +61,7 @@ def update_simulation_data(self, inputs: Inputs):
if "Time [s]" not in solution:
raise ValueError("The solution does not contain time data.")
- self.problem.time_data = solution["Time [s]"]
+ self.problem.domain_data = solution["Time [s]"]
self.problem.target = {key: solution[key] for key in self.problem.signal}
self.dt = solution["Time [s]"][1] - solution["Time [s]"][0]
diff --git a/pybop/costs/fitting_costs.py b/pybop/costs/fitting_costs.py
index 6204659ff..da40dd8c7 100644
--- a/pybop/costs/fitting_costs.py
+++ b/pybop/costs/fitting_costs.py
@@ -54,7 +54,7 @@ def compute(
# Calculate residuals and error
r = np.asarray([y[signal] - self._target[signal] for signal in self.signal])
- e = np.sqrt(np.mean(r**2, axis=1))
+ e = np.sqrt(np.mean(np.abs(r) ** 2, axis=1))
if calculate_grad:
de = np.mean((r * dy.T), axis=2) / (e + np.finfo(float).eps)
@@ -120,7 +120,7 @@ def compute(
# Calculate residuals and error
r = np.asarray([y[signal] - self._target[signal] for signal in self.signal])
- e = np.sum(np.sum(r**2, axis=0), axis=0)
+ e = np.sum(np.sum(np.abs(r) ** 2, axis=0), axis=0)
if calculate_grad is True:
de = 2 * np.sum(np.sum((r * dy.T), axis=2), axis=1)
@@ -335,6 +335,6 @@ def compute(
"""
inputs = self.parameters.as_dict()
log_likelihood = self._observer.log_likelihood(
- self._target, self._observer.time_data, inputs
+ self._target, self._observer.domain_data, inputs
)
return -log_likelihood
diff --git a/pybop/models/base_model.py b/pybop/models/base_model.py
index 47ae1e076..6f880bb8b 100644
--- a/pybop/models/base_model.py
+++ b/pybop/models/base_model.py
@@ -5,8 +5,10 @@
import casadi
import numpy as np
import pybamm
+from scipy.sparse import csc_matrix
+from scipy.sparse.linalg import spsolve
-from pybop import Dataset, Experiment, Parameters, ParameterSet
+from pybop import Dataset, Experiment, Parameters, ParameterSet, SymbolReplacer
from pybop.parameters.parameter import Inputs
@@ -65,7 +67,10 @@ class is designed to be subclassed for creating models with custom behaviour.
"""
def __init__(
- self, name: str = "Base Model", parameter_set: Optional[ParameterSet] = None
+ self,
+ name: str = "Base Model",
+ parameter_set: Optional[ParameterSet] = None,
+ eis=False,
):
"""
Initialise the BaseModel with an optional name and a parameter set.
@@ -91,6 +96,7 @@ def __init__(
(default: True).
"""
self.name = name
+ self.eis = eis
if parameter_set is None:
self._parameter_set = None
elif isinstance(parameter_set, dict):
@@ -144,20 +150,26 @@ def build(
# Clear the model if rebuild required (currently if any initial state)
self.set_initial_state(initial_state, inputs=inputs)
- if dataset is not None:
+ if not self.pybamm_model._built: # noqa: SLF001
+ self.pybamm_model.build_model()
+
+ if self.eis:
+ self.set_up_for_eis(self.pybamm_model)
+ self._parameter_set["Current function [A]"] = 0
+
+ V_scale = getattr(self.pybamm_model.variables["Voltage [V]"], "scale", 1)
+ I_scale = getattr(self.pybamm_model.variables["Current [A]"], "scale", 1)
+ self.z_scale = self._parameter_set.evaluate(V_scale / I_scale)
+
+ if dataset is not None and not self.eis:
self.set_current_function(dataset)
if self._built_model:
return
-
elif self.pybamm_model.is_discretised:
self._model_with_set_params = self.pybamm_model
self._built_model = self.pybamm_model
-
else:
- if not self.pybamm_model._built: # noqa: SLF001
- self.pybamm_model.build_model()
-
self.set_parameters()
self._mesh = pybamm.Mesh(self.geometry, self.submesh_types, self.var_pts)
self._disc = pybamm.Discretisation(
@@ -282,6 +294,55 @@ def set_parameters(self):
self._parameter_set.process_geometry(self._geometry)
self.pybamm_model = self._model_with_set_params
+ def set_up_for_eis(self, model):
+ """
+ Set up the model for electrochemical impedance spectroscopy (EIS) simulations.
+
+ This method sets up the model for EIS simulations by adding the necessary
+ algebraic equations and variables to the model.
+ Originally developed by pybamm-eis: https://github.com/pybamm-team/pybamm-eis
+
+ Parameters
+ ----------
+ model : pybamm.Model
+ The PyBaMM model to be used for EIS simulations.
+ """
+ V_cell = pybamm.Variable("Voltage variable [V]")
+ model.variables["Voltage variable [V]"] = V_cell
+ V = model.variables["Voltage [V]"]
+
+ # Add algebraic equation for the voltage
+ model.algebraic[V_cell] = V_cell - V
+ model.initial_conditions[V_cell] = model.param.ocv_init
+
+ # Create the FunctionControl submodel and extract variables
+ external_circuit_variables = pybamm.external_circuit.FunctionControl(
+ model.param, None, model.options, control="algebraic"
+ ).get_fundamental_variables()
+
+ # Perform the replacement
+ symbol_replacement_map = {
+ model.variables[name]: variable
+ for name, variable in external_circuit_variables.items()
+ }
+
+ # Don't replace initial conditions, as these should not contain
+ # Variable objects
+ replacer = SymbolReplacer(
+ symbol_replacement_map, process_initial_conditions=False
+ )
+ replacer.process_model(model, inplace=True)
+
+ # Add an algebraic equation for the current density variable
+ # External circuit submodels are always equations on the current
+ I_cell = model.variables["Current variable [A]"]
+ I = model.variables["Current [A]"]
+ I_applied = pybamm.FunctionParameter(
+ "Current function [A]", {"Time [s]": pybamm.t}
+ )
+ model.algebraic[I_cell] = I - I_applied
+ model.initial_conditions[I_cell] = 0
+
def clear(self):
"""
Clear any built PyBaMM model.
@@ -431,6 +492,110 @@ def simulate(
return self.solver.solve(self._built_model, inputs=inputs, t_eval=t_eval)
+ def simulateEIS(
+ self, inputs: Inputs, f_eval: list, initial_state: Optional[dict] = None
+ ) -> dict[str, np.ndarray]:
+ """
+ Compute the forward model simulation with electrochemical impedance spectroscopy
+ and return the result.
+
+ Parameters
+ ----------
+ inputs : dict or array-like
+ The input parameters for the simulation. If array-like, it will be
+ converted to a dictionary using the model's fit keys.
+ f_eval : array-like
+ An array of frequency points at which to evaluate the solution.
+
+ Returns
+ -------
+ array-like
+ The simulation result corresponding to the specified signal.
+
+ Raises
+ ------
+ ValueError
+ If the model has not been built before simulation.
+ """
+ inputs = self.parameters.verify(inputs)
+
+ # Build or rebuild if required
+ self.build(inputs=inputs, initial_state=initial_state)
+
+ if not self.check_params(
+ inputs=inputs,
+ allow_infeasible_solutions=self.allow_infeasible_solutions,
+ ):
+ raise ValueError("These parameter values are infeasible.")
+
+ self.initialise_eis_simulation(inputs)
+ zs = [self.calculate_impedance(frequency) for frequency in f_eval]
+
+ return {"Impedance": np.asarray(zs) * self.z_scale}
+
+ def initialise_eis_simulation(self, inputs: Optional[Inputs] = None):
+ """
+ Initialise the Electrochemical Impedance Spectroscopy (EIS) simulation.
+
+ This method sets up the mass matrix and solver, converts inputs to the appropriate format,
+ extracts necessary attributes from the model, and prepares matrices for the simulation.
+
+ Parameters
+ ----------
+ inputs : dict (optional)
+ The input parameters for the simulation.
+ """
+ # Setup mass matrix, solver
+ self.M = self._built_model.mass_matrix.entries
+ self._solver.set_up(self._built_model, inputs=inputs)
+
+ # Convert inputs to casadi format if needed
+ casadi_inputs = (
+ casadi.vertcat(*inputs.values())
+ if inputs is not None and self._built_model.convert_to_format == "casadi"
+ else inputs or []
+ )
+
+ # Extract necessary attributes from the model
+ self.y0 = self._built_model.concatenated_initial_conditions.evaluate(
+ 0, inputs=inputs
+ )
+ self.J = self._built_model.jac_rhs_algebraic_eval(
+ 0, self.y0, casadi_inputs
+ ).sparse()
+
+ # Convert to Compressed Sparse Column format
+ self.M = csc_matrix(self.M)
+ self.J = csc_matrix(self.J)
+
+ # Add forcing to the RHS on the current density
+ self.b = np.zeros(self.y0.shape)
+ self.b[-1] = -1
+
+ def calculate_impedance(self, frequency):
+ """
+ Calculate the impedance for a given frequency.
+
+ This method computes the system matrix, solves the linear system, and calculates
+ the impedance based on the solution.
+
+ Parameters
+ ----------
+ frequency (np.ndarray | list like): The frequency at which to calculate the impedance.
+
+ Returns
+ -------
+ The calculated impedance (complex np.ndarray).
+ """
+ # Compute the system matrix
+ A = 1.0j * 2 * np.pi * frequency * self.M - self.J
+
+ # Solve the system
+ x = spsolve(A, self.b)
+
+ # Calculate the impedance
+ return -x[-2] / x[-1]
+
def simulateS1(
self, inputs: Inputs, t_eval: np.array, initial_state: Optional[dict] = None
):
@@ -658,7 +823,7 @@ def new_copy(self):
"""
model_class = type(self)
if self.pybamm_model is None:
- model_args = {"parameter_set": self.parameter_set.copy()}
+ model_args = {"parameter_set": self._parameter_set.copy()}
else:
model_args = {
"options": self._unprocessed_model.options,
@@ -668,6 +833,7 @@ def new_copy(self):
"var_pts": self.pybamm_model.default_var_pts.copy(),
"spatial_methods": self.pybamm_model.default_spatial_methods.copy(),
"solver": self.pybamm_model.default_solver.copy(),
+ "eis": copy.copy(self.eis),
}
return model_class(**model_args)
diff --git a/pybop/models/empirical/base_ecm.py b/pybop/models/empirical/base_ecm.py
index d329cf384..882a220d4 100644
--- a/pybop/models/empirical/base_ecm.py
+++ b/pybop/models/empirical/base_ecm.py
@@ -45,6 +45,7 @@ def __init__(
var_pts=None,
spatial_methods=None,
solver=None,
+ eis=False,
**model_kwargs,
):
model_options = dict(build=False)
@@ -64,7 +65,7 @@ def __init__(
print("Setting open-circuit voltage to default function")
parameter_set["Open-circuit voltage [V]"] = default_ocp
- super().__init__(name=name, parameter_set=parameter_set)
+ super().__init__(name=name, parameter_set=parameter_set, eis=eis)
self.pybamm_model = pybamm_model
self._unprocessed_model = self.pybamm_model
diff --git a/pybop/models/empirical/ecm.py b/pybop/models/empirical/ecm.py
index f831aee83..95a1a170e 100644
--- a/pybop/models/empirical/ecm.py
+++ b/pybop/models/empirical/ecm.py
@@ -38,8 +38,12 @@ class Thevenin(ECircuitModel):
def __init__(
self,
name="Equivalent Circuit Thevenin Model",
+ eis=False,
**model_kwargs,
):
super().__init__(
- pybamm_model=pybamm_equivalent_circuit.Thevenin, name=name, **model_kwargs
+ pybamm_model=pybamm_equivalent_circuit.Thevenin,
+ name=name,
+ eis=eis,
+ **model_kwargs,
)
diff --git a/pybop/models/lithium_ion/base_echem.py b/pybop/models/lithium_ion/base_echem.py
index ae709825a..e919eb317 100644
--- a/pybop/models/lithium_ion/base_echem.py
+++ b/pybop/models/lithium_ion/base_echem.py
@@ -46,9 +46,10 @@ def __init__(
var_pts=None,
spatial_methods=None,
solver=None,
+ eis=False,
**model_kwargs,
):
- super().__init__(name=name, parameter_set=parameter_set)
+ super().__init__(name=name, parameter_set=parameter_set, eis=eis)
model_options = dict(build=False)
for key, value in model_kwargs.items():
diff --git a/pybop/models/lithium_ion/echem.py b/pybop/models/lithium_ion/echem.py
index 9fdd308e3..aac05b738 100644
--- a/pybop/models/lithium_ion/echem.py
+++ b/pybop/models/lithium_ion/echem.py
@@ -38,11 +38,13 @@ class SPM(EChemBaseModel):
def __init__(
self,
name="Single Particle Model",
+ eis=False,
**model_kwargs,
):
super().__init__(
pybamm_model=pybamm_lithium_ion.SPM,
name=name,
+ eis=eis,
**model_kwargs,
)
@@ -83,10 +85,11 @@ class SPMe(EChemBaseModel):
def __init__(
self,
name="Single Particle Model with Electrolyte",
+ eis=False,
**model_kwargs,
):
super().__init__(
- pybamm_model=pybamm_lithium_ion.SPMe, name=name, **model_kwargs
+ pybamm_model=pybamm_lithium_ion.SPMe, name=name, eis=eis, **model_kwargs
)
@@ -126,9 +129,12 @@ class DFN(EChemBaseModel):
def __init__(
self,
name="Doyle-Fuller-Newman",
+ eis=False,
**model_kwargs,
):
- super().__init__(pybamm_model=pybamm_lithium_ion.DFN, name=name, **model_kwargs)
+ super().__init__(
+ pybamm_model=pybamm_lithium_ion.DFN, name=name, eis=eis, **model_kwargs
+ )
class MPM(EChemBaseModel):
@@ -165,10 +171,12 @@ class MPM(EChemBaseModel):
def __init__(
self,
name="Many Particle Model",
+ eis=False,
**model_kwargs,
):
super().__init__(
pybamm_model=pybamm_lithium_ion.MPM,
+ eis=eis,
name=name,
**model_kwargs,
)
@@ -208,11 +216,13 @@ class MSMR(EChemBaseModel):
def __init__(
self,
name="Multi Species Multi Reactions Model",
+ eis=False,
**model_kwargs,
):
super().__init__(
pybamm_model=pybamm_lithium_ion.MSMR,
name=name,
+ eis=eis,
**model_kwargs,
)
@@ -241,5 +251,7 @@ class WeppnerHuggins(EChemBaseModel):
The solver to use for simulating the model. If None, the default solver from PyBaMM is used.
"""
- def __init__(self, name="Weppner & Huggins model", **model_kwargs):
- super().__init__(pybamm_model=BaseWeppnerHuggins, name=name, **model_kwargs)
+ def __init__(self, name="Weppner & Huggins model", eis=False, **model_kwargs):
+ super().__init__(
+ pybamm_model=BaseWeppnerHuggins, name=name, eis=eis, **model_kwargs
+ )
diff --git a/pybop/models/lithium_ion/weppner_huggins.py b/pybop/models/lithium_ion/weppner_huggins.py
index b8707cca9..9703269d3 100644
--- a/pybop/models/lithium_ion/weppner_huggins.py
+++ b/pybop/models/lithium_ion/weppner_huggins.py
@@ -99,6 +99,7 @@ def __init__(self, name="Weppner & Huggins model", **model_kwargs):
self.variables = {
"Voltage [V]": V,
"Time [s]": t,
+ "Current [A]": self.param.current_with_time,
}
# Set the built property on creation to prevent unnecessary model rebuilds
diff --git a/pybop/observers/observer.py b/pybop/observers/observer.py
index 81247c7c4..239f352a1 100644
--- a/pybop/observers/observer.py
+++ b/pybop/observers/observer.py
@@ -159,13 +159,13 @@ def evaluate(self, inputs: Inputs):
if self._dataset is not None:
for signal in self.signal:
ym = self._target[signal]
- for i, t in enumerate(self._time_data):
+ for i, t in enumerate(self._domain_data):
self.observe(t, ym[i])
ys.append(self.get_current_measure())
output[signal] = np.vstack(ys)
else:
for signal in self.signal:
- for t in self._time_data:
+ for t in self._domain_data:
self.observe(t)
ys.append(self.get_current_measure())
output[signal] = np.vstack(ys)
diff --git a/pybop/observers/unscented_kalman.py b/pybop/observers/unscented_kalman.py
index bdad6cc5d..69e019ee6 100644
--- a/pybop/observers/unscented_kalman.py
+++ b/pybop/observers/unscented_kalman.py
@@ -70,11 +70,11 @@ def __init__(
)
if dataset is not None:
# Check that the dataset contains necessary variables
- dataset.check([*self.signal, "Current function [A]"])
+ dataset.check(signal=[*self.signal, "Current function [A]"])
self._dataset = dataset.data
- self._time_data = self._dataset["Time [s]"]
- self.n_time_data = len(self._time_data)
+ self._domain_data = self._dataset["Time [s]"]
+ self.n_data = len(self._domain_data)
self._target = {signal: self._dataset[signal] for signal in self.signal}
# Observer initiation
diff --git a/pybop/plotting/nyquist.py b/pybop/plotting/nyquist.py
new file mode 100644
index 000000000..074c60b92
--- /dev/null
+++ b/pybop/plotting/nyquist.py
@@ -0,0 +1,145 @@
+import sys
+
+from pybop import StandardPlot
+from pybop.parameters.parameter import Inputs
+
+
+def nyquist(problem, problem_inputs: Inputs = None, show=True, **layout_kwargs):
+ """
+ Generates Nyquist plots for the given problem by evaluating the model's output and target values.
+
+ Parameters
+ ----------
+ problem : pybop.BaseProblem
+ An instance of a problem class (e.g., `pybop.EISProblem`) that contains the parameters and methods
+ for evaluation and target retrieval.
+ problem_inputs : Inputs, optional
+ Input parameters for the problem. If not provided, the default parameters from the problem
+ instance will be used. These parameters are verified before use (default is None).
+ show : bool, optional
+ If True, the plots will be displayed. If running in an IPython kernel (e.g., Jupyter Notebook),
+ the plots will be shown using SVG format for better quality (default is True).
+ **layout_kwargs : dict, optional
+ Additional keyword arguments for customising the plot layout. These arguments are passed to
+ `fig.update_layout()`.
+
+ Returns
+ -------
+ list
+ A list of plotly `Figure` objects, each representing a Nyquist plot for the model's output and target values.
+
+ Notes
+ -----
+ - The function extracts the real part of the impedance from the model's output and the real and imaginary parts
+ of the impedance from the target output.
+ - For each signal in the problem, a Nyquist plot is created with the model's impedance plotted as a scatter plot.
+ - An additional trace for the reference (target output) is added to the plot.
+ - The plot layout can be customised using `layout_kwargs`.
+
+ Example
+ -------
+ >>> problem = pybop.EISProblem()
+ >>> nyquist_figures = nyquist(problem, show=True, title="Nyquist Plot", xaxis_title="Real(Z)", yaxis_title="Imag(Z)")
+ >>> # The plots will be displayed and nyquist_figures will contain the list of figure objects.
+ """
+ if problem_inputs is None:
+ problem_inputs = problem.parameters.as_dict()
+ else:
+ problem_inputs = problem.parameters.verify(problem_inputs)
+
+ model_output = problem.evaluate(problem_inputs)
+ domain_data = model_output["Impedance"].real
+ target_output = problem.get_target()
+
+ figure_list = []
+ for i in problem.signal:
+ default_layout_options = dict(
+ title="Nyquist Plot",
+ font=dict(family="Arial", size=14),
+ plot_bgcolor="white",
+ paper_bgcolor="white",
+ xaxis=dict(
+ title=dict(text="Zre / Ω", font=dict(size=16), standoff=15),
+ showline=True,
+ linewidth=2,
+ linecolor="black",
+ mirror=True,
+ ticks="outside",
+ tickwidth=2,
+ tickcolor="black",
+ ticklen=5,
+ ),
+ yaxis=dict(
+ title=dict(text="-Zim / Ω", font=dict(size=16), standoff=15),
+ showline=True,
+ linewidth=2,
+ linecolor="black",
+ mirror=True,
+ ticks="outside",
+ tickwidth=2,
+ tickcolor="black",
+ ticklen=5,
+ scaleanchor="x",
+ scaleratio=1,
+ ),
+ legend=dict(
+ x=0.02,
+ y=0.98,
+ bgcolor="rgba(255, 255, 255, 0.5)",
+ bordercolor="black",
+ borderwidth=1,
+ ),
+ width=600,
+ height=600,
+ )
+
+ plot_dict = StandardPlot(
+ x=domain_data,
+ y=-model_output[i].imag,
+ layout_options=default_layout_options,
+ trace_names="Model",
+ )
+
+ plot_dict.traces[0].update(
+ mode="lines+markers",
+ line=dict(color="blue", width=2),
+ marker=dict(size=8, color="blue", symbol="circle"),
+ )
+
+ target_trace = plot_dict.create_trace(
+ x=target_output[i].real,
+ y=-target_output[i].imag,
+ name="Reference",
+ mode="markers",
+ marker=dict(size=8, color="red", symbol="circle-open"),
+ showlegend=True,
+ )
+ plot_dict.traces.append(target_trace)
+
+ fig = plot_dict(show=False)
+
+ # Add minor gridlines
+ fig.update_xaxes(
+ showgrid=True,
+ gridwidth=1,
+ gridcolor="lightgray",
+ minor=dict(showgrid=True, gridwidth=0.5, gridcolor="lightgray"),
+ )
+ fig.update_yaxes(
+ showgrid=True,
+ gridwidth=1,
+ gridcolor="lightgray",
+ minor=dict(showgrid=True, gridwidth=0.5, gridcolor="lightgray"),
+ )
+
+ # Overwrite with user-kwargs
+ fig.update_layout(**layout_kwargs)
+
+ if "ipykernel" in sys.modules and show:
+ fig.show("svg")
+ elif show:
+ fig.show()
+
+ figure_list.append(fig)
+
+ return figure_list
diff --git a/pybop/plotting/plot_dataset.py b/pybop/plotting/plot_dataset.py
index ecc84aa6e..5dd7bc246 100644
--- a/pybop/plotting/plot_dataset.py
+++ b/pybop/plotting/plot_dataset.py
@@ -31,7 +31,7 @@ def plot_dataset(dataset, signal=None, trace_names=None, show=True, **layout_kwa
# Get data dictionary
if signal is None:
signal = ["Voltage [V]"]
- dataset.check(signal)
+ dataset.check(signal=signal)
# Compile ydata and labels or legend
y = [dataset[s] for s in signal]
diff --git a/pybop/plotting/plot_problem.py b/pybop/plotting/plot_problem.py
index 32ef38b2c..90cd7d069 100644
--- a/pybop/plotting/plot_problem.py
+++ b/pybop/plotting/plot_problem.py
@@ -37,7 +37,7 @@ def quick_plot(problem, problem_inputs: Inputs = None, show=True, **layout_kwarg
problem_inputs = problem.parameters.verify(problem_inputs)
# Extract the time data and evaluate the model's output and target values
- xaxis_data = problem.time_data
+ xaxis_data = problem.domain_data
model_output = problem.evaluate(problem_inputs)
target_output = problem.get_target()
@@ -53,13 +53,13 @@ def quick_plot(problem, problem_inputs: Inputs = None, show=True, **layout_kwarg
# Create a plotting dictionary
if isinstance(problem, DesignProblem):
trace_name = "Optimised"
- opt_time_data = model_output["Time [s]"]
+ opt_domain_data = model_output["Time [s]"]
else:
trace_name = "Model"
- opt_time_data = xaxis_data
+ opt_domain_data = xaxis_data
plot_dict = StandardPlot(
- x=opt_time_data,
+ x=opt_domain_data,
y=model_output[i],
layout_options=default_layout_options,
trace_names=trace_name,
diff --git a/pybop/problems/base_problem.py b/pybop/problems/base_problem.py
index 5b9f4ec11..59c7ce286 100644
--- a/pybop/problems/base_problem.py
+++ b/pybop/problems/base_problem.py
@@ -63,15 +63,24 @@ def __init__(
self.parameters.reset_initial_value()
self._model = model
+ self.eis = False
+ self.domain = "Time [s]"
self.check_model = check_model
self.signal = signal or ["Voltage [V]"]
- self.additional_variables = additional_variables or []
self.set_initial_state(initial_state)
self._dataset = None
- self._time_data = None
+ self._domain_data = None
self._target = None
self.verbose = False
self.failure_output = np.asarray([np.inf])
+ if isinstance(self._model, BaseModel):
+ self.eis = self.model.eis
+ self.domain = "Frequency [Hz]" if self.eis else "Time [s]"
+
+ # 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)
@@ -106,7 +115,7 @@ def n_parameters(self):
def n_outputs(self):
return len(self.signal)
- def evaluate(self, inputs: Inputs):
+ def evaluate(self, inputs: Inputs, eis=False):
"""
Evaluate the model with the given parameters and return the signal.
@@ -179,12 +188,12 @@ def target(self, target):
self._target = target
@property
- def time_data(self):
- return self._time_data
+ def domain_data(self):
+ return self._domain_data
- @time_data.setter
- def time_data(self, time_data):
- self._time_data = time_data
+ @domain_data.setter
+ def domain_data(self, domain_data):
+ self._domain_data = domain_data
@property
def dataset(self):
diff --git a/pybop/problems/design_problem.py b/pybop/problems/design_problem.py
index 57c3597c1..15279b463 100644
--- a/pybop/problems/design_problem.py
+++ b/pybop/problems/design_problem.py
@@ -42,11 +42,6 @@ def __init__(
additional_variables: Optional[list[str]] = None,
initial_state: Optional[dict] = None,
):
- # Add time and current and remove duplicates
- additional_variables = additional_variables or []
- additional_variables.extend(["Time [s]", "Current [A]"])
- additional_variables = list(set(additional_variables))
-
super().__init__(
parameters, model, check_model, signal, additional_variables, initial_state
)
@@ -54,7 +49,7 @@ def __init__(
# Add an example dataset for plotting comparison
sol = self.evaluate(self.parameters.as_dict("initial"))
- self._time_data = sol["Time [s]"]
+ self._domain_data = sol["Time [s]"]
self._target = {key: sol[key] for key in self.signal}
self._dataset = None
self.warning_patterns = [
@@ -89,7 +84,7 @@ def set_initial_state(self, initial_state):
self.initial_state = initial_state
- def evaluate(self, inputs: Inputs, update_capacity=False):
+ def evaluate(self, inputs: Inputs, eis=False, update_capacity=False):
"""
Evaluate the model with the given parameters and return the signal.
diff --git a/pybop/problems/fitting_problem.py b/pybop/problems/fitting_problem.py
index 252a52d43..25609b6b1 100644
--- a/pybop/problems/fitting_problem.py
+++ b/pybop/problems/fitting_problem.py
@@ -34,10 +34,10 @@ class FittingProblem(BaseProblem):
---------------------
dataset : dictionary
The dictionary from a Dataset object containing the signal keys and values to fit the model to.
- time_data : np.ndarray
- The time points in the dataset.
- n_time_data : int
- The number of time points.
+ domain_data : np.ndarray
+ The domain points in the dataset.
+ n_domain_data : int
+ The number of domain points.
target : np.ndarray
The target values of the signals.
"""
@@ -52,11 +52,6 @@ def __init__(
additional_variables: Optional[list[str]] = None,
initial_state: Optional[dict] = None,
):
- # Add time and remove duplicates
- additional_variables = additional_variables or []
- additional_variables.extend(["Time [s]"])
- additional_variables = list(set(additional_variables))
-
super().__init__(
parameters, model, check_model, signal, additional_variables, initial_state
)
@@ -64,14 +59,14 @@ def __init__(
self._n_parameters = len(self.parameters)
# Check that the dataset contains necessary variables
- dataset.check([*self.signal, "Current function [A]"])
+ dataset.check(domain=self.domain, signal=[*self.signal, "Current function [A]"])
- # Unpack time and target data
- self._time_data = self._dataset["Time [s]"]
- self.n_time_data = len(self._time_data)
+ # Unpack domain and target data
+ self._domain_data = self._dataset[self.domain]
+ self.n_data = len(self._domain_data)
self.set_target(dataset)
- if model is not None:
+ if self._model is not None:
# Build the model from scratch
if self._model.built_model is not None:
self._model.clear()
@@ -107,7 +102,9 @@ def set_initial_state(self, initial_state: Optional[dict] = None):
self.initial_state = initial_state
def evaluate(
- self, inputs: Inputs, update_capacity=False
+ self,
+ inputs: Inputs,
+ update_capacity=False,
) -> dict[str, np.ndarray[np.float64]]:
"""
Evaluate the model with the given parameters and return the signal.
@@ -120,23 +117,56 @@ def evaluate(
Returns
-------
y : np.ndarray
- The model output y(t) simulated with given inputs.
+ The simulated model output y(t) for self.eis == False, and y(ω) for self.eis == True for the given inputs.
"""
inputs = self.parameters.verify(inputs)
+ if self.eis:
+ return self._evaluateEIS(inputs, update_capacity=update_capacity)
+ else:
+ try:
+ sol = self._model.simulate(
+ inputs=inputs,
+ t_eval=self._domain_data,
+ initial_state=self.initial_state,
+ )
+ except Exception as e:
+ if self.verbose:
+ print(f"Simulation error: {e}")
+ return {signal: self.failure_output for signal in self.signal}
+ return {
+ signal: sol[signal].data
+ for signal in (self.signal + self.additional_variables)
+ }
+
+ def _evaluateEIS(
+ self, inputs: Inputs, update_capacity=False
+ ) -> dict[str, np.ndarray[np.float64]]:
+ """
+ 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 simulated model output y(ω) for the given inputs.
+ """
try:
- sol = self._model.simulate(
- inputs=inputs, t_eval=self._time_data, initial_state=self.initial_state
+ sol = self._model.simulateEIS(
+ inputs=inputs,
+ f_eval=self._domain_data,
+ initial_state=self.initial_state,
)
except Exception as e:
if self.verbose:
print(f"Simulation error: {e}")
return {signal: self.failure_output for signal in self.signal}
- return {
- signal: sol[signal].data
- for signal in (self.signal + self.additional_variables)
- }
+ return sol
def evaluateS1(self, inputs: Inputs):
"""
@@ -158,7 +188,9 @@ def evaluateS1(self, inputs: Inputs):
try:
sol = self._model.simulateS1(
- inputs=inputs, t_eval=self._time_data, initial_state=self.initial_state
+ inputs=inputs,
+ t_eval=self._domain_data,
+ initial_state=self.initial_state,
)
except Exception as e:
print(f"Error: {e}")
@@ -186,4 +218,4 @@ def evaluateS1(self, inputs: Inputs):
axis=-1,
)
- return (y, np.asarray(dy))
+ return y, np.asarray(dy)
diff --git a/pybop/problems/multi_fitting_problem.py b/pybop/problems/multi_fitting_problem.py
index 263fc92dc..7b9f0dbb7 100644
--- a/pybop/problems/multi_fitting_problem.py
+++ b/pybop/problems/multi_fitting_problem.py
@@ -37,30 +37,31 @@ def __init__(self, *args):
combined_parameters.join(problem.parameters)
# Combine the target datasets
- combined_time_data = []
+ combined_domain_data = []
combined_signal = []
for problem in self.problems:
for signal in problem.signal:
- combined_time_data.extend(problem.time_data)
+ combined_domain_data.extend(problem.domain_data)
combined_signal.extend(problem.target[signal])
- combined_dataset = Dataset(
- {
- "Time [s]": np.asarray(combined_time_data),
- "Combined signal": np.asarray(combined_signal),
- }
- )
super().__init__(
parameters=combined_parameters,
model=None,
signal=["Combined signal"],
)
+
+ combined_dataset = Dataset(
+ {
+ self.domain: np.asarray(combined_domain_data),
+ "Combined signal": np.asarray(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)
+ # Unpack domain and target data
+ self._domain_data = self._dataset[self.domain]
+ self.n_domain_data = len(self._domain_data)
self.set_target(combined_dataset)
def set_initial_state(self, initial_state: Optional[dict] = None):
@@ -75,7 +76,7 @@ def set_initial_state(self, initial_state: Optional[dict] = None):
for problem in self.problems:
problem.set_initial_state(initial_state)
- def evaluate(self, inputs: Inputs, **kwargs):
+ def evaluate(self, inputs: Inputs, eis=False, **kwargs):
"""
Evaluate the model with the given parameters and return the signal.
diff --git a/tests/integration/test_eis_parameterisation.py b/tests/integration/test_eis_parameterisation.py
new file mode 100644
index 000000000..8022b67dc
--- /dev/null
+++ b/tests/integration/test_eis_parameterisation.py
@@ -0,0 +1,189 @@
+import numpy as np
+import pytest
+
+import pybop
+
+
+class TestEISParameterisation:
+ """
+ A class to test the eis parameterisation methods.
+ """
+
+ @pytest.fixture(autouse=True)
+ def setup(self):
+ self.sigma0 = 5e-4
+ self.ground_truth = np.clip(
+ np.asarray([0.55, 0.55]) + np.random.normal(loc=0.0, scale=0.05, size=2),
+ a_min=0.4,
+ a_max=0.75,
+ )
+
+ @pytest.fixture
+ def model(self):
+ parameter_set = pybop.ParameterSet.pybamm("Chen2020")
+ x = self.ground_truth
+ parameter_set.update(
+ {
+ "Negative electrode active material volume fraction": x[0],
+ "Positive electrode active material volume fraction": x[1],
+ }
+ )
+ return pybop.lithium_ion.SPM(
+ parameter_set=parameter_set,
+ eis=True,
+ options={"surface form": "differential"},
+ )
+
+ @pytest.fixture
+ def parameters(self):
+ return pybop.Parameters(
+ pybop.Parameter(
+ "Negative electrode active material volume fraction",
+ prior=pybop.Uniform(0.4, 0.75),
+ bounds=[0.375, 0.775],
+ ),
+ pybop.Parameter(
+ "Positive electrode active material volume fraction",
+ prior=pybop.Uniform(0.4, 0.75),
+ bounds=[0.375, 0.775],
+ ),
+ )
+
+ @pytest.fixture(params=[0.5])
+ def init_soc(self, request):
+ return request.param
+
+ @pytest.fixture(
+ params=[
+ pybop.GaussianLogLikelihoodKnownSigma,
+ pybop.GaussianLogLikelihood,
+ pybop.RootMeanSquaredError,
+ pybop.SumSquaredError,
+ pybop.SumofPower,
+ pybop.Minkowski,
+ pybop.MAP,
+ ]
+ )
+ def cost(self, request):
+ return request.param
+
+ def noise(self, 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
+
+ @pytest.fixture(
+ params=[
+ pybop.SciPyDifferentialEvolution,
+ pybop.CMAES,
+ pybop.CuckooSearch,
+ pybop.NelderMead,
+ pybop.SNES,
+ pybop.XNES,
+ ]
+ )
+ def optimiser(self, request):
+ return request.param
+
+ @pytest.fixture
+ def optim(self, optimiser, model, parameters, cost, init_soc):
+ n_frequency = 15
+ # Set frequency set
+ f_eval = np.logspace(-4, 5, n_frequency)
+
+ # Form dataset
+ solution = self.get_data(model, init_soc, f_eval)
+ dataset = pybop.Dataset(
+ {
+ "Frequency [Hz]": f_eval,
+ "Current function [A]": np.ones(n_frequency) * 0.0,
+ "Impedance": solution["Impedance"]
+ + self.noise(self.sigma0, len(solution["Impedance"])),
+ }
+ )
+
+ # Define the problem
+ signal = ["Impedance"]
+ problem = pybop.FittingProblem(model, parameters, dataset, signal=signal)
+
+ # Construct the cost
+ if cost in [pybop.GaussianLogLikelihoodKnownSigma]:
+ cost = cost(problem, sigma0=self.sigma0)
+ elif cost in [pybop.GaussianLogLikelihood]:
+ cost = cost(problem, sigma0=self.sigma0 * 4) # Initial sigma0 guess
+ elif cost in [pybop.MAP]:
+ cost = cost(
+ problem, pybop.GaussianLogLikelihoodKnownSigma, sigma0=self.sigma0
+ )
+ elif cost in [pybop.SumofPower, pybop.Minkowski]:
+ cost = cost(problem, p=2)
+ else:
+ cost = cost(problem)
+
+ # Construct optimisation object
+ common_args = {
+ "cost": cost,
+ "max_iterations": 250,
+ "absolute_tolerance": 1e-6,
+ "max_unchanged_iterations": 35,
+ }
+
+ if isinstance(cost, pybop.MAP):
+ for i in cost.parameters.keys():
+ cost.parameters[i].prior = pybop.Uniform(
+ 0.2, 2.0
+ ) # Increase range to avoid prior == np.inf
+
+ # Set sigma0 and create optimiser
+ sigma0 = 0.05 if isinstance(cost, pybop.MAP) else None
+ optim = optimiser(sigma0=sigma0, **common_args)
+
+ return optim
+
+ @pytest.mark.integration
+ def test_eis_optimisers(self, optim):
+ x0 = optim.parameters.initial_value()
+
+ # Add sigma0 to ground truth for GaussianLogLikelihood
+ if isinstance(optim.cost, pybop.GaussianLogLikelihood):
+ self.ground_truth = np.concatenate(
+ (self.ground_truth, np.asarray([self.sigma0]))
+ )
+
+ initial_cost = optim.cost(x0)
+ x, final_cost = optim.run()
+
+ # Assertions
+ if np.allclose(x0, self.ground_truth, atol=1e-5):
+ raise AssertionError("Initial guess is too close to ground truth")
+
+ # Assert on identified values, without sigma for GaussianLogLikelihood
+ # as the sigma values are small (5e-4), this is a difficult identification process
+ # and requires a high number of iterations, and parameter dependent step sizes.
+ if optim.minimising:
+ assert initial_cost > final_cost
+ else:
+ assert initial_cost < final_cost
+ np.testing.assert_allclose(x, self.ground_truth, atol=1.5e-2)
+
+ def get_data(self, model, init_soc, f_eval):
+ initial_state = {"Initial SoC": init_soc}
+ sim = model.simulateEIS(
+ inputs={
+ "Negative electrode active material volume fraction": self.ground_truth[
+ 0
+ ],
+ "Positive electrode active material volume fraction": self.ground_truth[
+ 1
+ ],
+ },
+ f_eval=f_eval,
+ initial_state=initial_state,
+ )
+
+ return sim
diff --git a/tests/integration/test_spm_parameterisations.py b/tests/integration/test_spm_parameterisations.py
index 482f59bbd..1392ca83d 100644
--- a/tests/integration/test_spm_parameterisations.py
+++ b/tests/integration/test_spm_parameterisations.py
@@ -13,8 +13,10 @@ class Test_SPM_Parameterisation:
@pytest.fixture(autouse=True)
def setup(self):
self.sigma0 = 0.002
- self.ground_truth = np.asarray([0.55, 0.55]) + np.random.normal(
- loc=0.0, scale=0.05, size=2
+ self.ground_truth = np.clip(
+ np.asarray([0.55, 0.55]) + np.random.normal(loc=0.0, scale=0.05, size=2),
+ a_min=0.4,
+ a_max=0.75,
)
@pytest.fixture
diff --git a/tests/unit/test_dataset.py b/tests/unit/test_dataset.py
index 2267f576c..bb881429c 100644
--- a/tests/unit/test_dataset.py
+++ b/tests/unit/test_dataset.py
@@ -55,3 +55,14 @@ def test_dataset(self):
# Test conversion of single signal to list
assert dataset.check()
+
+ # Form frequency dataset
+ data_dictionary = {
+ "Frequency [Hz]": np.linspace(-10, 0, 10),
+ "Current [A]": np.zeros(10),
+ "Impedance": np.zeros(10),
+ }
+ frequency_dataset = pybop.Dataset(data_dictionary)
+
+ with pytest.raises(ValueError, match="Frequencies cannot be negative."):
+ frequency_dataset.check(domain="Frequency [Hz]", signal="Impedance")
diff --git a/tests/unit/test_likelihoods.py b/tests/unit/test_likelihoods.py
index dfc36a3c3..cd1d676b5 100644
--- a/tests/unit/test_likelihoods.py
+++ b/tests/unit/test_likelihoods.py
@@ -70,7 +70,7 @@ def test_base_likelihood_init(self, problem_name, n_outputs, request):
likelihood = pybop.BaseLikelihood(problem)
assert likelihood.problem == problem
assert likelihood.n_outputs == n_outputs
- assert likelihood.n_time_data == problem.n_time_data
+ assert likelihood.n_data == problem.n_data
assert likelihood.n_parameters == 1
assert np.array_equal(likelihood.target, problem.target)
diff --git a/tests/unit/test_models.py b/tests/unit/test_models.py
index f1722366e..02d86f84a 100644
--- a/tests/unit/test_models.py
+++ b/tests/unit/test_models.py
@@ -313,6 +313,28 @@ def test_simulate(self):
with pytest.raises(ValueError):
ExponentialDecay(n_states=-1)
+ @pytest.mark.unit
+ def test_simulateEIS(self):
+ # Test EIS on SPM
+ model = pybop.lithium_ion.SPM(eis=True)
+
+ # Construct frequencies and solve
+ f_eval = np.linspace(100, 1000, 5)
+ sol = model.simulateEIS(inputs={}, f_eval=f_eval)
+ assert np.isfinite(sol["Impedance"]).all()
+
+ # Test infeasible parameter values
+ model.allow_infeasible_solutions = False
+ inputs = {
+ "Negative electrode active material volume fraction": 0.9,
+ "Positive electrode active material volume fraction": 0.9,
+ }
+ # Rebuild model
+ model.build(inputs=inputs)
+
+ with pytest.raises(ValueError, match="These parameter values are infeasible."):
+ model.simulateEIS(f_eval=f_eval, inputs=inputs)
+
@pytest.mark.unit
def test_basemodel(self):
base = pybop.BaseModel()
diff --git a/tests/unit/test_observers.py b/tests/unit/test_observers.py
index 5a784b3c2..1fe631ddf 100644
--- a/tests/unit/test_observers.py
+++ b/tests/unit/test_observers.py
@@ -73,7 +73,7 @@ def test_observer(self, model, parameters):
assert np.shape(covariance) == (n, n)
# Test evaluate with different inputs
- observer._time_data = t_eval
+ observer._domain_data = t_eval
observer.evaluate(parameters.as_dict())
observer.evaluate(parameters.current_value())
diff --git a/tests/unit/test_plots.py b/tests/unit/test_plots.py
index c110312dd..27f5b7941 100644
--- a/tests/unit/test_plots.py
+++ b/tests/unit/test_plots.py
@@ -218,3 +218,38 @@ def test_plot2d_prior_bounds(self, model, dataset):
):
warnings.simplefilter("always")
pybop.plot2d(cost)
+
+ @pytest.mark.unit
+ def test_nyquist(self):
+ # Define model
+ model = pybop.lithium_ion.SPM(
+ eis=True, options={"surface form": "differential"}
+ )
+
+ # Fitting parameters
+ parameters = pybop.Parameters(
+ pybop.Parameter(
+ "Positive electrode thickness [m]",
+ prior=pybop.Gaussian(60e-6, 1e-6),
+ bounds=[10e-6, 80e-6],
+ )
+ )
+
+ # Form dataset
+ dataset = pybop.Dataset(
+ {
+ "Frequency [Hz]": np.logspace(-4, 5, 10),
+ "Current function [A]": np.ones(10) * 0.0,
+ "Impedance": np.ones(10) * 0.0,
+ }
+ )
+
+ signal = ["Impedance"]
+ # Generate problem, cost function, and optimisation class
+ problem = pybop.FittingProblem(model, parameters, dataset, signal=signal)
+
+ # Plot the nyquist
+ pybop.nyquist(problem, problem_inputs=[60e-6], title="Optimised Comparison")
+
+ # Without inputs
+ pybop.nyquist(problem, title="Optimised Comparison")
diff --git a/tests/unit/test_problem.py b/tests/unit/test_problem.py
index b7fded8d6..3724cf847 100644
--- a/tests/unit/test_problem.py
+++ b/tests/unit/test_problem.py
@@ -150,6 +150,11 @@ def test_fitting_problem(self, parameters, dataset, model, signal):
# Test model.simulate with an initial state
problem.evaluate(inputs=[1e-5, 1e-5])
+ # Test try-except
+ problem.verbose = True
+ out = problem.evaluate(inputs=[0.0, 0.0])
+ assert not np.isfinite(out["Voltage [V]"])
+
# Test problem construction errors
for bad_dataset in [
pybop.Dataset({"Time [s]": np.array([0])}),
@@ -182,6 +187,34 @@ def test_fitting_problem(self, parameters, dataset, model, signal):
with pytest.raises(ValueError):
pybop.FittingProblem(model, parameters, bad_dataset, signal=two_signals)
+ @pytest.mark.unit
+ def test_fitting_problem_eis(self, parameters):
+ model = pybop.lithium_ion.SPM(eis=True)
+
+ dataset = pybop.Dataset(
+ {
+ "Frequency [Hz]": np.logspace(-4, 5, 30),
+ "Current function [A]": np.ones(30) * 0.0,
+ "Impedance": np.ones(30) * 0.0,
+ }
+ )
+
+ # Construct Problem
+ problem = pybop.FittingProblem(
+ model,
+ parameters,
+ dataset,
+ signal=["Impedance"],
+ initial_state={"Initial open-circuit voltage [V]": 4.0},
+ )
+ assert problem.eis == model.eis
+ assert problem.domain == "Frequency [Hz]"
+
+ # Test try-except
+ problem.verbose = True
+ out = problem.evaluate(inputs=[0.0, 0.0])
+ assert not np.isfinite(out["Impedance"])
+
@pytest.mark.unit
def test_multi_fitting_problem(self, model, parameters, dataset, signal):
problem_1 = pybop.FittingProblem(model, parameters, dataset, signal=signal)