Solve experiment with drive cycle

examples/scripts/experiment_drive_cycle.py

@@ -0,0 +1,92 @@
+#
+# Constant-current constant-voltage charge with US06 Drive Cycle using Experiment Class.
+#
+import pybamm
+import numpy as np
+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_ip_value, max_ip_value, min_op_value, max_op_value):
+    return (x - min_ip_value) * (max_op_value - min_op_value) / (
+        max_ip_value - min_ip_value
+    ) + min_op_value
+
+
+# Map Current to Voltage
+temp_volts = np.array([])
+min_ip_value = drive_cycle_current[:, 1].min()
+max_ip_value = drive_cycle_current[:, 1].max()
+min_Voltage = 3.5
+max_Voltage = 4.1
+for I in drive_cycle_current[:, 1]:
+    temp_volts = np.append(
+        temp_volts,
+        Map_Drive_Cycle(I, min_ip_value, max_ip_value, min_Voltage, max_Voltage),
+    )
+
+drive_cycle_voltage = drive_cycle_current
+drive_cycle_voltage = np.column_stack(
+    (np.delete(drive_cycle_voltage, 1, 1), np.array(temp_volts))
+)
+
+# Map Current to Power
+temp_volts = np.array([])
+min_ip_value = drive_cycle_current[:, 1].min()
+max_ip_value = drive_cycle_current[:, 1].max()
+min_Power = 2.5
+max_Power = 5.5
+for I in drive_cycle_current[:, 1]:
+    temp_volts = np.append(
+        temp_volts, Map_Drive_Cycle(I, min_ip_value, max_ip_value, min_Power, max_Power)
+    )
+
+drive_cycle_power = drive_cycle_current
+drive_cycle_power = np.column_stack(
+    (np.delete(drive_cycle_power, 1, 1), np.array(temp_volts))
+)
+experiment = pybamm.Experiment(
+    [
+        "Charge at 1 A until 4.1 V",
+        "Hold at 4.1 V until 50 mA",
+        "Rest for 1 hour",
+        "Run US06_A (A)",
+        "Rest for 1 hour",
+    ]
+    # + [
+    #     "Charge at 1 A until 4.1 V",
+    #     "Hold at 4.1 V until 50 mA",
+    #     "Rest for 1 hour",
+    #     "Run US06_V (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",
+        "Run US06_W (W)",
+        "Rest for 1 hour",
+    ],
+    drive_cycles={
+        "US06_A": drive_cycle_current,
+        "US06_V": drive_cycle_voltage,
+        "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,13 @@ def read_string(self, cond, drive_cycles):
                     )
                 )
 
-        return {"electric": electric, "time": time, "period": period}, events
+        return {
+            "electric": electric,
+            "time": time,
+            "period": period,
+            # "cond": cond,
+            "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

@@ -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 ndarray is recived from, create 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(
-                        {"Current switch": 1, "Current input [A]": I}
+                        {"Power switch": 1, "Power input [W]": drive_cycle_interpolant}
                     )
-            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})
+            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(
+                        {"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})
+
             # 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

@@ -1201,6 +1201,12 @@ def step(
         # Set up external variables and inputs
         external_variables = external_variables or {}
         inputs = inputs or {}
+        # if isinstance(inputs['Current input [A]'], pybamm.Interpolant):
+        #     del inputs['Current input [A]']
+        # elif isinstance(inputs['Voltage input [V]'], pybamm.Interpolant):
+        #     del inputs['Voltage input [V]']
+        # elif isinstance(inputs['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