Repository: pybamm-team/PyBaMM (python)

Issue:
Calculate summary variables only when they are called
As discussed in #3995 

- create a `SummaryVariable` object that saves `solution.last_state` and uses it to calculate the esoh variables only when they are called (the same way as `Solution` calculates the variables only when they are called)
- allow users to pass in functions that calculate summary variables based on other variables

Pull request:
Summary variables calculated only when called
# Description

Rather than calculating summary variables and storing them in every Solution class, a 'SummaryVariables' class is created which stores the first + last state of the solution, plus any relevant esoh parameters, so that the summary variables are only calculated and stored when they are called (similar to accessing standard variables). If an esoh variable is requested, they are all calculated and stored at once to reduce the number of calls to the esoh solver.

rather than viewing the list of summary variable names as
```
sol.summary_variables.keys()
```
you can now use
```
sol.summary_variables.all_variables
```
to view the names without calculating the variables.

Similarly, when plotting rather than getting the number of cycles as 
```
sol.summary_variables['Cycle number']
```
they are accessed as 
```
sol.summary_variables.cycle_number
```

Fixes #4058 

## Type of change

- [x] Optimization (back-end change that speeds up the code)

# Key checklist:

- [ ] No style issues: `$ pre-commit run` (or `$ nox -s pre-commit`) (see [CONTRIBUTING.md](https://github.com/pybamm-team/PyBaMM/blob/develop/CONTRIBUTING.md#installing-and-using-pre-commit) for how to set this up to run automatically when committing locally, in just two lines of code)
- [ ] All tests pass: `$ python run-tests.py --all` (or `$ nox -s tests`)
- [ ] The documentation builds: `$ python run-tests.py --doctest` (or `$ nox -s doctests`)

You can run integration tests, unit tests, and doctests together at once, using `$ python run-tests.py --quick` (or `$ nox -s quick`).

## Further checks:

- [ ] Code is commented, particularly in hard-to-understand areas
- [ ] Tests added that prove fix is effective or that feature works


Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index ebc582547..24c9be68e 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -8,6 +8,7 @@
 
 ## Breaking changes
 
+- Summary variables now calculated only when called, accessed via a class in the same manner as other variables rather than a dictionary. ([#4621](https://github.com/pybamm-team/PyBaMM/pull/4621))
 - The conda distribution (`pybamm`) now installs all optional dependencies available on conda-forge. Use the new `pybamm-base` conda
 package to install PyBaMM with only the required dependencies. ([conda-forge/pybamm-feedstock#70](https://github.com/conda-forge/pybamm-feedstock/pull/70))
 - Separated extrapolation options for `pybamm.BoundaryValue` and `pybamm.BoundaryGradient`, and updated the default to be "linear" for the value and "quadratic" for the gradient. ([#4614](https://github.com/pybamm-team/PyBaMM/pull/4614))
diff --git a/docs/source/api/solvers/index.rst b/docs/source/api/solvers/index.rst
index a9aa8ac1d..2bb15503b 100644
--- a/docs/source/api/solvers/index.rst
+++ b/docs/source/api/solvers/index.rst
@@ -13,3 +13,4 @@ Solvers
   algebraic_solvers
   solution
   processed_variable
+  summary_variables
diff --git a/docs/source/api/solvers/summary_variables.rst b/docs/source/api/solvers/summary_variables.rst
new file mode 100644
index 000000000..60344660c
--- /dev/null
+++ b/docs/source/api/solvers/summary_variables.rst
@@ -0,0 +1,5 @@
+Summary Variables
+======================
+
+.. autoclass:: pybamm.SummaryVariables
+  :members:
diff --git a/docs/source/examples/notebooks/simulations_and_experiments/simulating-long-experiments.ipynb b/docs/source/examples/notebooks/simulations_and_experiments/simulating-long-experiments.ipynb
index c7f1f0e63..3dcb46b1b 100644
--- a/docs/source/examples/notebooks/simulations_and_experiments/simulating-long-experiments.ipynb
+++ b/docs/source/examples/notebooks/simulations_and_experiments/simulating-long-experiments.ipynb
@@ -541,7 +541,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 7,
+   "execution_count": null,
    "id": "right-skiing",
    "metadata": {},
    "outputs": [
@@ -638,7 +638,7 @@
     }
    ],
    "source": [
-    "sorted(sol.summary_variables.keys())"
+    "sorted(sol.summary_variables.all_variables)"
    ]
   },
   {
@@ -1936,7 +1936,7 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Python 3 (ipykernel)",
+   "display_name": "venv",
    "language": "python",
    "name": "python3"
   },
@@ -1950,7 +1950,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.8.10"
+   "version": "3.11.10"
   },
   "toc": {
    "base_numbering": 1,
@@ -1964,11 +1964,6 @@
    "toc_position": {},
    "toc_section_display": true,
    "toc_window_display": true
-  },
-  "vscode": {
-   "interpreter": {
-    "hash": "612adcc456652826e82b485a1edaef831aa6d5abc680d008e93d513dd8724f14"
-   }
   }
  },
  "nbformat": 4,
diff --git a/src/pybamm/__init__.py b/src/pybamm/__init__.py
index b466c3896..df3d3d053 100644
--- a/src/pybamm/__init__.py
+++ b/src/pybamm/__init__.py
@@ -162,6 +162,7 @@ from .solvers.solution import Solution, EmptySolution, make_cycle_solution
 from .solvers.processed_variable_time_integral import ProcessedVariableTimeIntegral
 from .solvers.processed_variable import ProcessedVariable, process_variable
 from .solvers.processed_variable_computed import ProcessedVariableComputed
+from .solvers.summary_variable import SummaryVariables
 from .solvers.base_solver import BaseSolver
 from .solvers.dummy_solver import DummySolver
 from .solvers.algebraic_solver import AlgebraicSolver
diff --git a/src/pybamm/callbacks.py b/src/pybamm/callbacks.py
index 4e8c67c8b..29286ab5a 100644
--- a/src/pybamm/callbacks.py
+++ b/src/pybamm/callbacks.py
@@ -212,7 +212,7 @@ class LoggingCallback(Callback):
 
         voltage_stop = logs["stopping conditions"]["voltage"]
         if voltage_stop is not None:
-            min_voltage = logs["summary variables"]["Minimum voltage [V]"]
+            min_voltage = logs["Minimum voltage [V]"]
             if min_voltage > voltage_stop[0]:
                 self.logger.notice(
                     f"Minimum voltage is now {min_voltage:.3f} V "
diff --git a/src/pybamm/models/full_battery_models/lithium_ion/electrode_soh.py b/src/pybamm/models/full_battery_models/lithium_ion/electrode_soh.py
index a74391090..cdcb9a55f 100644
--- a/src/pybamm/models/full_battery_models/lithium_ion/electrode_soh.py
+++ b/src/pybamm/models/full_battery_models/lithium_ion/electrode_soh.py
@@ -312,6 +312,27 @@ class ElectrodeSOHSolver:
             self.__get_electrode_soh_sims_split
         )
 
+    def __getstate__(self):
+        """
+        Return dictionary of picklable items
+        """
+        result = self.__dict__.copy()
+        result["_get_electrode_soh_sims_full"] = None  # Exclude LRU cache
+        result["_get_electrode_soh_sims_split"] = None  # Exclude LRU cache
+        return result
+
+    def __setstate__(self, state):
+        """
+        Unpickle, restoring unpicklable relationships
+        """
+        self.__dict__ = state
+        self._get_electrode_soh_sims_full = lru_cache()(
+            self.__get_electrode_soh_sims_full
+        )
+        self._get_electrode_soh_sims_split = lru_cache()(
+            self.__get_electrode_soh_sims_split
+        )
+
     def _get_lims_ocp(self):
         parameter_values = self.parameter_values
 
diff --git a/src/pybamm/plotting/plot_summary_variables.py b/src/pybamm/plotting/plot_summary_variables.py
index bd4db0ee6..7181563fc 100644
--- a/src/pybamm/plotting/plot_summary_variables.py
+++ b/src/pybamm/plotting/plot_summary_variables.py
@@ -62,13 +62,13 @@ def plot_summary_variables(
         for solution in solutions:
             # plot summary variable v/s cycle number
             ax.plot(
-                solution.summary_variables["Cycle number"],
+                solution.summary_variables.cycle_number,
                 solution.summary_variables[var],
             )
         # label the axes
         ax.set_xlabel("Cycle number")
         ax.set_ylabel(var)
-        ax.set_xlim([1, solution.summary_variables["Cycle number"][-1]])
+        ax.set_xlim([1, solution.summary_variables.cycle_number[-1]])
 
     fig.tight_layout()
 
diff --git a/src/pybamm/simulation.py b/src/pybamm/simulation.py
index 9dafb23ce..d499f7740 100644
--- a/src/pybamm/simulation.py
+++ b/src/pybamm/simulation.py
@@ -920,7 +920,7 @@ class Simulation:
                 # See PR #3995
                 if voltage_stop is not None:
                     min_voltage = np.min(cycle_solution["Battery voltage [V]"].data)
-                    logs["summary variables"]["Minimum voltage [V]"] = min_voltage
+                    logs["Minimum voltage [V]"] = min_voltage
 
                 callbacks.on_cycle_end(logs)
 
@@ -941,7 +941,7 @@ class Simulation:
 
             if self._solution is not None and len(all_cycle_solutions) > 0:
                 self._solution.cycles = all_cycle_solutions
-                self._solution.set_summary_variables(all_summary_variables)
+                self._solution.update_summary_variables(all_summary_variables)
                 self._solution.all_first_states = all_first_states
 
             callbacks.on_experiment_end(logs)
diff --git a/src/pybamm/solvers/solution.py b/src/pybamm/solvers/solution.py
index 256d596fd..49e820983 100644
--- a/src/pybamm/solvers/solution.py
+++ b/src/pybamm/solvers/solution.py
@@ -563,16 +563,10 @@ class Solution:
         """Updates the initial start time of the experiment"""
         self._initial_start_time = value
 
-    def set_summary_variables(self, all_summary_variables):
-        summary_variables = {var: [] for var in all_summary_variables[0]}
-        for sum_vars in all_summary_variables:
-            for name, value in sum_vars.items():
-                summary_variables[name].append(value)
-
-        summary_variables["Cycle number"] = range(1, len(all_summary_variables) + 1)
+    def update_summary_variables(self, all_summary_variables):
         self.all_summary_variables = all_summary_variables
-        self._summary_variables = pybamm.FuzzyDict(
-            {name: np.array(value) for name, value in summary_variables.items()}
+        self._summary_variables = pybamm.SummaryVariables(
+            self, cycle_summary_variables=all_summary_variables
         )
 
     def update(self, variables):
@@ -1142,8 +1136,8 @@ def make_cycle_solution(
 
     cycle_solution.steps = step_solutions
 
-    cycle_summary_variables = _get_cycle_summary_variables(
-        cycle_solution, esoh_solver, user_inputs=inputs
+    cycle_summary_variables = pybamm.SummaryVariables(
+        cycle_solution, esoh_solver=esoh_solver, user_inputs=inputs
     )
 
     cycle_first_state = cycle_solution.first_state
@@ -1154,46 +1148,3 @@ def make_cycle_solution(
         cycle_solution = None
 
     return cycle_solution, cycle_summary_variables, cycle_first_state
-
-
-def _get_cycle_summary_variables(cycle_solution, esoh_solver, user_inputs=None):
-    user_inputs = user_inputs or {}
-    model = cycle_solution.all_models[0]
-    cycle_summary_variables = pybamm.FuzzyDict({})
-
-    # Summary variables
-    summary_variables = model.summary_variables
-    first_state = cycle_solution.first_state
-    last_state = cycle_solution.last_state
-    for var in summary_variables:
-        data_first = first_state[var].data
-        data_last = last_state[var].data
-        cycle_summary_variables[var] = data_last[0]
-        var_lowercase = var[0].lower() + var[1:]
-        cycle_summary_variables["Change in " + var_lowercase] = (
-            data_last[0] - data_first[0]
-        )
-
-    # eSOH variables (full-cell lithium-ion model only, for now)
-    if (
-        esoh_solver is not None
-        and isinstance(model, pybamm.lithium_ion.BaseModel)
-        and model.options.electrode_types["negative"] == "porous"
-        and "Negative electrode capacity [A.h]" in model.variables
-        and "Positive electrode capacity [A.h]" in model.variables
-    ):
-        Q_n = last_state["Negative electrode capacity [A.h]"].data[0]
-        Q_p = last_state["Positive electrode capacity [A.h]"].data[0]
-        Q_Li = last_state["Total lithium capacity in particles [A.h]"].data[0]
-        all_inputs = {**user_inputs, "Q_n": Q_n, "Q_p": Q_p, "Q_Li": Q_Li}
-        try:
-            esoh_sol = esoh_solver.solve(inputs=all_inputs)
-        except pybamm.SolverError as error:  # pragma: no cover
-            raise pybamm.SolverError(
-                "Could not solve for summary variables, run "
-                "`sim.solve(calc_esoh=False)` to skip this step"
-            ) from error
-
-        cycle_summary_variables.update(esoh_sol)
-
-    return cycle_summary_variables
diff --git a/src/pybamm/solvers/summary_variable.py b/src/pybamm/solvers/summary_variable.py
new file mode 100644
index 000000000..812dea72d
--- /dev/null
+++ b/src/pybamm/solvers/summary_variable.py
@@ -0,0 +1,188 @@
+#
+# Summary Variable class
+#
+from __future__ import annotations
+import pybamm
+import numpy as np
+from typing import Any
+
+
+class SummaryVariables:
+    """
+    Class for managing and calculating summary variables from a PyBaMM solution.
+    Summary variables are only calculated when simulations are run with PyBaMM
+    Experiments.
+
+    Parameters
+    ----------
+    solution : :class:`pybamm.Solution`
+        The solution object to be used for creating the processed variables.
+    cycle_summary_variables : list[pybamm.SummaryVariables], optional
+        A list of cycle summary variables.
+    esoh_solver : :class:`pybamm.lithium_ion.ElectrodeSOHSolver`, optional
+        Solver for electrode state-of-health (eSOH) calculations.
+    user_inputs : dict, optional
+        Additional user inputs for calculations.
+
+    Attributes
+    ----------
+    cycle_number : array[int]
+        Stores the cycle number for each saved cycle, for use when plotting.
+        Length is equal to the number of cycles in a solution.
+    """
+
+    def __init__(
+        self,
+        solution: pybamm.Solution,
+        cycle_summary_variables: list[SummaryVariables] | None = None,
+        esoh_solver: pybamm.lithium_ion.ElectrodeSOHSolver | None = None,
+        user_inputs: dict[str, Any] | None = None,
+    ):
+        self.user_inputs = user_inputs or {}
+        self.esoh_solver = esoh_solver
+        self._variables = {}  # Store computed variables
+        self.cycle_number = None
+
+        model = solution.all_models[0]
+        self._possible_variables = model.summary_variables  # minus esoh variables
+        self._esoh_variables = None  # Store eSOH variable names
+
+        # Flag if eSOH calculations are needed
+        self.calc_esoh = (
+            self.esoh_solver is not None
+            and isinstance(model, pybamm.lithium_ion.BaseModel)
+            and model.options.electrode_types["negative"] == "porous"
+            and "Negative electrode capacity [A.h]" in model.variables
+            and "Positive electrode capacity [A.h]" in model.variables
+        )
+
+        # Initialize based on cycle information
+        if cycle_summary_variables:
+            self._initialize_for_cycles(cycle_summary_variables)
+        else:
+            self.first_state = solution.first_state
+            self.last_state = solution.last_state
+            self.cycles = None
+
+    def _initialize_for_cycles(self, cycle_summary_variables: list[SummaryVariables]):
+        """Initialize attributes for when multiple cycles are provided."""
+        self.first_state = None
+        self.last_state = None
+        self.cycles = cycle_summary_variables
+        self.cycle_number = np.arange(1, len(self.cycles) + 1)
+        first_cycle = self.cycles[0]
+        self.calc_esoh = first_cycle.calc_esoh
+        self.esoh_solver = first_cycle.esoh_solver
+        self.user_inputs = first_cycle.user_inputs
+
+    @property
+    def all_variables(self) -> list[str]:
+        """
+        Return names of all possible summary variables, including eSOH variables
+         if appropriate.
+        """
+        try:
+            return self._all_variables
+        except AttributeError:
+            base_vars = self._possible_variables.copy()
+            base_vars.extend(
+                f"Change in {var[0].lower() + var[1:]}"
+                for var in self._possible_variables
+            )
+
+            if self.calc_esoh:
+                base_vars.extend(self.esoh_variables)
+
+            self._all_variables = base_vars
+            return self._all_variables
+
+    @property
+    def esoh_variables(self) -> list[str] | None:
+        """Return names of all eSOH variables."""
+        if self.calc_esoh and self._esoh_variables is None:
+            esoh_model = self.esoh_solver._get_electrode_soh_sims_full().model
+            esoh_vars = list(esoh_model.variables.keys())
+            self._esoh_variables = esoh_vars
+        return self._esoh_variables
+
+    def __getitem__(self, key: str) -> float | list[float]:
+        """
+        Access or compute a summary variable by its name.
+
+        Parameters
+        ----------
+        key : str
+            The name of the variable
+
+        Returns
+        -------
+        float or list[float]
+        """
+
+        if key in self._variables:
+            # return it if it exists
+            return self._variables[key]
+        elif key not in self.all_variables:
+            # check it's listed as a summary variable
+            raise KeyError(f"Variable '{key}' is not a summary variable.")
+        else:
+            # otherwise create it, save it and then return it
+            if self.calc_esoh and key in self._esoh_variables:
+                self.update_esoh()
+            else:
+                base_key = key.removeprefix("Change in ")
+                base_key = base_key[0].upper() + base_key[1:]
+                # this will create 'X' and 'Change in x' at the same time
+                self.update(base_key)
+            return self._variables[key]
+
+    def update(self, var: str):
+        """Compute and store a variable and its change."""
+        var_lowercase = var[0].lower() + var[1:]
+        if self.cycles:
+            self._update_multiple_cycles(var, var_lowercase)
+        else:
+            self._update(var, var_lowercase)
+
+    def _update_multiple_cycles(self, var: str, var_lowercase: str):
+        """Creates aggregated summary variables for where more than one cycle exists."""
+        var_cycle = [cycle[var] for cycle in self.cycles]
+        change_var_cycle = [
+            cycle[f"Change in {var_lowercase}"] for cycle in self.cycles
+        ]
+        self._variables[var] = var_cycle
+        self._variables[f"Change in {var_lowercase}"] = change_var_cycle
+
+    def _update(self, var: str, var_lowercase: str):
+        """Create variable `var` for a single cycle."""
+        data_first = self.first_state[var].data
+        data_last = self.last_state[var].data
+        self._variables[var] = data_last[0]
+        self._variables[f"Change in {var_lowercase}"] = data_last[0] - data_first[0]
+
+    def update_esoh(self):
+        """Create all aggregated eSOH variables"""
+        if self.cycles is not None:
+            var_cycle = [cycle._get_esoh_variables() for cycle in self.cycles]
+            aggregated_vars = {k: [] for k in var_cycle[0].keys()}
+            for cycle in var_cycle:
+                for k, v in cycle.items():
+                    aggregated_vars[k].append(v)
+            self._variables.update(aggregated_vars)
+        else:
+            self._variables.update(self._get_esoh_variables())
+
+    def _get_esoh_variables(self) -> dict[str, float]:
+        """Compute eSOH variables for a single solution."""
+        Q_n = self.last_state["Negative electrode capacity [A.h]"].data[0]
+        Q_p = self.last_state["Positive electrode capacity [A.h]"].data[0]
+        Q_Li = self.last_state["Total lithium capacity in particles [A.h]"].data[0]
+        all_inputs = {**self.user_inputs, "Q_n": Q_n, "Q_p": Q_p, "Q_Li": Q_Li}
+        try:
+            esoh_sol = self.esoh_solver.solve(inputs=all_inputs)
+        except pybamm.SolverError as error:  # pragma: no cover
+            raise pybamm.SolverError(
+                "Could not solve for eSOH summary variables"
+            ) from error
+
+        return esoh_sol


Test patch:
diff --git a/tests/integration/test_solvers/test_idaklu.py b/tests/integration/test_solvers/test_idaklu.py
index 3ee96a9cc..1766c1f2a 100644
--- a/tests/integration/test_solvers/test_idaklu.py
+++ b/tests/integration/test_solvers/test_idaklu.py
@@ -202,5 +202,5 @@ class TestIDAKLUSolver:
         )
 
         # check summary variables are the same if output variables are specified
-        for var in summary_vars[0].keys():
+        for var in model.summary_variables:
             assert summary_vars[0][var] == summary_vars[1][var]
diff --git a/tests/unit/test_plotting/test_plot_summary_variables.py b/tests/unit/test_plotting/test_plot_summary_variables.py
index 5f1a650ce..418eeaf71 100644
--- a/tests/unit/test_plotting/test_plot_summary_variables.py
+++ b/tests/unit/test_plotting/test_plot_summary_variables.py
@@ -45,9 +45,9 @@ class TestPlotSummaryVariables:
 
             cycle_number, var = ax.get_lines()[0].get_data()
             np.testing.assert_array_equal(
-                cycle_number, sol.summary_variables["Cycle number"]
+                cycle_number, sol.summary_variables.cycle_number
             )
-            np.testing.assert_array_equal(var, sol.summary_variables[output_var])
+            np.testing.assert_allclose(sol.summary_variables[output_var], var)
 
         axes = pybamm.plot_summary_variables(
             [sol, sol], labels=["SPM", "SPM"], show_plot=False
@@ -62,12 +62,12 @@ class TestPlotSummaryVariables:
 
             cycle_number, var = ax.get_lines()[0].get_data()
             np.testing.assert_array_equal(
-                cycle_number, sol.summary_variables["Cycle number"]
+                cycle_number, sol.summary_variables.cycle_number
             )
-            np.testing.assert_array_equal(var, sol.summary_variables[output_var])
+            np.testing.assert_allclose(sol.summary_variables[output_var], var)
 
             cycle_number, var = ax.get_lines()[1].get_data()
             np.testing.assert_array_equal(
-                cycle_number, sol.summary_variables["Cycle number"]
+                cycle_number, sol.summary_variables.cycle_number
             )
-            np.testing.assert_array_equal(var, sol.summary_variables[output_var])
+            np.testing.assert_allclose(sol.summary_variables[output_var], var)
diff --git a/tests/unit/test_solvers/test_summary_variables.py b/tests/unit/test_solvers/test_summary_variables.py
new file mode 100644
index 000000000..0132c42ea
--- /dev/null
+++ b/tests/unit/test_solvers/test_summary_variables.py
@@ -0,0 +1,171 @@
+#
+# Tests for the Summary Variables class
+#
+
+import pybamm
+import numpy as np
+import pytest
+
+
+class TestSummaryVariables:
+    @staticmethod
+    def create_sum_vars():
+        model = pybamm.BaseModel()
+        c = pybamm.Variable("c")
+        model.rhs = {c: -c}
+        model.initial_conditions = {c: 1}
+        model.variables["c"] = c
+        model.variables["2c"] = 2 * c
+        model.summary_variables = ["2c"]
+
+        solution = pybamm.ScipySolver().solve(model, np.linspace(0, 1))
+
+        sum_vars = pybamm.SummaryVariables(solution)
+
+        return sum_vars, solution
+
+    def test_init(self):
+        model = pybamm.BaseModel()
+        c = pybamm.Variable("c")
+        model.rhs = {c: -c}
+        model.initial_conditions = {c: 1}
+        model.variables["c"] = c
+        model.variables["2c"] = 2 * c
+        model.summary_variables = ["2c"]
+
+        solution = pybamm.ScipySolver().solve(model, np.linspace(0, 1))
+
+        sum_vars = pybamm.SummaryVariables(solution)
+
+        # no variables should have been calculated until called
+        assert sum_vars._variables == {}
+
+        assert sum_vars.first_state == solution.first_state
+        assert sum_vars.last_state == solution.last_state
+        assert sum_vars.cycles is None
+
+    def test_init_with_cycle_summary_variables(self):
+        model = pybamm.BaseModel()
+        c = pybamm.Variable("c")
+        model.rhs = {c: -c}
+        model.initial_conditions = {c: 1}
+        model.variables["c"] = c
+        model.variables["2c"] = 2 * c
+        model.summary_variables = ["2c"]
+
+        sol1 = pybamm.ScipySolver().solve(model, np.linspace(0, 1))
+        sol2 = pybamm.ScipySolver().solve(model, np.linspace(1, 2))
+        sol3 = pybamm.ScipySolver().solve(model, np.linspace(2, 3))
+
+        cycle_sol = sol1 + sol2 + sol3
+
+        all_sum_vars = [
+            pybamm.SummaryVariables(sol1),
+            pybamm.SummaryVariables(sol2),
+            pybamm.SummaryVariables(sol3),
+        ]
+
+        cycle_sum_vars = pybamm.SummaryVariables(
+            cycle_sol,
+            cycle_summary_variables=all_sum_vars,
+        )
+
+        assert cycle_sum_vars.first_state is None
+        assert cycle_sum_vars.last_state is None
+        assert cycle_sum_vars._variables == {}
+        assert cycle_sum_vars.cycles == all_sum_vars
+        np.testing.assert_array_equal(cycle_sum_vars.cycle_number, np.array([1, 2, 3]))
+
+    def test_get_variable(self):
+        sum_vars, solution = self.create_sum_vars()
+
+        summary_c = sum_vars["2c"]
+
+        assert summary_c == solution["2c"].data[-1]
+        assert list(sum_vars._variables.keys()) == ["2c", "Change in 2c"]
+
+    def test_get_esoh_variable(self):
+        model = pybamm.lithium_ion.SPM()
+        sim = pybamm.Simulation(model)
+        sol = sim.solve(np.linspace(0, 1))
+        esoh_solver = sim.get_esoh_solver(True)
+        sum_vars_esoh = pybamm.SummaryVariables(sol, esoh_solver=esoh_solver)
+
+        assert np.isclose(sum_vars_esoh["x_100"], 0.9493)
+
+        # all esoh vars should be calculated at the same time to reduce solver calls
+        assert "Practical NPR" in sum_vars_esoh._variables
+
+    def test_get_variable_error_not_summary_variable(self):
+        sum_vars, _ = self.create_sum_vars()
+
+        with pytest.raises(KeyError, match="Variable 'c' is not a summary variable"):
+            sum_vars["c"]
+
+    def test_summary_vars_all_variables(self):
+        # no esoh
+        sum_vars, _ = self.create_sum_vars()
+
+        assert sum_vars.all_variables == ["2c", "Change in 2c"]
+
+        # with esoh
+        model = pybamm.lithium_ion.SPM()
+        sim = pybamm.Simulation(model)
+        sol = sim.solve(np.linspace(0, 1))
+        esoh_solver = sim.get_esoh_solver(True)
+        sum_vars_esoh = pybamm.SummaryVariables(sol, esoh_solver=esoh_solver)
+
+        assert sum_vars_esoh.calc_esoh is True
+        assert "Total lithium lost [mol]" in sum_vars_esoh.all_variables
+        assert "x_100" in sum_vars_esoh.all_variables
+
+        assert "x_100" in sum_vars_esoh._esoh_variables
+
+    def test_get_with_cycle_summary_variables(self):
+        model = pybamm.BaseModel()
+        c = pybamm.Variable("c")
+        model.rhs = {c: -c}
+        model.initial_conditions = {c: 1}
+        model.variables["c"] = c
+        model.variables["2c"] = 2 * c
+        model.summary_variables = ["2c"]
+
+        sol1 = pybamm.ScipySolver().solve(model, np.linspace(0, 1))
+        sol2 = pybamm.ScipySolver().solve(model, np.linspace(1, 2))
+        sol3 = pybamm.ScipySolver().solve(model, np.linspace(2, 3))
+
+        cycle_sol = sol1 + sol2 + sol3
+
+        all_sum_vars = [
+            pybamm.SummaryVariables(sol1),
+            pybamm.SummaryVariables(sol2),
+            pybamm.SummaryVariables(sol3),
+        ]
+
+        cycle_sum_vars = pybamm.SummaryVariables(
+            cycle_sol,
+            cycle_summary_variables=all_sum_vars,
+        )
+
+        np.testing.assert_array_equal(cycle_sum_vars.cycle_number, np.array([1, 2, 3]))
+        np.testing.assert_allclose(
+            cycle_sum_vars["2c"], np.array([0.735758, 0.735758, 0.735758])
+        )
+
+    def test_get_esoh_cycle_summary_vars(self):
+        experiment = pybamm.Experiment(
+            [
+                (
+                    "Discharge at 1C for 1 sec",
+                    "Charge at 1C for 1 sec",
+                ),
+            ]
+            * 10,
+        )
+        model = pybamm.lithium_ion.SPM()
+        sim = pybamm.Simulation(model, experiment=experiment)
+        sol = sim.solve()
+
+        assert len(sol.summary_variables.cycles) == 10
+        assert len(sol.summary_variables["x_100"]) == 10
+        assert np.isclose(sol.summary_variables["x_100"][0], 0.9493)