Feature: Add Positive and Negative EIS calculation

CHANGELOG.md

@@ -2,6 +2,7 @@
 
 ## Features
 
+- Added three-electrode EIS support to `pybamm.EISSimulation`, including automatic default reference-electrode insertion, named positive/negative electrode impedance outputs, and component-aware Nyquist plotting. ([#5648](https://github.com/pybamm-team/PyBaMM/pull/5648))
 - Added unstructured mesh support (`UnstructuredSubMesh`, generators, and interface coupling) for arbitrary 2D/3D domains. Hexahedra must have planar faces (warped hexes raise a `GeometryError`), and `UserSuppliedUnstructuredMesh` accepts tetrahedral, triangular, and quadrilateral cells only. ([#5687](https://github.com/pybamm-team/PyBaMM/pull/5687))
 - Generalised `VectorField` to N components and added `Component`/`Norm` operators for multi-dimensional vector fields. ([#5686](https://github.com/pybamm-team/PyBaMM/pull/5686))
 - Removed the left sidebar from the documentation home page for a cleaner landing experience. ([#5699](https://github.com/pybamm-team/PyBaMM/pull/5699))

examples/scripts/run_eis_three_electrode_simulation.py

@@ -0,0 +1,14 @@
+import numpy as np
+
+import pybamm
+
+model = pybamm.lithium_ion.SPM(options={"surface form": "differential"})
+
+frequencies = np.logspace(-4, 4, 50)
+
+eis_sim = pybamm.EISSimulation(
+    model,
+    three_electrodes=True,
+)
+eis_sim.solve(frequencies)
+eis_sim.nyquist_plot()

packages/pybamm/src/pybamm/simulation/eis_simulation.py

@@ -29,6 +29,10 @@ class EISSimulation(BaseSimulation):
         A dictionary of the types of spatial method to use on each domain.
     skip_surface_form_check : bool, optional
         If True, skip the 'surface form' model option validation. Defaults to False.
+    three_electrodes: bool, optional
+        Whether to use a three-electrode model. If True, the model must have a
+        reference electrode inserted. Default is False. If True, returns the impedance
+        of both the working and counter electrodes, as well as the cell impedance.
     """
 
     def __init__(
@@ -40,17 +44,28 @@ def __init__(
         var_pts=None,
         spatial_methods=None,
         skip_surface_form_check=False,
+        three_electrodes=False,
     ):
         timer = pybamm.Timer()
         model_name = model.name
 
         parameter_values = parameter_values or model.default_parameter_values
 
+        self.three_electrodes = three_electrodes
+        if three_electrodes and (
+            "Positive electrode 3E potential [V]" not in model.variables
+            or "Negative electrode 3E potential [V]" not in model.variables
+        ):
+            reference_position = (
+                parameter_values["Negative electrode thickness [m]"]
+                + parameter_values["Separator thickness [m]"] / 2
+            )
+            model.insert_reference_electrode(reference_position)
         # Validate required variables and surface form before any processing
-        self._validate_model_for_eis(model, skip_surface_form_check)
+        self._validate_model_for_eis(model, skip_surface_form_check, three_electrodes)
 
         pybamm.logger.info(f"Setting up {model_name} for EIS")
-        eis_model = self._set_up_model_for_eis(model)
+        eis_model = self._set_up_model_for_eis(model, three_electrodes)
 
         # Compute impedance scale factor after model transformation: the
         # external-circuit FunctionControl swap rescales "Current [A]" by
@@ -78,7 +93,9 @@ def __init__(
         pybamm.citations.register("Hallemans2025")
 
     @staticmethod
-    def _validate_model_for_eis(model, skip_surface_form_check=False):
+    def _validate_model_for_eis(
+        model, skip_surface_form_check=False, three_electrodes=False
+    ):
         """Validate that a model is suitable for frequency-domain EIS.
 
         Raises
@@ -87,6 +104,14 @@ def _validate_model_for_eis(model, skip_surface_form_check=False):
             If the model is missing required variables or options.
         """
         required_vars = ["Voltage [V]", "Current [A]"]
+        if three_electrodes:
+            required_vars.extend(
+                [
+                    "Positive electrode 3E potential [V]",
+                    "Negative electrode 3E potential [V]",
+                ]
+            )
+
         for var in required_vars:
             if var not in model.variables:
                 raise ValueError(
@@ -105,7 +130,7 @@ def _validate_model_for_eis(model, skip_surface_form_check=False):
             )
 
     @staticmethod
-    def _set_up_model_for_eis(model):
+    def _set_up_model_for_eis(model, three_electrodes=False):
         """Prepare a model for frequency-domain EIS.
 
         Creates a copy of the model with voltage and current as algebraic
@@ -132,6 +157,27 @@ def _set_up_model_for_eis(model):
         new_model.algebraic[V_cell] = V_cell - V
         new_model.initial_conditions[V_cell] = new_model.param.ocv_init
 
+        if three_electrodes:
+            V_pos_3e = new_model.variables["Positive electrode 3E potential [V]"]
+            V_pos_3e_var = pybamm.Variable(
+                "Positive electrode 3E potential variable [V]"
+            )
+            new_model.variables["Positive electrode 3E potential variable [V]"] = (
+                V_pos_3e_var
+            )
+            new_model.algebraic[V_pos_3e_var] = V_pos_3e_var - V_pos_3e
+            new_model.initial_conditions[V_pos_3e_var] = new_model.param.p.prim.U_init
+
+            V_neg_3e = new_model.variables["Negative electrode 3E potential [V]"]
+            V_neg_3e_var = pybamm.Variable(
+                "Negative electrode 3E potential variable [V]"
+            )
+            new_model.variables["Negative electrode 3E potential variable [V]"] = (
+                V_neg_3e_var
+            )
+            new_model.algebraic[V_neg_3e_var] = V_neg_3e_var - V_neg_3e
+            new_model.initial_conditions[V_neg_3e_var] = new_model.param.n.prim.U_init
+
         # Replace current with a FunctionControl variable
         external_circuit_variables = pybamm.external_circuit.FunctionControl(
             model.param, None, model.options, control="algebraic"
@@ -212,7 +258,7 @@ def _build_matrix_problem(self, inputs_dict=None):
         return self._cached_M, neg_J, self._cached_b
 
     @staticmethod
-    def _calculate_impedance(frequency, M, neg_J, b):
+    def _calculate_impedance(frequency, M, neg_J, b, three_electrodes=False):
         """Calculate impedance at a single frequency.
 
         Parameters
@@ -233,9 +279,15 @@ def _calculate_impedance(frequency, M, neg_J, b):
         """
         A = 1.0j * 2 * np.pi * frequency * M + neg_J
         x = spsolve(A, b)
-        # Voltage is penultimate, current is last (by construction in
-        # _set_up_model_for_eis)
-        return -x[-2] / x[-1]
+        if three_electrodes:
+            z_cell = -x[-4] / x[-1]
+            z_pos = -x[-3] / x[-1]
+            z_neg = -x[-2] / x[-1]
+            return z_cell, z_pos, z_neg
+        else:
+            # Voltage is penultimate, current is last (by construction in
+            # _set_up_model_for_eis)
+            return -x[-2] / x[-1]
 
     def solve(self, frequencies, inputs=None, initial_soc=None):
         """Compute impedance at the given frequencies.
@@ -269,9 +321,20 @@ def solve(self, frequencies, inputs=None, initial_soc=None):
 
         M, neg_J, b = self._build_matrix_problem(inputs_dict=inputs)
 
-        zs = [self._calculate_impedance(f, M, neg_J, b) for f in frequencies]
+        zs = [
+            self._calculate_impedance(f, M, neg_J, b, self.three_electrodes)
+            for f in frequencies
+        ]
         impedance = np.array(zs) * self._z_scale
-        self._solution = pybamm.EISSolution(frequencies, impedance)
+        if self.three_electrodes:
+            self._solution = pybamm.EISSolution(frequencies, impedance[:, 0])
+            self._solution._data["Cell impedance [Ohm]"] = impedance[:, 0]
+            self._solution._data["Positive electrode impedance [Ohm]"] = impedance[:, 1]
+            self._solution._data["Negative electrode impedance [Ohm]"] = -impedance[
+                :, 2
+            ]
+        else:
+            self._solution = pybamm.EISSolution(frequencies, impedance)
         self._solution.set_up_time = self.set_up_time
 
         self.solve_time = timer.time()

packages/pybamm/src/pybamm/solvers/solution.py

@@ -173,7 +173,35 @@ def nyquist_plot(self, **kwargs):
         """
         from pybamm.plotting.nyquist_plot import nyquist_plot
 
-        return nyquist_plot(self.impedance, **kwargs)
+        component_keys = {
+            "Cell": "Cell impedance [Ohm]",
+            "Positive electrode": "Positive electrode impedance [Ohm]",
+            "Negative electrode": "Negative electrode impedance [Ohm]",
+        }
+        if not all(key in self._data for key in component_keys.values()):
+            return nyquist_plot(self.impedance, **kwargs)
+
+        plot_kwargs = dict(kwargs)
+        ax = plot_kwargs.pop("ax", None)
+        show_plot = plot_kwargs.pop("show_plot", True)
+        plot_kwargs.pop("label", None)
+        fig = None
+        for label, key in component_keys.items():
+            fig_i, ax = nyquist_plot(
+                self._data[key],
+                ax=ax,
+                show_plot=False,
+                label=label,
+                **plot_kwargs,
+            )
+            fig = fig or fig_i
+        ax.legend()
+
+        if show_plot:  # pragma: no cover
+            plt = pybamm.import_optional_dependency("matplotlib.pyplot")
+            plt.show()
+
+        return fig, ax
 
 
 _DEFAULT_SOLUTION_OPTIONS = {

packages/pybamm/tests/unit/test_eis_simulation.py

@@ -29,6 +29,12 @@ def test_skip_surface_form_check(self):
             model, skip_surface_form_check=True
         )
 
+    def test_three_electrode_inserts_reference_electrode(self):
+        model = pybamm.lithium_ion.SPM(options={"surface form": "differential"})
+        pybamm.EISSimulation(model, three_electrodes=True)
+        assert "Positive electrode 3E potential [V]" in model.variables
+        assert "Negative electrode 3E potential [V]" in model.variables
+
 
 class TestEISSolution:
     """Tests for the EISSolution class."""
@@ -251,6 +257,23 @@ def test_high_freq_intercept_matches_contact_resistance(self):
         assert high_freq_z.real == pytest.approx(1.0, abs=1e-3)
         assert abs(high_freq_z.imag) < 1e-3
 
+    @pytest.mark.parametrize(
+        "model_class",
+        [pybamm.lithium_ion.SPM, pybamm.lithium_ion.SPMe, pybamm.lithium_ion.DFN],
+    )
+    def test_three_electrode_impedances_sum_to_cell_impedance(self, model_class):
+        model = model_class(options={"surface form": "differential"})
+        eis_sim = pybamm.EISSimulation(model, three_electrodes=True)
+        frequencies = np.logspace(-2, 2, 5)
+
+        result = eis_sim.solve(frequencies)
+        z_cell = result["Cell impedance [Ohm]"]
+        z_pos = result["Positive electrode impedance [Ohm]"]
+        z_neg = result["Negative electrode impedance [Ohm]"]
+
+        np.testing.assert_allclose(z_pos + z_neg, z_cell, rtol=1e-5, atol=1e-5)
+        np.testing.assert_allclose(result.impedance, z_cell)
+
 
 class TestNyquistPlot:
     """Tests for Nyquist plotting."""
@@ -290,6 +313,26 @@ def test_eis_nyquist_plot(self):
         assert fig is not None
         assert ax is not None
 
+    def test_eis_nyquist_plot_components(self):
+        import matplotlib
+
+        matplotlib.use("Agg")
+
+        impedance = np.array([1 + 0.5j, 2 + 1j, 3 + 1.5j])
+        solution = pybamm.EISSolution(np.array([1.0, 10.0, 100.0]), impedance)
+        solution._data["Cell impedance [Ohm]"] = impedance
+        solution._data["Positive electrode impedance [Ohm]"] = 0.4 * impedance
+        solution._data["Negative electrode impedance [Ohm]"] = 0.6 * impedance
+
+        fig, ax = solution.nyquist_plot(show_plot=False)
+
+        assert fig is not None
+        assert [line.get_label() for line in ax.get_lines()] == [
+            "Cell",
+            "Positive electrode",
+            "Negative electrode",
+        ]
+
     def test_nyquist_plot_before_solve_raises(self):
         model = pybamm.lithium_ion.SPM(options={"surface form": "differential"})
         eis_sim = pybamm.EISSimulation(model)