diff --git a/CHANGELOG.md b/CHANGELOG.md index f683a96db..b90957881 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,7 @@ ## Features +- [#975](https://github.com/pybop-team/PyBOP/pull/975) - Adds coupled GITT-EIS parameterisation: `pybop.pybamm.EISSimulator` accepts a `protocol` and computes an impedance spectrum about the state reached at each acquisition time, solving the time-domain trajectory once. Adds `pybop.get_impedance_variables` and `pybop.parse_impedance_variables` to name the impedance variables of a `Dataset`, plots the acquired spectra as Nyquist subplots in `pybop.plot.problem`, and adds an example script `gitt_eis.py`. - [#969](https://github.com/pybop-team/PyBOP/pull/969) - Updates synthetic data and adds example script for thermal parameterisation. - [#965](https://github.com/pybop-team/PyBOP/pull/965) - Adds synthetic data and example scripts for OCV parameterisation. - [#963](https://github.com/pybop-team/PyBOP/pull/963) - Adds an example for generating synthetic data from a specification and exporting it to a PyProBE-compatible parquet file. @@ -12,18 +13,21 @@ ## Optimisations +- [#975](https://github.com/pybop-team/PyBOP/pull/975) - Set up the EIS solver, mass matrix and forcing vector once rather than on every evaluation. - [#967](https://github.com/pybop-team/PyBOP/pull/967) - Add `Dataset.get_discontinuities` and update the `pybop.pybamm.RecommendedSolver` options. - [#946](https://github.com/pybop-team/PyBOP/pull/946) - Use `vectorized` evaluation for SciPy differential evolution by default instead of multiprocessing `workers`. - [#925](https://github.com/pybop-team/PyBOP/pull/925) - Add `UnboundedDistribution` and the `get_transformed_distribution` functionality. ## Bug Fixes +- [#975](https://github.com/pybop-team/PyBOP/pull/975) - Fixes the ordering of the input parameters when evaluating the EIS Jacobian, which swapped their values for two or more parameters, and ensures that each cost within a `WeightedCost` keeps its own target. - [#915](https://github.com/pybop-team/PyBOP/pull/915) - Fixes axis labels for non-standard domain names, adds `Dataset` length property and adds `kind` property to `Interpolant`. - [#911](https://github.com/pybop-team/PyBOP/pull/911) - Fixes the passing of the cost log to the Voronoi surface plot. - [#905](https://github.com/pybop-team/PyBOP/pull/905) - Remove restriction on numpy. ## Breaking Changes +- [#975](https://github.com/pybop-team/PyBOP/pull/975) - `pybop.pybamm.EISSimulator` now requires the `surface form` model option and no longer modifies the model passed to it. Adds `equal_aspect` to the plotting backends, removes the unused `Simulator.time_data` property and removes the superseded `operando_eis.py` example. - [#928](https://github.com/pybop-team/PyBOP/pull/928) - Deprecates `StandardPlot` and `StandardSubplot` in favour of new standardised backend functionality. - [#960](https://github.com/pybop-team/PyBOP/pull/960) - Remove `asv` benchmarking. - [#938](https://github.com/pybop-team/PyBOP/pull/938) - Make SALib an optional dependency and remove `sensitivity_analysis` in favour of using SALib directly. diff --git a/examples/scripts/battery_parameterisation/gitt_eis.py b/examples/scripts/battery_parameterisation/gitt_eis.py new file mode 100644 index 000000000..cca09f6c6 --- /dev/null +++ b/examples/scripts/battery_parameterisation/gitt_eis.py @@ -0,0 +1,162 @@ +import numpy as np +import pybamm + +import pybop + +""" +Example demonstrating coupled GITT-EIS parameterisation: a GITT experiment in which an +EIS spectrum is also acquired at the end of each pulse. + +A synthetic dataset is built in two stages: a GITT experiment is simulated to give the +time-domain voltage, then an impedance spectrum is computed about the state reached at +the end of each pulse. Both data types are stored in a single dataset on the time +domain: the voltage is recorded at every time, and the impedance variables hold the real +and imaginary components of each spectrum at the times of acquisition and zero +everywhere else. + +Diffusivity and the exchange-current density are fitted together, which is the pairing +the two measurements are meant to separate: the relaxation after each pulse constrains +transport, while the charge-transfer semicircle of the spectrum constrains kinetics. +""" + +# Define the model +model = pybamm.lithium_ion.SPM( + options={"surface form": "differential", "contact resistance": "true"}, +) +parameter_values = pybamm.ParameterValues("Chen2020") +parameter_values["Contact resistance [Ohm]"] = 0.01 + + +# The exchange-current density of Chen2020 hard-codes its prefactor, so redefine it with +# the prefactor exposed as a parameter which can then be fitted +def positive_exchange_current_density(c_e, c_s_surf, c_s_max, T): + m_ref = pybamm.Parameter( + "Positive electrode reference exchange-current density [A.m-2]" + ) + E_r = 17800 + arrhenius = pybamm.exp(E_r / pybamm.constants.R * (1 / 298.15 - 1 / T)) + return m_ref * arrhenius * c_e**0.5 * c_s_surf**0.5 * (c_s_max - c_s_surf) ** 0.5 + + +parameter_values.update( + {"Positive electrode reference exchange-current density [A.m-2]": 3.42e-6}, + check_already_exists=False, +) +parameter_values["Positive electrode exchange-current density [A.m-2]"] = ( + positive_exchange_current_density +) +parameter_values.set_initial_state(0.9) + +# Simulate a GITT experiment: repeated pulses at 1C, each followed by a rest +n_pulses = 5 +pulse_duration = 300 # s +rest_duration = 2400 # s +period = 1.0 # s +experiment = pybamm.Experiment( + [ + ( + f"Discharge at 1C for {pulse_duration} seconds", + f"Rest for {rest_duration} seconds", + ) + ] + * n_pulses, + period=f"{period} seconds", +) +gitt = pybamm.Simulation( + model, parameter_values=parameter_values, experiment=experiment +).solve() + +# Resample onto a uniform grid, avoiding the repeated times at each step change +time = np.arange(0, gitt.t[-1], period) +current = gitt["Current [A]"](time) +voltage = gitt["Voltage [V]"](time) + +# Acquire a spectrum at the end of each pulse, meaning at the end of each rest where the +# cell has relaxed and the current is zero +f_eval = np.logspace(-3, 4, 30) +impedance_variables = pybop.get_impedance_variables(f_eval) +eis_times = [ + (i + 1) * (pulse_duration + rest_duration) - period for i in range(n_pulses) +] +eis_rows = [int(np.argmin(np.abs(time - t))) for t in eis_times] + +# Assemble the dataset. The impedance variables start as a marker of which times were +# acquired, which is how the simulator learns where to compute a spectrum; the measured +# values replace the markers once they have been simulated below. +acquired = np.isin(np.arange(len(time)), eis_rows) +dataset = pybop.Dataset( + { + "Time [s]": time, + "Current [A]": current, + "Voltage [V]": voltage, + **{name: acquired.astype(float) for name in impedance_variables}, + }, + domain="Time [s]", +) + +# Simulate the impedance about the state reached at the end of each pulse. This solves +# the protocol above once and linearises the model at each of the acquisition times +sigma_v = 1e-3 # V +sigma_z = 1e-4 # Ohm +solution = pybop.pybamm.EISSimulator( + model, parameter_values=parameter_values, protocol=dataset, f_eval=f_eval +).solve() + +# Complete the synthetic dataset with the simulated spectra, adding noise to both data +# types. Only the acquired spectra carry noise; the remaining entries stay at zero +dataset["Voltage [V]"] = pybop.add_noise(voltage, sigma_v) +for name in impedance_variables: + dataset[name] = np.where( + acquired, pybop.add_noise(solution[name].data, sigma_z), 0.0 + ) + +# Save the true values +true_values = [ + parameter_values[p] + for p in [ + "Positive particle diffusivity [m2.s-1]", + "Positive electrode reference exchange-current density [A.m-2]", + ] +] + +# Fitting parameters, each searched over an order of magnitude around the true value +parameter_values.update( + { + "Positive particle diffusivity [m2.s-1]": pybop.Parameter( + pybop.Uniform(1e-16, 1e-13) + ), + "Positive electrode reference exchange-current density [A.m-2]": pybop.Parameter( + pybop.Uniform(1e-6, 1e-5) + ), + } +) + +# Build the problem. The two error measures share the dataset and the time domain, so +# they can be combined with a weight setting their relative importance. Sum-based +# measures are used because the impedance variables are zero at most times, which would +# otherwise dilute the impedance term. +simulator = pybop.pybamm.EISSimulator( + model, parameter_values=parameter_values, protocol=dataset, f_eval=f_eval +) +voltage_cost = pybop.SumSquaredError(dataset, target=["Voltage [V]"]) +impedance_cost = pybop.SumSquaredError(dataset, target=impedance_variables) +cost = pybop.WeightedCost(voltage_cost, impedance_cost, weights=[1.0, 1e3]) +problem = pybop.Problem(simulator, cost) + +# Set up the optimiser +options = pybop.PintsOptions(max_iterations=100, max_unchanged_iterations=25) +optim = pybop.XNES(problem, options=options) + +# Run the optimisation +result = optim.run() +print(result) + +# Compare identified to true parameter values +print("True parameters:", true_values) +print("Identified parameters:", result.x) + + +# Plot the optimisation result +pybop.plot.problem(problem, inputs=result.best_inputs, title="Optimised Comparison") +result.plot_convergence() +result.plot_parameters() diff --git a/pybop/__init__.py b/pybop/__init__.py index 228d1e7f2..0e8645da5 100644 --- a/pybop/__init__.py +++ b/pybop/__init__.py @@ -30,7 +30,7 @@ # # Dataset class # -from .processing.dataset import Dataset, import_pybamm_solution, import_pyprobe_result +from .processing.dataset import Dataset, get_impedance_variables, parse_impedance_variables, import_pybamm_solution, import_pyprobe_result from .processing.interpolate_current import generate_consistent_current, downsample_constant_current # diff --git a/pybop/costs/weighted_cost.py b/pybop/costs/weighted_cost.py index 87a9f88d1..b3f1446d7 100644 --- a/pybop/costs/weighted_cost.py +++ b/pybop/costs/weighted_cost.py @@ -121,8 +121,13 @@ def set_target( dataset: Dataset | None = None, ): """Set the target variable for all costs. Expecting a list of list[str] the same length as self.costs.""" - target = [target] if isinstance(target, str) else target or self._target - if isinstance(target[0], str): + if target is None: + # Keep the target of each cost, which may differ between them. Broadcasting + # self._target here would instead give every cost the union of the targets. + target = [cost.target for cost in self.costs] + elif isinstance(target, str): + target = [[target]] * len(self.costs) + elif isinstance(target[0], str): target = [target] * len(self.costs) self._target = [] diff --git a/pybop/plot/backends/base.py b/pybop/plot/backends/base.py index 4059db937..c576720cd 100644 --- a/pybop/plot/backends/base.py +++ b/pybop/plot/backends/base.py @@ -221,6 +221,21 @@ def update_plot_titles(self, figs, axes, titles, pad): """ raise NotImplementedError + @abstractmethod + def equal_aspect(self, fig, ax=None): + """ + Constrain the axes of a subplot to an equal aspect ratio, so that one unit is + the same length on both axes. Used for e.g. Nyquist plots. + + Parameters + ---------- + fig : Figure + Figure containing the axes to update. + ax : tuple, optional + Subplot location. + """ + raise NotImplementedError + @abstractmethod def update_axes_ranges(self, fig, ax, xaxis_range, yaxis_range): """ diff --git a/pybop/plot/backends/matplotlib.py b/pybop/plot/backends/matplotlib.py index 4c10ae075..619a7c432 100644 --- a/pybop/plot/backends/matplotlib.py +++ b/pybop/plot/backends/matplotlib.py @@ -309,6 +309,20 @@ def update_plot_titles(self, figures, axes, titles, max_text_width=40, pad=0): title = wrap_text(title, width=max_text_width) ax.set_title(title, pad=pad) + def equal_aspect(self, fig, ax=None): + """ + Constrain the axes of a subplot to an equal aspect ratio. + + Parameters + ---------- + fig : matplotlib.figure.Figure + Figure containing the axes to update. + ax : matplotlib.axes.Axes, optional + The axes to update. Defaults to the current axes of the figure. + """ + ax = ax if ax is not None else fig.gca() + ax.set_aspect("equal", adjustable="datalim") + def update_axes_ranges(self, fig, ax=None, xaxis_range=None, yaxis_range=None): """ Update the ranges of the axes in the provided figure. diff --git a/pybop/plot/backends/plotly.py b/pybop/plot/backends/plotly.py index 5e2df966d..3813e753c 100644 --- a/pybop/plot/backends/plotly.py +++ b/pybop/plot/backends/plotly.py @@ -387,6 +387,29 @@ def update_plot_titles(self, figures, axes, titles, max_text_width=40, pad=0): font=dict(size=14), ) + def equal_aspect(self, fig, ax=None): + """ + Constrain the axes of a subplot to an equal aspect ratio. + + Parameters + ---------- + fig : plotly.graph_objects.Figure + Figure containing the axes to update. + ax : tuple, optional + Subplot location. + """ + if ax is None: + fig.update_yaxes(scaleanchor="x", scaleratio=1) + else: + self._check_axis_input(ax) + # Each subplot anchors to its own x-axis, e.g. "x5" + fig.update_yaxes( + scaleanchor=fig.get_subplot(*ax).yaxis.anchor, + scaleratio=1, + row=ax[0], + col=ax[1], + ) + def update_axes_ranges(self, fig, ax, xaxis_range, yaxis_range): """ Update the ranges of the axes in the provided figure. diff --git a/pybop/plot/problem.py b/pybop/plot/problem.py index 61adb604e..a4367ee43 100644 --- a/pybop/plot/problem.py +++ b/pybop/plot/problem.py @@ -1,3 +1,5 @@ +import math + import numpy as np from pybop.costs.design_cost import DesignCost @@ -6,6 +8,7 @@ from pybop.plot.util import get_backend_from_figure, remove_brackets from pybop.problems.meta_problem import MetaProblem from pybop.problems.problem import Problem +from pybop.processing.dataset import parse_impedance_variables from pybop.simulators.solution import Solution @@ -77,15 +80,35 @@ def problem( model_output = problem.simulate(inputs) model_domain = target_domain[: len(model_output[target].data)] - # Create a plot for each output # Import plotting backend backend = get_backend_from_figure(backend, figures) + # Impedance variables are plotted as a Nyquist plot per acquisition, not as a series + # over the domain, so separate them from the remaining targets + frequencies, real_variables, imaginary_variables = parse_impedance_variables( + problem.target + ) + targets = [ + var + for var in problem.target + if var not in set(real_variables) | set(imaginary_variables) + ] + acquisitions = ( + _acquisition_indices(target_output, real_variables, imaginary_variables) + if len(frequencies) > 0 + else [] + ) + # Process input figures, axes, create_figure, _ = backend.parse_input_axes( - figures, axes, num_plots=len(problem.target), allow_single_axis=False + figures, + axes, + num_plots=len(targets), + allow_single_axis=False, ) - for i, var in enumerate(problem.target): + + # Create a plot for each output + for i, var in enumerate(targets): ax = axes[i % len(axes)] if create_figure: fig = backend.create_figure( @@ -146,5 +169,83 @@ def problem( if show: backend.show_figure(fig) + # Collect the acquired spectra into a single figure of Nyquist subplots + if len(acquisitions) > 0: + num_cols = int(math.ceil(math.sqrt(len(acquisitions)))) + num_rows = int(math.ceil(len(acquisitions) / num_cols)) + impedance_figure, impedance_axes = backend.make_subplots( + num_rows=num_rows, + num_cols=num_cols, + num_plots=len(acquisitions), + title=title, + style={ + "bg_color": "white", + "width": 400 * num_cols, + "height": 400 * num_rows, + }, + ) + backend.update_axes_titles( + impedance_figure, + impedance_axes, + r"$Z_{re} / \Omega$", + r"$-Z_{im} / \Omega$", + ) + backend.update_plot_titles( + impedance_figure, + impedance_axes, + [ + f"{remove_brackets(domain)} = {target_domain[row]:g}" + for row in acquisitions + ], + ) + + # Fixed styles, matching plot.nyquist, so that the colours mean the same thing + # in every subplot rather than following the shared colour cycle + impedance_styles = ( + ( + target_output, + "Reference", + { + "linestyle": "none", + "marker": "o", + "fillstyle": "none", + "markeredgecolor": "#636EFA", + "color": "#636EFA", + }, + ), + (model_output, "Model", {"linestyle": "solid", "color": "#00CC96"}), + ) + + for ax, row in zip(impedance_axes, acquisitions, strict=True): + for output, label, style in impedance_styles: + backend.plot_trace( + backend.line( + x=[output[name].data[row] for name in real_variables], + y=[-output[name].data[row] for name in imaginary_variables], + label=label, + style=style, + ), + impedance_figure, + ax=ax, + ) + + # A Nyquist plot is only readable with an equal aspect ratio + backend.equal_aspect(impedance_figure, ax=ax) + backend.legend(impedance_figure, axes=ax) + + figures = np.append(figures, impedance_figure) + if show: + backend.show_figure(impedance_figure) + if not show: return figures[0] if len(figures) == 1 else figures + + +def _acquisition_indices( + target_output, real_variables: list[str], imaginary_variables: list[str] +) -> np.ndarray: + """Return the indices of the domain points at which a spectrum was acquired.""" + measured = np.asarray( + [target_output[name].data for name in (*real_variables, *imaginary_variables)] + ) + return np.flatnonzero(np.any(measured != 0.0, axis=0)) diff --git a/pybop/processing/dataset.py b/pybop/processing/dataset.py index f2b2b45cd..75041bdd9 100644 --- a/pybop/processing/dataset.py +++ b/pybop/processing/dataset.py @@ -1,3 +1,4 @@ +import re import warnings from typing import Protocol @@ -66,6 +67,16 @@ def __getitem__(self, key): return self.data[key] + def keys(self): + """ + Return the keys of the data dictionary. + + Returns + ------- + dict_keys + """ + return self.data.keys() + def __len__(self) -> int: """Return the length of the data, based on the length of the domain data.""" return len(self.data[self.domain]) @@ -131,7 +142,7 @@ def _check_data_consistency( ) -> None: n_domain_data = len(domain_data) for s in signals: - if len(self.data[s]) != n_domain_data: + if np.shape(self.data[s])[-1] != n_domain_data: raise ValueError( f"{self.domain} data and {s} data must be the same length." ) @@ -180,6 +191,67 @@ def get_interpolant(self, control: str = "Current [A]") -> Interpolant: return Interpolant(self.data["Time [s]"], self.data[control], pybamm_t) +def get_impedance_variables(frequencies: np.ndarray | list[float]) -> list[str]: + """ + Return the names of the dataset variables which hold an impedance spectrum, being + two real-valued variables, the real and the imaginary component, per frequency. + + Complex data is split into real components because the error measures square the + residual, and `r**2` is not `abs(r)**2` for a complex array. + + Parameters + ---------- + frequencies : np.ndarray or list[float] + The frequencies at which the impedance is evaluated. + + Returns + ------- + list[str] + Variable names, ordered (real, imaginary) for each frequency in turn. + """ + return [ + f"Impedance {component} [Ohm] ({frequency:.6g} Hz)" + for frequency in frequencies + for component in ("real", "imaginary") + ] + + +def parse_impedance_variables( + variables: list[str], +) -> tuple[np.ndarray, list[str], list[str]]: + """ + Pick out the impedance variables from a list of variable names, inverting + :func:`get_impedance_variables`. + + Parameters + ---------- + variables : list[str] + Variable names, which may or may not include impedance variables. + + Returns + ------- + tuple[np.ndarray, list[str], list[str]] + The frequencies in increasing order, and the corresponding real and imaginary + variable names. All three are empty if no impedance variables are present. + """ + pattern = re.compile(r"^Impedance (real|imaginary) \[Ohm\] \((\S+) Hz\)$") + + frequencies = {} + for variable in variables: + match = pattern.match(variable) + if match: + component, frequency = match.group(1), float(match.group(2)) + frequencies.setdefault(frequency, {})[component] = variable + + # Ignore any frequency which is missing one of its two components + ordered = sorted(f for f, components in frequencies.items() if len(components) == 2) + return ( + np.asarray(ordered), + [frequencies[f]["real"] for f in ordered], + [frequencies[f]["imaginary"] for f in ordered], + ) + + def import_pybamm_solution( solution: Solution, variables: list[str] | None = None, diff --git a/pybop/pybamm/eis_simulator.py b/pybop/pybamm/eis_simulator.py index 8445b0319..03dfa7e95 100644 --- a/pybop/pybamm/eis_simulator.py +++ b/pybop/pybamm/eis_simulator.py @@ -10,6 +10,7 @@ if TYPE_CHECKING: from pybop.parameters.parameter import Inputs +from pybop.processing.dataset import Dataset, get_impedance_variables from pybop.pybamm.simulator import Simulator from pybop.pybamm.utils import SymbolReplacer from pybop.simulators.base_simulator import BaseSimulator, Solution @@ -38,6 +39,13 @@ class EISSimulator(BaseSimulator): The frequencies at which to evaluate the impedance. parameter_values : pybamm.ParameterValues, optional The parameter values to be used in the model. + protocol : pybop.Dataset, optional + A dataset defining a time-domain protocol, containing the domain data, a control + variable (e.g. "Current [A]") and the impedance variables given + by `pybop.get_impedance_variables(f_eval)`. These are non-zero at the times at + which a spectrum was measured, and zero elsewhere; the simulator computes a + spectrum at exactly those times, coupling the time-domain simulation to the + EIS. If None, a single spectrum is computed about the initial state. initial_state : dict, optional A valid initial state, e.g. `"Initial open-circuit voltage [V]"` or ``"Initial SoC"`. Defaults to None, indicating that the existing initial state of charge (for an ECM) @@ -68,6 +76,7 @@ def __init__( model: pybamm.BaseModel, f_eval: np.ndarray | list[float], parameter_values: pybamm.ParameterValues | None = None, + protocol: Dataset | None = None, initial_state: float | str | None = None, solver: pybamm.BaseSolver | None = None, geometry: pybamm.Geometry | None = None, @@ -86,10 +95,17 @@ def __init__( super().__init__(parameters=parameter_values) - # Set up a simulation + # Locate the times at which to compute a spectrum, if any + self._impedance_variables = get_impedance_variables(f_eval) + self._acquisition_indices = self._locate_acquisition_times(protocol) + + # Set up a simulation. When a protocol is given, the Simulator installs the + # control interpolant and sets t_interp to the domain data, so the entries of + # the solution align with the rows of the dataset. self._simulator = Simulator( model, parameter_values=parameter_values, + protocol=protocol, initial_state=initial_state, solver=solver, geometry=geometry, @@ -104,13 +120,81 @@ def __init__( # Initialise self.M = None - self._jac = None self.b = None v_scale = getattr(model.variables["Voltage [V]"], "scale", 1) i_scale = getattr(model.variables["Current [A]"], "scale", 1) self.z_scale = self.parameter_values.evaluate(v_scale / i_scale) + def _locate_acquisition_times(self, protocol: Dataset | None) -> np.ndarray | None: + """ + Return the indices of the rows at which a spectrum was measured, or None for a + stationary simulation. + """ + if protocol is None: + return None + + missing = set(self._impedance_variables) - set(protocol.keys()) + if missing: + raise ValueError( + "The protocol dataset is missing impedance variables, e.g. " + f"'{sorted(missing)[0]}'. Name them with pybop.get_impedance_variables(f_eval), " + "using zeros at the times where no spectrum was measured." + ) + + measured = np.asarray([protocol[name] for name in self._impedance_variables]) + indices = np.flatnonzero(np.any(measured != 0.0, axis=0)) + if indices.size == 0: + raise ValueError( + "The impedance variables are zero everywhere, so there is nothing to fit. " + "Set them to the measured spectra at the times of acquisition." + ) + return indices + + def set_output_variables(self, target: list[str]): + """ + Deliberately a no-op. Restricting the solver to a list of output variables stops + PyBaMM from returning the state vector, which is required to linearise the model + about the state at each acquisition time. + """ + return None + + def _set_up_matrices(self, inputs: "Inputs") -> None: + """ + Set up the solver and the parts of the linear system which do not depend on the + operating point: the mass matrix and the forcing vector. Called once, unless the + model has to be rebuilt. + """ + built_model = self.simulation.built_model + self.simulation.solver.set_up(built_model, inputs=inputs) + + self.M = csc_matrix(built_model.mass_matrix.entries) + + # Add forcing to the RHS on the current density + self.b = np.zeros((self.M.shape[0], 1)) + self.b[-1] = -1 + + def _jacobian(self, t: float, y: np.ndarray, inputs: "Inputs") -> csc_matrix: + """Evaluate the Jacobian of the built model about the state y at time t.""" + built_model = self.simulation.built_model + + # PyBaMM orders the input parameters of the built model alphabetically, so the + # values must be stacked in that order and not in the order given by the caller + casadi_inputs = ( + casadi.vertcat(*[inputs[name] for name in sorted(inputs)]) + if inputs is not None and built_model.convert_to_format == "casadi" + else inputs or [] + ) + jac = built_model.jac_rhs_algebraic_eval(t, y, casadi_inputs).sparse() + return csc_matrix(jac) + + def _calculate_spectrum(self, jac: csc_matrix) -> np.ndarray: + """Compute the impedance at every frequency for a given Jacobian.""" + return ( + np.asarray([self.calculate_impedance(f, jac) for f in self._f_eval]) + * self.z_scale + ) + def set_up_for_eis(self, model: pybamm.BaseModel) -> pybamm.BaseModel: """ Set up the model for electrochemical impedance spectroscopy (EIS) simulations. @@ -125,12 +209,12 @@ def set_up_for_eis(self, model: pybamm.BaseModel) -> pybamm.BaseModel: Returns ------- pybamm.BaseModel - The modified model ready for EIS simulations. + A modified copy of the model, ready for EIS simulations. Raises ------ ValueError - If the model is missing required variables. + If the model is missing required variables or options. """ # Verify model has required variables required_vars = ["Voltage [V]", "Current [A]"] @@ -140,6 +224,18 @@ def set_up_for_eis(self, model: pybamm.BaseModel) -> pybamm.BaseModel: f"Model must contain variable '{var}' for EIS simulation" ) + # Without a surface form, the double layer is absent from the model and the + # computed impedance is silently meaningless + surface_form = model.options.get("surface form", "false") + if surface_form != "differential": + raise ValueError( + "EIS simulation requires the 'surface form' model option to be " + f"'differential', got '{surface_form}'." + ) + + # Work on a copy, so that the model given by the user is left untouched + model = model.new_copy() + V_cell = pybamm.Variable("Voltage variable [V]") model.variables["Voltage variable [V]"] = V_cell V = model.variables["Voltage [V]"] @@ -182,46 +278,18 @@ def set_up_for_eis(self, model: pybamm.BaseModel) -> pybamm.BaseModel: return model def _model_rebuild(self, inputs: "Inputs") -> None: - """Update the parameter values and rebuild the EIS model.""" + """ + Rebuild the EIS model if required, and set up the operating-point-independent + matrices. Mirroring the Simulator, the model and these matrices are set up once + unless a rebuild is required on every evaluation. + """ if self._simulator.requires_model_rebuild: self.parameter_values.update(inputs) self._simulator.create_simulation() self.simulation.build(initial_soc=self._simulator.initial_state) - self._initialise_eis_matrices(inputs=inputs) - - def _initialise_eis_matrices(self, inputs: "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 the necessary attributes from the model, and prepares matrices for the simulation. - - Raises - ------ - RuntimeError - If the model hasn't been built yet. - """ - built_model = self.simulation.built_model - M = built_model.mass_matrix.entries - self.simulation.solver.set_up(built_model, inputs=inputs) - - # Convert inputs to casadi format if needed - casadi_inputs = ( - casadi.vertcat(*inputs.values()) - if inputs is not None and built_model.convert_to_format == "casadi" - else inputs or [] - ) - - # Extract the necessary attributes from the model - y0 = built_model.concatenated_initial_conditions.evaluate(0, inputs=inputs) - jac = built_model.jac_rhs_algebraic_eval(0, y0, casadi_inputs).sparse() - - # Convert to Compressed Sparse Column format - self.M = csc_matrix(M) - self._jac = csc_matrix(jac) - - # Add forcing to the RHS on the current density - self.b = np.zeros(y0.shape) - self.b[-1] = -1 + self._set_up_matrices(inputs=inputs) + elif self.M is None: + self._set_up_matrices(inputs=inputs) def solve( self, @@ -289,7 +357,9 @@ def _catch_errors(self, inputs: "list[Inputs]") -> list[Solution | FailedSolutio try: simulations.append(self._solve(x)) except (ZeroDivisionError, RuntimeError, ValueError): - simulations.append(FailedSolution(["Impedance"], x.keys())) + simulations.append( + FailedSolution(self.solution_variables, x.keys()) + ) return simulations simulations = [] @@ -301,29 +371,78 @@ def _solve(self, inputs: "Inputs") -> Solution: """ Run the EIS simulation to calculate impedance at all specified frequencies. + For a stationary simulation, one spectrum is computed about the initial state. + When coupled to a protocol, the time-domain trajectory is solved once and each + spectrum is computed by linearising about the state at the requested time. + Parameters ---------- inputs : Inputs Input parameters. - calculate_sensitivities : bool - Whether to calculate sensitivities (default: False). - Currently not implemented for EIS. Returns ------- Solution - Complex impedance results. + Complex impedance results, or the voltage and the real and imaginary + impedance components over the time domain when coupled to a protocol. """ - # Always run initialise_eis_matrices, after rebuilding the model if necessary + # Rebuild the model only if necessary, then set up the constant matrices self._model_rebuild(inputs) - zs = [self.calculate_impedance(frequency) for frequency in self._f_eval] + if self._acquisition_indices is None: + y0 = self.simulation.built_model.concatenated_initial_conditions.evaluate( + 0, inputs=inputs + ) + solution = Solution() + solution.set_solution_variable( + "Impedance", + data=self._calculate_spectrum(self._jacobian(0, y0, inputs)), + ) + return solution + + return self._solve_along_protocol(inputs) + + def _solve_along_protocol(self, inputs: "Inputs") -> Solution: + """ + Solve the time-domain protocol once, then compute a spectrum about the state at + each of the requested times. + """ + sim_solution = self._simulator.solve(inputs) + if isinstance(sim_solution, FailedSolution): + raise ValueError("The time-domain simulation failed.") + + t, y = sim_solution.t, sim_solution.y + if self._acquisition_indices[-1] >= len(t): + raise ValueError( + "The time-domain simulation terminated before the last EIS time." + ) solution = Solution() - solution.set_solution_variable("Impedance", data=np.asarray(zs) * self.z_scale) + solution.set_solution_variable("Time [s]", data=t) + solution.set_solution_variable( + "Voltage [V]", data=sim_solution["Voltage [V]"].data + ) + + # Zero away from the times of acquisition, matching the dataset convention + impedance = {name: np.zeros(len(t)) for name in self._impedance_variables} + for i in self._acquisition_indices: + y_i = np.asarray(y[:, i]).reshape(-1, 1) + if np.abs(y_i[-1]) > 1e-10: + warnings.warn( + f"The current is not zero at the acquisition time t={t[i]} s, " + "so the impedance is linearised about a non-zero operating point.", + stacklevel=2, + ) + zs = self._calculate_spectrum(self._jacobian(t[i], y_i, inputs)) + for j, z in enumerate(zs): + impedance[self._impedance_variables[2 * j]][i] = z.real + impedance[self._impedance_variables[2 * j + 1]][i] = z.imag + + for name, data in impedance.items(): + solution.set_solution_variable(name, data=data) return solution - def calculate_impedance(self, frequency): + def calculate_impedance(self, frequency: float, jac: csc_matrix) -> complex: """ Calculate the impedance for a given frequency. @@ -334,6 +453,8 @@ def calculate_impedance(self, frequency): ---------- frequency : float The frequency at which to calculate the impedance in Hz. + jac : csc_matrix + The Jacobian of the built model about the operating point. Returns ------- @@ -342,7 +463,7 @@ def calculate_impedance(self, frequency): """ # Compute the system matrix - A = 1.0j * 2 * np.pi * frequency * self.M - self._jac + A = 1.0j * 2 * np.pi * frequency * self.M - jac # Solve the system x = spsolve(A, self.b) @@ -362,6 +483,13 @@ def parameter_values(self): def input_parameter_names(self): return self._simulator.input_parameter_names + @property + def solution_variables(self) -> list[str]: + """The names of the variables set by a solve.""" + if self._acquisition_indices is None: + return ["Impedance"] + return ["Time [s]", "Voltage [V]", *self._impedance_variables] + @property def has_sensitivities(self): return False diff --git a/tests/integration/models/test_grouped_models.py b/tests/integration/models/test_grouped_models.py index deda9f635..f655b8813 100644 --- a/tests/integration/models/test_grouped_models.py +++ b/tests/integration/models/test_grouped_models.py @@ -139,6 +139,9 @@ def test_voltage_fitting(self, dataset, model_config, parameters): ) def test_eis_fitting(self, eis_dataset, model_config, parameters): + if model_config["model"].options.get("surface form") != "differential": + pytest.skip("EIS simulation requires a differential surface form") + parameter_values = model_config["parameter_values"] parameter_values.update(parameters) simulator = pybop.pybamm.EISSimulator( diff --git a/tests/unit/test_cost.py b/tests/unit/test_cost.py index 443881c76..72015c6d3 100644 --- a/tests/unit/test_cost.py +++ b/tests/unit/test_cost.py @@ -335,6 +335,18 @@ def noisy_problem(self, ground_truth, parameters, experiment): model, parameter_values=parameter_values, protocol=noisy_dataset ) + def test_weighted_cost_keeps_distinct_targets(self, noisy_problem): + """Building a Problem must not give every cost the union of the targets.""" + dataset, simulator = noisy_problem + cost1 = pybop.SumSquaredError(dataset, target=["Voltage [V]"]) + cost2 = pybop.SumSquaredError(dataset, target=["Current [A]"]) + weighted_cost = pybop.WeightedCost(cost1, cost2) + + pybop.Problem(simulator, weighted_cost) + + assert cost1.target == ["Voltage [V]"] + assert cost2.target == ["Current [A]"] + def test_weighted_fitting_cost(self, noisy_problem, parameters, dataset): dataset, simulator = noisy_problem cost1 = pybop.SumSquaredError(dataset) diff --git a/tests/unit/test_problem.py b/tests/unit/test_problem.py index 7a9eb5987..574029d73 100644 --- a/tests/unit/test_problem.py +++ b/tests/unit/test_problem.py @@ -117,7 +117,7 @@ def test_fitting_problem(self, simulator, dataset): assert not np.isfinite(out["Voltage [V]"].data) def test_fitting_problem_eis(self, parameters): - model = pybamm.lithium_ion.SPM() + model = pybamm.lithium_ion.SPM(options={"surface form": "differential"}) dataset = pybop.Dataset( { "Frequency [Hz]": np.logspace(-4, 5, 30), diff --git a/tests/unit/test_pybamm_utils.py b/tests/unit/test_pybamm_utils.py index ecf8df674..7908b0eda 100644 --- a/tests/unit/test_pybamm_utils.py +++ b/tests/unit/test_pybamm_utils.py @@ -25,7 +25,8 @@ class TestPybammUtils: def test_simulate_procedure(self, tmp_path): import pyprobe - model = pybamm.lithium_ion.SPM() + # The procedures include an EIS sweep, which requires a surface form + model = pybamm.lithium_ion.SPM(options={"surface form": "differential"}) full_cell_parameters = pybamm.ParameterValues("Chen2020") cell_info = { "Cell type": "LG M50 Synthetic", diff --git a/tests/unit/test_simulator.py b/tests/unit/test_simulator.py index 748617d17..df71f9e35 100644 --- a/tests/unit/test_simulator.py +++ b/tests/unit/test_simulator.py @@ -1,5 +1,6 @@ from copy import copy +import numpy as np import pybamm import pytest @@ -78,3 +79,188 @@ def test_set_output_variables(self): ValueError, match="Not a variable is not a variable in the model." ): simulator.set_output_variables(["Not a variable"]) + + +class TestCoupledEISSimulator: + """ + A class to test the pybamm.EISSimulator class when coupled to a protocol. + """ + + pytestmark = pytest.mark.unit + + @pytest.fixture + def model(self): + return pybamm.lithium_ion.SPM(options={"surface form": "differential"}) + + @pytest.fixture + def parameter_values(self): + parameter_values = pybamm.ParameterValues("Chen2020") + parameter_values.set_initial_state(0.9) + return parameter_values + + @pytest.fixture + def frequencies(self): + return np.logspace(-1, 3, 4) + + @pytest.fixture + def impedance_variables(self, frequencies): + return pybop.get_impedance_variables(frequencies) + + @pytest.fixture + def acquisitions(self): + """The rows at which a spectrum is acquired: one at the start, one at rest.""" + return [0, 40] + + @pytest.fixture + def dataset(self, impedance_variables, acquisitions): + time = np.arange(0, 601, 10.0) + data = { + "Time [s]": time, + "Current [A]": np.zeros_like(time), + "Voltage [V]": np.zeros_like(time), + } + for name in impedance_variables: + variable = np.zeros(len(time)) + variable[acquisitions] = 1.0 + data[name] = variable + + return pybop.Dataset(data, domain="Time [s]") + + def test_output_shape_and_zero_padding( + self, + model, + parameter_values, + dataset, + frequencies, + impedance_variables, + acquisitions, + ): + simulator = pybop.pybamm.EISSimulator( + model, + parameter_values=parameter_values, + protocol=dataset, + f_eval=frequencies, + ) + solution = simulator.solve() + + n_time = len(dataset["Time [s]"]) + assert len(solution["Voltage [V]"].data) == n_time + for name in impedance_variables: + data = solution[name].data + assert len(data) == n_time + # Non-zero only at the times of acquisition + np.testing.assert_array_equal(np.flatnonzero(data), acquisitions) + + def test_matches_stationary_at_initial_state( + self, model, parameter_values, dataset, frequencies, impedance_variables + ): + """The spectrum at t=0 must match a stationary simulation of the same state.""" + coupled = pybop.pybamm.EISSimulator( + model, + parameter_values=parameter_values, + protocol=dataset, + f_eval=frequencies, + ).solve() + stationary = pybop.pybamm.EISSimulator( + model, parameter_values=parameter_values, f_eval=frequencies + ).solve() + + impedance = np.asarray( + [ + coupled[impedance_variables[2 * j]].data[0] + + 1j * coupled[impedance_variables[2 * j + 1]].data[0] + for j in range(len(frequencies)) + ] + ) + np.testing.assert_allclose(impedance, stationary["Impedance"].data, rtol=1e-10) + + def test_builds_once(self, model, parameter_values, dataset, frequencies): + """The model and the constant matrices are set up once, not per evaluation.""" + parameter_values["Negative electrode active material volume fraction"] = ( + pybop.Parameter(pybop.Uniform(0.4, 0.75)) + ) + simulator = pybop.pybamm.EISSimulator( + model, + parameter_values=parameter_values, + protocol=dataset, + f_eval=frequencies, + ) + + calls = {"create": 0, "set_up": 0} + original_create = simulator._simulator.create_simulation + original_set_up = simulator._set_up_matrices + + def counted_create(*args, **kwargs): + calls["create"] += 1 + return original_create(*args, **kwargs) + + def counted_set_up(*args, **kwargs): + calls["set_up"] += 1 + return original_set_up(*args, **kwargs) + + simulator._simulator.create_simulation = counted_create + simulator._set_up_matrices = counted_set_up + + inputs = {"Negative electrode active material volume fraction": 0.6} + for _ in range(3): + simulator.solve(inputs) + + assert calls["create"] == 0 # built during construction + assert calls["set_up"] == 1 + + def test_model_is_not_modified(self, model, parameter_values, dataset, frequencies): + """Setting up for EIS must copy the model, not modify the caller's.""" + n_algebraic = len(model.algebraic) + pybop.pybamm.EISSimulator( + model, + parameter_values=parameter_values, + protocol=dataset, + f_eval=frequencies, + ) + assert len(model.algebraic) == n_algebraic + + def test_surface_form_required(self, parameter_values, dataset, frequencies): + with pytest.raises(ValueError, match="surface form"): + pybop.pybamm.EISSimulator( + pybamm.lithium_ion.SPM(), + parameter_values=parameter_values, + protocol=dataset, + f_eval=frequencies, + ) + + def test_impedance_variables_round_trip(self, frequencies, impedance_variables): + # Variables may reach the parser in any order, e.g. via a set + parsed, real, imaginary = pybop.parse_impedance_variables( + ["Voltage [V]", *reversed(impedance_variables)] + ) + # The names carry six significant figures, which bounds the round trip + np.testing.assert_allclose(parsed, frequencies, rtol=1e-5) + assert real == impedance_variables[::2] + assert imaginary == impedance_variables[1::2] + + # No impedance variables present + parsed, real, imaginary = pybop.parse_impedance_variables( + ["Voltage [V]", "Current [A]"] + ) + assert len(parsed) == 0 and real == [] and imaginary == [] + + def test_dataset_errors( + self, model, parameter_values, dataset, frequencies, impedance_variables + ): + with pytest.raises(ValueError, match="missing impedance variables"): + pybop.pybamm.EISSimulator( + model, + parameter_values=parameter_values, + protocol=dataset, + f_eval=np.append(frequencies, 1e4), + ) + + for name in impedance_variables: + dataset[name] = np.zeros(len(dataset["Time [s]"])) + with pytest.raises(ValueError, match="zero everywhere"): + pybop.pybamm.EISSimulator( + model, + parameter_values=parameter_values, + protocol=dataset, + f_eval=frequencies, + )