Repository: pybamm-team/PyBaMM

Issue:
[Bug]: `starting_solution` does not work if last cycle is not saved
### PyBaMM Version

23.5

### Python Version

3.9.17

### Describe the bug

When using `starting_solution` in solving a simulation, if the last cycle of the `starting_solution` is not saved, then the initialisation doesn't work: the new solution starts at `t=0` and does not contain the previous states (not even the first cycle which should be always saved).

My guess is that when checking the solution for the last cycle it finds it is None and thus it treats it as if no `starting_solution` has been provided.

### Steps to Reproduce

```python3
        model = pybamm.lithium_ion.SPM()

        experiment = pybamm.Experiment([
            # pybamm.step.string("Discharge at C/2 for 10 minutes", start_time=datetime(2023, 1, 1, 8, 0, 0)),
            pybamm.step.string("Discharge at C/2 for 10 minutes"),
            pybamm.step.string("Rest for 10 minutes"),
        ])

        sim = pybamm.Simulation(model, experiment=experiment)
        solution = sim.solve(save_at_cycles=[1])

        experiment = pybamm.Experiment([
            # pybamm.step.string("Discharge at C/2 for 20 minutes", start_time=datetime(2023, 1, 1, 10, 0, 0)),
            pybamm.step.string("Discharge at C/2 for 20 minutes"),
            pybamm.step.string("Rest for 20 minutes"),
        ])

        sim = pybamm.Simulation(model, experiment=experiment)
        new_solution = sim.solve(calc_esoh=False, starting_solution=solution)
```

Plotting `new_solution` yields:
![image](https://github.com/pybamm-team/PyBaMM/assets/28443643/38b12725-90a5-40ac-9daf-1aeea4cc1291)

If, instead, I do `save_at_cycles=None` in the first solve, I get the expected result:
![image](https://github.com/pybamm-team/PyBaMM/assets/28443643/b61ce54f-7157-415f-b81c-03b75fd181d1)


### Relevant log output

_No response_

Discussion:

Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 6488dbef4..ccd396050 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -13,10 +13,14 @@
 - Fixed a bug where the "basic" lithium-ion models gave incorrect results when using nonlinear particle diffusivity ([#3207](https://github.com/pybamm-team/PyBaMM/pull/3207))
 - Particle size distributions now work with SPMe and NewmanTobias models ([#3207](https://github.com/pybamm-team/PyBaMM/pull/3207))
 - Fix to simulate c_rate steps with drive cycles ([#3186](https://github.com/pybamm-team/PyBaMM/pull/3186))
-- Parameters in `Prada2013` have been updated to better match those given in the paper, which is a 2.3 Ah cell, instead of the mix-and-match with the 1.1 Ah cell from Lain2019.
-- Error generated when invalid parameter values are passed. ([#3132](https://github.com/pybamm-team/PyBaMM/pull/3132))
-- Thevenin() model is now constructed with standard variables: `Time [s], Time [min], Time [h]` ([#3143](https://github.com/pybamm-team/PyBaMM/pull/3143))
+- Always save last cycle in experiment, to fix issues with `starting_solution` and `last_state` ([#3177](https://github.com/pybamm-team/PyBaMM/pull/3177))
+- Fix simulations with `starting_solution` to work with `start_time` experiments ([#3177](https://github.com/pybamm-team/PyBaMM/pull/3177))
 - Fix SEI Example Notebook ([#3166](https://github.com/pybamm-team/PyBaMM/pull/3166))
+- Thevenin() model is now constructed with standard variables: `Time [s]`, `Time [min]`, `Time [h]` ([#3143](https://github.com/pybamm-team/PyBaMM/pull/3143))
+- Error generated when invalid parameter values are passed ([#3132](https://github.com/pybamm-team/PyBaMM/pull/3132))
+- Parameters in `Prada2013` have been updated to better match those given in the paper, which is a 2.3 Ah cell, instead of the mix-and-match with the 1.1 Ah cell from Lain2019 ([#3096](https://github.com/pybamm-team/PyBaMM/pull/3096))
+
+
 
 ## Breaking changes
 
diff --git a/pybamm/simulation.py b/pybamm/simulation.py
index 52c192254..8511d01af 100644
--- a/pybamm/simulation.py
+++ b/pybamm/simulation.py
@@ -641,6 +641,21 @@ class Simulation:
                     starting_solution.all_first_states.copy()
                 )
 
+            # set simulation initial_start_time
+            if starting_solution is None:
+                initial_start_time = self.experiment.initial_start_time
+            else:
+                initial_start_time = starting_solution.initial_start_time
+
+            if (
+                initial_start_time is None
+                and self.experiment.initial_start_time is not None
+            ):
+                raise ValueError(
+                    "When using experiments with `start_time`, the starting_solution "
+                    "must have a `start_time` too."
+                )
+
             cycle_offset = len(starting_solution_cycles)
             all_cycle_solutions = starting_solution_cycles
             all_summary_variables = starting_solution_summary_variables
@@ -653,6 +668,49 @@ class Simulation:
             idx = 0
             num_cycles = len(self.experiment.cycle_lengths)
             feasible = True  # simulation will stop if experiment is infeasible
+
+            # Add initial padding rest if current time is earlier than first start time
+            # This could be the case when using a starting solution
+            if starting_solution is not None:
+                op_conds = self.experiment.operating_conditions_steps[0]
+                if op_conds.start_time is not None:
+                    rest_time = (
+                        op_conds.start_time
+                        - (
+                            initial_start_time
+                            + timedelta(seconds=float(current_solution.t[-1]))
+                        )
+                    ).total_seconds()
+                    if rest_time > pybamm.settings.step_start_offset:
+                        # logs["step operating conditions"] = "Initial rest for padding"
+                        # callbacks.on_step_start(logs)
+
+                        kwargs["inputs"] = {
+                            **user_inputs,
+                            "Ambient temperature [K]": (
+                                op_conds.temperature or self._original_temperature
+                            ),
+                            "start time": current_solution.t[-1],
+                        }
+                        steps = current_solution.cycles[-1].steps
+                        step_solution = current_solution.cycles[-1].steps[-1]
+
+                        step_solution_with_rest = self.run_padding_rest(
+                            kwargs, rest_time, step_solution
+                        )
+                        steps[-1] = step_solution + step_solution_with_rest
+
+                        cycle_solution, _, _ = pybamm.make_cycle_solution(
+                            steps, esoh_solver=esoh_solver, save_this_cycle=True
+                        )
+                        old_cycles = current_solution.cycles.copy()
+                        old_cycles[-1] = cycle_solution
+                        current_solution += step_solution_with_rest
+                        current_solution.cycles = old_cycles
+
+                        # Update _solution
+                        self._solution = current_solution
+
             for cycle_num, cycle_length in enumerate(
                 # tqdm is the progress bar.
                 tqdm.tqdm(
@@ -676,6 +734,8 @@ class Simulation:
                 save_this_cycle = (
                     # always save cycle 1
                     cycle_num == 1
+                    # always save last cycle
+                    or cycle_num == num_cycles
                     # None: save all cycles
                     or save_at_cycles is None
                     # list: save all cycles in the list
@@ -703,7 +763,7 @@ class Simulation:
                             (
                                 op_conds.end_time
                                 - (
-                                    self.experiment.initial_start_time
+                                    initial_start_time
                                     + timedelta(seconds=float(start_time))
                                 )
                             ).total_seconds(),
@@ -758,41 +818,25 @@ class Simulation:
                         rest_time = (
                             op_conds.next_start_time
                             - (
-                                self.experiment.initial_start_time
+                                initial_start_time
                                 + timedelta(seconds=float(step_solution.t[-1]))
                             )
                         ).total_seconds()
                         if rest_time > pybamm.settings.step_start_offset:
-                            start_time = step_solution.t[-1]
-                            # Let me know if you have a better name
-                            op_conds_str = "Rest for padding"
-                            model = self.op_conds_to_built_models[op_conds_str]
-                            solver = self.op_conds_to_built_solvers[op_conds_str]
-
                             logs["step number"] = (step_num, cycle_length)
-                            logs["step operating conditions"] = op_conds_str
+                            logs["step operating conditions"] = "Rest for padding"
                             callbacks.on_step_start(logs)
 
-                            ambient_temp = (
-                                op_conds.temperature or self._original_temperature
-                            )
                             kwargs["inputs"] = {
                                 **user_inputs,
-                                "Ambient temperature [K]": ambient_temp,
-                                "start time": start_time,
+                                "Ambient temperature [K]": (
+                                    op_conds.temperature or self._original_temperature
+                                ),
+                                "start time": step_solution.t[-1],
                             }
-                            # Make sure we take at least 2 timesteps
-                            # The period is hardcoded to 10 minutes, the user can
-                            # always override it by adding a rest step
-                            npts = max(int(round(rest_time / 600)) + 1, 2)
-
-                            step_solution_with_rest = solver.step(
-                                step_solution,
-                                model,
-                                rest_time,
-                                npts=npts,
-                                save=False,
-                                **kwargs,
+
+                            step_solution_with_rest = self.run_padding_rest(
+                                kwargs, rest_time, step_solution
                             )
                             step_solution += step_solution_with_rest
 
@@ -896,8 +940,30 @@ class Simulation:
 
             callbacks.on_experiment_end(logs)
 
+            # record initial_start_time of the solution
+            self.solution.initial_start_time = initial_start_time
+
         return self.solution
 
+    def run_padding_rest(self, kwargs, rest_time, step_solution):
+        model = self.op_conds_to_built_models["Rest for padding"]
+        solver = self.op_conds_to_built_solvers["Rest for padding"]
+
+        # Make sure we take at least 2 timesteps. The period is hardcoded to 10
+        # minutes,the user can always override it by adding a rest step
+        npts = max(int(round(rest_time / 600)) + 1, 2)
+
+        step_solution_with_rest = solver.step(
+            step_solution,
+            model,
+            rest_time,
+            npts=npts,
+            save=False,
+            **kwargs,
+        )
+
+        return step_solution_with_rest
+
     def step(
         self, dt, solver=None, npts=2, save=True, starting_solution=None, **kwargs
     ):
diff --git a/pybamm/solvers/solution.py b/pybamm/solvers/solution.py
index 411341f88..9f51bd6c2 100644
--- a/pybamm/solvers/solution.py
+++ b/pybamm/solvers/solution.py
@@ -131,6 +131,9 @@ class Solution(object):
         # Initialize empty summary variables
         self._summary_variables = None
 
+        # Initialise initial start time
+        self.initial_start_time = None
+
         # Solution now uses CasADi
         pybamm.citations.register("Andersson2019")
 
@@ -421,6 +424,15 @@ class Solution(object):
     def summary_variables(self):
         return self._summary_variables
 
+    @property
+    def initial_start_time(self):
+        return self._initial_start_time
+
+    @initial_start_time.setter
+    def initial_start_time(self, value):
+        """Updates the reason for termination"""
+        self._initial_start_time = value
+
     def set_summary_variables(self, all_summary_variables):
         summary_variables = {var: [] for var in all_summary_variables[0]}
         for sum_vars in all_summary_variables:


Test patch:
diff --git a/tests/unit/test_experiments/test_simulation_with_experiment.py b/tests/unit/test_experiments/test_simulation_with_experiment.py
index 460d82c8c..13531281d 100644
--- a/tests/unit/test_experiments/test_simulation_with_experiment.py
+++ b/tests/unit/test_experiments/test_simulation_with_experiment.py
@@ -382,9 +382,9 @@ class TestSimulationExperiment(TestCase):
             solver=pybamm.CasadiSolver("fast with events"), save_at_cycles=[3, 4, 5, 9]
         )
         # Note offset by 1 (0th cycle is cycle 1)
-        for cycle_num in [1, 5, 6, 7, 9]:
+        for cycle_num in [1, 5, 6, 7]:
             self.assertIsNone(sol.cycles[cycle_num])
-        for cycle_num in [0, 2, 3, 4, 8]:
+        for cycle_num in [0, 2, 3, 4, 8, 9]:  # first & last cycle always saved
             self.assertIsNotNone(sol.cycles[cycle_num])
         # Summary variables are not None
         self.assertIsNotNone(sol.summary_variables["Capacity [A.h]"])
@@ -605,7 +605,7 @@ class TestSimulationExperiment(TestCase):
             pybamm.lithium_ion.SPM,
         )
 
-    def test_run_time_stamped_experiment(self):
+    def test_run_start_time_experiment(self):
         model = pybamm.lithium_ion.SPM()
 
         # Test experiment is cut short if next_start_time is early
@@ -640,6 +640,95 @@ class TestSimulationExperiment(TestCase):
         sol = sim.solve(calc_esoh=False)
         self.assertEqual(sol["Time [s]"].entries[-1], 10800)
 
+    def test_starting_solution(self):
+        model = pybamm.lithium_ion.SPM()
+
+        experiment = pybamm.Experiment(
+            [
+                pybamm.step.string("Discharge at C/2 for 10 minutes"),
+                pybamm.step.string("Rest for 5 minutes"),
+                pybamm.step.string("Rest for 5 minutes"),
+            ]
+        )
+
+        sim = pybamm.Simulation(model, experiment=experiment)
+        solution = sim.solve(save_at_cycles=[1])
+
+        # test that the last state is correct (i.e. final cycle is saved)
+        self.assertEqual(solution.last_state.t[-1], 1200)
+
+        experiment = pybamm.Experiment(
+            [
+                pybamm.step.string("Discharge at C/2 for 20 minutes"),
+                pybamm.step.string("Rest for 20 minutes"),
+            ]
+        )
+
+        sim = pybamm.Simulation(model, experiment=experiment)
+        new_solution = sim.solve(calc_esoh=False, starting_solution=solution)
+
+        # test that the final time is correct (i.e. starting solution correctly set)
+        self.assertEqual(new_solution["Time [s]"].entries[-1], 3600)
+
+    def test_experiment_start_time_starting_solution(self):
+        model = pybamm.lithium_ion.SPM()
+
+        # Test error raised if starting_solution does not have start_time
+        experiment = pybamm.Experiment(
+            [pybamm.step.string("Discharge at C/2 for 10 minutes")]
+        )
+        sim = pybamm.Simulation(model, experiment=experiment)
+        solution = sim.solve()
+
+        experiment = pybamm.Experiment(
+            [
+                pybamm.step.string(
+                    "Discharge at C/2 for 10 minutes",
+                    start_time=datetime(1, 1, 1, 9, 0, 0),
+                )
+            ]
+        )
+
+        sim = pybamm.Simulation(model, experiment=experiment)
+        with self.assertRaisesRegex(ValueError, "experiments with `start_time`"):
+            sim.solve(starting_solution=solution)
+
+        # Test starting_solution works well with start_time
+        experiment = pybamm.Experiment(
+            [
+                pybamm.step.string(
+                    "Discharge at C/2 for 10 minutes",
+                    start_time=datetime(1, 1, 1, 8, 0, 0),
+                ),
+                pybamm.step.string(
+                    "Discharge at C/2 for 10 minutes",
+                    start_time=datetime(1, 1, 1, 8, 20, 0),
+                ),
+            ]
+        )
+
+        sim = pybamm.Simulation(model, experiment=experiment)
+        solution = sim.solve()
+
+        experiment = pybamm.Experiment(
+            [
+                pybamm.step.string(
+                    "Discharge at C/2 for 10 minutes",
+                    start_time=datetime(1, 1, 1, 9, 0, 0),
+                ),
+                pybamm.step.string(
+                    "Discharge at C/2 for 10 minutes",
+                    start_time=datetime(1, 1, 1, 9, 20, 0),
+                ),
+            ]
+        )
+
+        sim = pybamm.Simulation(model, experiment=experiment)
+        new_solution = sim.solve(starting_solution=solution)
+
+        # test that the final time is correct (i.e. starting solution correctly set)
+        self.assertEqual(new_solution["Time [s]"].entries[-1], 5400)
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")