Repository: pybamm-team/PyBaMM

Issue:
Add an Experiment class 
**Summary**
Add an `Experiment` class where the user can define all thing related to the specific set up of an experiment, including:
- initial conditions
- cycling profile
- SOC to concentrations calculator
- cooling set up? 

Feel free to suggest other things that should be in this

 **Motivation**
Performing different experiments on models will be one of the main use cases of PyBaMM so we should have a clear and simple method for setting this up. 


Discussion:

Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 4ab9fea93..94e4d1707 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -2,6 +2,8 @@
 
 ## Features
 
+-   Added functionality to solve DAE models with non-smooth current inputs ([#808](https://github.com/pybamm-team/PyBaMM/pull/808))
+-   Added functionality to simulate experiments and testing protocols ([#807](https://github.com/pybamm-team/PyBaMM/pull/807))
 -   Added fuzzy string matching for parameters and variables ([#796](https://github.com/pybamm-team/PyBaMM/pull/796))
 -   Changed ParameterValues to raise an error when a parameter that wasn't previously defined is updated ([#796](https://github.com/pybamm-team/PyBaMM/pull/796))
 -   Added some basic models (BasicSPM and BasicDFN) in order to clearly demonstrate the PyBaMM model structure for battery models ([#795](https://github.com/pybamm-team/PyBaMM/pull/795))
@@ -63,6 +65,7 @@
 
 ## Breaking changes
 
+-   Removed `ParameterValues.update_model`, whose functionality is now replaced by `InputParameter` ([#801](https://github.com/pybamm-team/PyBaMM/pull/801))
 -   Removed `Outer` and `Kron` nodes as no longer used ([#777](https://github.com/pybamm-team/PyBaMM/pull/777))
 -   Moved `results` to separate repositories ([#761](https://github.com/pybamm-team/PyBaMM/pull/761))
 -   The parameters "Bruggeman coefficient" must now be specified separately as "Bruggeman coefficient (electrolyte)" and "Bruggeman coefficient (electrode)"
diff --git a/docs/index.rst b/docs/index.rst
index d991331d6..7b59f96c4 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -30,6 +30,7 @@ Contents
     source/meshes/index
     source/spatial_methods/index
     source/solvers/index
+    source/experiments/index
     source/processed_variable
     source/util
     source/simulation
diff --git a/docs/source/experiments/experiment.rst b/docs/source/experiments/experiment.rst
new file mode 100644
index 000000000..52db069f2
--- /dev/null
+++ b/docs/source/experiments/experiment.rst
@@ -0,0 +1,5 @@
+Base Experiment Class
+=====================
+
+.. autoclass:: pybamm.Experiment
+  :members:
diff --git a/docs/source/experiments/index.rst b/docs/source/experiments/index.rst
new file mode 100644
index 000000000..05dafa486
--- /dev/null
+++ b/docs/source/experiments/index.rst
@@ -0,0 +1,8 @@
+Experiments
+===========
+
+Classes to help set operating conditions for some standard battery modelling experiments
+
+.. toctree::
+
+  experiment
\ No newline at end of file
diff --git a/docs/source/solvers/solution.rst b/docs/source/solvers/solution.rst
index 3f6dc3040..25a9977ec 100644
--- a/docs/source/solvers/solution.rst
+++ b/docs/source/solvers/solution.rst
@@ -1,5 +1,8 @@
 Solution
 ========
 
+.. autoclass:: pybamm._BaseSolution
+  :members:
+
 .. autoclass:: pybamm.Solution
   :members:
diff --git a/examples/scripts/DFN.py b/examples/scripts/DFN.py
index 4811df8a1..7c4e8539e 100644
--- a/examples/scripts/DFN.py
+++ b/examples/scripts/DFN.py
@@ -7,15 +7,15 @@ import numpy as np
 
 pybamm.set_logging_level("INFO")
 
+
 # load model
-model = pybamm.lithium_ion.DFN({"operating mode": "voltage"})
+model = pybamm.lithium_ion.DFN()
 
 # create geometry
 geometry = model.default_geometry
 
 # load parameter values and process model and geometry
 param = model.default_parameter_values
-param.update({"Voltage function [V]": 4.1}, check_already_exists=False)
 param.process_model(model)
 param.process_geometry(geometry)
 
diff --git a/examples/scripts/experimental_protocols/cccv.py b/examples/scripts/experimental_protocols/cccv.py
new file mode 100644
index 000000000..d287e0b91
--- /dev/null
+++ b/examples/scripts/experimental_protocols/cccv.py
@@ -0,0 +1,39 @@
+#
+# Constant-current constant-voltage charge
+#
+import pybamm
+import matplotlib.pyplot as plt
+
+pybamm.set_logging_level("INFO")
+experiment = pybamm.Experiment(
+    [
+        "Discharge at C/10 for 13 hours or until 3.3 V",
+        "Rest for 1 hour",
+        "Charge at 1 A until 4.1 V",
+        "Hold at 4.1 V until 50 mA",
+        "Rest for 1 hour",
+    ]
+    * 3,
+    period="2 minutes",
+)
+model = pybamm.lithium_ion.DFN()
+sim = pybamm.Simulation(model, experiment=experiment, solver=pybamm.CasadiSolver())
+sim.solve()
+
+# Plot voltages from the discharge segments only
+fig, ax = plt.subplots()
+for i in range(3):
+    # Extract sub solutions
+    sol = sim.solution.sub_solutions[i * 5]
+    # Extract variables
+    t = sol["Time [h]"].entries
+    V = sol["Terminal voltage [V]"].entries
+    # Plot
+    ax.plot(t - t[0], V, label="Discharge {}".format(i + 1))
+    ax.set_xlabel("Time [h]")
+    ax.set_ylabel("Voltage [V]")
+    ax.set_xlim([0, 13])
+ax.legend()
+
+# Show all plots
+sim.plot()
diff --git a/examples/scripts/experimental_protocols/cccv_lead_acid.py b/examples/scripts/experimental_protocols/cccv_lead_acid.py
new file mode 100644
index 000000000..4614fd514
--- /dev/null
+++ b/examples/scripts/experimental_protocols/cccv_lead_acid.py
@@ -0,0 +1,19 @@
+#
+# Constant-current constant-voltage charge
+#
+import pybamm
+
+pybamm.set_logging_level("INFO")
+experiment = pybamm.Experiment(
+    [
+        "Discharge at C/2 until 11 V",
+        "Rest for 1 hour",
+        "Charge at C/2 until 14.5 V",
+        "Hold at 14.5 V until 200 mA",
+        "Rest for 1 hour",
+    ]
+)
+model = pybamm.lead_acid.Full()
+sim = pybamm.Simulation(model, experiment=experiment, solver=pybamm.CasadiSolver())
+sim.solve()
+sim.plot()
diff --git a/examples/scripts/experimental_protocols/gitt.py b/examples/scripts/experimental_protocols/gitt.py
new file mode 100644
index 000000000..dfebc7096
--- /dev/null
+++ b/examples/scripts/experimental_protocols/gitt.py
@@ -0,0 +1,13 @@
+#
+# GITT discharge
+#
+import pybamm
+
+pybamm.set_logging_level("INFO")
+experiment = pybamm.Experiment(
+    ["Discharge at C/20 for 1 hour", "Rest for 1 hour"] * 20,
+)
+model = pybamm.lithium_ion.DFN()
+sim = pybamm.Simulation(model, experiment=experiment, solver=pybamm.CasadiSolver())
+sim.solve()
+sim.plot()
diff --git a/pybamm/__init__.py b/pybamm/__init__.py
index 9f560939e..4b5dbcbb6 100644
--- a/pybamm/__init__.py
+++ b/pybamm/__init__.py
@@ -237,7 +237,7 @@ from .spatial_methods.scikit_finite_element import ScikitFiniteElement
 #
 # Solver classes
 #
-from .solvers.solution import Solution
+from .solvers.solution import Solution, _BaseSolution
 from .solvers.base_solver import BaseSolver
 from .solvers.algebraic_solver import AlgebraicSolver
 from .solvers.casadi_solver import CasadiSolver
@@ -246,6 +246,12 @@ from .solvers.scikits_ode_solver import ScikitsOdeSolver, have_scikits_odes
 from .solvers.scipy_solver import ScipySolver
 from .solvers.idaklu_solver import IDAKLUSolver, have_idaklu
 
+#
+# Experiments
+#
+from .experiments.experiment import Experiment
+from . import experiments
+
 #
 # other
 #
diff --git a/pybamm/experiments/__init__.py b/pybamm/experiments/__init__.py
new file mode 100644
index 000000000..e69de29bb
diff --git a/pybamm/experiments/experiment.py b/pybamm/experiments/experiment.py
new file mode 100644
index 000000000..fb7bbef64
--- /dev/null
+++ b/pybamm/experiments/experiment.py
@@ -0,0 +1,251 @@
+#
+# Experiment class
+#
+
+examples = """
+
+    Discharge at 1C for 0.5 hours,
+    Discharge at C/20 for 0.5 hours,
+    Charge at 0.5 C for 45 minutes,
+    Discharge at 1 A for 90 seconds,
+    Charge at 200mA for 45 minutes (1 minute period),
+    Discharge at 1 W for 0.5 hours,
+    Charge at 200 mW for 45 minutes,
+    Rest for 10 minutes (5 minute period),
+    Hold at 1 V for 20 seconds,
+    Charge at 1 C until 4.1V,
+    Hold at 4.1 V until 50 mA,
+    Hold at 3V until C/50,
+    """
+
+
+class Experiment:
+    """
+    Base class for experimental conditions under which to run the model. In general, a
+    list of operating conditions should be passed in. Each operating condition should
+    be of the form "Do this for this long" or "Do this until this happens". For example,
+    "Charge at 1 C for 1 hour", or "Charge at 1 C until 4.2 V", or "Charge at 1 C for 1
+    hour or until 4.2 V". The instructions can be of the form "(Dis)charge at x A/C/W",
+    "Rest", or "Hold at x V". The running time should be a time in seconds, minutes or
+    hours, e.g. "10 seconds", "3 minutes" or "1 hour". The stopping conditions should be
+    a circuit state, e.g. "1 A", "C/50" or "3 V".
+
+    Parameters
+    ----------
+    operating_conditions : list
+        List of operating conditions
+    parameters : dict
+        Dictionary of parameters to use for this experiment, replacing default
+        parameters as appropriate
+    period : string, optional
+        Period (1/frequency) at which to record outputs. Default is 1 minute. Can be
+        overwritten by individual operating conditions.
+
+    """
+
+    def __init__(self, operating_conditions, parameters=None, period="1 minute"):
+        self.period = self.convert_time_to_seconds(period.split())
+        self.operating_conditions_strings = operating_conditions
+        self.operating_conditions, self.events = self.read_operating_conditions(
+            operating_conditions
+        )
+        parameters = parameters or {}
+        if isinstance(parameters, dict):
+            self.parameters = parameters
+        else:
+            raise TypeError("experimental parameters should be a dictionary")
+
+    def __str__(self):
+        return str(self.operating_conditions_strings)
+
+    def __repr__(self):
+        return "pybamm.Experiment({!s})".format(self)
+
+    def read_operating_conditions(self, operating_conditions):
+        """
+        Convert operating conditions to the appropriate format
+
+        Parameters
+        ----------
+        operating_conditions : list
+            List of operating conditions
+
+        Returns
+        -------
+        operating_conditions : list
+            Operating conditions in the tuple format
+        """
+        converted_operating_conditions = []
+        events = []
+        for cond in operating_conditions:
+            if isinstance(cond, str):
+                next_op, next_event = self.read_string(cond)
+                converted_operating_conditions.append(next_op)
+                events.append(next_event)
+            else:
+                raise TypeError(
+                    """Operating conditions should be strings, not {}. For example: {}
+                    """.format(
+                        type(cond), examples
+                    )
+                )
+
+        return converted_operating_conditions, events
+
+    def read_string(self, cond):
+        """
+        Convert a string to a tuple of the right format
+
+        Parameters
+        ----------
+        cond : str
+            String of appropriate form for example "Charge at x C for y hours". x and y
+            must be numbers, 'C' denotes the unit of the external circuit (can be A for
+            current, C for C-rate, V for voltage or W for power), and 'hours' denotes
+            the unit of time (can be second(s), minute(s) or hour(s))
+        """
+        # Read period
+        if " period)" in cond:
+            cond, time_period = cond.split("(")
+            time, _ = time_period.split(" period)")
+            period = self.convert_time_to_seconds(time.split())
+        else:
+            period = self.period
+        # Read instructions
+        if "for" in cond and "or until" in cond:
+            # e.g. for 3 hours or until 4.2 V
+            cond_list = cond.split()
+            idx_for = cond_list.index("for")
+            idx_until = cond_list.index("or")
+            electric = self.convert_electric(cond_list[:idx_for])
+            time = self.convert_time_to_seconds(cond_list[idx_for + 1 : idx_until])
+            events = self.convert_electric(cond_list[idx_until + 2 :])
+        elif "for" in cond:
+            # e.g. for 3 hours
+            cond_list = cond.split()
+            idx = cond_list.index("for")
+            electric = self.convert_electric(cond_list[:idx])
+            time = self.convert_time_to_seconds(cond_list[idx + 1 :])
+            events = None
+        elif "until" in cond:
+            # e.g. until 4.2 V
+            cond_list = cond.split()
+            idx = cond_list.index("until")
+            electric = self.convert_electric(cond_list[:idx])
+            time = None
+            events = self.convert_electric(cond_list[idx + 1 :])
+        else:
+            raise ValueError(
+                """Operating conditions must contain keyword 'for' or 'until'.
+                For example: {}""".format(
+                    examples
+                )
+            )
+        return electric + (time,) + (period,), events
+
+    def convert_electric(self, electric):
+        "Convert electrical instructions to consistent output"
+        # Rest == zero current
+        if electric[0].lower() == "rest":
+            return (0, "A")
+        else:
+            if len(electric) in [3, 4]:
+                if len(electric) == 4:
+                    # e.g. Charge at 4 A, Hold at 3 V
+                    instruction, _, value, unit = electric
+                elif len(electric) == 3:
+                    # e.g. Discharge at C/2, Charge at 1A
+                    instruction, _, value_unit = electric
+                    if value_unit[0] == "C":
+                        # e.g. C/2
+                        unit = value_unit[0]
+                        value = 1 / float(value_unit[2:])
+                    else:
+                        # e.g. 1A
+                        if "m" in value_unit:
+                            # e.g. 1mA
+                            unit = value_unit[-2:]
+                            value = float(value_unit[:-2])
+                        else:
+                            # e.g. 1A
+                            unit = value_unit[-1]
+                            value = float(value_unit[:-1])
+                # Read instruction
+                if instruction.lower() in ["discharge", "hold"]:
+                    sign = 1
+                elif instruction.lower() == "charge":
+                    sign = -1
+                else:
+                    raise ValueError(
+                        """instruction must be 'discharge', 'charge', 'rest' or 'hold'.
+                        For example: {}""".format(
+                            examples
+                        )
+                    )
+            elif len(electric) == 2:
+                # e.g. 3 A, 4.1 V
+                value, unit = electric
+                sign = 1
+            elif len(electric) == 1:
+                # e.g. C/2, 1A
+                value_unit = electric[0]
+                if value_unit[0] == "C":
+                    # e.g. C/2
+                    unit = value_unit[0]
+                    value = 1 / float(value_unit[2:])
+                else:
+                    if "m" in value_unit:
+                        # e.g. 1mA
+                        unit = value_unit[-2:]
+                        value = float(value_unit[:-2])
+                    else:
+                        # e.g. 1A
+                        unit = value_unit[-1]
+                        value = float(value_unit[:-1])
+                sign = 1
+            else:
+                raise ValueError(
+                    """Instruction '{}' not recognized. Some acceptable examples are: {}
+                    """.format(
+                        " ".join(electric), examples
+                    )
+                )
+            # Read value and units
+            if unit == "C":
+                return (sign * float(value), "C")
+            elif unit == "A":
+                return (sign * float(value), "A")
+            elif unit == "mA":
+                return (sign * float(value) / 1000, "A")
+            elif unit == "V":
+                return (float(value), "V")
+            elif unit == "W":
+                return (sign * float(value), "W")
+            elif unit == "mW":
+                return (sign * float(value) / 1000, "W")
+            else:
+                raise ValueError(
+                    """units must be 'C', 'A', 'mA', 'V', 'W' or 'mW', not '{}'.
+                    For example: {}
+                    """.format(
+                        unit, examples
+                    )
+                )
+
+    def convert_time_to_seconds(self, time_and_units):
+        "Convert a time in seconds, minutes or hours to a time in seconds"
+        time, units = time_and_units
+        if units in ["second", "seconds", "s", "sec"]:
+            time_in_seconds = float(time)
+        elif units in ["minute", "minutes", "m", "min"]:
+            time_in_seconds = float(time) * 60
+        elif units in ["hour", "hours", "h", "hr"]:
+            time_in_seconds = float(time) * 3600
+        else:
+            raise ValueError(
+                """time units must be 'seconds', 'minutes' or 'hours'. For example: {}
+                """.format(
+                    examples
+                )
+            )
+        return time_in_seconds
diff --git a/pybamm/expression_tree/exceptions.py b/pybamm/expression_tree/exceptions.py
index cbeca72b5..7097c475a 100644
--- a/pybamm/expression_tree/exceptions.py
+++ b/pybamm/expression_tree/exceptions.py
@@ -67,3 +67,4 @@ class InputError(Exception):
     """
 
     pass
+
diff --git a/pybamm/models/base_model.py b/pybamm/models/base_model.py
index 519ea9bee..d09d43046 100644
--- a/pybamm/models/base_model.py
+++ b/pybamm/models/base_model.py
@@ -304,6 +304,11 @@ class BaseModel(object):
     def options(self, options):
         self._options = options
 
+    @property
+    def timescale(self):
+        "Default timescale for a model is 1 second"
+        return pybamm.Scalar(1)
+
     def __getitem__(self, key):
         return self.rhs[key]
 
diff --git a/pybamm/models/full_battery_models/base_battery_model.py b/pybamm/models/full_battery_models/base_battery_model.py
index a8e74f218..c179a8ec0 100644
--- a/pybamm/models/full_battery_models/base_battery_model.py
+++ b/pybamm/models/full_battery_models/base_battery_model.py
@@ -253,6 +253,11 @@ class BaseBatteryModel(pybamm.BaseModel):
 
         self._options = options
 
+    @property
+    def timescale(self):
+        "Default timescale for a battery model is the discharge timescale"
+        return self.param.tau_discharge
+
     def set_standard_output_variables(self):
         # Standard output variables
 
diff --git a/pybamm/models/submodels/external_circuit/function_control_external_circuit.py b/pybamm/models/submodels/external_circuit/function_control_external_circuit.py
index fb1dd39f6..6e3812f99 100644
--- a/pybamm/models/submodels/external_circuit/function_control_external_circuit.py
+++ b/pybamm/models/submodels/external_circuit/function_control_external_circuit.py
@@ -8,9 +8,9 @@ from .base_external_circuit import BaseModel, LeadingOrderBaseModel
 class FunctionControl(BaseModel):
     """External circuit with an arbitrary function. """
 
-    def __init__(self, param, external_circuit_class):
+    def __init__(self, param, external_circuit_function):
         super().__init__(param)
-        self.external_circuit_class = external_circuit_class
+        self.external_circuit_function = external_circuit_function
 
     def _get_current_variable(self):
         return pybamm.Variable("Total current density")
@@ -23,14 +23,6 @@ class FunctionControl(BaseModel):
         # Add discharge capacity variable
         variables.update(super().get_fundamental_variables())
 
-        # Add switches
-        # These are not implemented yet but can be used later with the Experiment class
-        # to simulate different external circuit conditions sequentially within a
-        # single model (for example Constant Current - Constant Voltage)
-        # for i in range(self.external_circuit_class.num_switches):
-        #     s = pybamm.Parameter("Switch {}".format(i + 1))
-        #     variables["Switch {}".format(i + 1)] = s
-
         return variables
 
     def set_initial_conditions(self, variables):
@@ -44,7 +36,7 @@ class FunctionControl(BaseModel):
         # The external circuit function should fix either the current, or the voltage,
         # or a combination (e.g. I*V for power control)
         i_cell = variables["Total current density"]
-        self.algebraic[i_cell] = self.external_circuit_class(variables)
+        self.algebraic[i_cell] = self.external_circuit_function(variables)
 
 
 class VoltageFunctionControl(FunctionControl):
@@ -53,31 +45,25 @@ class VoltageFunctionControl(FunctionControl):
     """
 
     def __init__(self, param):
-        super().__init__(param, ConstantVoltage())
+        super().__init__(param, constant_voltage)
 
 
-class ConstantVoltage:
-    num_switches = 0
-
-    def __call__(self, variables):
-        V = variables["Terminal voltage [V]"]
-        return V - pybamm.FunctionParameter("Voltage function [V]", pybamm.t)
+def constant_voltage(variables):
+    V = variables["Terminal voltage [V]"]
+    return V - pybamm.FunctionParameter("Voltage function [V]", pybamm.t)
 
 
 class PowerFunctionControl(FunctionControl):
     """External circuit with power control. """
 
     def __init__(self, param):
-        super().__init__(param, ConstantPower())
-
+        super().__init__(param, constant_power)
 
-class ConstantPower:
-    num_switches = 0
 
-    def __call__(self, variables):
-        I = variables["Current [A]"]
-        V = variables["Terminal voltage [V]"]
-        return I * V - pybamm.FunctionParameter("Power function [W]", pybamm.t)
+def constant_power(variables):
+    I = variables["Current [A]"]
+    V = variables["Terminal voltage [V]"]
+    return I * V - pybamm.FunctionParameter("Power function [W]", pybamm.t)
 
 
 class LeadingOrderFunctionControl(FunctionControl, LeadingOrderBaseModel):
@@ -97,12 +83,12 @@ class LeadingOrderVoltageFunctionControl(LeadingOrderFunctionControl):
     """
 
     def __init__(self, param):
-        super().__init__(param, ConstantVoltage())
+        super().__init__(param, constant_voltage)
 
 
 class LeadingOrderPowerFunctionControl(LeadingOrderFunctionControl):
     """External circuit with power control, at leading order. """
 
     def __init__(self, param):
-        super().__init__(param, ConstantPower())
+        super().__init__(param, constant_power)
 
diff --git a/pybamm/parameters/electrical_parameters.py b/pybamm/parameters/electrical_parameters.py
index 7e9d1ecfa..c83405ca5 100644
--- a/pybamm/parameters/electrical_parameters.py
+++ b/pybamm/parameters/electrical_parameters.py
@@ -13,6 +13,7 @@ C_rate = abs(I_typ / Q)
 n_electrodes_parallel = pybamm.Parameter(
     "Number of electrodes connected in parallel to make a cell"
 )
+n_cells = pybamm.Parameter("Number of cells connected in series to make a battery")
 i_typ = pybamm.Function(
     np.abs, I_typ / (n_electrodes_parallel * pybamm.geometric_parameters.A_cc)
 )
diff --git a/pybamm/simulation.py b/pybamm/simulation.py
index d94e47e7c..b16c83f46 100644
--- a/pybamm/simulation.py
+++ b/pybamm/simulation.py
@@ -20,6 +20,26 @@ def isnotebook():
         return False  # Probably standard Python interpreter
 
 
+def constant_current_constant_voltage_constant_power(variables):
+    I = variables["Current [A]"]
+    V = variables["Terminal voltage [V]"]
+    s_I = pybamm.InputParameter("Current switch")
+    s_V = pybamm.InputParameter("Voltage switch")
+    s_P = pybamm.InputParameter("Power switch")
+    n_electrodes_parallel = pybamm.electrical_parameters.n_electrodes_parallel
+    n_cells = pybamm.electrical_parameters.n_cells
+    return (
+        s_I * (I - pybamm.InputParameter("Current input [A]") / n_electrodes_parallel)
+        + s_V * (V - pybamm.InputParameter("Voltage input [V]") / n_cells)
+        + s_P
+        * (
+            V * I
+            - pybamm.InputParameter("Power input [W]")
+            / (n_cells * n_electrodes_parallel)
+        )
+    )
+
+
 class Simulation:
     """A Simulation class for easy building and running of PyBaMM simulations.
 
@@ -27,6 +47,8 @@ class Simulation:
     ----------
     model : :class:`pybamm.BaseModel`
         The model to be simulated
+    experiment : : class:`pybamm.Experiment` (optional)
+        The experimental conditions under which to solve the model
     geometry: :class:`pybamm.Geometry` (optional)
         The geometry upon which to solve the model
     parameter_values: dict (optional)
@@ -52,6 +74,7 @@ class Simulation:
     def __init__(
         self,
         model,
+        experiment=None,
         geometry=None,
         parameter_values=None,
         submesh_types=None,
@@ -61,19 +84,23 @@ class Simulation:
         quick_plot_vars=None,
         C_rate=None,
     ):
-        self.model = model
-
-        self.geometry = geometry or model.default_geometry
         self._parameter_values = parameter_values or model.default_parameter_values
-        self._submesh_types = submesh_types or model.default_submesh_types
-        self._var_pts = var_pts or model.default_var_pts
-        self._spatial_methods = spatial_methods or model.default_spatial_methods
-        self._solver = solver or self._model.default_solver
-        self._quick_plot_vars = quick_plot_vars
 
-        self.C_rate = C_rate
-        if self.C_rate:
-            self._parameter_values.update({"C-rate": self.C_rate})
+        if experiment is None:
+            self.operating_mode = "without experiment"
+            self.C_rate = C_rate
+            if self.C_rate:
+                self._parameter_values.update({"C-rate": self.C_rate})
+            self.model = model
+        else:
+            self.set_up_experiment(model, experiment)
+
+        self.geometry = geometry or self.model.default_geometry
+        self._submesh_types = submesh_types or self.model.default_submesh_types
+        self._var_pts = var_pts or self.model.default_var_pts
+        self._spatial_methods = spatial_methods or self.model.default_spatial_methods
+        self._solver = solver or self.model.default_solver
+        self._quick_plot_vars = quick_plot_vars
 
         self.reset(update_model=False)
 
@@ -83,6 +110,132 @@ class Simulation:
 
             warnings.filterwarnings("ignore")
 
+    def set_up_experiment(self, model, experiment):
+        """
+        Set up a simulation to run with an experiment. This creates a dictionary of
+        inputs (current/voltage/power, running time, stopping condition) for each
+        operating condition in the experiment. The model will then be solved by
+        integrating the model successively with each group of inputs, one group at a
+        time.
+        """
+        self.operating_mode = "with experiment"
+        self.model = model.new_copy(
+            options={
+                **model.options,
+                "operating mode": constant_current_constant_voltage_constant_power,
+            }
+        )
+        if not isinstance(experiment, pybamm.Experiment):
+            raise TypeError("experiment must be a pybamm `Experiment` instance")
+        # Save the experiment
+        self.experiment = experiment
+        # Update parameter values with experiment parameters
+        self._parameter_values.update(experiment.parameters)
+        # Create a new submodel for each set of operating conditions and update
+        # parameters and events accordingly
+        self._experiment_inputs = []
+        self._experiment_times = []
+        for op, events in zip(experiment.operating_conditions, experiment.events):
+            if op[1] in ["A", "C"]:
+                # Update inputs for constant current
+                if op[1] == "A":
+                    I = op[0]
+                else:
+                    # Scale C-rate with capacity to obtain current
+                    capacity = self._parameter_values["Cell capacity [A.h]"]
+                    I = op[0] * capacity
+                operating_inputs = {
+                    "Current switch": 1,
+                    "Voltage switch": 0,
+                    "Power switch": 0,
+                    "Current input [A]": I,
+                    "Voltage input [V]": 0,  # doesn't matter
+                    "Power input [W]": 0,  # doesn't matter
+                }
+            elif op[1] == "V":
+                # Update inputs for constant voltage
+                V = op[0]
+                operating_inputs = {
+                    "Current switch": 0,
+                    "Voltage switch": 1,
+                    "Power switch": 0,
+                    "Current input [A]": 0,  # doesn't matter
+                    "Voltage input [V]": V,
+                    "Power input [W]": 0,  # doesn't matter
+                }
+            elif op[1] == "W":
+                # Update inputs for constant power
+                P = op[0]
+                operating_inputs = {
+                    "Current switch": 0,
+                    "Voltage switch": 0,
+                    "Power switch": 1,
+                    "Current input [A]": 0,  # doesn't matter
+                    "Voltage input [V]": 0,  # doesn't matter
+                    "Power input [W]": P,
+                }
+            # Update period
+            operating_inputs["period"] = op[3]
+            # Update events
+            if events is None:
+                # make current and voltage values that won't be hit
+                operating_inputs.update(
+                    {"Current cut-off [A]": -1e10, "Voltage cut-off [V]": -1e10}
+                )
+            elif events[1] in ["A", "C"]:
+                # update current cut-off, make voltage a value that won't be hit
+                if events[1] == "A":
+                    I = events[0]
+                else:
+                    # Scale C-rate with capacity to obtain current
+                    capacity = self._parameter_values["Cell capacity [A.h]"]
+                    I = events[0] * capacity
+                operating_inputs.update(
+                    {"Current cut-off [A]": I, "Voltage cut-off [V]": -1e10}
+                )
+            elif events[1] == "V":
+                # update voltage cut-off, make current a value that won't be hit
+                V = events[0]
+                operating_inputs.update(
+                    {"Current cut-off [A]": -1e10, "Voltage cut-off [V]": V}
+                )
+
+            self._experiment_inputs.append(operating_inputs)
+            # Convert time to dimensionless
+            dt_dimensional = op[2]
+            if dt_dimensional is None:
+                # max simulation time: 1 week
+                dt_dimensional = 7 * 24 * 3600
+            tau = self._parameter_values.evaluate(self.model.timescale)
+            dt_dimensionless = dt_dimensional / tau
+            self._experiment_times.append(dt_dimensionless)
+
+        # add current and voltage events to the model
+        # current events both negative and positive to catch specification
+        n_electrodes_parallel = pybamm.electrical_parameters.n_electrodes_parallel
+        n_cells = pybamm.electrical_parameters.n_cells
+        self.model.events.extend(
+            [
+                pybamm.Event(
+                    "Current cut-off (positive) [A] [experiment]",
+                    self.model.variables["Current [A]"]
+                    - abs(pybamm.InputParameter("Current cut-off [A]"))
+                    / n_electrodes_parallel,
+                ),
+                pybamm.Event(
+                    "Current cut-off (negative) [A] [experiment]",
+                    self.model.variables["Current [A]"]
+                    + abs(pybamm.InputParameter("Current cut-off [A]"))
+                    / n_electrodes_parallel,
+                ),
+                pybamm.Event(
+                    "Voltage cut-off [V] [experiment]",
+                    self.model.variables["Terminal voltage [V]"]
+                    - pybamm.InputParameter("Voltage cut-off [V]") / n_cells,
+                ),
+            ]
+        )
+
     def set_defaults(self):
         """
         A method to set all the simulation specs to default values for the
@@ -181,28 +334,68 @@ class Simulation:
             If True, model checks are performed after discretisation (see
             :meth:`pybamm.Discretisation.process_model`). Default is True.
         """
+        # Setup
         self.build(check_model=check_model)
-
-        if t_eval is None:
-            try:
-                # Try to compute discharge time
-                tau = self._parameter_values.evaluate(self.model.param.tau_discharge)
-                C_rate = self._parameter_values["C-rate"]
-                t_end = 3600 / tau / C_rate
-                t_eval = np.linspace(0, t_end, 100)
-            except AttributeError:
-                t_eval = np.linspace(0, 1, 100)
-
         if solver is None:
             solver = self.solver
 
-        self.t_eval = t_eval
-        self._solution = solver.solve(
-            self.built_model,
-            t_eval,
-            external_variables=external_variables,
-            inputs=inputs,
-        )
+        if self.operating_mode == "without experiment":
+            # Solve the normal way, with a single solve
+            if t_eval is None:
+                try:
+                    # Try to compute discharge time
+                    tau = self._parameter_values.evaluate(self.model.param.timescale)
+                    C_rate = self._parameter_values["C-rate"]
+                    t_end = 3600 / tau / C_rate
+                    t_eval = np.linspace(0, t_end, 100)
+                except AttributeError:
+                    t_eval = np.linspace(0, 1, 100)
+
+            self.t_eval = t_eval
+            self._solution = solver.solve(self.built_model, t_eval, inputs=inputs)
+        elif self.operating_mode == "with experiment":
+            if t_eval is not None:
+                pybamm.logger.warning(
+                    "Ignoring t_eval as solution times are specified by the experiment"
+                )
+            # Step through all experimental conditions
+            inputs = inputs or {}
+            pybamm.logger.info("Start running experiment")
+            timer = pybamm.Timer()
+            for idx, (exp_inputs, dt) in enumerate(
+                zip(self._experiment_inputs, self._experiment_times)
+            ):
+                pybamm.logger.info(self.experiment.operating_conditions_strings[idx])
+                inputs.update(exp_inputs)
+                # Non-dimensionalise period
+                tau = self._parameter_values.evaluate(self.model.timescale)
+                freq = exp_inputs["period"] / tau
+                # Make sure we take at least 2 timesteps
+                npts = max(int(round(dt / freq)) + 1, 2)
+                self.step(
+                    dt, npts=npts, external_variables=external_variables, inputs=inputs
+                )
+                # Only allow events specified by experiment
+                if not (
+                    self._solution.termination == "final time"
+                    or "[experiment]" in self._solution.termination
+                ):
+                    pybamm.logger.warning(
+                        """
+                        Experiment is infeasible: '{}' was triggered during '{}'. Try
+                        reducing current, shortening the time interval, or reducing
+                        the period.
+                        """.format(
+                            self._solution.termination,
+                            self.experiment.operating_conditions_strings[idx],
+                        )
+                    )
+                    break
+            pybamm.logger.info(
+                "Finish experiment simulation, took {}".format(
+                    timer.format(timer.time())
+                )
+            )
 
     def step(
         self, dt, solver=None, npts=2, external_variables=None, inputs=None, save=True
@@ -469,7 +662,7 @@ class Simulation:
     def save(self, filename):
         """Save simulation using pickle"""
         if self.model.convert_to_format == "python":
-            # We currently cannot save models in the 'python'
+            # We currently cannot save models in the 'python' format
             raise NotImplementedError(
                 """
                 Cannot save simulation if model format is python.
diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py
index 38e53903c..450735eee 100644
--- a/pybamm/solvers/base_solver.py
+++ b/pybamm/solvers/base_solver.py
@@ -47,6 +47,7 @@ class BaseSolver(object):
         self.root_tol = root_tol
         self.max_steps = max_steps
 
+        self.models_set_up = set()
         self.model_step_times = {}
 
         # Defaults, can be overwritten by specific solver
@@ -432,9 +433,13 @@ class BaseSolver(object):
         inputs = inputs or {}
         ext_and_inputs = {**external_variables, **inputs}
 
-        self.set_up(model, ext_and_inputs)
-        set_up_time = timer.time()
-
+        # Set up (if not done already)
+        if model not in self.models_set_up:
+            self.set_up(model, ext_and_inputs)
+            set_up_time = timer.time()
+            self.models_set_up.add(model)
+        else:
+            set_up_time = 0
         # Solve
         # Set inputs and external
         self.set_inputs(model, ext_and_inputs)
@@ -568,8 +573,13 @@ class BaseSolver(object):
             If an empty model is passed (`model.rhs = {}` and `model.algebraic={}`)
 
         """
-        if old_solution is not None and old_solution.termination != "final time":
+
+        if old_solution is not None and not (
+            old_solution.termination == "final time"
+            or "[experiment]" in old_solution.termination
+        ):
             # Return same solution as an event has already been triggered
+            # With hack to allow stepping past experiment current / voltage cut-off
             return old_solution
 
         # Make sure model isn't empty
diff --git a/pybamm/solvers/casadi_solver.py b/pybamm/solvers/casadi_solver.py
index 01e031c1c..e29c78904 100644
--- a/pybamm/solvers/casadi_solver.py
+++ b/pybamm/solvers/casadi_solver.py
@@ -107,9 +107,7 @@ class CasadiSolver(pybamm.BaseSolver):
                 np.concatenate([event(0, model.y0)
                                 for event in model.terminate_events_eval])
             )
-            pybamm.logger.info(
-                "Start solving {} with {} in 'safe' mode".format(model.name, self.name)
-            )
+            pybamm.logger.info("Start solving {} with {}".format(model.name, self.name))
             t = t_eval[0]
             y0 = model.y0
             # Initialize solution
diff --git a/pybamm/solvers/solution.py b/pybamm/solvers/solution.py
index 31f9ef387..1763997ef 100644
--- a/pybamm/solvers/solution.py
+++ b/pybamm/solvers/solution.py
@@ -1,6 +1,7 @@
 #
 # Solution class
 #
+import copy
 import numbers
 import numpy as np
 import pickle
@@ -8,10 +9,11 @@ import pybamm
 from collections import defaultdict
 
 
-class Solution(object):
+class _BaseSolution(object):
     """
-    Class containing the solution of, and various attributes associated with, a PyBaMM
-    model.
+    (Semi-private) class containing the solution of, and various attributes associated
+    with, a PyBaMM model. This class is automatically created by the `Solution` class,
+    and should never be called from outside the `Solution` class.
 
     Parameters
     ----------
@@ -27,48 +29,56 @@ class Solution(object):
         the event happens.
     termination : str
         String to indicate why the solution terminated
+    copy_this : :class:`pybamm.Solution`, optional
+        A solution to copy, if provided. Default is None.
 
     """
 
-    def __init__(self, t, y, t_event=None, y_event=None, termination="final time"):
-        self.t = t
-        self.y = y
-        self.t_event = t_event
-        self.y_event = y_event
-        self.termination = termination
+    def __init__(
+        self,
+        t,
+        y,
+        t_event=None,
+        y_event=None,
+        termination="final time",
+        copy_this=None,
+    ):
+        self._t = t
+        self._y = y
+        self._t_event = t_event
+        self._y_event = y_event
+        self._termination = termination
         # initialize empty inputs and model, to be populated later
-        self.inputs = {}
-        self._model = None
+        if copy_this is None:
+            self._inputs = pybamm.FuzzyDict()
+            self._model = None
+            self.set_up_time = None
+            self.solve_time = None
+        else:
+            self._inputs = copy.copy(copy_this.inputs)
+            self._model = copy_this.model
+            self.set_up_time = copy_this.set_up_time
+            self.solve_time = copy_this.solve_time
 
         # initiaize empty variables and data
-        self._variables = {}
-        self.data = {}
+        self._variables = pybamm.FuzzyDict()
+        self.data = pybamm.FuzzyDict()
 
         # initialize empty known evals
-        self.known_evals = defaultdict(dict)
+        self._known_evals = defaultdict(dict)
         for time in t:
-            self.known_evals[time] = {}
+            self._known_evals[time] = {}
 
     @property
     def t(self):
         "Times at which the solution is evaluated"
         return self._t
 
-    @t.setter
-    def t(self, value):
-        "Updates the solution times"
-        self._t = value
-
     @property
     def y(self):
         "Values of the solution"
         return self._y
 
-    @y.setter
-    def y(self, value):
-        "Updates the solution values"
-        self._y = value
-
     @property
     def inputs(self):
         "Values of the inputs"
@@ -124,44 +134,6 @@ class Solution(object):
         "Updates the reason for termination"
         self._termination = value
 
-    def __add__(self, other):
-        "See :meth:`Solution.append`"
-        self.append(other)
-        return self
-
-    def append(self, solution, start_index=1):
-        """
-
-        Appends solution.t and solution.y onto self.t and self.y.
-
-        Note: by default this process removes the initial time and state of solution to
-        avoid duplicate times and states being stored (self.t[-1] is equal to
-        solution.t[0], and self.y[:, -1] is equal to solution.y[:, 0]). Set the optional
-        argument ``start_index`` to override this behavior
-
-        """
-        # Update t, y and inputs
-        self.t = np.concatenate((self.t, solution.t[start_index:]))
-        self.y = np.concatenate((self.y, solution.y[:, start_index:]), axis=1)
-        for name, inp in self.inputs.items():
-            solution_inp = solution.inputs[name]
-            if isinstance(solution_inp, numbers.Number):
-                solution_inp = solution_inp * np.ones_like(solution.t)
-            self.inputs[name] = np.concatenate((inp, solution_inp[start_index:]))
-        # Update solution time
-        self.solve_time += solution.solve_time
-        # Update termination
-        self.termination = solution.termination
-        self.t_event = solution.t_event
-        self.y_event = solution.y_event
-
-        # Update known_evals
-        for t, evals in solution.known_evals.items():
-            self.known_evals[t].update(evals)
-        # Recompute existing variables
-        for var in self._variables.keys():
-            self.update(var)
-
     @property
     def total_time(self):
         return self.set_up_time + self.solve_time
@@ -175,12 +147,12 @@ class Solution(object):
         for key in variables:
             pybamm.logger.debug("Post-processing {}".format(key))
             var = pybamm.ProcessedVariable(
-                self.model.variables[key], self, self.known_evals
+                self.model.variables[key], self, self._known_evals
             )
 
             # Update known_evals in order to process any other variables faster
             for t in var.known_evals:
-                self.known_evals[t].update(var.known_evals[t])
+                self._known_evals[t].update(var.known_evals[t])
 
             # Save variable and data
             self._variables[key] = var
@@ -230,3 +202,93 @@ class Solution(object):
         with open(filename, "wb") as f:
             pickle.dump(self.data, f, pickle.HIGHEST_PROTOCOL)
 
+
+class Solution(_BaseSolution):
+    """
+    Class extending the base solution, with additional functionality for concatenating
+    different solutions together
+
+    **Extends**: :class:`_BaseSolution`
+
+    """
+
+    def __init__(
+        self, t, y, t_event=None, y_event=None, termination="final time",
+    ):
+        super().__init__(t, y, t_event, y_event, termination)
+
+    @property
+    def sub_solutions(self):
+        "List of sub solutions that have been concatenated to form the full solution"
+        try:
+            return self._sub_solutions
+        except AttributeError:
+            raise AttributeError(
+                "sub solutions are only created once other solutions have been appended"
+            )
+
+    def __add__(self, other):
+        "See :meth:`Solution.append`"
+        self.append(other, create_sub_solutions=True)
+        return self
+
+    def append(self, solution, start_index=1, create_sub_solutions=False):
+        """
+        Appends solution.t and solution.y onto self.t and self.y.
+
+        Note: by default this process removes the initial time and state of solution to
+        avoid duplicate times and states being stored (self.t[-1] is equal to
+        solution.t[0], and self.y[:, -1] is equal to solution.y[:, 0]). Set the optional
+        argument ``start_index`` to override this behavior
+        """
+        # Create sub-solutions if necessary
+        # sub-solutions are 'BaseSolution' objects, which have slightly reduced
+        # functionality compared to normal solutions (can't append other solutions)
+        if create_sub_solutions and not hasattr(self, "_sub_solutions"):
+            self._sub_solutions = [
+                _BaseSolution(
+                    self.t,
+                    self.y,
+                    self.t_event,
+                    self.y_event,
+                    self.termination,
+                    copy_this=self,
+                )
+            ]
+
+        # (Create and) update sub-solutions
+        # Create a list of sub-solutions, which are simpler BaseSolution classes
+
+        # Update t, y and inputs
+        self._t = np.concatenate((self._t, solution.t[start_index:]))
+        self._y = np.concatenate((self._y, solution.y[:, start_index:]), axis=1)
+        for name, inp in self.inputs.items():
+            solution_inp = solution.inputs[name]
+            self.inputs[name] = np.concatenate((inp, solution_inp[start_index:]))
+        # Update solution time
+        self.solve_time += solution.solve_time
+        # Update termination
+        self._termination = solution.termination
+        self._t_event = solution._t_event
+        self._y_event = solution._y_event
+
+        # Update known_evals
+        for t, evals in solution._known_evals.items():
+            self._known_evals[t].update(evals)
+        # Recompute existing variables
+        for var in self._variables.keys():
+            self.update(var)
+
+        # Append sub_solutions
+        if create_sub_solutions:
+            self._sub_solutions.append(
+                _BaseSolution(
+                    solution.t,
+                    solution.y,
+                    solution.t_event,
+                    solution.y_event,
+                    solution.termination,
+                    copy_this=solution,
+                )
+            )
+


Test patch:
diff --git a/tests/integration/test_experiments.py b/tests/integration/test_experiments.py
new file mode 100644
index 000000000..8f3a3b688
--- /dev/null
+++ b/tests/integration/test_experiments.py
@@ -0,0 +1,60 @@
+#
+# Test some experiments
+#
+import pybamm
+import numpy as np
+import unittest
+
+
+class TestExperiments(unittest.TestCase):
+    def test_discharge_rest_charge(self):
+        experiment = pybamm.Experiment(
+            [
+                "Discharge at C/2 for 1 hour",
+                "Rest for 1 hour",
+                "Charge at C/2 for 1 hour",
+            ],
+            period="0.25 hours",
+        )
+        model = pybamm.lithium_ion.SPM()
+        sim = pybamm.Simulation(
+            model, experiment=experiment, solver=pybamm.CasadiSolver()
+        )
+        sim.solve()
+        np.testing.assert_array_almost_equal(
+            sim._solution["Time [h]"].entries, np.linspace(0, 3, 13)
+        )
+        cap = model.default_parameter_values["Cell capacity [A.h]"]
+        np.testing.assert_array_almost_equal(
+            sim._solution["Current [A]"].entries,
+            [cap / 2] * 5 + [0] * 4 + [-cap / 2] * 4,
+        )
+
+    def test_gitt(self):
+        experiment = pybamm.Experiment(
+            ["Discharge at C/20 for 1 hour", "Rest for 1 hour"] * 10,
+            period="6 minutes",
+        )
+        model = pybamm.lithium_ion.SPM()
+        sim = pybamm.Simulation(
+            model, experiment=experiment, solver=pybamm.CasadiSolver()
+        )
+        sim.solve()
+        np.testing.assert_array_almost_equal(
+            sim._solution["Time [h]"].entries, np.arange(0, 20.01, 0.1)
+        )
+        cap = model.default_parameter_values["Cell capacity [A.h]"]
+        np.testing.assert_array_almost_equal(
+            sim._solution["Current [A]"].entries,
+            [cap / 20] * 11 + [0] * 10 + ([cap / 20] * 10 + [0] * 10) * 9,
+        )
+
+
+if __name__ == "__main__":
+    print("Add -v for more debug output")
+    import sys
+
+    if "-v" in sys.argv:
+        debug = True
+    pybamm.settings.debug_mode = True
+    unittest.main()
diff --git a/tests/integration/test_models/test_submodels/test_external_circuit/test_function_control.py b/tests/integration/test_models/test_submodels/test_external_circuit/test_function_control.py
index 538846c37..1dec8e3e8 100644
--- a/tests/integration/test_models/test_submodels/test_external_circuit/test_function_control.py
+++ b/tests/integration/test_models/test_submodels/test_external_circuit/test_function_control.py
@@ -8,17 +8,14 @@ import unittest
 
 class TestFunctionControl(unittest.TestCase):
     def test_constant_current(self):
-        class ConstantCurrent:
-            num_switches = 0
-
-            def __call__(self, variables):
-                I = variables["Current [A]"]
-                return I + 1
+        def constant_current(variables):
+            I = variables["Current [A]"]
+            return I + 1
 
         # load models
         models = [
             pybamm.lithium_ion.SPM(),
-            pybamm.lithium_ion.SPM({"operating mode": ConstantCurrent()}),
+            pybamm.lithium_ion.SPM({"operating mode": constant_current}),
         ]
 
         # load parameter values and process models and geometry
@@ -56,17 +53,14 @@ class TestFunctionControl(unittest.TestCase):
         )
 
     def test_constant_voltage(self):
-        class ConstantVoltage:
-            num_switches = 0
-
-            def __call__(self, variables):
-                V = variables["Terminal voltage [V]"]
-                return V - 4.1
+        def constant_voltage(variables):
+            V = variables["Terminal voltage [V]"]
+            return V - 4.1
 
         # load models
         models = [
             pybamm.lithium_ion.SPM({"operating mode": "voltage"}),
-            pybamm.lithium_ion.SPM({"operating mode": ConstantVoltage()}),
+            pybamm.lithium_ion.SPM({"operating mode": constant_voltage}),
         ]
 
         # load parameter values and process models and geometry
@@ -102,18 +96,15 @@ class TestFunctionControl(unittest.TestCase):
         np.testing.assert_array_almost_equal(abs((I1 - I0) / I0), 0, decimal=1)
 
     def test_constant_power(self):
-        class ConstantPower:
-            num_switches = 0
-
-            def __call__(self, variables):
-                I = variables["Current [A]"]
-                V = variables["Terminal voltage [V]"]
-                return I * V - 4
+        def constant_power(variables):
+            I = variables["Current [A]"]
+            V = variables["Terminal voltage [V]"]
+            return I * V - 4
 
         # load models
         models = [
             pybamm.lithium_ion.SPM({"operating mode": "power"}),
-            pybamm.lithium_ion.SPM({"operating mode": ConstantPower()}),
+            pybamm.lithium_ion.SPM({"operating mode": constant_power}),
         ]
 
         # load parameter values and process models and geometry
diff --git a/tests/unit/test_experiments/__init__.py b/tests/unit/test_experiments/__init__.py
new file mode 100644
index 000000000..e69de29bb
diff --git a/tests/unit/test_experiments/test_experiment.py b/tests/unit/test_experiments/test_experiment.py
new file mode 100644
index 000000000..dbde2faa6
--- /dev/null
+++ b/tests/unit/test_experiments/test_experiment.py
@@ -0,0 +1,125 @@
+#
+# Test the base experiment class
+#
+import pybamm
+import unittest
+
+
+class TestExperiment(unittest.TestCase):
+    def test_read_strings(self):
+        experiment = pybamm.Experiment(
+            [
+                "Discharge at 1C for 0.5 hours",
+                "Discharge at C/20 for 0.5 hours",
+                "Charge at 0.5 C for 45 minutes",
+                "Discharge at 1 A for 0.5 hours",
+                "Charge at 200 mA for 45 minutes (1 minute period)",
+                "Discharge at 1W for 0.5 hours",
+                "Charge at 200 mW for 45 minutes",
+                "Rest for 10 minutes (5 minute period)",
+                "Hold at 1V for 20 seconds",
+                "Charge at 1 C until 4.1V",
+                "Hold at 4.1 V until 50mA",
+                "Hold at 3V until C/50",
+                "Discharge at C/3 for 2 hours or until 2.5 V",
+            ],
+            {"test": "test"},
+            period="20 seconds",
+        )
+        self.assertEqual(
+            experiment.operating_conditions,
+            [
+                (1, "C", 1800.0, 20.0),
+                (0.05, "C", 1800.0, 20.0),
+                (-0.5, "C", 2700.0, 20.0),
+                (1, "A", 1800.0, 20.0),
+                (-0.2, "A", 2700.0, 60.0),
+                (1, "W", 1800.0, 20.0),
+                (-0.2, "W", 2700.0, 20.0),
+                (0, "A", 600.0, 300.0),
+                (1, "V", 20.0, 20.0),
+                (-1, "C", None, 20.0),
+                (4.1, "V", None, 20.0),
+                (3, "V", None, 20.0),
+                (1 / 3, "C", 7200.0, 20.0),
+            ],
+        )
+        self.assertEqual(
+            experiment.events,
+            [
+                None,
+                None,
+                None,
+                None,
+                None,
+                None,
+                None,
+                None,
+                None,
+                (4.1, "V"),
+                (0.05, "A"),
+                (0.02, "C"),
+                (2.5, "V"),
+            ],
+        )
+        self.assertEqual(experiment.parameters, {"test": "test"})
+        self.assertEqual(experiment.period, 20)
+
+    def test_read_strings_repeat(self):
+        experiment = pybamm.Experiment(
+            ["Discharge at 10 mA for 0.5 hours"]
+            + ["Charge at 0.5 C for 45 minutes", "Hold at 1 V for 20 seconds"] * 2,
+        )
+        self.assertEqual(
+            experiment.operating_conditions,
+            [
+                (0.01, "A", 1800.0, 60),
+                (-0.5, "C", 2700.0, 60),
+                (1, "V", 20.0, 60),
+                (-0.5, "C", 2700.0, 60),
+                (1, "V", 20.0, 60),
+            ],
+        )
+        self.assertEqual(experiment.period, 60)
+
+    def test_str_repr(self):
+        conds = ["Discharge at 1 C for 20 seconds", "Charge at 0.5 W for 10 minutes"]
+        experiment = pybamm.Experiment(conds)
+        self.assertEqual(str(experiment), str(conds))
+        self.assertEqual(
+            repr(experiment),
+            "pybamm.Experiment(['Discharge at 1 C for 20 seconds'"
+            + ", 'Charge at 0.5 W for 10 minutes'])",
+        )
+
+    def test_bad_strings(self):
+        with self.assertRaisesRegex(
+            TypeError, "Operating conditions should be strings"
+        ):
+            pybamm.Experiment([1, 2, 3])
+        with self.assertRaisesRegex(ValueError, "Operating conditions must contain"):
+            pybamm.Experiment(["Discharge at 1 A at 2 hours"])
+        with self.assertRaisesRegex(ValueError, "instruction must be"):
+            pybamm.Experiment(["Run at 1 A for 2 hours"])
+        with self.assertRaisesRegex(
+            ValueError, "Instruction 'Run at at 1 A' not recognized"
+        ):
+            pybamm.Experiment(["Run at at 1 A for 2 hours"])
+        with self.assertRaisesRegex(ValueError, "units must be"):
+            pybamm.Experiment(["Discharge at 1 B for 2 hours"])
+        with self.assertRaisesRegex(ValueError, "time units must be"):
+            pybamm.Experiment(["Discharge at 1 A for 2 years"])
+        with self.assertRaisesRegex(
+            TypeError, "experimental parameters should be a dictionary"
+        ):
+            pybamm.Experiment([], "not a dictionary")
+
+
+if __name__ == "__main__":
+    print("Add -v for more debug output")
+    import sys
+
+    if "-v" in sys.argv:
+        debug = True
+    pybamm.settings.debug_mode = True
+    unittest.main()
diff --git a/tests/unit/test_experiments/test_simulation_with_experiment.py b/tests/unit/test_experiments/test_simulation_with_experiment.py
new file mode 100644
index 000000000..08d7bc50f
--- /dev/null
+++ b/tests/unit/test_experiments/test_simulation_with_experiment.py
@@ -0,0 +1,106 @@
+#
+# Test setting up a simulation with an experiment
+#
+import pybamm
+import unittest
+
+
+class TestSimulationExperiment(unittest.TestCase):
+    def test_set_up(self):
+        experiment = pybamm.Experiment(
+            [
+                "Discharge at C/20 for 1 hour",
+                "Charge at 1 A until 4.1 V",
+                "Hold at 4.1 V until 50 mA",
+                "Discharge at 2 W for 1 hour",
+            ],
+        )
+        model = pybamm.lithium_ion.DFN()
+        sim = pybamm.Simulation(model, experiment=experiment)
+
+        self.assertEqual(sim.experiment, experiment)
+        self.assertEqual(
+            sim._experiment_inputs[0]["Current input [A]"],
+            1 / 20 * model.default_parameter_values["Cell capacity [A.h]"],
+        )
+        self.assertEqual(sim._experiment_inputs[0]["Current switch"], 1)
+        self.assertEqual(sim._experiment_inputs[0]["Voltage switch"], 0)
+        self.assertEqual(sim._experiment_inputs[0]["Power switch"], 0)
+        self.assertEqual(sim._experiment_inputs[0]["Current cut-off [A]"], -1e10)
+        self.assertEqual(sim._experiment_inputs[0]["Voltage cut-off [V]"], -1e10)
+        self.assertEqual(sim._experiment_inputs[1]["Current input [A]"], -1)
+        self.assertEqual(sim._experiment_inputs[1]["Current switch"], 1)
+        self.assertEqual(sim._experiment_inputs[1]["Voltage switch"], 0)
+        self.assertEqual(sim._experiment_inputs[1]["Power switch"], 0)
+        self.assertEqual(sim._experiment_inputs[1]["Current cut-off [A]"], -1e10)
+        self.assertEqual(sim._experiment_inputs[1]["Voltage cut-off [V]"], 4.1)
+        self.assertEqual(sim._experiment_inputs[2]["Current switch"], 0)
+        self.assertEqual(sim._experiment_inputs[2]["Voltage switch"], 1)
+        self.assertEqual(sim._experiment_inputs[2]["Power switch"], 0)
+        self.assertEqual(sim._experiment_inputs[2]["Voltage input [V]"], 4.1)
+        self.assertEqual(sim._experiment_inputs[2]["Current cut-off [A]"], 0.05)
+        self.assertEqual(sim._experiment_inputs[2]["Voltage cut-off [V]"], -1e10)
+        self.assertEqual(sim._experiment_inputs[3]["Current switch"], 0)
+        self.assertEqual(sim._experiment_inputs[3]["Voltage switch"], 0)
+        self.assertEqual(sim._experiment_inputs[3]["Power switch"], 1)
+        self.assertEqual(sim._experiment_inputs[3]["Power input [W]"], 2)
+        self.assertEqual(sim._experiment_inputs[3]["Current cut-off [A]"], -1e10)
+        self.assertEqual(sim._experiment_inputs[3]["Voltage cut-off [V]"], -1e10)
+
+        tau = sim._parameter_values.evaluate(model.timescale)
+        self.assertEqual(
+            sim._experiment_times,
+            [t / tau for t in [3600, 7 * 24 * 3600, 7 * 24 * 3600, 3600]],
+        )
+
+        self.assertIn(
+            "Current cut-off (positive) [A] [experiment]",
+            [event.name for event in sim.model.events],
+        )
+        self.assertIn(
+            "Current cut-off (negative) [A] [experiment]",
+            [event.name for event in sim.model.events],
+        )
+        self.assertIn(
+            "Voltage cut-off [V] [experiment]",
+            [event.name for event in sim.model.events],
+        )
+
+        # fails if trying to set up with something that isn't an experiment
+        with self.assertRaisesRegex(TypeError, "experiment must be"):
+            pybamm.Simulation(model, experiment=0)
+
+    def test_run_experiment(self):
+        experiment = pybamm.Experiment(
+            [
+                "Discharge at C/20 for 1 hour",
+                "Charge at 1 A until 4.1 V",
+                "Hold at 4.1 V until C/2",
+                "Discharge at 2 W for 1 hour",
+            ],
+        )
+        model = pybamm.lithium_ion.SPM()
+        sim = pybamm.Simulation(model, experiment=experiment)
+        sim.solve()
+        self.assertEqual(sim._solution.termination, "final time")
+
+    def test_run_experiment_breaks_early(self):
+        experiment = pybamm.Experiment(["Discharge at 2 C for 1 hour"])
+        model = pybamm.lithium_ion.SPM()
+        sim = pybamm.Simulation(model, experiment=experiment)
+        pybamm.set_logging_level("ERROR")
+        # giving the time, should get ignored
+        t_eval = [0, 1]
+        sim.solve(t_eval)
+        pybamm.set_logging_level("WARNING")
+        self.assertIn("event", sim._solution.termination)
+
+
+if __name__ == "__main__":
+    print("Add -v for more debug output")
+    import sys
+
+    if "-v" in sys.argv:
+        debug = True
+    pybamm.settings.debug_mode = True
+    unittest.main()
diff --git a/tests/unit/test_models/test_base_model.py b/tests/unit/test_models/test_base_model.py
index 7feb27350..b9c5e6c70 100644
--- a/tests/unit/test_models/test_base_model.py
+++ b/tests/unit/test_models/test_base_model.py
@@ -401,6 +401,10 @@ class TestStandardBatteryBaseModel(unittest.TestCase):
         )
         os.chdir(cwd)
 
+    def test_timescale(self):
+        model = pybamm.BaseModel()
+        self.assertEqual(model.timescale.evaluate(), 1)
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")
diff --git a/tests/unit/test_models/test_full_battery_models/test_lead_acid/test_loqs.py b/tests/unit/test_models/test_full_battery_models/test_lead_acid/test_loqs.py
index 59697e2b3..0c7757d17 100644
--- a/tests/unit/test_models/test_full_battery_models/test_lead_acid/test_loqs.py
+++ b/tests/unit/test_models/test_full_battery_models/test_lead_acid/test_loqs.py
@@ -167,15 +167,12 @@ class TestLeadAcidLOQSExternalCircuits(unittest.TestCase):
         model.check_well_posedness()
 
     def test_well_posed_function(self):
-        class ExternalCircuitFunction:
-            num_switches = 0
+        def external_circuit_function(variables):
+            I = variables["Current [A]"]
+            V = variables["Terminal voltage [V]"]
+            return V + I - pybamm.FunctionParameter("Function", pybamm.t)
 
-            def __call__(self, variables):
-                I = variables["Current [A]"]
-                V = variables["Terminal voltage [V]"]
-                return V + I - pybamm.FunctionParameter("Function", pybamm.t)
-
-        options = {"operating mode": ExternalCircuitFunction()}
+        options = {"operating mode": external_circuit_function}
         model = pybamm.lead_acid.LOQS(options)
         model.check_well_posedness()
 
diff --git a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_spm.py b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_spm.py
index 3df0760bb..87da7b98b 100644
--- a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_spm.py
+++ b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_spm.py
@@ -190,15 +190,12 @@ class TestSPMExternalCircuits(unittest.TestCase):
         model.check_well_posedness()
 
     def test_well_posed_function(self):
-        class ExternalCircuitFunction:
-            num_switches = 0
+        def external_circuit_function(variables):
+            I = variables["Current [A]"]
+            V = variables["Terminal voltage [V]"]
+            return V + I - pybamm.FunctionParameter("Function", pybamm.t)
 
-            def __call__(self, variables):
-                I = variables["Current [A]"]
-                V = variables["Terminal voltage [V]"]
-                return V + I - pybamm.FunctionParameter("Function", pybamm.t)
-
-        options = {"operating mode": ExternalCircuitFunction()}
+        options = {"operating mode": external_circuit_function}
         model = pybamm.lithium_ion.SPM(options)
         model.check_well_posedness()
 
diff --git a/tests/unit/test_models/test_submodels/test_external_circuit/test_function_control.py b/tests/unit/test_models/test_submodels/test_external_circuit/test_function_control.py
index ec961b431..4a06002ba 100644
--- a/tests/unit/test_models/test_submodels/test_external_circuit/test_function_control.py
+++ b/tests/unit/test_models/test_submodels/test_external_circuit/test_function_control.py
@@ -6,22 +6,17 @@ import tests
 import unittest
 
 
-class ExternalCircuitFunction:
-    num_switches = 0
-
-    def __call__(self, variables):
-        I = variables["Current [A]"]
-        V = variables["Terminal voltage [V]"]
-        return (
-            V + I - pybamm.FunctionParameter("Current plus voltage function", pybamm.t)
-        )
+def external_circuit_function(variables):
+    I = variables["Current [A]"]
+    V = variables["Terminal voltage [V]"]
+    return V + I - pybamm.FunctionParameter("Current plus voltage function", pybamm.t)
 
 
 class TestFunctionControl(unittest.TestCase):
     def test_public_functions(self):
         param = pybamm.standard_parameters_lithium_ion
         submodel = pybamm.external_circuit.FunctionControl(
-            param, ExternalCircuitFunction()
+            param, external_circuit_function
         )
         variables = {"Terminal voltage [V]": pybamm.Scalar(0)}
         std_tests = tests.StandardSubModelTests(submodel, variables)
diff --git a/tests/unit/test_solvers/test_solution.py b/tests/unit/test_solvers/test_solution.py
index 7aa4db34f..dc97165e7 100644
--- a/tests/unit/test_solvers/test_solution.py
+++ b/tests/unit/test_solvers/test_solution.py
@@ -19,25 +19,47 @@ class TestSolution(unittest.TestCase):
         self.assertEqual(sol.inputs, {})
         self.assertEqual(sol.model, None)
 
+        with self.assertRaisesRegex(AttributeError, "sub solutions"):
+            print(sol.sub_solutions)
+
     def test_append(self):
         # Set up first solution
         t1 = np.linspace(0, 1)
         y1 = np.tile(t1, (20, 1))
         sol1 = pybamm.Solution(t1, y1)
         sol1.solve_time = 1.5
-        sol1.inputs = {}
+        sol1.model = pybamm.BaseModel()
+        sol1.inputs = {"a": 1}
 
         # Set up second solution
         t2 = np.linspace(1, 2)
         y2 = np.tile(t2, (20, 1))
         sol2 = pybamm.Solution(t2, y2)
         sol2.solve_time = 1
-        sol1.append(sol2)
+        sol2.inputs = {"a": 2}
+        sol1.append(sol2, create_sub_solutions=True)
 
         # Test
         self.assertEqual(sol1.solve_time, 2.5)
         np.testing.assert_array_equal(sol1.t, np.concatenate([t1, t2[1:]]))
         np.testing.assert_array_equal(sol1.y, np.concatenate([y1, y2[:, 1:]], axis=1))
+        np.testing.assert_array_equal(
+            sol1.inputs["a"],
+            np.concatenate([1 * np.ones_like(t1), 2 * np.ones_like(t2[1:])]),
+        )
+
+        # Test sub-solutions
+        self.assertEqual(len(sol1.sub_solutions), 2)
+        np.testing.assert_array_equal(sol1.sub_solutions[0].t, t1)
+        np.testing.assert_array_equal(sol1.sub_solutions[1].t, t2)
+        self.assertEqual(sol1.sub_solutions[0].model, sol1.model)
+        np.testing.assert_array_equal(
+            sol1.sub_solutions[0].inputs["a"], 1 * np.ones_like(t1)
+        )
+        self.assertEqual(sol1.sub_solutions[1].model, sol2.model)
+        np.testing.assert_array_equal(
+            sol1.sub_solutions[1].inputs["a"], 2 * np.ones_like(t2)
+        )
 
     def test_total_time(self):
         sol = pybamm.Solution([], None)