Repository: pybamm-team/PyBaMM

Issue:
[Bug]: Wrong starting_solution for cycle-experiments
### PyBaMM Version

22.5

### Python Version

3.9.12

### Describe the bug

After initialising a cycle-simulation, solving it multiple times in a row leads to wrong starting_solution/initial values for the 2nd/3rd/... solution. In detail this occured by simulating a cycle-experiment for SEI-Growth by varying SEI-parameters for each simulation. The very first simulation contains the correct intial values, however the second simulation starts with a SEI-thickness close to the final solution of the previous simulation. This continues for each following simulation.

I think it is supposed to start with the same initial inner/outer SEI-thickness for each simulation. Note that the starting SEI-thickness of the (i+1)-th simulation is not exactly the final solution of the i-th simulation, but always a bit below that value. I`ve attached a MWE below for this behaviour, it is straightforward to reprocude.
Is this how it is supposed to work, or how can I set the starting_solution to the initial values of the simulation without re-initialising the simulation again?

### Steps to Reproduce

`
import pybamm
import matplotlib.pyplot as plt
import numpy as np
from pybamm.models.full_battery_models.lithium_ion.electrode_soh import get_initial_stoichiometries
from copy import deepcopy


seigrowth= "interstitial-diffusion limited"

model = pybamm.lithium_ion.SPMe({"SEI": seigrowth})


parameter_values = model.default_parameter_values
param = pybamm.lithium_ion.SPMe().param

c_n_max = parameter_values.evaluate(param.n.c_max)
c_p_max = parameter_values.evaluate(param.p.c_max)

V = 3.9
Vmin= 3.0
Vmax= 4.1

experiment = pybamm.Experiment(
    [ "Discharge at 0.1C until 3 V",
      f"Charge at 0.1C until {V}V",
     "Discharge at 0C for 10 hours",
     ]*4,)

parameter_values.update({'Lower voltage cut-off [V]': Vmin,'Upper voltage cut-off [V]': Vmax,})
parameter_values.update({'Initial inner SEI thickness [m]': 1e-12,
                            'Initial outer SEI thickness [m]': 1e-12,})
initial_soc = 0.9
x,y = get_initial_stoichiometries(initial_soc=initial_soc, parameter_values=parameter_values)

parameter_values.update({
    "Initial concentration in negative electrode [mol.m-3]": x * c_n_max,
    "Initial concentration in positive electrode [mol.m-3]": y * c_p_max,})

free_parameter0={'Inner SEI lithium interstitial diffusivity [m2.s-1]': 1e-16}
free_parameter1={'Inner SEI lithium interstitial diffusivity [m2.s-1]': 2e-16}
free_parameter2={'Inner SEI lithium interstitial diffusivity [m2.s-1]': 3e-16}
free_parameter3={'Inner SEI lithium interstitial diffusivity [m2.s-1]': 4e-16}
free_parameter_values=[free_parameter0,free_parameter1,free_parameter2,free_parameter3]
parameter_values.update({name: "[input]" for name in free_parameter0})


solver = pybamm.Simulation(model, experiment=experiment, parameter_values=parameter_values)

for i in range(0,4):
    solution = solver.solve(inputs=deepcopy(free_parameter_values[i]))
    plt.plot(solution.t,solution["X-averaged total SEI thickness [m]"].entries, label=model.name, color="red")

plt.show()
`

### Relevant log output


Discussion:

Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 69894e5a5..dd3159739 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -21,6 +21,7 @@
 
 ## Bug fixes
 
+-   Fixed a bug in `Simulation` that caused initial conditions to change when solving an experiment multiple times ([#2204](https://github.com/pybamm-team/PyBaMM/pull/2204))
 -   Fixed labels and ylims in `plot_voltage_components`([#2183](https://github.com/pybamm-team/PyBaMM/pull/2183))
 -   Fixed 2D interpolant ([#2180](https://github.com/pybamm-team/PyBaMM/pull/2180))
 -   Fixes a bug where the SPMe always builds even when `build=False` ([#2169](https://github.com/pybamm-team/PyBaMM/pull/2169))
diff --git a/pybamm/models/base_model.py b/pybamm/models/base_model.py
index 0cb0b22cb..d9355b524 100644
--- a/pybamm/models/base_model.py
+++ b/pybamm/models/base_model.py
@@ -579,7 +579,7 @@ class BaseModel:
 
         self.build_model_equations()
 
-    def set_initial_conditions_from(self, solution, inplace=True):
+    def set_initial_conditions_from(self, solution, inplace=True, return_type="model"):
         """
         Update initial conditions with the final states from a Solution object or from
         a dictionary.
@@ -590,14 +590,12 @@ class BaseModel:
         ----------
         solution : :class:`pybamm.Solution`, or dict
             The solution to use to initialize the model
-        inplace : bool
-            Whether to modify the model inplace or create a new model
+        inplace : bool, optional
+            Whether to modify the model inplace or create a new model (default True)
+        return_type : str, optional
+            Whether to return the model (default) or initial conditions ("ics")
         """
-        if inplace is True:
-            model = self
-        else:
-            model = self.new_copy()
-
+        initial_conditions = {}
         if isinstance(solution, pybamm.Solution):
             solution = solution.last_state
         for var in self.initial_conditions:
@@ -623,7 +621,7 @@ class BaseModel:
                     final_state_eval = final_state[:, :, -1].flatten(order="F")
                 else:
                     raise NotImplementedError("Variable must be 0D, 1D, or 2D")
-                model.initial_conditions[var] = pybamm.Vector(final_state_eval)
+                initial_conditions[var] = pybamm.Vector(final_state_eval)
             elif isinstance(var, pybamm.Concatenation):
                 children = []
                 for child in var.orphans:
@@ -645,7 +643,7 @@ class BaseModel:
                             "Variable in concatenation must be 1D"
                         )
                     children.append(final_state_eval)
-                model.initial_conditions[var] = pybamm.Vector(np.concatenate(children))
+                initial_conditions[var] = pybamm.Vector(np.concatenate(children))
 
             else:
                 raise NotImplementedError(
@@ -670,14 +668,26 @@ class BaseModel:
                     )
                 else:
                     slices.append(y_slices[symbol][0])
-            equations = list(model.initial_conditions.values())
+            equations = list(initial_conditions.values())
             # sort equations according to slices
             sorted_equations = [eq for _, eq in sorted(zip(slices, equations))]
-            model.concatenated_initial_conditions = pybamm.NumpyConcatenation(
+            concatenated_initial_conditions = pybamm.NumpyConcatenation(
                 *sorted_equations
             )
+        else:
+            concatenated_initial_conditions = None
+
+        if return_type == "model":
+            if inplace is True:
+                model = self
+            else:
+                model = self.new_copy()
 
-        return model
+            model.initial_conditions = initial_conditions
+            model.concatenated_initial_conditions = concatenated_initial_conditions
+            return model
+        elif return_type == "ics":
+            return initial_conditions, concatenated_initial_conditions
 
     def check_and_combine_dict(self, dict1, dict2):
         # check that the key ids are distinct
diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py
index 6fd62ff1e..15374a5ad 100644
--- a/pybamm/solvers/base_solver.py
+++ b/pybamm/solvers/base_solver.py
@@ -1296,9 +1296,12 @@ class BaseSolver(object):
                 # initialize with old solution
                 model.y0 = old_solution.all_ys[-1][:, -1]
             else:
-                model.y0 = model.set_initial_conditions_from(
-                    old_solution
-                ).concatenated_initial_conditions.evaluate(0, inputs=ext_and_inputs)
+                _, concatenated_initial_conditions = model.set_initial_conditions_from(
+                    old_solution, return_type="ics"
+                )
+                model.y0 = concatenated_initial_conditions.evaluate(
+                    0, inputs=ext_and_inputs
+                )
         set_up_time = timer.time()
 
         # (Re-)calculate consistent initial conditions


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 a8a01d4c6..b07ba90e5 100644
--- a/tests/unit/test_experiments/test_simulation_with_experiment.py
+++ b/tests/unit/test_experiments/test_simulation_with_experiment.py
@@ -124,6 +124,26 @@ class TestSimulationExperiment(unittest.TestCase):
         self.assertEqual(len(sol3.cycles), 2)
         os.remove("test_experiment.sav")
 
+    def test_run_experiment_multiple_times(self):
+        experiment = pybamm.Experiment(
+            [
+                (
+                    "Discharge at C/20 for 1 hour",
+                    "Charge at C/20 until 4.1 V",
+                )
+            ]
+            * 3
+        )
+        model = pybamm.lithium_ion.DFN()
+        sim = pybamm.Simulation(model, experiment=experiment)
+
+        # Test that solving twice gives the same solution (see #2193)
+        sol1 = sim.solve()
+        sol2 = sim.solve()
+        np.testing.assert_array_equal(
+            sol1["Terminal voltage [V]"].data, sol2["Terminal voltage [V]"].data
+        )
+
     def test_run_experiment_cccv_ode(self):
         experiment_2step = pybamm.Experiment(
             [