Repository: pybamm-team/PyBaMM

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

Discussion:
codecov[bot]: ## [Codecov](https://app.codecov.io/gh/pybamm-team/PyBaMM/pull/4621?dropdown=coverage&src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pybamm-team) Report
Attention: Patch coverage is `98.01980%` with `2 lines` in your changes missing coverage. Please review.
> Project coverage is 99.17%. Comparing base [(`403e1d3`)](https://app.codecov.io/gh/pybamm-team/PyBaMM/commit/403e1d31148770cde7bb55f5fc74463465a0edf2?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pybamm-team) to head [(`b02ea14`)](https://app.codecov.io/gh/pybamm-team/PyBaMM/commit/b02ea14dcd6de1fd110c8f8bfd8dbf71f01d129f?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pybamm-team).
> Report is 1 commits behind head on develop.

| [Files with missing lines](https://app.codecov.io/gh/pybamm-team/PyBaMM/pull/4621?dropdown=coverage&src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pybamm-team) | Patch % | Lines |
|---|---|---|
| [src/pybamm/solvers/summary\_variable.py](https://app.codecov.io/gh/pybamm-team/PyBaMM/pull/4621?src=pr&el=tree&filepath=src%2Fpybamm%2Fsolvers%2Fsummary_variable.py&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pybamm-team#diff-c3JjL3B5YmFtbS9zb2x2ZXJzL3N1bW1hcnlfdmFyaWFibGUucHk=) | 97.64% | [2 Missing :warning: ](https://app.codecov.io/gh/pybamm-team/PyBaMM/pull/4621?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pybamm-team) |

<details><summary>Additional details and impacted files</summary>


```diff
@@             Coverage Diff             @@
##           develop    #4621      +/-   ##
===========================================
- Coverage    99.18%   99.17%   -0.01%     
===========================================
  Files          302      303       +1     
  Lines        22840    22908      +68     
===========================================
+ Hits         22654    22720      +66     
- Misses         186      188       +2     
```

</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/pybamm-team/PyBaMM/pull/4621?dropdown=coverage&src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pybamm-team).   
:loudspeaker: Have feedback on the report? [Share it here](https://about.codecov.io/codecov-pr-comment-feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pybamm-team).


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)