From 9bd42e8f81238401248f91bc5da058610de576ef Mon Sep 17 00:00:00 2001
From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com>
Date: Mon, 24 Mar 2025 10:33:18 +0000
Subject: [PATCH 01/20] Enable non-stationary EIS

---
 pybop/models/base_model.py | 83 +++++++++++++++++++++++++++++---------
 1 file changed, 65 insertions(+), 18 deletions(-)

diff --git a/pybop/models/base_model.py b/pybop/models/base_model.py
index 3dc87adc5..0b90e151a 100644
--- a/pybop/models/base_model.py
+++ b/pybop/models/base_model.py
@@ -165,17 +165,16 @@ def build(
         if not self.pybamm_model._built:  # noqa: SLF001
             self.pybamm_model.build_model()
 
+        if dataset is not None:
+            self.set_current_function(dataset)
+
         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:
@@ -329,7 +328,9 @@ def set_up_for_eis(self, model):
             "Current function [A]", {"Time [s]": pybamm.t}
         )
         model.algebraic[I_cell] = I - I_applied
-        model.initial_conditions[I_cell] = 0
+        model.initial_conditions[I_cell] = self._parameter_set.evaluate(
+            I_applied, {"Time [s]": 0}
+        )
 
     def clear(self):
         """
@@ -491,7 +492,11 @@ def simulate(
         return self._pybamm_solution
 
     def simulateEIS(
-        self, inputs: Inputs, f_eval: list, initial_state: Optional[dict] = None
+        self,
+        inputs: Inputs,
+        f_eval: list,
+        initial_state: Optional[dict] = None,
+        t_eval=None,
     ) -> dict[str, np.ndarray]:
         """
         Compute the forward model simulation with electrochemical impedance spectroscopy
@@ -504,10 +509,13 @@ def simulateEIS(
             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.
+        t_eval : array-like, optional
+            An array of time points at which to simulate operando EIS. Defaults to None,
+            indicating that stationary EIS should be simulated at time t=0, with I=0.
 
         Returns
         -------
-        array-like
+        dict
             The simulation result corresponding to the specified signal.
 
         Raises
@@ -517,6 +525,10 @@ def simulateEIS(
         """
         inputs = self.parameters.verify(inputs)
 
+        # Perform stationary EIS by default
+        if t_eval is None:
+            self._parameter_set["Current function [A]"] = 0
+
         # Build or rebuild if required
         self.build(inputs=inputs, initial_state=initial_state)
 
@@ -526,12 +538,26 @@ def simulateEIS(
         ):
             raise ValueError("These parameter values are infeasible.")
 
-        self.initialise_eis_simulation(inputs)
-        zs = [self.calculate_impedance(frequency) for frequency in f_eval]
+        self.initialise_eis_simulation(inputs, t_eval=t_eval)
 
-        return {"Impedance": np.asarray(zs) * self.z_scale}
+        if t_eval is None:
+            ## Stationary EIS
+            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):
+        else:
+            ## Operando EIS
+            zs_at_t_eval = []
+            for i in range(len(t_eval)):
+                self.J = self.J_at_t_eval[i]
+                zs = [self.calculate_impedance(frequency) for frequency in f_eval]
+                zs_at_t_eval.append(zs)
+            return {
+                "Time [s]": np.asarray(t_eval),
+                "Impedance": np.asarray(zs_at_t_eval) * self.z_scale,
+            }
+
+    def initialise_eis_simulation(self, inputs: Optional[Inputs] = None, t_eval=None):
         """
         Initialise the Electrochemical Impedance Spectroscopy (EIS) simulation.
 
@@ -540,26 +566,28 @@ def initialise_eis_simulation(self, inputs: Optional[Inputs] = None):
 
         Parameters
         ----------
-        inputs : dict (optional)
+        inputs : dict, optional
             The input parameters for the simulation.
+        t_eval : array-like, optional
+            An array of time points at which to simulate operando EIS. Defaults to None,
+            indicating that EIS should be simulated at the initial time point (t=0).
         """
         # 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 = (
+        self._casadi_inputs = (
             casadi.vertcat(*inputs.values())
             if inputs is not None and self._built_model.convert_to_format == "casadi"
             else inputs or []
         )
 
+        ## Stationary EIS
         # Extract necessary attributes from the model
-        self.y0 = self._built_model.concatenated_initial_conditions.evaluate(
-            0, inputs=inputs
-        )
+        y = self._built_model.concatenated_initial_conditions.evaluate(0, inputs=inputs)
         self.J = self._built_model.jac_rhs_algebraic_eval(
-            0, self.y0, casadi_inputs
+            0, y, self._casadi_inputs
         ).sparse()
 
         # Convert to Compressed Sparse Column format
@@ -567,9 +595,28 @@ def initialise_eis_simulation(self, inputs: Optional[Inputs] = None):
         self.J = csc_matrix(self.J)
 
         # Add forcing to the RHS on the current density
-        self.b = np.zeros(self.y0.shape)
+        self.b = np.zeros(y.shape)
         self.b[-1] = -1
 
+        ## Operando EIS
+        if t_eval is not None:
+            # Initial state
+            state = self.get_state(inputs or {}, 0, y)
+
+            self.J_at_t_eval = [self.J]
+            for t in t_eval[1:]:
+                # Step forwards in time
+                state = self.step(state, t)
+
+                # Extract necessary attributes from the model
+                y = state.as_ndarray()
+                J = self._built_model.jac_rhs_algebraic_eval(
+                    t, y, self._casadi_inputs
+                ).sparse()
+
+                # Convert to Compressed Sparse Column format
+                self.J_at_t_eval.append(csc_matrix(J))
+
     def calculate_impedance(self, frequency):
         """
         Calculate the impedance for a given frequency.

From 050fc2ab9875b51634aadd2965f8f6d66af10b51 Mon Sep 17 00:00:00 2001
From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com>
Date: Mon, 24 Mar 2025 10:58:18 +0000
Subject: [PATCH 02/20] Update current when given time

---
 pybop/models/base_model.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pybop/models/base_model.py b/pybop/models/base_model.py
index 0b90e151a..f9c2cd054 100644
--- a/pybop/models/base_model.py
+++ b/pybop/models/base_model.py
@@ -165,7 +165,7 @@ def build(
         if not self.pybamm_model._built:  # noqa: SLF001
             self.pybamm_model.build_model()
 
-        if dataset is not None:
+        if dataset is not None and "Time [s]" in dataset.keys():
             self.set_current_function(dataset)
 
         if self.eis:

From 42ccc791b7bd89cf9e45d2eba92c01e165b73caf Mon Sep 17 00:00:00 2001
From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com>
Date: Mon, 24 Mar 2025 13:35:57 +0000
Subject: [PATCH 03/20] Add dataset.keys

---
 pybop/_dataset.py | 10 ++++++++++
 1 file changed, 10 insertions(+)

diff --git a/pybop/_dataset.py b/pybop/_dataset.py
index 35fc9c611..f2b0e4634 100644
--- a/pybop/_dataset.py
+++ b/pybop/_dataset.py
@@ -80,6 +80,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 check(self, domain: str = None, signal: Union[str, list[str]] = None) -> bool:
         """
         Check the consistency of a PyBOP Dataset against the expected format.

From 1f88c8d0302455eaa8df7b448ed164645da48740 Mon Sep 17 00:00:00 2001
From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com>
Date: Fri, 28 Nov 2025 21:30:27 +0000
Subject: [PATCH 04/20] Move changes to EISSimulator

---
 pybop/pybamm/eis_simulator.py | 131 +++++++++++++++++++++++++++++-----
 pybop/pybamm/simulator.py     |   4 ++
 2 files changed, 119 insertions(+), 16 deletions(-)

diff --git a/pybop/pybamm/eis_simulator.py b/pybop/pybamm/eis_simulator.py
index fc58211ec..b980935bc 100644
--- a/pybop/pybamm/eis_simulator.py
+++ b/pybop/pybamm/eis_simulator.py
@@ -1,5 +1,6 @@
 import warnings
 from copy import copy
+from dataclasses import dataclass
 from typing import TYPE_CHECKING
 
 import casadi
@@ -10,12 +11,37 @@
 
 if TYPE_CHECKING:
     from pybop.parameters.parameter import Inputs
+from pybop._dataset import Dataset
 from pybop._utils import FailedSolution, SymbolReplacer
 from pybop.parameters.parameter import Parameter, Parameters
 from pybop.pybamm.simulator import Simulator
 from pybop.simulators.base_simulator import BaseSimulator, Solution
 
 
+@dataclass
+class TimeSeriesState:
+    """
+    The current state of a time series model that is a PyBaMM model.
+    """
+
+    sol: pybamm.Solution
+    inputs: "Inputs"
+    t: float = 0.0
+
+    def as_ndarray(self) -> np.ndarray:
+        ncol = self.sol.y.shape[1]
+        if ncol > 1:
+            y = self.sol.y[:, -1]
+        else:
+            y = self.sol.y
+        if isinstance(y, casadi.DM):
+            y = y.full()
+        return y
+
+    def __len__(self):
+        return self.sol.y.shape[0]
+
+
 class EISSimulator(BaseSimulator):
     """
     A class to extend a PyBaMM model for EIS, automatically build/rebuild a pybamm.Simulation to obtain
@@ -42,6 +68,9 @@ class EISSimulator(BaseSimulator):
         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)
         or initial concentrations (for an EChem model) will be used.
+    protocol : Dataset | np.ndarray, optional
+        A 1D array of values or dataset containing the time points at which to simulate
+        operando EIS. Defaults to None, corresponding to stationary EIS at time t=0, with I=0.
     solver : pybamm.BaseSolver, optional
         The solver to simulate the composed Simulator. If None, uses `pybop.RecommendedSolver`.
     geometry : pybamm.Geometry, optional
@@ -66,6 +95,7 @@ def __init__(
         f_eval: np.ndarray | list[float],
         parameter_values: pybamm.ParameterValues | None = None,
         initial_state: float | str | None = None,
+        protocol: Dataset | np.ndarray | None = None,
         solver: pybamm.BaseSolver | None = None,
         geometry: pybamm.Geometry | None = None,
         submesh_types: dict | None = None,
@@ -76,9 +106,21 @@ def __init__(
     ):
         # Set-up model for EIS
         self._f_eval = f_eval
-        model = self.set_up_for_eis(model)
         parameter_values = parameter_values or model.default_parameter_values
-        parameter_values["Current function [A]"] = 0
+        if protocol is None:  # perform stationary EIS by default
+            parameter_values["Current function [A]"] = 0
+            initial_current = 0
+        elif isinstance(protocol, pybamm.Experiment):
+            raise ValueError("EISSimulator cannot simulate a pybamm.Experiment.")
+        elif (
+            isinstance(protocol, Dataset)
+            and "Current function [A]" in protocol.data.keys()
+        ):
+            parameter_values["Current function [A]"] = pybamm.Interpolant(
+                protocol["Time [s]"], protocol["Current function [A]"], pybamm.t
+            )
+            initial_current = protocol["Current function [A]"][0]
+        model = self.set_up_for_eis(model, initial_current=float(initial_current))
 
         # Unpack the uncertain parameters from the parameter values
         parameters = Parameters()
@@ -92,6 +134,7 @@ def __init__(
             model,
             parameter_values=parameter_values,
             initial_state=initial_state,
+            protocol=protocol,
             solver=solver,
             geometry=geometry,
             submesh_types=submesh_types,
@@ -104,7 +147,7 @@ def __init__(
         self.debug_mode = False
 
         # Initialise
-        self.M = None
+        self._mass = None
         self._jac = None
         self.b = None
 
@@ -112,7 +155,9 @@ def __init__(
         i_scale = getattr(model.variables["Current [A]"], "scale", 1)
         self.z_scale = self._simulation.parameter_values.evaluate(v_scale / i_scale)
 
-    def set_up_for_eis(self, model: pybamm.BaseModel) -> pybamm.BaseModel:
+    def set_up_for_eis(
+        self, model: pybamm.BaseModel, initial_current: float
+    ) -> pybamm.BaseModel:
         """
         Set up the model for electrochemical impedance spectroscopy (EIS) simulations.
         This method adds the necessary algebraic equations and variables to the model.
@@ -178,7 +223,7 @@ def set_up_for_eis(self, model: pybamm.BaseModel) -> pybamm.BaseModel:
             "Current function [A]", {"Time [s]": pybamm.t}
         )
         model.algebraic[I_cell] = I - I_applied
-        model.initial_conditions[I_cell] = 0
+        model.initial_conditions[I_cell] = initial_current
 
         return model
 
@@ -200,7 +245,7 @@ def _initialise_eis_matrices(self, inputs: "Inputs") -> None:
             If the model hasn't been built yet.
         """
         built_model = self._simulation.built_model
-        M = self._simulation.built_model.mass_matrix.entries
+        M = built_model.mass_matrix.entries
         self._simulation.solver.set_up(built_model, inputs=inputs)
 
         # Convert inputs to casadi format if needed
@@ -210,18 +255,49 @@ def _initialise_eis_matrices(self, inputs: "Inputs") -> None:
             else inputs or []
         )
 
+        ## Stationary EIS
         # 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()
+        y = built_model.concatenated_initial_conditions.evaluate(0, inputs=inputs)
+        J = built_model.jac_rhs_algebraic_eval(0, y, casadi_inputs).sparse()
 
         # Convert to Compressed Sparse Column format
-        self.M = csc_matrix(M)
-        self._jac = csc_matrix(jac)
+        self._mass = csc_matrix(M)
+        self._jac = csc_matrix(J)
 
         # Add forcing to the RHS on the current density
-        self.b = np.zeros(y0.shape)
+        self.b = np.zeros(y.shape)
         self.b[-1] = -1
 
+        ## Operando EIS
+        if self.time_data is not None:
+            # Initial state
+            t = np.asarray([0])
+            inputs = inputs or {}
+            sol = pybamm.Solution([t], [y], built_model, inputs)
+            state = TimeSeriesState(sol=sol, inputs=inputs, t=t)
+
+            self._jac_at_time_t = [self._jac]
+            for t in self.time_data[1:]:
+                # Step forwards in time
+                dt = (t - state.t).item()
+                new_sol = self._simulation.solver.step(
+                    state.sol, built_model, dt, inputs=state.inputs, save=False
+                )
+                state = TimeSeriesState(sol=new_sol, inputs=state.inputs, t=t)
+
+                # Extract necessary attributes from the model
+                y = state.as_ndarray()
+                J = built_model.jac_rhs_algebraic_eval(t, y, casadi_inputs).sparse()
+
+                if np.abs(y[-1]) > 1e-10:
+                    warnings.warn(
+                        f"The current is not zero at the requested EIS point at V={y[-2]} V.",
+                        stacklevel=2,
+                    )
+
+                # Convert to Compressed Sparse Column format
+                self._jac_at_time_t.append(csc_matrix(J))
+
     def solve(
         self,
         inputs: "Inputs | list[Inputs] | None" = None,
@@ -288,7 +364,10 @@ def _catch_errors(self, inputs: "list[Inputs]") -> list[Solution | FailedSolutio
                 try:
                     simulations.append(self._solve(x))
                 except (ZeroDivisionError, RuntimeError, ValueError) as e:
-                    if isinstance(e, ValueError) and str(e) not in self.exception:
+                    if (
+                        isinstance(e, ValueError)
+                        and str(e) not in self._simulation.exception
+                    ):
                         raise  # Raise the error if it doesn't match the expected list
                     simulations.append(
                         FailedSolution(["Impedance"], [k for k in x.keys()])
@@ -320,10 +399,26 @@ def _solve(self, inputs: "Inputs") -> Solution:
         # Always run initialise_eis_matrices, after rebuilding the model if necessary
         self._model_rebuild(inputs)
 
-        zs = [self.calculate_impedance(frequency) for frequency in self._f_eval]
-
         solution = Solution()
-        solution.set_solution_variable("Impedance", data=np.asarray(zs) * self.z_scale)
+        if self.time_data is None:
+            ## Stationary EIS
+            zs = [self.calculate_impedance(frequency) for frequency in self._f_eval]
+            solution.set_solution_variable(
+                "Impedance", data=np.asarray(zs) * self.z_scale
+            )
+
+        else:
+            ## Operando EIS
+            zs_at_time_t = []
+            for i in range(len(self.time_data)):
+                self._jac = self._jac_at_time_t[i]
+                zs = [self.calculate_impedance(frequency) for frequency in self._f_eval]
+                zs_at_time_t.append(zs)
+            solution.set_solution_variable("Time [s]", data=np.asarray(self.time_data))
+            solution.set_solution_variable(
+                "Impedance", data=np.asarray(zs_at_time_t) * self.z_scale
+            )
+
         return solution
 
     def calculate_impedance(self, frequency):
@@ -345,7 +440,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._mass - self._jac
 
         # Solve the system
         x = spsolve(A, self.b)
@@ -365,6 +460,10 @@ def parameter_values(self):
     def input_parameter_names(self):
         return self._simulation.input_parameter_names
 
+    @property
+    def time_data(self):
+        return self._simulation.time_data
+
     @property
     def has_sensitivities(self):
         return False
diff --git a/pybop/pybamm/simulator.py b/pybop/pybamm/simulator.py
index 0140dcad8..3827837af 100644
--- a/pybop/pybamm/simulator.py
+++ b/pybop/pybamm/simulator.py
@@ -547,6 +547,10 @@ def initial_state(self):
     def experiment(self):
         return self._experiment
 
+    @property
+    def time_data(self):
+        return self._t_interp if self._t_interp is not None else self._t_eval
+
     @property
     def solver(self):
         return self._solver

From 3a75778baf3079cfa3aea9a5273d5efc1551b5de Mon Sep 17 00:00:00 2001
From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com>
Date: Thu, 4 Dec 2025 12:20:38 +0000
Subject: [PATCH 05/20] Create operando_eis.py

---
 .../battery_parameterisation/operando_eis.py  | 61 +++++++++++++++++++
 1 file changed, 61 insertions(+)
 create mode 100644 examples/scripts/battery_parameterisation/operando_eis.py

diff --git a/examples/scripts/battery_parameterisation/operando_eis.py b/examples/scripts/battery_parameterisation/operando_eis.py
new file mode 100644
index 000000000..7ea3455ff
--- /dev/null
+++ b/examples/scripts/battery_parameterisation/operando_eis.py
@@ -0,0 +1,61 @@
+import matplotlib.pyplot as plt
+import numpy as np
+import pybamm
+
+import pybop
+
+"""
+Example demonstrating EIS applied during operation (slow dis/charge) of the cell.
+"""
+
+# Define model and parameter values
+model = pybamm.lithium_ion.SPMe(
+    options={"surface form": "differential", "contact resistance": "true"}
+)
+parameter_values = pybamm.ParameterValues("Chen2020")
+parameter_values["Contact resistance [Ohm]"] = 0.02
+parameter_values.set_initial_state("2.85 V", options=model.options)
+
+# Set up and run a charge/discharge experiment
+C_rate = parameter_values["Nominal cell capacity [A.h]"]
+dataset = pybop.Dataset(
+    {
+        "Time [s]": np.asarray(
+            [0, 1, 1001, 2001, 3001, 3002, 3003, 4003, 5003, 6003, 6004]
+        ),
+        "Current function [A]": np.asarray([0, -1, -1, -1, -1, 0, 1, 1, 1, 1, 0])
+        * C_rate
+        / 3,
+    }
+)
+
+sim = pybop.pybamm.Simulator(
+    model, parameter_values=parameter_values, protocol=dataset
+)
+solution = sim.solve()
+solution.plot()
+
+# Set up and run the simulation
+n_frequency = 60
+solution = pybop.pybamm.EISSimulator(
+    model,
+    parameter_values=parameter_values,
+    f_eval=np.logspace(-4, 5, n_frequency),
+    protocol=dataset,
+).solve()
+
+fig, ax = plt.subplots()
+n_time_steps = len(solution["Time [s]"].data)
+for i in range(n_time_steps):
+    impedance = solution["Impedance"].data[i, :]
+    ax.plot(
+        np.real(impedance),
+        -np.imag(impedance),
+        "-" if i < n_time_steps / 2 else "--",
+        label=f"t={solution['Time [s]'].data[i]}s",
+    )
+ax.set(xlabel=r"$Z_r(\omega)$ [$\Omega$]", ylabel=r"$-Z_j(\omega)$ [$\Omega$]")
+ax.set_aspect("equal", "box")
+ax.legend()
+ax.set_ylim([0, ax.get_xlim()[1]])
+plt.show()

From fa83b9593e8be21b3b93e54a75da79b291cbfe0c Mon Sep 17 00:00:00 2001
From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com>
Date: Mon, 12 Jan 2026 12:35:51 +0000
Subject: [PATCH 06/20] Fix merge

---
 pybop/pybamm/eis_simulator.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/pybop/pybamm/eis_simulator.py b/pybop/pybamm/eis_simulator.py
index 776849bb0..2ee1dc586 100644
--- a/pybop/pybamm/eis_simulator.py
+++ b/pybop/pybamm/eis_simulator.py
@@ -281,7 +281,7 @@ def _initialise_eis_matrices(self, inputs: "Inputs") -> None:
             for t in self.time_data[1:]:
                 # Step forwards in time
                 dt = (t - state.t).item()
-                new_sol = self._simulation.solver.step(
+                new_sol = self.simulation.solver.step(
                     state.sol, built_model, dt, inputs=state.inputs, save=False
                 )
                 state = TimeSeriesState(sol=new_sol, inputs=state.inputs, t=t)
@@ -458,7 +458,7 @@ def input_parameter_names(self):
 
     @property
     def time_data(self):
-        return self._simulation.time_data
+        return self._simulator.time_data
 
     @property
     def has_sensitivities(self):

From f7380a06844adddb72b40f1e9fa9e43e2344dbba Mon Sep 17 00:00:00 2001
From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com>
Date: Thu, 12 Feb 2026 17:20:47 +0000
Subject: [PATCH 07/20] Finish merging

---
 .../battery_parameterisation/operando_eis.py  |  8 ++--
 pybop/pybamm/eis_simulator.py                 | 38 +++----------------
 2 files changed, 8 insertions(+), 38 deletions(-)

diff --git a/examples/scripts/battery_parameterisation/operando_eis.py b/examples/scripts/battery_parameterisation/operando_eis.py
index 7ea3455ff..9c2d33a7b 100644
--- a/examples/scripts/battery_parameterisation/operando_eis.py
+++ b/examples/scripts/battery_parameterisation/operando_eis.py
@@ -23,15 +23,13 @@
         "Time [s]": np.asarray(
             [0, 1, 1001, 2001, 3001, 3002, 3003, 4003, 5003, 6003, 6004]
         ),
-        "Current function [A]": np.asarray([0, -1, -1, -1, -1, 0, 1, 1, 1, 1, 0])
+        "Current [A]": np.asarray([0, -1, -1, -1, -1, 0, 1, 1, 1, 1, 0])
         * C_rate
         / 3,
     }
 )
 
-sim = pybop.pybamm.Simulator(
-    model, parameter_values=parameter_values, protocol=dataset
-)
+sim = pybop.pybamm.Simulator(model, parameter_values=parameter_values, protocol=dataset)
 solution = sim.solve()
 solution.plot()
 
@@ -47,7 +45,7 @@
 fig, ax = plt.subplots()
 n_time_steps = len(solution["Time [s]"].data)
 for i in range(n_time_steps):
-    impedance = solution["Impedance"].data[i, :]
+    impedance = solution["Impedance"].data[:, i]
     ax.plot(
         np.real(impedance),
         -np.imag(impedance),
diff --git a/pybop/pybamm/eis_simulator.py b/pybop/pybamm/eis_simulator.py
index 4936d1130..96cda731c 100644
--- a/pybop/pybamm/eis_simulator.py
+++ b/pybop/pybamm/eis_simulator.py
@@ -12,38 +12,13 @@
 if TYPE_CHECKING:
     from pybop.parameters.parameter import Inputs
 from pybop.parameters.parameter import Parameter, Parameters
+from pybop.processing.dataset import Dataset
 from pybop.pybamm.simulator import Simulator
 from pybop.pybamm.utils import SymbolReplacer
 from pybop.simulators.base_simulator import BaseSimulator, Solution
-from pybop.processing.dataset import Dataset
-from pybop.pybamm.utils import SymbolReplacer
 from pybop.simulators.failed_solution import FailedSolution
 
 
-@dataclass
-class TimeSeriesState:
-    """
-    The current state of a time series model that is a PyBaMM model.
-    """
-
-    sol: pybamm.Solution
-    inputs: "Inputs"
-    t: float = 0.0
-
-    def as_ndarray(self) -> np.ndarray:
-        ncol = self.sol.y.shape[1]
-        if ncol > 1:
-            y = self.sol.y[:, -1]
-        else:
-            y = self.sol.y
-        if isinstance(y, casadi.DM):
-            y = y.full()
-        return y
-
-    def __len__(self):
-        return self.sol.y.shape[0]
-
-
 @dataclass
 class TimeSeriesState:
     """
@@ -134,18 +109,15 @@ def __init__(
         self._f_eval = f_eval
         parameter_values = parameter_values or model.default_parameter_values
         if protocol is None:  # perform stationary EIS by default
-            parameter_values["Current function [A]"] = 0
+            parameter_values["Current [A]"] = 0
             initial_current = 0
         elif isinstance(protocol, pybamm.Experiment):
             raise ValueError("EISSimulator cannot simulate a pybamm.Experiment.")
-        elif (
-            isinstance(protocol, Dataset)
-            and "Current function [A]" in protocol.data.keys()
-        ):
+        elif isinstance(protocol, Dataset) and "Current [A]" in protocol.data.keys():
             parameter_values["Current function [A]"] = pybamm.Interpolant(
-                protocol["Time [s]"], protocol["Current function [A]"], pybamm.t
+                protocol["Time [s]"], protocol["Current [A]"], pybamm.t
             )
-            initial_current = protocol["Current function [A]"][0]
+            initial_current = protocol["Current [A]"][0]
         model = self.set_up_for_eis(model, initial_current=float(initial_current))
 
         # Unpack the uncertain parameters from the parameter values

From a5a39dbb6c7444d2134a736d33e983ea634295c6 Mon Sep 17 00:00:00 2001
From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com>
Date: Thu, 12 Feb 2026 17:19:47 +0000
Subject: [PATCH 08/20] Catch solver errors

---
 pybop/pybamm/eis_simulator.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/pybop/pybamm/eis_simulator.py b/pybop/pybamm/eis_simulator.py
index 96cda731c..a4049b2e1 100644
--- a/pybop/pybamm/eis_simulator.py
+++ b/pybop/pybamm/eis_simulator.py
@@ -6,6 +6,7 @@
 import casadi
 import numpy as np
 import pybamm
+from pybamm import SolverError
 from scipy.sparse import csc_matrix
 from scipy.sparse.linalg import spsolve
 
@@ -362,7 +363,7 @@ def _catch_errors(self, inputs: "list[Inputs]") -> list[Solution | FailedSolutio
             for x in inputs:
                 try:
                     simulations.append(self._solve(x))
-                except (ZeroDivisionError, RuntimeError, ValueError):
+                except (SolverError, ZeroDivisionError, RuntimeError, ValueError):
                     simulations.append(
                         FailedSolution(["Impedance"], [k for k in x.keys()])
                     )

From 79e58ee0b04f7f702edff397508fb5f38bc6deae Mon Sep 17 00:00:00 2001
From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com>
Date: Fri, 13 Feb 2026 11:18:06 +0000
Subject: [PATCH 09/20] Allow solve at time zero

---
 pybop/pybamm/eis_simulator.py | 50 +++++++++++++++--------------------
 1 file changed, 22 insertions(+), 28 deletions(-)

diff --git a/pybop/pybamm/eis_simulator.py b/pybop/pybamm/eis_simulator.py
index a4049b2e1..25c771a5d 100644
--- a/pybop/pybamm/eis_simulator.py
+++ b/pybop/pybamm/eis_simulator.py
@@ -277,13 +277,14 @@ def _initialise_eis_matrices(self, inputs: "Inputs") -> None:
             state = TimeSeriesState(sol=sol, inputs=inputs, t=t)
 
             self._jac_at_time_t = [self._jac]
-            for t in self.time_data[1:]:
+            for t in self.time_data:
                 # Step forwards in time
                 dt = (t - state.t).item()
-                new_sol = self.simulation.solver.step(
-                    state.sol, built_model, dt, inputs=state.inputs, save=False
-                )
-                state = TimeSeriesState(sol=new_sol, inputs=state.inputs, t=t)
+                if dt > 0:
+                    new_sol = self.simulation.solver.step(
+                        state.sol, built_model, dt, inputs=state.inputs, save=False
+                    )
+                    state = TimeSeriesState(sol=new_sol, inputs=state.inputs, t=t)
 
                 # Extract necessary attributes from the model
                 y = state.as_ndarray()
@@ -319,16 +320,11 @@ def solve(
         Solution | list[Solution]
             Complex impedance results.
         """
-        if calculate_sensitivities:
-            warnings.warn(
-                "Sensitivity calculation not implemented for EIS simulations",
-                stacklevel=2,
-            )
-
+        inputs = inputs or {}
         if not isinstance(inputs, list):
-            return self._catch_errors([inputs])[0]
+            return self.solve_batch([inputs], calculate_sensitivities)[0]
 
-        return self._catch_errors(inputs)
+        return self.solve_batch(inputs, calculate_sensitivities)
 
     def solve_batch(
         self, inputs: "list[Inputs]" = None, calculate_sensitivities: bool = False
@@ -355,25 +351,23 @@ def solve_batch(
                 stacklevel=2,
             )
 
-        return self._catch_errors(inputs)
-
-    def _catch_errors(self, inputs: "list[Inputs]") -> list[Solution | FailedSolution]:
-        if not self.debug_mode:
-            simulations = []
-            for x in inputs:
-                try:
-                    simulations.append(self._solve(x))
-                except (SolverError, ZeroDivisionError, RuntimeError, ValueError):
-                    simulations.append(
-                        FailedSolution(["Impedance"], [k for k in x.keys()])
-                    )
-            return simulations
+        if len(inputs) == 1:
+            return [self._catch_errors(inputs[0])]
 
         simulations = []
         for x in inputs:
-            simulations.append(self._solve(x))
+            simulations.append(self._catch_errors(x))
         return simulations
 
+    def _catch_errors(self, inputs: "Inputs") -> Solution | FailedSolution:
+        if not self.debug_mode:
+            try:
+                return self._solve(inputs)
+            except (SolverError, ZeroDivisionError, RuntimeError, ValueError):
+                return FailedSolution(["Impedance"], [k for k in inputs.keys()])
+
+        return self._solve(inputs)
+
     def _solve(self, inputs: "Inputs") -> Solution:
         """
         Run the EIS simulation to calculate impedance at all specified frequencies.
@@ -411,7 +405,7 @@ def _solve(self, inputs: "Inputs") -> Solution:
                 zs_at_time_t.append(zs)
             solution.set_solution_variable("Time [s]", data=np.asarray(self.time_data))
             solution.set_solution_variable(
-                "Impedance", data=np.asarray(zs_at_time_t) * self.z_scale
+                "Impedance", data=np.asarray(zs_at_time_t).T * self.z_scale
             )
 
         return solution

From 942f7b1dafdd7f1f9fd40561ee399bfb2962bbfa Mon Sep 17 00:00:00 2001
From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com>
Date: Fri, 13 Feb 2026 14:38:13 +0000
Subject: [PATCH 10/20] Allow time points as protocol

---
 pybop/pybamm/eis_simulator.py | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/pybop/pybamm/eis_simulator.py b/pybop/pybamm/eis_simulator.py
index 25c771a5d..6c360a16a 100644
--- a/pybop/pybamm/eis_simulator.py
+++ b/pybop/pybamm/eis_simulator.py
@@ -119,7 +119,9 @@ def __init__(
                 protocol["Time [s]"], protocol["Current [A]"], pybamm.t
             )
             initial_current = protocol["Current [A]"][0]
-        model = self.set_up_for_eis(model, initial_current=float(initial_current))
+        elif isinstance(protocol, np.ndarray):
+            initial_current = 0  # assumption
+        model = self.set_up_for_eis(model.new_copy(), initial_current=float(initial_current))
 
         # Unpack the uncertain parameters from the parameter values
         parameters = Parameters()
@@ -276,7 +278,7 @@ def _initialise_eis_matrices(self, inputs: "Inputs") -> None:
             sol = pybamm.Solution([t], [y], built_model, inputs)
             state = TimeSeriesState(sol=sol, inputs=inputs, t=t)
 
-            self._jac_at_time_t = [self._jac]
+            self._jac_at_time_t = []
             for t in self.time_data:
                 # Step forwards in time
                 dt = (t - state.t).item()

From e0346c8a092f44efa4d546caddda1a08b6db8041 Mon Sep 17 00:00:00 2001
From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com>
Date: Fri, 13 Feb 2026 15:34:41 +0000
Subject: [PATCH 11/20] Switch dataset shape back

---
 examples/scripts/battery_parameterisation/operando_eis.py | 2 +-
 pybop/processing/dataset.py                               | 2 +-
 pybop/pybamm/eis_simulator.py                             | 2 +-
 3 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/examples/scripts/battery_parameterisation/operando_eis.py b/examples/scripts/battery_parameterisation/operando_eis.py
index 9c2d33a7b..c8be85d99 100644
--- a/examples/scripts/battery_parameterisation/operando_eis.py
+++ b/examples/scripts/battery_parameterisation/operando_eis.py
@@ -45,7 +45,7 @@
 fig, ax = plt.subplots()
 n_time_steps = len(solution["Time [s]"].data)
 for i in range(n_time_steps):
-    impedance = solution["Impedance"].data[:, i]
+    impedance = solution["Impedance"].data[i, :]
     ax.plot(
         np.real(impedance),
         -np.imag(impedance),
diff --git a/pybop/processing/dataset.py b/pybop/processing/dataset.py
index cf5a91830..b2aa2fa19 100644
--- a/pybop/processing/dataset.py
+++ b/pybop/processing/dataset.py
@@ -137,7 +137,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."
                 )
diff --git a/pybop/pybamm/eis_simulator.py b/pybop/pybamm/eis_simulator.py
index 6c360a16a..3e358955d 100644
--- a/pybop/pybamm/eis_simulator.py
+++ b/pybop/pybamm/eis_simulator.py
@@ -407,7 +407,7 @@ def _solve(self, inputs: "Inputs") -> Solution:
                 zs_at_time_t.append(zs)
             solution.set_solution_variable("Time [s]", data=np.asarray(self.time_data))
             solution.set_solution_variable(
-                "Impedance", data=np.asarray(zs_at_time_t).T * self.z_scale
+                "Impedance", data=np.asarray(zs_at_time_t) * self.z_scale
             )
 
         return solution

From f7d494c65b62edd3b66d17a62e6803456a07b26e Mon Sep 17 00:00:00 2001
From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com>
Date: Fri, 13 Feb 2026 17:04:46 +0000
Subject: [PATCH 12/20] Try with an extra initialise

---
 pybop/pybamm/eis_simulator.py | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/pybop/pybamm/eis_simulator.py b/pybop/pybamm/eis_simulator.py
index 3e358955d..ce88f081f 100644
--- a/pybop/pybamm/eis_simulator.py
+++ b/pybop/pybamm/eis_simulator.py
@@ -155,6 +155,10 @@ def __init__(
         i_scale = getattr(model.variables["Current [A]"], "scale", 1)
         self.z_scale = self.parameter_values.evaluate(v_scale / i_scale)
 
+        self._initialise_eis_matrices(
+            inputs=self.parameters.to_dict(self.parameters.get_initial_values())
+        )  # not sure why this extra initialise is required before the first solve...
+
     def set_up_for_eis(
         self, model: pybamm.BaseModel, initial_current: float
     ) -> pybamm.BaseModel:

From 85649fe482bd0da311f077fced33ad4f671604e8 Mon Sep 17 00:00:00 2001
From: Ombrini <91598680+Ombrini@users.noreply.github.com>
Date: Wed, 12 Aug 2026 14:20:11 +0200
Subject: [PATCH 13/20] First attempt: add GITT parameter estimation example
 and EIS functionality

- Introduced a new example script for parameter estimation from GITT experiments with operando EIS.
- Enhanced the EISSimulator class to support operando simulations and added methods for handling impedance data.
- Updated the WeightedCost class to maintain distinct targets for different cost functions.
- Added unit tests for the operando EIS simulator to ensure correct functionality and output.
---
 .../battery_parameterisation/gitt_eis.py      | 162 +++++++++++
 pybop/costs/weighted_cost.py                  |   9 +-
 pybop/plot/problem.py                         |  71 ++++-
 pybop/pybamm/__init__.py                      |   2 +-
 pybop/pybamm/eis_simulator.py                 | 260 ++++++++++++++----
 tests/unit/test_cost.py                       |  12 +
 tests/unit/test_simulator.py                  | 141 ++++++++++
 7 files changed, 603 insertions(+), 54 deletions(-)
 create mode 100644 examples/scripts/battery_parameterisation/gitt_eis.py

diff --git a/examples/scripts/battery_parameterisation/gitt_eis.py b/examples/scripts/battery_parameterisation/gitt_eis.py
new file mode 100644
index 000000000..a82116e2e
--- /dev/null
+++ b/examples/scripts/battery_parameterisation/gitt_eis.py
@@ -0,0 +1,162 @@
+import numpy as np
+import pybamm
+
+import pybop
+
+"""
+Example demonstrating parameter estimation from a GITT experiment in which an EIS
+spectrum is acquired at the end of each pulse ("operando" EIS).
+
+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 columns 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
+operando EIS is 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.0
+
+
+# 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(-2, 4, 50)
+columns = pybop.pybamm.eis_column_names(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 columns 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 columns},
+    },
+    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 columns:
+    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-7, 1e-4)
+        ),
+    }
+)
+
+# 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 columns 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=columns)
+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/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/problem.py b/pybop/plot/problem.py
index 61adb604e..6f0d289a1 100644
--- a/pybop/plot/problem.py
+++ b/pybop/plot/problem.py
@@ -6,6 +6,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.pybamm.eis_simulator import parse_eis_column_names
 from pybop.simulators.solution import Solution
 
 
@@ -77,15 +78,33 @@ 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 columns are plotted as a Nyquist plot per acquisition, not as a series
+    # over the domain, so separate them from the remaining targets
+    frequencies, real_names, imaginary_names = parse_eis_column_names(problem.target)
+    variables = [
+        var
+        for var in problem.target
+        if var not in set(real_names) | set(imaginary_names)
+    ]
+    acquisitions = (
+        _eis_acquisitions(target_output, real_names, imaginary_names)
+        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(variables) + len(acquisitions),
+        allow_single_axis=False,
     )
-    for i, var in enumerate(problem.target):
+
+    # Create a plot for each output
+    for i, var in enumerate(variables):
         ax = axes[i % len(axes)]
         if create_figure:
             fig = backend.create_figure(
@@ -146,5 +165,51 @@ def problem(
         if show:
             backend.show_figure(fig)
 
+    # Add a Nyquist plot for each acquired spectrum
+    for j, row in enumerate(acquisitions):
+        i = len(variables) + j
+        ax = axes[i % len(axes)]
+        if create_figure:
+            fig = backend.create_figure(
+                style={"bg_color": "white", "width": 600, "height": 600},
+            )
+            figures = np.append(figures, fig)
+        else:
+            fig = figures[i % len(figures)]
+
+        backend.update_axes_titles(fig, ax, r"$Z_{re} / \Omega$", r"$-Z_{im} / \Omega$")
+        backend.update_plot_titles(
+            fig, ax, f"{title}: {remove_brackets(domain)} = {target_domain[row]:g}"
+        )
+
+        for output, label, style in (
+            (target_output, "Reference", {"linestyle": "none", "marker": "."}),
+            (model_output, "Model", {"linestyle": "solid", "marker": "none"}),
+        ):
+            backend.plot_trace(
+                backend.line(
+                    x=[output[name].data[row] for name in real_names],
+                    y=[-output[name].data[row] for name in imaginary_names],
+                    label=label,
+                    style=style,
+                ),
+                fig,
+                ax=ax,
+            )
+
+        backend.legend(fig, axes=ax)
+        if show:
+            backend.show_figure(fig)
+
     if not show:
         return figures[0] if len(figures) == 1 else figures
+
+
+def _eis_acquisitions(
+    target_output, real_names: list[str], imaginary_names: 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_names, *imaginary_names)]
+    )
+    return np.flatnonzero(np.any(measured != 0.0, axis=0))
diff --git a/pybop/pybamm/__init__.py b/pybop/pybamm/__init__.py
index a1cfcce9a..4648d5513 100644
--- a/pybop/pybamm/__init__.py
+++ b/pybop/pybamm/__init__.py
@@ -1,5 +1,5 @@
 from .simulator import Simulator
-from .eis_simulator import EISSimulator
+from .eis_simulator import EISSimulator, eis_column_names, parse_eis_column_names
 from .parameter_utils import set_formation_concentrations, cell_mass, cell_volume
 from .design_variables import add_variable_to_model
 from .utils import RecommendedSolver, SymbolReplacer
diff --git a/pybop/pybamm/eis_simulator.py b/pybop/pybamm/eis_simulator.py
index 8445b0319..6f9304bf0 100644
--- a/pybop/pybamm/eis_simulator.py
+++ b/pybop/pybamm/eis_simulator.py
@@ -1,3 +1,4 @@
+import re
 import warnings
 from copy import copy
 from typing import TYPE_CHECKING
@@ -10,12 +11,71 @@
 
 if TYPE_CHECKING:
     from pybop.parameters.parameter import Inputs
+from pybop.processing.dataset import Dataset
 from pybop.pybamm.simulator import Simulator
 from pybop.pybamm.utils import SymbolReplacer
 from pybop.simulators.base_simulator import BaseSimulator, Solution
 from pybop.simulators.failed_solution import FailedSolution
 
 
+def eis_column_names(f_eval: np.ndarray | list[float]) -> list[str]:
+    """
+    Return the dataset column names for an impedance spectrum, two real-valued columns
+    (real and imaginary) per frequency.
+
+    Complex data is split into real columns because the error measures square the
+    residual, and `r**2` is not `abs(r)**2` for a complex array.
+
+    Parameters
+    ----------
+    f_eval : np.ndarray | list[float]
+        The frequencies at which the impedance is evaluated.
+
+    Returns
+    -------
+    list[str]
+        Column names, ordered (real, imaginary) for each frequency in turn.
+    """
+    return [
+        f"Impedance {part} [Ohm] ({f:.6g} Hz)"
+        for f in f_eval
+        for part in ("real", "imaginary")
+    ]
+
+
+def parse_eis_column_names(names: list[str]) -> tuple[np.ndarray, list[str], list[str]]:
+    """
+    Pick out the impedance columns from a list of variable names, inverting
+    `eis_column_names`.
+
+    Parameters
+    ----------
+    names : list[str]
+        Variable names, which may or may not include impedance columns.
+
+    Returns
+    -------
+    tuple[np.ndarray, list[str], list[str]]
+        The frequencies in increasing order, and the corresponding real and imaginary
+        column names. All three are empty if no impedance columns are present.
+    """
+    pattern = re.compile(r"^Impedance (real|imaginary) \[Ohm\] \((\S+) Hz\)$")
+
+    frequencies = {}
+    for name in names:
+        match = pattern.match(name)
+        if match:
+            part, frequency = match.group(1), float(match.group(2))
+            frequencies.setdefault(frequency, {})[part] = name
+
+    ordered = sorted(f for f, parts in frequencies.items() if len(parts) == 2)
+    return (
+        np.asarray(ordered),
+        [frequencies[f]["real"] for f in ordered],
+        [frequencies[f]["imaginary"] for f in ordered],
+    )
+
+
 class EISSimulator(BaseSimulator):
     """
     A class to extend a PyBaMM model for EIS, automatically build/rebuild a pybamm.Simulation to obtain
@@ -38,6 +98,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 columns given by
+        `eis_column_names(f_eval)`. The impedance columns are non-zero at the times at
+        which a spectrum was measured, and zero elsewhere; the simulator computes a
+        spectrum at exactly those times ("operando" EIS). If None, a single spectrum is
+        computed about the initial state ("stationary" EIS).
     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 +135,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 +154,17 @@ def __init__(
 
         super().__init__(parameters=parameter_values)
 
-        # Set up a simulation
+        # Locate the times at which to compute a spectrum, if any
+        self._eis_columns = eis_column_names(f_eval)
+        self._eis_indices = self._locate_eis_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 columns 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 +179,73 @@ 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_eis_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._eis_columns) - set(protocol.keys())
+        if missing:
+            raise ValueError(
+                "The protocol dataset is missing impedance columns, e.g. "
+                f"'{sorted(missing)[0]}'. Name them with pybop.pybamm.eis_column_names(f_eval), "
+                "using zeros at the times where no spectrum was measured."
+            )
+
+        measured = np.asarray([protocol[name] for name in self._eis_columns])
+        indices = np.flatnonzero(np.any(measured != 0.0, axis=0))
+        if indices.size == 0:
+            raise ValueError(
+                "The impedance columns are zero everywhere, so there is nothing to fit. "
+                "Set them to the measured spectra at the times of acquisition."
+            )
+        return indices
+
+    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 _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.
@@ -182,46 +317,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 +396,7 @@ 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.output_names, x.keys()))
             return simulations
 
         simulations = []
@@ -301,29 +408,77 @@ 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.
+        For an operando simulation, 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 for an operando simulation.
         """
-        # 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._eis_indices is None:
+            y0 = self.simulation.built_model.concatenated_initial_conditions.evaluate(
+                0, inputs=inputs
+            )
+            solution = Solution()
+            solution.set_solution_variable(
+                "Impedance", data=self._spectrum(self._jacobian(0, y0, inputs))
+            )
+            return solution
+
+        return self._solve_operando(inputs)
+
+    def _solve_operando(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._eis_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
+        columns = {name: np.zeros(len(t)) for name in self._eis_columns}
+        for i in self._eis_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 EIS point at t={t[i]} s, "
+                    "so the impedance is linearised about a non-zero operating point.",
+                    stacklevel=2,
+                )
+            zs = self._spectrum(self._jacobian(t[i], y_i, inputs))
+            for j, z in enumerate(zs):
+                columns[self._eis_columns[2 * j]][i] = z.real
+                columns[self._eis_columns[2 * j + 1]][i] = z.imag
+
+        for name, data in columns.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 +489,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 +499,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 +519,13 @@ def parameter_values(self):
     def input_parameter_names(self):
         return self._simulator.input_parameter_names
 
+    @property
+    def output_names(self) -> list[str]:
+        """The names of the variables returned by a solve."""
+        if self._eis_indices is None:
+            return ["Impedance"]
+        return ["Time [s]", "Voltage [V]", *self._eis_columns]
+
     @property
     def has_sensitivities(self):
         return False
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_simulator.py b/tests/unit/test_simulator.py
index 748617d17..375ca01d6 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,143 @@ 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 TestOperandoEISSimulator:
+    """
+    A class to test the operando mode of the pybamm.EISSimulator class.
+    """
+
+    pytestmark = pytest.mark.unit
+
+    @pytest.fixture
+    def setup(self):
+        model = pybamm.lithium_ion.SPM(
+            options={"surface form": "differential", "contact resistance": "true"}
+        )
+        parameter_values = pybamm.ParameterValues("Chen2020")
+        parameter_values["Contact resistance [Ohm]"] = 0.0
+        parameter_values.set_initial_state(0.9)
+
+        time = np.arange(0, 601, 10.0)
+        f_eval = np.logspace(-1, 3, 4)
+        columns = pybop.pybamm.eis_column_names(f_eval)
+
+        # Two spectra: one at the start, one part-way through the rest
+        eis_rows = [0, 40]
+        data = {
+            "Time [s]": time,
+            "Current [A]": np.zeros_like(time),
+            "Voltage [V]": np.zeros_like(time),
+        }
+        for name in columns:
+            column = np.zeros(len(time))
+            column[eis_rows] = 1.0
+            data[name] = column
+
+        dataset = pybop.Dataset(data, domain="Time [s]")
+        return model, parameter_values, dataset, f_eval, columns, eis_rows
+
+    def test_output_shape_and_zero_padding(self, setup):
+        model, parameter_values, dataset, f_eval, columns, eis_rows = setup
+        simulator = pybop.pybamm.EISSimulator(
+            model, parameter_values=parameter_values, protocol=dataset, f_eval=f_eval
+        )
+        solution = simulator.solve()
+
+        n_time = len(dataset["Time [s]"])
+        assert len(solution["Voltage [V]"].data) == n_time
+        for name in columns:
+            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), eis_rows)
+
+    def test_matches_stationary_at_initial_state(self, setup):
+        """The spectrum at t=0 must match a stationary simulation of the same state."""
+        model, parameter_values, dataset, f_eval, columns, _ = setup
+        operando = pybop.pybamm.EISSimulator(
+            model, parameter_values=parameter_values, protocol=dataset, f_eval=f_eval
+        ).solve()
+        stationary = pybop.pybamm.EISSimulator(
+            model, parameter_values=parameter_values, f_eval=f_eval
+        ).solve()
+
+        impedance = np.asarray(
+            [
+                operando[columns[2 * j]].data[0]
+                + 1j * operando[columns[2 * j + 1]].data[0]
+                for j in range(len(f_eval))
+            ]
+        )
+        np.testing.assert_allclose(impedance, stationary["Impedance"].data, rtol=1e-10)
+
+    def test_builds_once(self, setup):
+        """The model and the constant matrices are set up once, not per evaluation."""
+        model, parameter_values, dataset, f_eval, _, _ = setup
+        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=f_eval
+        )
+
+        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_column_names_round_trip(self):
+        f_eval = np.logspace(-2, 4, 5)
+        columns = pybop.pybamm.eis_column_names(f_eval)
+        # Names may reach the parser in any order, e.g. via a set
+        frequencies, real, imaginary = pybop.pybamm.parse_eis_column_names(
+            ["Voltage [V]", *reversed(columns)]
+        )
+        np.testing.assert_allclose(frequencies, f_eval, rtol=1e-6)
+        assert real == columns[::2]
+        assert imaginary == columns[1::2]
+
+        # No impedance columns present
+        frequencies, real, imaginary = pybop.pybamm.parse_eis_column_names(
+            ["Voltage [V]", "Current [A]"]
+        )
+        assert len(frequencies) == 0 and real == [] and imaginary == []
+
+    def test_dataset_errors(self, setup):
+        model, parameter_values, dataset, f_eval, columns, _ = setup
+
+        with pytest.raises(ValueError, match="missing impedance columns"):
+            pybop.pybamm.EISSimulator(
+                model,
+                parameter_values=parameter_values,
+                protocol=dataset,
+                f_eval=np.append(f_eval, 1e4),
+            )
+
+        for name in columns:
+            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=f_eval,
+            )

From 45782608a388c1a372d02b7db3d8931185f74ddb Mon Sep 17 00:00:00 2001
From: Ombrini <91598680+Ombrini@users.noreply.github.com>
Date: Wed, 12 Aug 2026 17:32:23 +0200
Subject: [PATCH 14/20] Refactor EIS functionality: update variable names, add
 impedance variable handling, and enhance plotting capabilities

---
 .../battery_parameterisation/gitt_eis.py      |  16 +-
 pybop/__init__.py                             |   2 +-
 pybop/plot/backends/base.py                   |  15 ++
 pybop/plot/backends/matplotlib.py             |  14 ++
 pybop/plot/backends/plotly.py                 |  23 +++
 pybop/plot/problem.py                         | 118 +++++++++-----
 pybop/processing/dataset.py                   |  62 ++++++++
 pybop/pybamm/__init__.py                      |   2 +-
 pybop/pybamm/eis_simulator.py                 | 142 +++++++----------
 tests/unit/test_problem.py                    |   2 +-
 tests/unit/test_simulator.py                  | 149 ++++++++++++------
 11 files changed, 353 insertions(+), 192 deletions(-)

diff --git a/examples/scripts/battery_parameterisation/gitt_eis.py b/examples/scripts/battery_parameterisation/gitt_eis.py
index a82116e2e..4922902db 100644
--- a/examples/scripts/battery_parameterisation/gitt_eis.py
+++ b/examples/scripts/battery_parameterisation/gitt_eis.py
@@ -10,7 +10,7 @@
 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 columns hold the real
+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.
 
@@ -20,7 +20,7 @@
 """
 
 # Define the model
-model = pybamm.lithium_ion.SPM(
+model = pybamm.lithium_ion.SPMe(
     options={"surface form": "differential", "contact resistance": "true"},
 )
 parameter_values = pybamm.ParameterValues("Chen2020")
@@ -74,13 +74,13 @@ def positive_exchange_current_density(c_e, c_s_surf, c_s_max, T):
 # 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(-2, 4, 50)
-columns = pybop.pybamm.eis_column_names(f_eval)
+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 columns start as a marker of which times were
+# 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)
@@ -89,7 +89,7 @@ def positive_exchange_current_density(c_e, c_s_surf, c_s_max, T):
         "Time [s]": time,
         "Current [A]": current,
         "Voltage [V]": voltage,
-        **{name: acquired.astype(float) for name in columns},
+        **{name: acquired.astype(float) for name in impedance_variables},
     },
     domain="Time [s]",
 )
@@ -105,7 +105,7 @@ def positive_exchange_current_density(c_e, c_s_surf, c_s_max, T):
 # 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 columns:
+for name in impedance_variables:
     dataset[name] = np.where(
         acquired, pybop.add_noise(solution[name].data, sigma_z), 0.0
     )
@@ -133,13 +133,13 @@ def positive_exchange_current_density(c_e, c_s_surf, c_s_max, T):
 
 # 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 columns are zero at most times, which would
+# 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=columns)
+impedance_cost = pybop.SumSquaredError(dataset, target=impedance_variables)
 cost = pybop.WeightedCost(voltage_cost, impedance_cost, weights=[1.0, 1e3])
 problem = pybop.Problem(simulator, cost)
 
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/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 6f0d289a1..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,7 +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.pybamm.eis_simulator import parse_eis_column_names
+from pybop.processing.dataset import parse_impedance_variables
 from pybop.simulators.solution import Solution
 
 
@@ -81,16 +83,18 @@ def problem(
     # Import plotting backend
     backend = get_backend_from_figure(backend, figures)
 
-    # Impedance columns are plotted as a Nyquist plot per acquisition, not as a series
+    # 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_names, imaginary_names = parse_eis_column_names(problem.target)
-    variables = [
+    frequencies, real_variables, imaginary_variables = parse_impedance_variables(
+        problem.target
+    )
+    targets = [
         var
         for var in problem.target
-        if var not in set(real_names) | set(imaginary_names)
+        if var not in set(real_variables) | set(imaginary_variables)
     ]
     acquisitions = (
-        _eis_acquisitions(target_output, real_names, imaginary_names)
+        _acquisition_indices(target_output, real_variables, imaginary_variables)
         if len(frequencies) > 0
         else []
     )
@@ -99,12 +103,12 @@ def problem(
     figures, axes, create_figure, _ = backend.parse_input_axes(
         figures,
         axes,
-        num_plots=len(variables) + len(acquisitions),
+        num_plots=len(targets),
         allow_single_axis=False,
     )
 
     # Create a plot for each output
-    for i, var in enumerate(variables):
+    for i, var in enumerate(targets):
         ax = axes[i % len(axes)]
         if create_figure:
             fig = backend.create_figure(
@@ -165,51 +169,83 @@ def problem(
         if show:
             backend.show_figure(fig)
 
-    # Add a Nyquist plot for each acquired spectrum
-    for j, row in enumerate(acquisitions):
-        i = len(variables) + j
-        ax = axes[i % len(axes)]
-        if create_figure:
-            fig = backend.create_figure(
-                style={"bg_color": "white", "width": 600, "height": 600},
-            )
-            figures = np.append(figures, fig)
-        else:
-            fig = figures[i % len(figures)]
-
-        backend.update_axes_titles(fig, ax, r"$Z_{re} / \Omega$", r"$-Z_{im} / \Omega$")
+    # 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(
-            fig, ax, f"{title}: {remove_brackets(domain)} = {target_domain[row]:g}"
+            impedance_figure,
+            impedance_axes,
+            [
+                f"{remove_brackets(domain)} = {target_domain[row]:g}"
+                for row in acquisitions
+            ],
         )
 
-        for output, label, style in (
-            (target_output, "Reference", {"linestyle": "none", "marker": "."}),
-            (model_output, "Model", {"linestyle": "solid", "marker": "none"}),
-        ):
-            backend.plot_trace(
-                backend.line(
-                    x=[output[name].data[row] for name in real_names],
-                    y=[-output[name].data[row] for name in imaginary_names],
-                    label=label,
-                    style=style,
-                ),
-                fig,
-                ax=ax,
-            )
+        # 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"}),
+        )
 
-        backend.legend(fig, axes=ax)
+        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(fig)
+            backend.show_figure(impedance_figure)
 
     if not show:
         return figures[0] if len(figures) == 1 else figures
 
 
-def _eis_acquisitions(
-    target_output, real_names: list[str], imaginary_names: list[str]
+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_names, *imaginary_names)]
+        [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 15caabea9..75041bdd9 100644
--- a/pybop/processing/dataset.py
+++ b/pybop/processing/dataset.py
@@ -1,3 +1,4 @@
+import re
 import warnings
 from typing import Protocol
 
@@ -190,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/__init__.py b/pybop/pybamm/__init__.py
index 4648d5513..a1cfcce9a 100644
--- a/pybop/pybamm/__init__.py
+++ b/pybop/pybamm/__init__.py
@@ -1,5 +1,5 @@
 from .simulator import Simulator
-from .eis_simulator import EISSimulator, eis_column_names, parse_eis_column_names
+from .eis_simulator import EISSimulator
 from .parameter_utils import set_formation_concentrations, cell_mass, cell_volume
 from .design_variables import add_variable_to_model
 from .utils import RecommendedSolver, SymbolReplacer
diff --git a/pybop/pybamm/eis_simulator.py b/pybop/pybamm/eis_simulator.py
index 6f9304bf0..49453cba6 100644
--- a/pybop/pybamm/eis_simulator.py
+++ b/pybop/pybamm/eis_simulator.py
@@ -1,4 +1,3 @@
-import re
 import warnings
 from copy import copy
 from typing import TYPE_CHECKING
@@ -11,71 +10,13 @@
 
 if TYPE_CHECKING:
     from pybop.parameters.parameter import Inputs
-from pybop.processing.dataset import Dataset
+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
 from pybop.simulators.failed_solution import FailedSolution
 
 
-def eis_column_names(f_eval: np.ndarray | list[float]) -> list[str]:
-    """
-    Return the dataset column names for an impedance spectrum, two real-valued columns
-    (real and imaginary) per frequency.
-
-    Complex data is split into real columns because the error measures square the
-    residual, and `r**2` is not `abs(r)**2` for a complex array.
-
-    Parameters
-    ----------
-    f_eval : np.ndarray | list[float]
-        The frequencies at which the impedance is evaluated.
-
-    Returns
-    -------
-    list[str]
-        Column names, ordered (real, imaginary) for each frequency in turn.
-    """
-    return [
-        f"Impedance {part} [Ohm] ({f:.6g} Hz)"
-        for f in f_eval
-        for part in ("real", "imaginary")
-    ]
-
-
-def parse_eis_column_names(names: list[str]) -> tuple[np.ndarray, list[str], list[str]]:
-    """
-    Pick out the impedance columns from a list of variable names, inverting
-    `eis_column_names`.
-
-    Parameters
-    ----------
-    names : list[str]
-        Variable names, which may or may not include impedance columns.
-
-    Returns
-    -------
-    tuple[np.ndarray, list[str], list[str]]
-        The frequencies in increasing order, and the corresponding real and imaginary
-        column names. All three are empty if no impedance columns are present.
-    """
-    pattern = re.compile(r"^Impedance (real|imaginary) \[Ohm\] \((\S+) Hz\)$")
-
-    frequencies = {}
-    for name in names:
-        match = pattern.match(name)
-        if match:
-            part, frequency = match.group(1), float(match.group(2))
-            frequencies.setdefault(frequency, {})[part] = name
-
-    ordered = sorted(f for f, parts in frequencies.items() if len(parts) == 2)
-    return (
-        np.asarray(ordered),
-        [frequencies[f]["real"] for f in ordered],
-        [frequencies[f]["imaginary"] for f in ordered],
-    )
-
-
 class EISSimulator(BaseSimulator):
     """
     A class to extend a PyBaMM model for EIS, automatically build/rebuild a pybamm.Simulation to obtain
@@ -100,8 +41,8 @@ class EISSimulator(BaseSimulator):
         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 columns given by
-        `eis_column_names(f_eval)`. The impedance columns are non-zero at the times at
+        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 ("operando" EIS). If None, a single spectrum is
         computed about the initial state ("stationary" EIS).
@@ -155,11 +96,11 @@ def __init__(
         super().__init__(parameters=parameter_values)
 
         # Locate the times at which to compute a spectrum, if any
-        self._eis_columns = eis_column_names(f_eval)
-        self._eis_indices = self._locate_eis_times(protocol)
+        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 columns of
+        # 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,
@@ -185,7 +126,7 @@ def __init__(
         i_scale = getattr(model.variables["Current [A]"], "scale", 1)
         self.z_scale = self.parameter_values.evaluate(v_scale / i_scale)
 
-    def _locate_eis_times(self, protocol: Dataset | None) -> np.ndarray | None:
+    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.
@@ -193,23 +134,31 @@ def _locate_eis_times(self, protocol: Dataset | None) -> np.ndarray | None:
         if protocol is None:
             return None
 
-        missing = set(self._eis_columns) - set(protocol.keys())
+        missing = set(self._impedance_variables) - set(protocol.keys())
         if missing:
             raise ValueError(
-                "The protocol dataset is missing impedance columns, e.g. "
-                f"'{sorted(missing)[0]}'. Name them with pybop.pybamm.eis_column_names(f_eval), "
+                "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._eis_columns])
+        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 columns are zero everywhere, so there is nothing to fit. "
+                "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
@@ -239,7 +188,7 @@ def _jacobian(self, t: float, y: np.ndarray, inputs: "Inputs") -> csc_matrix:
         jac = built_model.jac_rhs_algebraic_eval(t, y, casadi_inputs).sparse()
         return csc_matrix(jac)
 
-    def _spectrum(self, jac: csc_matrix) -> np.ndarray:
+    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])
@@ -260,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]"]
@@ -275,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]"]
@@ -396,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(self.output_names, x.keys()))
+                    simulations.append(
+                        FailedSolution(self.solution_variables, x.keys())
+                    )
             return simulations
 
         simulations = []
@@ -426,13 +389,14 @@ def _solve(self, inputs: "Inputs") -> Solution:
         # Rebuild the model only if necessary, then set up the constant matrices
         self._model_rebuild(inputs)
 
-        if self._eis_indices is None:
+        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._spectrum(self._jacobian(0, y0, inputs))
+                "Impedance",
+                data=self._calculate_spectrum(self._jacobian(0, y0, inputs)),
             )
             return solution
 
@@ -448,7 +412,7 @@ def _solve_operando(self, inputs: "Inputs") -> Solution:
             raise ValueError("The time-domain simulation failed.")
 
         t, y = sim_solution.t, sim_solution.y
-        if self._eis_indices[-1] >= len(t):
+        if self._acquisition_indices[-1] >= len(t):
             raise ValueError(
                 "The time-domain simulation terminated before the last EIS time."
             )
@@ -460,21 +424,21 @@ def _solve_operando(self, inputs: "Inputs") -> Solution:
         )
 
         # Zero away from the times of acquisition, matching the dataset convention
-        columns = {name: np.zeros(len(t)) for name in self._eis_columns}
-        for i in self._eis_indices:
+        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 EIS point at t={t[i]} s, "
+                    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._spectrum(self._jacobian(t[i], y_i, inputs))
+            zs = self._calculate_spectrum(self._jacobian(t[i], y_i, inputs))
             for j, z in enumerate(zs):
-                columns[self._eis_columns[2 * j]][i] = z.real
-                columns[self._eis_columns[2 * j + 1]][i] = z.imag
+                impedance[self._impedance_variables[2 * j]][i] = z.real
+                impedance[self._impedance_variables[2 * j + 1]][i] = z.imag
 
-        for name, data in columns.items():
+        for name, data in impedance.items():
             solution.set_solution_variable(name, data=data)
         return solution
 
@@ -520,11 +484,11 @@ def input_parameter_names(self):
         return self._simulator.input_parameter_names
 
     @property
-    def output_names(self) -> list[str]:
-        """The names of the variables returned by a solve."""
-        if self._eis_indices is None:
+    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._eis_columns]
+        return ["Time [s]", "Voltage [V]", *self._impedance_variables]
 
     @property
     def has_sensitivities(self):
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_simulator.py b/tests/unit/test_simulator.py
index 375ca01d6..0ac3de1be 100644
--- a/tests/unit/test_simulator.py
+++ b/tests/unit/test_simulator.py
@@ -89,75 +89,103 @@ class TestOperandoEISSimulator:
     pytestmark = pytest.mark.unit
 
     @pytest.fixture
-    def setup(self):
-        model = pybamm.lithium_ion.SPM(
-            options={"surface form": "differential", "contact resistance": "true"}
+    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["Contact resistance [Ohm]"] = 0.0
         parameter_values.set_initial_state(0.9)
+        return parameter_values
 
-        time = np.arange(0, 601, 10.0)
-        f_eval = np.logspace(-1, 3, 4)
-        columns = pybop.pybamm.eis_column_names(f_eval)
+    @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]
 
-        # Two spectra: one at the start, one part-way through the rest
-        eis_rows = [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 columns:
-            column = np.zeros(len(time))
-            column[eis_rows] = 1.0
-            data[name] = column
-
-        dataset = pybop.Dataset(data, domain="Time [s]")
-        return model, parameter_values, dataset, f_eval, columns, eis_rows
-
-    def test_output_shape_and_zero_padding(self, setup):
-        model, parameter_values, dataset, f_eval, columns, eis_rows = setup
+        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=f_eval
+            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 columns:
+        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), eis_rows)
+            np.testing.assert_array_equal(np.flatnonzero(data), acquisitions)
 
-    def test_matches_stationary_at_initial_state(self, setup):
+    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."""
-        model, parameter_values, dataset, f_eval, columns, _ = setup
         operando = pybop.pybamm.EISSimulator(
-            model, parameter_values=parameter_values, protocol=dataset, f_eval=f_eval
+            model,
+            parameter_values=parameter_values,
+            protocol=dataset,
+            f_eval=frequencies,
         ).solve()
         stationary = pybop.pybamm.EISSimulator(
-            model, parameter_values=parameter_values, f_eval=f_eval
+            model, parameter_values=parameter_values, f_eval=frequencies
         ).solve()
 
         impedance = np.asarray(
             [
-                operando[columns[2 * j]].data[0]
-                + 1j * operando[columns[2 * j + 1]].data[0]
-                for j in range(len(f_eval))
+                operando[impedance_variables[2 * j]].data[0]
+                + 1j * operando[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, setup):
+    def test_builds_once(self, model, parameter_values, dataset, frequencies):
         """The model and the constant matrices are set up once, not per evaluation."""
-        model, parameter_values, dataset, f_eval, _, _ = setup
         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=f_eval
+            model,
+            parameter_values=parameter_values,
+            protocol=dataset,
+            f_eval=frequencies,
         )
 
         calls = {"create": 0, "set_up": 0}
@@ -182,40 +210,59 @@ def counted_set_up(*args, **kwargs):
         assert calls["create"] == 0  # built during construction
         assert calls["set_up"] == 1
 
-    def test_column_names_round_trip(self):
-        f_eval = np.logspace(-2, 4, 5)
-        columns = pybop.pybamm.eis_column_names(f_eval)
-        # Names may reach the parser in any order, e.g. via a set
-        frequencies, real, imaginary = pybop.pybamm.parse_eis_column_names(
-            ["Voltage [V]", *reversed(columns)]
+    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,
         )
-        np.testing.assert_allclose(frequencies, f_eval, rtol=1e-6)
-        assert real == columns[::2]
-        assert imaginary == columns[1::2]
+        assert len(model.algebraic) == n_algebraic
 
-        # No impedance columns present
-        frequencies, real, imaginary = pybop.pybamm.parse_eis_column_names(
-            ["Voltage [V]", "Current [A]"]
+    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)]
         )
-        assert len(frequencies) == 0 and real == [] and imaginary == []
+        # 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]
 
-    def test_dataset_errors(self, setup):
-        model, parameter_values, dataset, f_eval, columns, _ = setup
+        # No impedance variables present
+        parsed, real, imaginary = pybop.parse_impedance_variables(
+            ["Voltage [V]", "Current [A]"]
+        )
+        assert len(parsed) == 0 and real == [] and imaginary == []
 
-        with pytest.raises(ValueError, match="missing impedance columns"):
+    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(f_eval, 1e4),
+                f_eval=np.append(frequencies, 1e4),
             )
 
-        for name in columns:
+        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=f_eval,
+                f_eval=frequencies,
             )

From 8d8eec9fbfb0cdfdbb48ee61874f48b0cfc88deb Mon Sep 17 00:00:00 2001
From: Ombrini <91598680+Ombrini@users.noreply.github.com>
Date: Wed, 12 Aug 2026 17:38:41 +0200
Subject: [PATCH 15/20] Update example

---
 examples/scripts/battery_parameterisation/gitt_eis.py | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/examples/scripts/battery_parameterisation/gitt_eis.py b/examples/scripts/battery_parameterisation/gitt_eis.py
index 4922902db..c15563aba 100644
--- a/examples/scripts/battery_parameterisation/gitt_eis.py
+++ b/examples/scripts/battery_parameterisation/gitt_eis.py
@@ -20,11 +20,11 @@
 """
 
 # Define the model
-model = pybamm.lithium_ion.SPMe(
+model = pybamm.lithium_ion.SPM(
     options={"surface form": "differential", "contact resistance": "true"},
 )
 parameter_values = pybamm.ParameterValues("Chen2020")
-parameter_values["Contact resistance [Ohm]"] = 0.0
+parameter_values["Contact resistance [Ohm]"] = 0.01
 
 
 # The exchange-current density of Chen2020 hard-codes its prefactor, so redefine it with
@@ -73,7 +73,7 @@ def positive_exchange_current_density(c_e, c_s_surf, c_s_max, T):
 
 # 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(-2, 4, 50)
+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)
@@ -126,7 +126,7 @@ def positive_exchange_current_density(c_e, c_s_surf, c_s_max, T):
             pybop.Uniform(1e-16, 1e-13)
         ),
         "Positive electrode reference exchange-current density [A.m-2]": pybop.Parameter(
-            pybop.Uniform(1e-7, 1e-4)
+            pybop.Uniform(1e-6, 1e-5)
         ),
     }
 )

From 0d957e06bef17aecd675a625739bc76fe0aaf6ed Mon Sep 17 00:00:00 2001
From: Ombrini <91598680+Ombrini@users.noreply.github.com>
Date: Wed, 12 Aug 2026 18:09:49 +0200
Subject: [PATCH 16/20] Refactor GITT-EIS example and simulator: enhance
 documentation, remove operando EIS file, and update method names for clarity

---
 .../battery_parameterisation/gitt_eis.py      |  8 +--
 .../battery_parameterisation/operando_eis.py  | 59 -------------------
 pybop/pybamm/eis_simulator.py                 | 12 ++--
 pybop/pybamm/simulator.py                     |  4 --
 .../integration/models/test_grouped_models.py |  3 +
 tests/unit/test_pybamm_utils.py               |  3 +-
 tests/unit/test_simulator.py                  | 10 ++--
 7 files changed, 20 insertions(+), 79 deletions(-)
 delete mode 100644 examples/scripts/battery_parameterisation/operando_eis.py

diff --git a/examples/scripts/battery_parameterisation/gitt_eis.py b/examples/scripts/battery_parameterisation/gitt_eis.py
index c15563aba..cca09f6c6 100644
--- a/examples/scripts/battery_parameterisation/gitt_eis.py
+++ b/examples/scripts/battery_parameterisation/gitt_eis.py
@@ -4,8 +4,8 @@
 import pybop
 
 """
-Example demonstrating parameter estimation from a GITT experiment in which an EIS
-spectrum is acquired at the end of each pulse ("operando" EIS).
+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
@@ -15,8 +15,8 @@
 everywhere else.
 
 Diffusivity and the exchange-current density are fitted together, which is the pairing
-operando EIS is meant to separate: the relaxation after each pulse constrains transport,
-while the charge-transfer semicircle of the spectrum constrains kinetics.
+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
diff --git a/examples/scripts/battery_parameterisation/operando_eis.py b/examples/scripts/battery_parameterisation/operando_eis.py
deleted file mode 100644
index c8be85d99..000000000
--- a/examples/scripts/battery_parameterisation/operando_eis.py
+++ /dev/null
@@ -1,59 +0,0 @@
-import matplotlib.pyplot as plt
-import numpy as np
-import pybamm
-
-import pybop
-
-"""
-Example demonstrating EIS applied during operation (slow dis/charge) of the cell.
-"""
-
-# Define model and parameter values
-model = pybamm.lithium_ion.SPMe(
-    options={"surface form": "differential", "contact resistance": "true"}
-)
-parameter_values = pybamm.ParameterValues("Chen2020")
-parameter_values["Contact resistance [Ohm]"] = 0.02
-parameter_values.set_initial_state("2.85 V", options=model.options)
-
-# Set up and run a charge/discharge experiment
-C_rate = parameter_values["Nominal cell capacity [A.h]"]
-dataset = pybop.Dataset(
-    {
-        "Time [s]": np.asarray(
-            [0, 1, 1001, 2001, 3001, 3002, 3003, 4003, 5003, 6003, 6004]
-        ),
-        "Current [A]": np.asarray([0, -1, -1, -1, -1, 0, 1, 1, 1, 1, 0])
-        * C_rate
-        / 3,
-    }
-)
-
-sim = pybop.pybamm.Simulator(model, parameter_values=parameter_values, protocol=dataset)
-solution = sim.solve()
-solution.plot()
-
-# Set up and run the simulation
-n_frequency = 60
-solution = pybop.pybamm.EISSimulator(
-    model,
-    parameter_values=parameter_values,
-    f_eval=np.logspace(-4, 5, n_frequency),
-    protocol=dataset,
-).solve()
-
-fig, ax = plt.subplots()
-n_time_steps = len(solution["Time [s]"].data)
-for i in range(n_time_steps):
-    impedance = solution["Impedance"].data[i, :]
-    ax.plot(
-        np.real(impedance),
-        -np.imag(impedance),
-        "-" if i < n_time_steps / 2 else "--",
-        label=f"t={solution['Time [s]'].data[i]}s",
-    )
-ax.set(xlabel=r"$Z_r(\omega)$ [$\Omega$]", ylabel=r"$-Z_j(\omega)$ [$\Omega$]")
-ax.set_aspect("equal", "box")
-ax.legend()
-ax.set_ylim([0, ax.get_xlim()[1]])
-plt.show()
diff --git a/pybop/pybamm/eis_simulator.py b/pybop/pybamm/eis_simulator.py
index 49453cba6..03dfa7e95 100644
--- a/pybop/pybamm/eis_simulator.py
+++ b/pybop/pybamm/eis_simulator.py
@@ -44,8 +44,8 @@ class EISSimulator(BaseSimulator):
         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 ("operando" EIS). If None, a single spectrum is
-        computed about the initial state ("stationary" EIS).
+        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)
@@ -372,7 +372,7 @@ 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.
-        For an operando simulation, the time-domain trajectory is solved once and each
+        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
@@ -384,7 +384,7 @@ def _solve(self, inputs: "Inputs") -> Solution:
         -------
         Solution
             Complex impedance results, or the voltage and the real and imaginary
-            impedance components over the time domain for an operando simulation.
+            impedance components over the time domain when coupled to a protocol.
         """
         # Rebuild the model only if necessary, then set up the constant matrices
         self._model_rebuild(inputs)
@@ -400,9 +400,9 @@ def _solve(self, inputs: "Inputs") -> Solution:
             )
             return solution
 
-        return self._solve_operando(inputs)
+        return self._solve_along_protocol(inputs)
 
-    def _solve_operando(self, inputs: "Inputs") -> Solution:
+    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.
diff --git a/pybop/pybamm/simulator.py b/pybop/pybamm/simulator.py
index 728d73775..582a57225 100644
--- a/pybop/pybamm/simulator.py
+++ b/pybop/pybamm/simulator.py
@@ -463,10 +463,6 @@ def initial_state(self):
     def experiment(self):
         return self._experiment
 
-    @property
-    def time_data(self):
-        return self._t_interp if self._t_interp is not None else self._t_eval
-
     @property
     def solver(self):
         return self._solver
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_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 0ac3de1be..b60793bb2 100644
--- a/tests/unit/test_simulator.py
+++ b/tests/unit/test_simulator.py
@@ -81,9 +81,9 @@ def test_set_output_variables(self):
             simulator.set_output_variables(["Not a variable"])
 
 
-class TestOperandoEISSimulator:
+class TestCoupledEISSimulator:
     """
-    A class to test the operando mode of the pybamm.EISSimulator class.
+    A class to test the pybamm.EISSimulator class when coupled to a protocol.
     """
 
     pytestmark = pytest.mark.unit
@@ -157,7 +157,7 @@ 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."""
-        operando = pybop.pybamm.EISSimulator(
+        coupled = pybop.pybamm.EISSimulator(
             model,
             parameter_values=parameter_values,
             protocol=dataset,
@@ -169,8 +169,8 @@ def test_matches_stationary_at_initial_state(
 
         impedance = np.asarray(
             [
-                operando[impedance_variables[2 * j]].data[0]
-                + 1j * operando[impedance_variables[2 * j + 1]].data[0]
+                coupled[impedance_variables[2 * j]].data[0]
+                + 1j * coupled[impedance_variables[2 * j + 1]].data[0]
                 for j in range(len(frequencies))
             ]
         )

From 662c6e2b5be6207721498df1af3b7c0f6d7e3fbe Mon Sep 17 00:00:00 2001
From: "pre-commit-ci[bot]"
 <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Date: Wed, 12 Aug 2026 16:22:02 +0000
Subject: [PATCH 17/20] style: pre-commit fixes

---
 tests/unit/test_simulator.py | 4 +---
 1 file changed, 1 insertion(+), 3 deletions(-)

diff --git a/tests/unit/test_simulator.py b/tests/unit/test_simulator.py
index b60793bb2..df71f9e35 100644
--- a/tests/unit/test_simulator.py
+++ b/tests/unit/test_simulator.py
@@ -90,9 +90,7 @@ class TestCoupledEISSimulator:
 
     @pytest.fixture
     def model(self):
-        return pybamm.lithium_ion.SPM(
-            options={"surface form": "differential"}
-        )
+        return pybamm.lithium_ion.SPM(options={"surface form": "differential"})
 
     @pytest.fixture
     def parameter_values(self):

From 74af9cc03fbbd1fdcf244de4820087070ef39581 Mon Sep 17 00:00:00 2001
From: Ombrini <91598680+Ombrini@users.noreply.github.com>
Date: Wed, 12 Aug 2026 18:22:57 +0200
Subject: [PATCH 18/20] Update CHANGELOG.md

---
 CHANGELOG.md | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index f683a96db..c84cdf3da 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -2,6 +2,7 @@
 
 ## Features
 
+- [#973](https://github.com/pybop-team/PyBOP/pull/973) - 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
 
+- [#973](https://github.com/pybop-team/PyBOP/pull/973) - 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
 
+- [#973](https://github.com/pybop-team/PyBOP/pull/973) - 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
 
+- [#973](https://github.com/pybop-team/PyBOP/pull/973) - `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.

From 27441db39db1cf1e3b5cc9633962b7e52b639b23 Mon Sep 17 00:00:00 2001
From: Ombrini <91598680+Ombrini@users.noreply.github.com>
Date: Wed, 12 Aug 2026 18:23:10 +0200
Subject: [PATCH 19/20] Add CLAUDE.md

---
 CHANGELOG.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index c84cdf3da..d45614f8c 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -2,7 +2,7 @@
 
 ## Features
 
-- [#973](https://github.com/pybop-team/PyBOP/pull/973) - 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`.
+- [#975](https://github.com/pybop-team/PyBOP/pull/973) - 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.

From fffaabae03a4c98085665618b406f4398fc8af04 Mon Sep 17 00:00:00 2001
From: Ombrini <91598680+Ombrini@users.noreply.github.com>
Date: Wed, 12 Aug 2026 18:36:03 +0200
Subject: [PATCH 20/20] fix: correct pull request references in CHANGELOG.md

---
 CHANGELOG.md | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index d45614f8c..b90957881 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -2,7 +2,7 @@
 
 ## Features
 
-- [#975](https://github.com/pybop-team/PyBOP/pull/973) - 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`.
+- [#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.
@@ -13,21 +13,21 @@
 
 ## Optimisations
 
-- [#973](https://github.com/pybop-team/PyBOP/pull/973) - Set up the EIS solver, mass matrix and forcing vector once rather than on every evaluation.
+- [#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
 
-- [#973](https://github.com/pybop-team/PyBOP/pull/973) - 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.
+- [#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
 
-- [#973](https://github.com/pybop-team/PyBOP/pull/973) - `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.
+- [#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.