Merge remote-tracking branch 'alibh95/develop' into issue-1788-solve-…

…experiment-drive-cycle

examples/scripts/experiment_drive_cycle.py

@@ -0,0 +1,87 @@
+#
+# 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_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

@@ -273,7 +273,12 @@ def read_string(self, cond, drive_cycles):
                     )
                 )
 
-        return {"electric": electric, "time": time, "period": period}, events
+        return {
+            "electric": electric,
+            "time": time,
+            "period": period,
+            "cond": cond,
+        }, events
 
     def extend_drive_cycle(self, drive_cycle, end_time):
         "Extends the drive cycle to enable for event"

pybamm/simulation.py

@@ -174,41 +174,62 @@ def set_up_experiment(self, model, experiment):
                 "Power input [W]": 0,  # doesn't matter
             }
             op_control = op["electric"][1]
-            if op_control in ["A", "C"]:
-                capacity = self._parameter_values["Nominal cell capacity [A.h]"]
+            if isinstance(op[0], np.ndarray):
+                # If ndarray is recived from, create interpolant
+                timescale = self._parameter_values.evaluate(model.timescale)
+                drive_cycle_interpolant = pybamm.Interpolant(
+                    op[0][:, 0], op[0][:, 1], timescale * pybamm.t
+                )
                 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[1] == "V":
+                    operating_inputs.update(
+                        {"Voltage switch": 1, "Voltage input [V]": drive_cycle_interpolant}
+                    )
+                if op[1] == "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
@@ -1177,12 +1198,10 @@ def save(self, filename):
             and self._solver.integrator_specs != {}
         ):
             self._solver.integrator_specs = {}
-        if self.solution is not None:
-            self.solution.clear_casadi_attributes()
         with open(filename, "wb") as f:
             pickle.dump(self, f, pickle.HIGHEST_PROTOCOL)
 
 
 def load_sim(filename):
     """Load a saved simulation"""
-    return pybamm.load(filename)
+    return pybamm.load(filename)

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