Solve experiment with drive cycle

CHANGELOG.md

@@ -2,6 +2,7 @@
 
 ## Features
 
+-   `Experiment`s with drive cycles can be solved ([#1793](https://github.com/pybamm-team/PyBaMM/pull/1793))
 -   Added surface area to volume ratio as a factor to the SEI equations ([#1790](https://github.com/pybamm-team/PyBaMM/pull/1790))
 -   Half-cell SPM and SPMe have been implemented ([#1731](https://github.com/pybamm-team/PyBaMM/pull/1731))
 

examples/scripts/experiment_drive_cycle.py

@@ -0,0 +1,52 @@
+#
+# Constant-current constant-voltage charge with US06 Drive Cycle using Experiment Class.
+#
+import pybamm
+import pandas as pd
+import os
+
+os.chdir(pybamm.__path__[0] + "/..")
+
+pybamm.set_logging_level("INFO")
+
+# import drive cycle from file
+drive_cycle_current = pd.read_csv(
+    "pybamm/input/drive_cycles/US06.csv", comment="#", header=None
+).to_numpy()
+
+
+# Map Drive Cycle
+def map_drive_cycle(x, min_op_value, max_op_value):
+    min_ip_value = x[:, 1].min()
+    max_ip_value = x[:, 1].max()
+    x[:, 1] = (x[:, 1] - min_ip_value) / (max_ip_value - min_ip_value) * (
+        max_op_value - min_op_value
+    ) + min_op_value
+    return x
+
+
+# Map current drive cycle to voltage and power
+drive_cycle_power = map_drive_cycle(drive_cycle_current, 1.5, 3.5)
+
+experiment = pybamm.Experiment(
+    [
+        "Charge at 1 A until 4.0 V",
+        "Hold at 4.0 V until 50 mA",
+        "Rest for 30 minutes",
+        "Run US06_A (A)",
+        "Rest for 30 minutes",
+        "Run US06_W (W)",
+        "Rest for 30 minutes",
+    ],
+    drive_cycles={
+        "US06_A": drive_cycle_current,
+        "US06_W": drive_cycle_power,
+    },
+)
+
+model = pybamm.lithium_ion.DFN()
+sim = pybamm.Simulation(model, experiment=experiment, solver=pybamm.CasadiSolver())
+sim.solve()
+
+# Show all plots
+sim.plot()

pybamm/experiments/experiment.py

@@ -195,6 +195,7 @@ def read_string(self, cond, drive_cycles):
                 "electric": op_CC["electric"] + op_CV["electric"],
                 "time": op_CV["time"],
                 "period": op_CV["period"],
+                "dc_data": None,
             }, event_CV
         # Read period
         if " period)" in cond:
@@ -223,26 +224,29 @@ def read_string(self, cond, drive_cycles):
                     drive_cycles[cond_list[1]], end_time
                 )
                 # Drive cycle as numpy array
+                dc_name = cond_list[1] + "_ext_{}".format(end_time)
                 dc_data = ext_drive_cycle
                 # Find the type of drive cycle ("A", "V", or "W")
                 typ = cond_list[2][1]
-                electric = (dc_data, typ)
+                electric = (dc_name, typ)
                 time = ext_drive_cycle[:, 0][-1]
                 period = np.min(np.diff(ext_drive_cycle[:, 0]))
                 events = None
             else:
                 # e.g. Run US06
                 # Drive cycle as numpy array
+                dc_name = cond_list[1]
                 dc_data = drive_cycles[cond_list[1]]
                 # Find the type of drive cycle ("A", "V", or "W")
                 typ = cond_list[2][1]
-                electric = (dc_data, typ)
+                electric = (dc_name, typ)
                 # Set time and period to 1 second for first step and
                 # then calculate the difference in consecutive time steps
                 time = drive_cycles[cond_list[1]][:, 0][-1]
                 period = np.min(np.diff(drive_cycles[cond_list[1]][:, 0]))
                 events = None
         else:
+            dc_data = None
             if "for" in cond and "or until" in cond:
                 # e.g. for 3 hours or until 4.2 V
                 cond_list = cond.split()
@@ -273,7 +277,12 @@ def read_string(self, cond, drive_cycles):
                     )
                 )
 
-        return {"electric": electric, "time": time, "period": period}, events
+        return {
+            "electric": electric,
+            "time": time,
+            "period": period,
+            "dc_data": dc_data,
+        }, events
 
     def extend_drive_cycle(self, drive_cycle, end_time):
         "Extends the drive cycle to enable for event"

pybamm/simulation.py

@@ -21,7 +21,7 @@ def is_notebook():
             return False  # Terminal running IPython
         elif shell == "Shell":  # pragma: no cover
             return True  # Google Colab notebook
-        else:
+        else:  # pragma: no cover
             return False  # Other type (?)
     except NameError:
         return False  # Probably standard Python interpreter
@@ -170,45 +170,76 @@ def set_up_experiment(self, model, experiment):
                 "Power switch": 0,
                 "CCCV switch": 0,
                 "Current input [A]": 0,
-                "Voltage input [V]": 0,  # doesn't matter
-                "Power input [W]": 0,  # doesn't matter
+                "Voltage input [V]": 0,
+                "Power input [W]": 0,
             }
             op_control = op["electric"][1]
-            if op_control in ["A", "C"]:
-                capacity = self._parameter_values["Nominal cell capacity [A.h]"]
+            if op["dc_data"] is not None:
+                # If operating condition includes a drive cycle, define the interpolant
+                timescale = self._parameter_values.evaluate(model.timescale)
+                drive_cycle_interpolant = pybamm.Interpolant(
+                    op["dc_data"][:, 0],
+                    op["dc_data"][:, 1],
+                    timescale * (pybamm.t - pybamm.InputParameter("start time")),
+                )
                 if op_control == "A":
-                    I = op["electric"][0]
-                    Crate = I / capacity
-                else:
-                    # Scale C-rate with capacity to obtain current
-                    Crate = op["electric"][0]
-                    I = Crate * capacity
-                if len(op["electric"]) == 4:
-                    # Update inputs for CCCV
-                    op_control = "CCCV"  # change to CCCV
-                    V = op["electric"][2]
                     operating_inputs.update(
                         {
-                            "CCCV switch": 1,
-                            "Current input [A]": I,
-                            "Voltage input [V]": V,
+                            "Current switch": 1,
+                            "Current input [A]": drive_cycle_interpolant,
                         }
                     )
-                else:
-                    # Update inputs for constant current
+                if op_control == "V":
+                    operating_inputs.update(
+                        {
+                            "Voltage switch": 1,
+                            "Voltage input [V]": drive_cycle_interpolant,
+                        }
+                    )
+                if op_control == "W":
+                    operating_inputs.update(
+                        {"Power switch": 1, "Power input [W]": drive_cycle_interpolant}
+                    )
+            else:
+                if op_control in ["A", "C"]:
+                    capacity = self._parameter_values["Nominal cell capacity [A.h]"]
+                    if op_control == "A":
+                        I = op["electric"][0]
+                        Crate = I / capacity
+                    else:
+                        # Scale C-rate with capacity to obtain current
+                        Crate = op["electric"][0]
+                        I = Crate * capacity
+                    if len(op["electric"]) == 4:
+                        # Update inputs for CCCV
+                        op_control = "CCCV"  # change to CCCV
+                        V = op["electric"][2]
+                        operating_inputs.update(
+                            {
+                                "CCCV switch": 1,
+                                "Current input [A]": I,
+                                "Voltage input [V]": V,
+                            }
+                        )
+                    else:
+                        # Update inputs for constant current
+                        operating_inputs.update(
+                            {"Current switch": 1, "Current input [A]": I}
+                        )
+                elif op_control == "V":
+                    # Update inputs for constant voltage
+                    V = op["electric"][0]
                     operating_inputs.update(
-                        {"Current switch": 1, "Current input [A]": I}
+                        {"Voltage switch": 1, "Voltage input [V]": V}
                     )
-            elif op_control == "V":
-                # Update inputs for constant voltage
-                V = op["electric"][0]
-                operating_inputs.update({"Voltage switch": 1, "Voltage input [V]": V})
-            elif op_control == "W":
-                # Update inputs for constant power
-                P = op["electric"][0]
-                operating_inputs.update({"Power switch": 1, "Power input [W]": P})
+                elif op_control == "W":
+                    # Update inputs for constant power
+                    P = op["electric"][0]
+                    operating_inputs.update({"Power switch": 1, "Power input [W]": P})
+
             # Update period
             operating_inputs["period"] = op["period"]
+
             # Update events
             if events is None:
                 # make current and voltage values that won't be hit
@@ -851,6 +882,11 @@ def solve(
                         f"step {step_num}/{cycle_length}: {op_conds_str}"
                     )
                     inputs.update(exp_inputs)
+                    if current_solution is None:
+                        start_time = 0
+                    else:
+                        start_time = current_solution.t[-1]
+                    inputs.update({"start time": start_time})
                     kwargs["inputs"] = inputs
                     # Make sure we take at least 2 timesteps
                     npts = max(int(round(dt / exp_inputs["period"])) + 1, 2)

pybamm/solvers/base_solver.py

@@ -1160,7 +1160,7 @@ def step(
         external_variables : dict
             A dictionary of external variables and their corresponding
             values at the current time
-        inputs_dict : dict, optional
+        inputs : dict, optional
             Any input parameters to pass to the model when solving
         save : bool
             Turn on to store the solution of all previous timesteps
@@ -1201,6 +1201,14 @@ def step(
         # Set up external variables and inputs
         external_variables = external_variables or {}
         inputs = inputs or {}
+
+        # Remove interpolant inputs as Casadi can't handle them
+        if isinstance(inputs.get("Current input [A]"), pybamm.Interpolant):
+            del inputs["Current input [A]"]
+        elif isinstance(inputs.get("Voltage input [V]"), pybamm.Interpolant):
+            del inputs["Voltage input [V]"]
+        elif isinstance(inputs.get("Power input [W]"), pybamm.Interpolant):
+            del inputs["Power input [W]"]
         ext_and_inputs = {**external_variables, **inputs}
 
         # Check that any inputs that may affect the scaling have not changed