Repository: pybamm-team/PyBaMM (python)

Issue:
add `DiscreteTimeSum` subclass of `pybamm.Symbol`
### Description

This would add a new unary operator in a pybamm expression tree that would represent a discrete sum over time. This would be a discrete version of the already existing `pybamm.ExplicitTimeIntegral` class.

### Motivation

This would give model developers and users the ability to add "sum of squares" type variables to a pybamm model that would calculate the difference between, for example, a solution variable and a user-provided dataset. This would be useful for implementing parameter inference using pybamm



### Possible Implementation

Similar to the `ExplicitTimeIntegral`, this would only be evaluated in the `Solution` class (see `_update_variable` function). I would propose that the sum is done over the list of time points in the solution (rather than have a separate list of points in the expression tree node), that way a user can specify the time-points via `t_interp`.

I'm a bit unsure how to allow users to provide their data to the expression. Perhaps they would wrap it in an Interpolant:

```python
data = pybamm.Interpolant(data_t, data_y, pybamm.t)
model.variable["data_comparison"] = pybamm.DiscreteTimeSum((model.variable["Voltage"] - data)**2)
```

### Additional context

see https://github.com/pybop-team/PyBOP/issues/513

Pull request description:
feat: add discrete time sum expression tree node
# Description

Adds two new subclasses to `pybamm.Symbol`:
1. `pybamm.DiscreteTimeData` which is an interpolant for holding time series data
2. `pybamm.DiscreteTimeSum` which is a unary operator representing a sum of its child evaluation over a set of data time points. Its child must have exactly one `pybamm.DiscreteTimeData` object. 

When `pybamm.DiscreteTimeSum` is evaluated normally it returns the evaluation of its child. If it is evaluated as part of a `pybamm.Solution` then it interpolates this solution onto the data time points and sums these over the time axis

This also updates the evaluation of `pybamm.ExplicitTimeIntegral`. Previously this performed a trapezoidal integral over all the time points passed to the `__call__` function of `pybamm.ProcessedVariable`. Now it performs the integral over the solution time points, which is more consistent with its definition and always returns the same value irrespective of the time values that the user passes to `__call__`. In the future this integral should be updated to use the gradient information of the solution if available.

Fixes #4485

## Type of change

Please add a line in the relevant section of [CHANGELOG.md](https://github.com/pybamm-team/PyBaMM/blob/develop/CHANGELOG.md) to document the change (include PR #) - note reverse order of PR #s. If necessary, also add to the list of breaking changes.

- [x] New feature (non-breaking change which adds functionality)

# Key checklist:

- [x] 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)
- [x] All tests pass: `$ python run-tests.py --all` (or `$ nox -s tests`)
- [x] 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:

- [x] Code is commented, particularly in hard-to-understand areas
- [x] Tests added that prove fix is effective or that feature works


Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 8ffc210a7..fd69942fd 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -9,6 +9,7 @@
 - Added phase-dependent particle options to LAM
   ([#4369](https://github.com/pybamm-team/PyBaMM/pull/4369))
 - Added a lithium ion equivalent circuit model with split open circuit voltages for each electrode (`SplitOCVR`). ([#4330](https://github.com/pybamm-team/PyBaMM/pull/4330))
+- Added the `pybamm.DiscreteTimeSum` expression node to sum an expression over a sequence of data times, and accompanying `pybamm.DiscreteTimeData` class to store the data times and values ([#4501](https://github.com/pybamm-team/PyBaMM/pull/4501))
 
 ## Optimizations
 
diff --git a/src/pybamm/__init__.py b/src/pybamm/__init__.py
index f51ba05d0..2deb30f30 100644
--- a/src/pybamm/__init__.py
+++ b/src/pybamm/__init__.py
@@ -36,6 +36,7 @@ from .expression_tree.averages import _BaseAverage
 from .expression_tree.broadcasts import *
 from .expression_tree.functions import *
 from .expression_tree.interpolant import Interpolant
+from .expression_tree.discrete_time_sum import *
 from .expression_tree.input_parameter import InputParameter
 from .expression_tree.parameter import Parameter, FunctionParameter
 from .expression_tree.scalar import Scalar
@@ -158,6 +159,7 @@ from .spatial_methods.scikit_finite_element import ScikitFiniteElement
 
 # Solver classes
 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.base_solver import BaseSolver
diff --git a/src/pybamm/expression_tree/__init__.py b/src/pybamm/expression_tree/__init__.py
index 0b06746e6..7ac80e535 100644
--- a/src/pybamm/expression_tree/__init__.py
+++ b/src/pybamm/expression_tree/__init__.py
@@ -2,4 +2,4 @@ __all__ = ['array', 'averages', 'binary_operators', 'broadcasts',
            'concatenations', 'exceptions', 'functions', 'independent_variable',
            'input_parameter', 'interpolant', 'matrix', 'operations',
            'parameter', 'printing', 'scalar', 'state_vector', 'symbol',
-           'unary_operators', 'variable', 'vector']
+           'unary_operators', 'variable', 'vector', 'discrete_time_sum' ]
diff --git a/src/pybamm/expression_tree/discrete_time_sum.py b/src/pybamm/expression_tree/discrete_time_sum.py
new file mode 100644
index 000000000..41cd14960
--- /dev/null
+++ b/src/pybamm/expression_tree/discrete_time_sum.py
@@ -0,0 +1,88 @@
+import pybamm
+import numpy as np
+
+
+class DiscreteTimeData(pybamm.Interpolant):
+    """
+    A class for representing data that is only defined at discrete points in time.
+    This is implemented as a 1D interpolant with the time points as the nodes.
+
+    Parameters
+    ----------
+
+    time_points : :class:`numpy.ndarray`
+        The time points at which the data is defined
+    data : :class:`numpy.ndarray`
+        The data to be interpolated
+    name : str
+        The name of the data
+
+    """
+
+    def __init__(self, time_points: np.ndarray, data: np.ndarray, name: str):
+        super().__init__(time_points, data, pybamm.t, name)
+
+    def create_copy(self, new_children=None, perform_simplifications=True):
+        """See :meth:`pybamm.Symbol.new_copy()`."""
+        return pybamm.DiscreteTimeData(self.x[0], self.y, self.name)
+
+
+class DiscreteTimeSum(pybamm.UnaryOperator):
+    """
+    A node in the expression tree representing a discrete time sum operator.
+
+    .. math::
+        \\sum_{i=0}^{N} f(y(t_i), t_i)
+
+    where f is the expression given by the child, and the sum is over the discrete
+    time points t_i. The set of time points is given by the :class:`pybamm.DiscreteTimeData` node,
+    which must be somewhere in the expression tree given by the child. If the child
+    does not contain a :class:`pybamm.DiscreteTimeData` node, then an error will be raised when
+    the node is created. If the child contains multiple :class:`pybamm.DiscreteTimeData` nodes,
+    an error will be raised when the node is created.
+
+
+    Parameters
+    ----------
+    child: :class:`pybamm.Symbol`
+        The symbol to be summed
+
+    Attributes
+    ----------
+    data: :class:`pybamm.DiscreteTimeData`
+        The discrete time data node in the child
+
+    Raises
+    ------
+    :class:`pybamm.ModelError`
+        If the child does not contain a :class:`pybamm.DiscreteTimeData` node, or if the child
+        contains multiple :class:`pybamm.DiscreteTimeData` nodes.
+    """
+
+    def __init__(self, child: pybamm.Symbol):
+        self.data = None
+        for node in child.pre_order():
+            if isinstance(node, DiscreteTimeData):
+                # Check that there is exactly one DiscreteTimeData node in the child
+                if self.data is not None:
+                    raise pybamm.ModelError(
+                        "DiscreteTimeSum can only have one DiscreteTimeData node in the child"
+                    )
+                self.data = node
+        if self.data is None:
+            raise pybamm.ModelError(
+                "DiscreteTimeSum must contain a DiscreteTimeData node"
+            )
+        super().__init__("discrete time sum", child)
+
+    @property
+    def sum_values(self):
+        return self.data.y
+
+    @property
+    def sum_times(self):
+        return self.data.x[0]
+
+    def _unary_evaluate(self, child):
+        # return result of evaluating the child, we'll only implement the sum once the model is solved (in pybamm.ProcessedVariable)
+        return child
diff --git a/src/pybamm/solvers/__init__.py b/src/pybamm/solvers/__init__.py
index fc8be7e2f..e9d9d306f 100644
--- a/src/pybamm/solvers/__init__.py
+++ b/src/pybamm/solvers/__init__.py
@@ -2,4 +2,4 @@ __all__ = ['algebraic_solver', 'base_solver', 'c_solvers',
            'casadi_algebraic_solver', 'casadi_solver', 'dummy_solver',
            'idaklu_jax', 'idaklu_solver', 'jax_bdf_solver', 'jax_solver',
            'lrudict', 'processed_variable', 'processed_variable_computed',
-           'scipy_solver', 'solution']
+           'scipy_solver', 'solution', 'processed_variable_time_integral']
diff --git a/src/pybamm/solvers/processed_variable.py b/src/pybamm/solvers/processed_variable.py
index 5cf928ca7..12cf140b3 100644
--- a/src/pybamm/solvers/processed_variable.py
+++ b/src/pybamm/solvers/processed_variable.py
@@ -1,6 +1,7 @@
 #
 # Processed Variable class
 #
+from typing import Optional
 import casadi
 import numpy as np
 import pybamm
@@ -29,6 +30,8 @@ class ProcessedVariable:
         `base_Variable.evaluate` (but more efficiently).
     solution : :class:`pybamm.Solution`
         The solution object to be used to create the processed variables
+    time_integral : :class:`pybamm.ProcessedVariableTimeIntegral`, optional
+        Not none if the variable is to be time-integrated (default is None)
     """
 
     def __init__(
@@ -36,7 +39,7 @@ class ProcessedVariable:
         base_variables,
         base_variables_casadi,
         solution,
-        cumtrapz_ic=None,
+        time_integral: Optional[pybamm.ProcessedVariableTimeIntegral] = None,
     ):
         self.base_variables = base_variables
         self.base_variables_casadi = base_variables_casadi
@@ -50,7 +53,7 @@ class ProcessedVariable:
         self.mesh = base_variables[0].mesh
         self.domain = base_variables[0].domain
         self.domains = base_variables[0].domains
-        self.cumtrapz_ic = cumtrapz_ic
+        self.time_integral = time_integral
 
         # Process spatial variables
         geometry = solution.all_models[0].geometry
@@ -271,18 +274,21 @@ class ProcessedVariable:
                     self._coords_raw,
                 )
 
-            processed_entries = self._xr_interpolate(
-                entries_for_interp,
-                coords,
-                observe_raw,
-                t,
-                x,
-                r,
-                y,
-                z,
-                R,
-                fill_value,
-            )
+            if self.time_integral is None:
+                processed_entries = self._xr_interpolate(
+                    entries_for_interp,
+                    coords,
+                    observe_raw,
+                    t,
+                    x,
+                    r,
+                    y,
+                    z,
+                    R,
+                    fill_value,
+                )
+            else:
+                processed_entries = entries_for_interp
         else:
             processed_entries = entries
 
@@ -343,6 +349,16 @@ class ProcessedVariable:
             t_observe (np.ndarray): time points to observe
             observe_raw (bool): True if observing the raw data
         """
+        # if this is a time integral variable, t must be None and we observe either the
+        # data times (for a discrete sum) or the solution times (for a continuous sum)
+        if self.time_integral is not None:
+            if self.time_integral.method == "discrete":
+                # discrete sum should be observed at the discrete times
+                t = self.time_integral.discrete_times
+            else:
+                # assume we can do a sufficiently accurate trapezoidal integration at t_pts
+                t = self.t_pts
+
         observe_raw = (t is None) or (
             np.asarray(t).size == len(self.t_pts) and np.all(t == self.t_pts)
         )
@@ -483,14 +499,14 @@ class ProcessedVariable0D(ProcessedVariable):
         base_variables,
         base_variables_casadi,
         solution,
-        cumtrapz_ic=None,
+        time_integral: Optional[pybamm.ProcessedVariableTimeIntegral] = None,
     ):
         self.dimensions = 0
         super().__init__(
             base_variables,
             base_variables_casadi,
             solution,
-            cumtrapz_ic=cumtrapz_ic,
+            time_integral=time_integral,
         )
 
     def _observe_raw_python(self):
@@ -510,13 +526,19 @@ class ProcessedVariable0D(ProcessedVariable):
                 idx += 1
         return entries
 
-    def _observe_postfix(self, entries, _):
-        if self.cumtrapz_ic is None:
+    def _observe_postfix(self, entries, t):
+        if self.time_integral is None:
             return entries
-
-        return cumulative_trapezoid(
-            entries, self.t_pts, initial=float(self.cumtrapz_ic)
-        )
+        if self.time_integral.method == "discrete":
+            return np.sum(entries, axis=0, initial=self.time_integral.initial_condition)
+        elif self.time_integral.method == "continuous":
+            return cumulative_trapezoid(
+                entries, self.t_pts, initial=float(self.time_integral.initial_condition)
+            )
+        else:
+            raise ValueError(
+                "time_integral method must be 'discrete' or 'continuous'"
+            )  # pragma: no cover
 
     def _interp_setup(self, entries, t):
         # save attributes for interpolation
@@ -556,14 +578,14 @@ class ProcessedVariable1D(ProcessedVariable):
         base_variables,
         base_variables_casadi,
         solution,
-        cumtrapz_ic=None,
+        time_integral: Optional[pybamm.ProcessedVariableTimeIntegral] = None,
     ):
         self.dimensions = 1
         super().__init__(
             base_variables,
             base_variables_casadi,
             solution,
-            cumtrapz_ic=cumtrapz_ic,
+            time_integral=time_integral,
         )
 
     def _observe_raw_python(self):
@@ -653,14 +675,14 @@ class ProcessedVariable2D(ProcessedVariable):
         base_variables,
         base_variables_casadi,
         solution,
-        cumtrapz_ic=None,
+        time_integral: Optional[pybamm.ProcessedVariableTimeIntegral] = None,
     ):
         self.dimensions = 2
         super().__init__(
             base_variables,
             base_variables_casadi,
             solution,
-            cumtrapz_ic=cumtrapz_ic,
+            time_integral=time_integral,
         )
         first_dim_nodes = self.mesh.nodes
         first_dim_edges = self.mesh.edges
@@ -819,14 +841,14 @@ class ProcessedVariable2DSciKitFEM(ProcessedVariable2D):
         base_variables,
         base_variables_casadi,
         solution,
-        cumtrapz_ic=None,
+        time_integral: Optional[pybamm.ProcessedVariableTimeIntegral] = None,
     ):
         self.dimensions = 2
         super(ProcessedVariable2D, self).__init__(
             base_variables,
             base_variables_casadi,
             solution,
-            cumtrapz_ic=cumtrapz_ic,
+            time_integral=time_integral,
         )
         y_sol = self.mesh.edges["y"]
         z_sol = self.mesh.edges["z"]
diff --git a/src/pybamm/solvers/processed_variable_time_integral.py b/src/pybamm/solvers/processed_variable_time_integral.py
new file mode 100644
index 000000000..4fcdfb56b
--- /dev/null
+++ b/src/pybamm/solvers/processed_variable_time_integral.py
@@ -0,0 +1,28 @@
+from dataclasses import dataclass
+from typing import Literal, Optional, Union
+import numpy as np
+import pybamm
+
+
+@dataclass
+class ProcessedVariableTimeIntegral:
+    method: Literal["discrete", "continuous"]
+    initial_condition: np.ndarray
+    discrete_times: Optional[np.ndarray]
+
+    @staticmethod
+    def from_pybamm_var(
+        var: Union[pybamm.DiscreteTimeSum, pybamm.ExplicitTimeIntegral],
+    ) -> "ProcessedVariableTimeIntegral":
+        if isinstance(var, pybamm.DiscreteTimeSum):
+            return ProcessedVariableTimeIntegral(
+                method="discrete", initial_condition=0.0, discrete_times=var.sum_times
+            )
+        elif isinstance(var, pybamm.ExplicitTimeIntegral):
+            return ProcessedVariableTimeIntegral(
+                method="continuous",
+                initial_condition=var.initial_condition.evaluate(),
+                discrete_times=None,
+            )
+        else:
+            raise ValueError("Unsupported variable type")  # pragma: no cover
diff --git a/src/pybamm/solvers/solution.py b/src/pybamm/solvers/solution.py
index c884e79e3..1aa540ab3 100644
--- a/src/pybamm/solvers/solution.py
+++ b/src/pybamm/solvers/solution.py
@@ -571,7 +571,7 @@ class Solution:
             self._update_variable(variable)
 
     def _update_variable(self, variable):
-        cumtrapz_ic = None
+        time_integral = None
         pybamm.logger.debug(f"Post-processing {variable}")
         vars_pybamm = [
             model.variables_and_events[variable] for model in self.all_models
@@ -591,16 +591,22 @@ class Solution:
                     "solve. Please re-run the solve with `output_variables` set to "
                     "include this variable."
                 )
-            elif isinstance(var_pybamm, pybamm.ExplicitTimeIntegral):
-                cumtrapz_ic = var_pybamm.initial_condition
-                cumtrapz_ic = cumtrapz_ic.evaluate()
-                var_pybamm = var_pybamm.child
-                var_casadi = self.process_casadi_var(
-                    var_pybamm,
-                    inputs,
-                    ys.shape,
+            elif isinstance(
+                var_pybamm, (pybamm.ExplicitTimeIntegral, pybamm.DiscreteTimeSum)
+            ):
+                time_integral = pybamm.ProcessedVariableTimeIntegral.from_pybamm_var(
+                    var_pybamm
                 )
-                model._variables_casadi[variable] = var_casadi
+                var_pybamm = var_pybamm.child
+                if variable in model._variables_casadi:
+                    var_casadi = model._variables_casadi[variable]
+                else:
+                    var_casadi = self.process_casadi_var(
+                        var_pybamm,
+                        inputs,
+                        ys.shape,
+                    )
+                    model._variables_casadi[variable] = var_casadi
                 vars_pybamm[i] = var_pybamm
             elif variable in model._variables_casadi:
                 var_casadi = model._variables_casadi[variable]
@@ -613,7 +619,7 @@ class Solution:
                 model._variables_casadi[variable] = var_casadi
             vars_casadi.append(var_casadi)
         var = pybamm.process_variable(
-            vars_pybamm, vars_casadi, self, cumtrapz_ic=cumtrapz_ic
+            vars_pybamm, vars_casadi, self, time_integral=time_integral
         )
 
         self._variables[variable] = var


Test patch:
diff --git a/tests/unit/test_expression_tree/test_unary_operators.py b/tests/unit/test_expression_tree/test_unary_operators.py
index 812cbd8f6..0dbafa38c 100644
--- a/tests/unit/test_expression_tree/test_unary_operators.py
+++ b/tests/unit/test_expression_tree/test_unary_operators.py
@@ -724,6 +724,36 @@ class TestUnaryOperators:
         assert expr.create_copy() == expr
         assert not expr.is_constant()
 
+    def test_discrete_time_sum(self):
+        times = np.array([1, 2, 3, 4, 5])
+        values = np.array([2, 2, 3, 3, 1])
+        data = pybamm.DiscreteTimeData(times, values, "test")
+        assert data.name == "test"
+        np.testing.assert_array_equal(data.x[0], times)
+        np.testing.assert_array_equal(data.y, values)
+
+        y = pybamm.StateVector(slice(0, 1))
+
+        # check that raises error if data is not present
+        with pytest.raises(pybamm.ModelError, match="must contain a DiscreteTimeData"):
+            pybamm.DiscreteTimeSum(2 * y)
+
+        # check that raises error if two data are present
+        data2 = pybamm.DiscreteTimeData(values, times, "test2")
+        with pytest.raises(pybamm.ModelError, match="only have one DiscreteTimeData"):
+            pybamm.DiscreteTimeSum(data + data2)
+
+        sum = pybamm.DiscreteTimeSum(2 * data - y)
+        np.testing.assert_array_equal(sum.sum_times, times)
+        np.testing.assert_array_equal(sum.sum_values, values)
+        y = np.array([1])
+
+        # evaluate should return the values to sum up
+        for i in range(len(times)):
+            eval = sum.evaluate(y=y, t=times[i])
+            expect = 2 * values[i] - y
+            np.testing.assert_array_equal(eval[0], expect)
+
     def test_to_from_json(self, mocker):
         # UnaryOperator
         a = pybamm.Symbol("a", domain=["test"])
diff --git a/tests/unit/test_solvers/test_processed_variable.py b/tests/unit/test_solvers/test_processed_variable.py
index 04de88963..b3ab85060 100644
--- a/tests/unit/test_solvers/test_processed_variable.py
+++ b/tests/unit/test_solvers/test_processed_variable.py
@@ -160,6 +160,67 @@ class TestProcessedVariable:
                 processed_var._observe_raw_cpp(), processed_var._observe_raw_python()
             )
 
+    @pytest.mark.parametrize("hermite_interp", _hermite_args)
+    def test_processed_variable_0D_discrete_data(self, hermite_interp):
+        y = pybamm.StateVector(slice(0, 1))
+        t_sol = np.linspace(0, 1)
+        y_sol = np.array([np.linspace(0, 5)])
+        data_const = 3.6
+        if hermite_interp:
+            yp_sol = 5 * np.ones_like(y_sol)
+        else:
+            yp_sol = None
+
+        # hermite interpolation can do order 2 interpolation, otherwise make sure result is linear
+        order = 2 if hermite_interp else 1
+
+        # data is same timepoints as solution
+        data_t = t_sol
+        data_v = -data_const * data_t
+        data = pybamm.DiscreteTimeData(data_t, data_v, "test_data")
+        var = (y - data) ** order
+        expected_entries = (y_sol - data_v) ** order
+        var.mesh = None
+        model = pybamm.BaseModel()
+        var_casadi = to_casadi(var, y_sol)
+        processed_var = pybamm.process_variable(
+            [var],
+            [var_casadi],
+            self._sol_default(t_sol, y_sol, yp_sol, model),
+        )
+        np.testing.assert_array_equal(processed_var.entries, expected_entries.flatten())
+        np.testing.assert_array_equal(processed_var(data_t), expected_entries.flatten())
+
+        # data is different timepoints as solution
+        data_t = np.linspace(0, 1, 7)
+        data_v = -data_const * data_t
+        y_sol_interp = (np.interp(data_t, t_sol, y_sol[0]),)
+        data_v_interp = np.interp(t_sol, data_t, data_v)
+        data = pybamm.DiscreteTimeData(data_t, data_v, "test_data")
+
+        # check data interp
+        np.testing.assert_array_almost_equal(
+            data.evaluate(t=t_sol).flatten(), data_v_interp
+        )
+
+        var = (y - data) ** order
+        expected = (y_sol_interp - data_v) ** order
+        expected_entries = (y_sol - data_v_interp) ** order
+        var.mesh = None
+        model = pybamm.BaseModel()
+        var_casadi = to_casadi(var, y_sol)
+        processed_var = pybamm.process_variable(
+            [var],
+            [var_casadi],
+            self._sol_default(t_sol, y_sol, yp_sol, model),
+        )
+        np.testing.assert_array_almost_equal(
+            processed_var.entries, expected_entries.flatten(), decimal=10
+        )
+        np.testing.assert_array_almost_equal(
+            processed_var(t=data_t), expected.flatten(), decimal=10
+        )
+
     @pytest.mark.parametrize("hermite_interp", _hermite_args)
     def test_processed_variable_0D_no_sensitivity(self, hermite_interp):
         # without space
diff --git a/tests/unit/test_solvers/test_solution.py b/tests/unit/test_solvers/test_solution.py
index 4ac931253..2fb25f79a 100644
--- a/tests/unit/test_solvers/test_solution.py
+++ b/tests/unit/test_solvers/test_solution.py
@@ -420,3 +420,43 @@ class TestSolution:
         sim.solve(t_eval=np.linspace(0, 10, 10), inputs=inputs)
         time = sim.solution["Time [h]"](sim.solution.t)
         assert len(time) == 10
+
+    _solver_classes = [pybamm.CasadiSolver]
+    if pybamm.has_idaklu():
+        _solver_classes.append(pybamm.IDAKLUSolver)
+
+    @pytest.mark.parametrize("solver_class", _solver_classes)
+    def test_discrete_data_sum(self, solver_class):
+        model = pybamm.BaseModel(name="test_model")
+        c = pybamm.Variable("c")
+        model.rhs = {c: -c}
+        model.initial_conditions = {c: 1}
+        model.variables["c"] = c
+
+        data_times = np.linspace(0, 1, 10)
+        if solver_class == pybamm.IDAKLUSolver:
+            t_eval = [data_times[0], data_times[-1]]
+            t_interp = data_times
+        else:
+            t_eval = data_times
+            t_interp = None
+        solver = solver_class()
+        data_values = solver.solve(model, t_eval=t_eval, t_interp=t_interp)["c"].entries
+
+        data = pybamm.DiscreteTimeData(data_times, data_values, "test_data")
+        data_comparison = pybamm.DiscreteTimeSum((c - data) ** 2)
+
+        model = pybamm.BaseModel(name="test_model2")
+        a = pybamm.InputParameter("a")
+        model.rhs = {c: -a * c}
+        model.initial_conditions = {c: 1}
+        model.variables["data_comparison"] = data_comparison
+
+        solver = solver_class()
+        for a in [0.5, 1.0, 2.0]:
+            sol = solver.solve(model, t_eval=t_eval, inputs={"a": a})
+            y_sol = np.exp(-a * data_times)
+            expected = np.sum((y_sol - data_values) ** 2)
+            np.testing.assert_array_almost_equal(
+                sol["data_comparison"](), expected, decimal=2
+            )