diff --git a/CHANGELOG.md b/CHANGELOG.md index a264dfc6b2..0b93524c14 100644 --- a/CHANGELOG.md +++ b/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)) diff --git a/examples/scripts/experiment_drive_cycle.py b/examples/scripts/experiment_drive_cycle.py new file mode 100644 index 0000000000..2ead334dad --- /dev/null +++ b/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() diff --git a/pybamm/experiments/experiment.py b/pybamm/experiments/experiment.py index ab2f3d709d..c175a94815 100644 --- a/pybamm/experiments/experiment.py +++ b/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" diff --git a/pybamm/simulation.py b/pybamm/simulation.py index 58f435868e..589b762cd9 100644 --- a/pybamm/simulation.py +++ b/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) diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py index ca2b43dc1a..98fa02dd59 100644 --- a/pybamm/solvers/base_solver.py +++ b/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 diff --git a/tests/unit/test_experiments/test_experiment.py b/tests/unit/test_experiments/test_experiment.py index db6f6196e2..b91645460e 100644 --- a/tests/unit/test_experiments/test_experiment.py +++ b/tests/unit/test_experiments/test_experiment.py @@ -46,19 +46,44 @@ def test_read_strings(self): self.assertEqual( experiment.operating_conditions[:-3], [ - {"electric": (1, "C"), "time": 1800.0, "period": 20.0}, - {"electric": (0.05, "C"), "time": 1800.0, "period": 20.0}, - {"electric": (-0.5, "C"), "time": 2700.0, "period": 20.0}, - {"electric": (1, "A"), "time": 1800.0, "period": 20.0}, - {"electric": (-0.2, "A"), "time": 2700.0, "period": 60.0}, - {"electric": (1, "W"), "time": 1800.0, "period": 20.0}, - {"electric": (-0.2, "W"), "time": 2700.0, "period": 20.0}, - {"electric": (0, "A"), "time": 600.0, "period": 300.0}, - {"electric": (1, "V"), "time": 20.0, "period": 20.0}, - {"electric": (-1, "C"), "time": None, "period": 20.0}, - {"electric": (4.1, "V"), "time": None, "period": 20.0}, - {"electric": (3, "V"), "time": None, "period": 20.0}, - {"electric": (1 / 3, "C"), "time": 7200.0, "period": 20.0}, + {"electric": (1, "C"), "time": 1800.0, "period": 20.0, "dc_data": None}, + { + "electric": (0.05, "C"), + "time": 1800.0, + "period": 20.0, + "dc_data": None, + }, + { + "electric": (-0.5, "C"), + "time": 2700.0, + "period": 20.0, + "dc_data": None, + }, + {"electric": (1, "A"), "time": 1800.0, "period": 20.0, "dc_data": None}, + { + "electric": (-0.2, "A"), + "time": 2700.0, + "period": 60.0, + "dc_data": None, + }, + {"electric": (1, "W"), "time": 1800.0, "period": 20.0, "dc_data": None}, + { + "electric": (-0.2, "W"), + "time": 2700.0, + "period": 20.0, + "dc_data": None, + }, + {"electric": (0, "A"), "time": 600.0, "period": 300.0, "dc_data": None}, + {"electric": (1, "V"), "time": 20.0, "period": 20.0, "dc_data": None}, + {"electric": (-1, "C"), "time": None, "period": 20.0, "dc_data": None}, + {"electric": (4.1, "V"), "time": None, "period": 20.0, "dc_data": None}, + {"electric": (3, "V"), "time": None, "period": 20.0, "dc_data": None}, + { + "electric": (1 / 3, "C"), + "time": 7200.0, + "period": 20.0, + "dc_data": None, + }, ], ) # Calculation for operating conditions of drive cycle @@ -72,19 +97,19 @@ def test_read_strings(self): period_2 = np.min(np.diff(drive_cycle_2[:, 0])) # Check drive cycle operating conditions np.testing.assert_array_equal( - experiment.operating_conditions[-3]["electric"][0], drive_cycle + experiment.operating_conditions[-3]["dc_data"], drive_cycle ) self.assertEqual(experiment.operating_conditions[-3]["electric"][1], "A") self.assertEqual(experiment.operating_conditions[-3]["time"], time_0) self.assertEqual(experiment.operating_conditions[-3]["period"], period_0) np.testing.assert_array_equal( - experiment.operating_conditions[-2]["electric"][0], drive_cycle_1 + experiment.operating_conditions[-2]["dc_data"], drive_cycle_1 ) self.assertEqual(experiment.operating_conditions[-2]["electric"][1], "V") self.assertEqual(experiment.operating_conditions[-2]["time"], time_1) self.assertEqual(experiment.operating_conditions[-2]["period"], period_1) np.testing.assert_array_equal( - experiment.operating_conditions[-1]["electric"][0], drive_cycle_2 + experiment.operating_conditions[-1]["dc_data"], drive_cycle_2 ) self.assertEqual(experiment.operating_conditions[-1]["electric"][1], "W") self.assertEqual(experiment.operating_conditions[-1]["time"], time_2) @@ -128,9 +153,24 @@ def test_read_strings_cccv_combined(self): self.assertEqual( experiment.operating_conditions, [ - {"electric": (0.05, "C"), "time": 1800.0, "period": 60.0}, - {"electric": (-0.5, "C", 1, "V"), "time": None, "period": 60.0}, - {"electric": (0.05, "C"), "time": 1800.0, "period": 60.0}, + { + "electric": (0.05, "C"), + "time": 1800.0, + "period": 60.0, + "dc_data": None, + }, + { + "electric": (-0.5, "C", 1, "V"), + "time": None, + "period": 60.0, + "dc_data": None, + }, + { + "electric": (0.05, "C"), + "time": 1800.0, + "period": 60.0, + "dc_data": None, + }, ], ) self.assertEqual(experiment.events, [None, (0.02, "C"), None]) @@ -146,8 +186,13 @@ def test_read_strings_cccv_combined(self): self.assertEqual( experiment.operating_conditions, [ - {"electric": (-0.5, "C"), "time": None, "period": 60.0}, - {"electric": (1, "V"), "time": None, "period": 60.0}, + { + "electric": (-0.5, "C"), + "time": None, + "period": 60.0, + "dc_data": None, + }, + {"electric": (1, "V"), "time": None, "period": 60.0, "dc_data": None}, ], ) experiment = pybamm.Experiment( @@ -160,8 +205,13 @@ def test_read_strings_cccv_combined(self): self.assertEqual( experiment.operating_conditions, [ - {"electric": (-0.5, "C"), "time": 120.0, "period": 60.0}, - {"electric": (1, "V"), "time": None, "period": 60.0}, + { + "electric": (-0.5, "C"), + "time": 120.0, + "period": 60.0, + "dc_data": None, + }, + {"electric": (1, "V"), "time": None, "period": 60.0, "dc_data": None}, ], ) @@ -173,11 +223,26 @@ def test_read_strings_repeat(self): self.assertEqual( experiment.operating_conditions, [ - {"electric": (0.01, "A"), "time": 1800.0, "period": 60}, - {"electric": (-0.5, "C"), "time": 2700.0, "period": 60}, - {"electric": (1, "V"), "time": 20.0, "period": 60}, - {"electric": (-0.5, "C"), "time": 2700.0, "period": 60}, - {"electric": (1, "V"), "time": 20.0, "period": 60}, + { + "electric": (0.01, "A"), + "time": 1800.0, + "period": 60, + "dc_data": None, + }, + { + "electric": (-0.5, "C"), + "time": 2700.0, + "period": 60, + "dc_data": None, + }, + {"electric": (1, "V"), "time": 20.0, "period": 60, "dc_data": None}, + { + "electric": (-0.5, "C"), + "time": 2700.0, + "period": 60, + "dc_data": None, + }, + {"electric": (1, "V"), "time": 20.0, "period": 60, "dc_data": None}, ], ) self.assertEqual(experiment.period, 60) @@ -193,10 +258,30 @@ def test_cycle_unpacking(self): self.assertEqual( experiment.operating_conditions, [ - {"electric": (0.05, "C"), "time": 1800.0, "period": 60.0}, - {"electric": (-0.2, "C"), "time": 2700.0, "period": 60.0}, - {"electric": (0.05, "C"), "time": 1800.0, "period": 60.0}, - {"electric": (-0.2, "C"), "time": 2700.0, "period": 60.0}, + { + "electric": (0.05, "C"), + "time": 1800.0, + "period": 60.0, + "dc_data": None, + }, + { + "electric": (-0.2, "C"), + "time": 2700.0, + "period": 60.0, + "dc_data": None, + }, + { + "electric": (0.05, "C"), + "time": 1800.0, + "period": 60.0, + "dc_data": None, + }, + { + "electric": (-0.2, "C"), + "time": 2700.0, + "period": 60.0, + "dc_data": None, + }, ], ) self.assertEqual(experiment.cycle_lengths, [2, 1, 1]) diff --git a/tests/unit/test_experiments/test_simulation_with_experiment.py b/tests/unit/test_experiments/test_simulation_with_experiment.py index 31eb795df9..5b47b672b8 100644 --- a/tests/unit/test_experiments/test_simulation_with_experiment.py +++ b/tests/unit/test_experiments/test_simulation_with_experiment.py @@ -158,6 +158,24 @@ def test_run_experiment_cccv_ode(self): ) self.assertEqual(solutions[1].termination, "final time") + def test_run_experiment_drive_cycle(self): + drive_cycle = np.array([np.arange(10), np.arange(10)]).T + experiment = pybamm.Experiment( + [ + ( + "Run drive_cycle (A)", + "Run drive_cycle (V)", + "Run drive_cycle (W)", + ) + ], + drive_cycles={"drive_cycle": drive_cycle} + ) + model = pybamm.lithium_ion.DFN() + sim = pybamm.Simulation(model, experiment=experiment) + self.assertIn(('drive_cycle', 'A'), sim.op_conds_to_model_and_param) + self.assertIn(('drive_cycle', 'V'), sim.op_conds_to_model_and_param) + self.assertIn(('drive_cycle', 'W'), sim.op_conds_to_model_and_param) + def test_run_experiment_old_setup_type(self): experiment = pybamm.Experiment( [ diff --git a/tests/unit/test_solvers/test_base_solver.py b/tests/unit/test_solvers/test_base_solver.py index 7967c96949..180946765a 100644 --- a/tests/unit/test_solvers/test_base_solver.py +++ b/tests/unit/test_solvers/test_base_solver.py @@ -296,6 +296,21 @@ def test_timescale_input_fail(self): with self.assertRaisesRegex(pybamm.SolverError, "The model timescale"): sol = solver.step(old_solution=sol, model=model, dt=1.0, inputs={"a": 20}) + def test_inputs_step(self): + # Make sure interpolant inputs are dropped + model = pybamm.BaseModel() + v = pybamm.Variable("v") + model.rhs = {v: -1} + model.initial_conditions = {v: 1} + x = np.array([0, 1]) + interp = pybamm.Interpolant(x, x, pybamm.t) + solver = pybamm.CasadiSolver() + for input_key in ["Current input [A]", "Voltage input [V]", "Power input [W]"]: + sol = solver.step( + old_solution=None, model=model, dt=1.0, inputs={input_key: interp} + ) + self.assertFalse(input_key in sol.all_inputs[0]) + def test_extrapolation_warnings(self): # Make sure the extrapolation warnings work model = pybamm.BaseModel() @@ -327,29 +342,25 @@ def test_extrapolation_warnings(self): @unittest.skipIf(not pybamm.have_idaklu(), "idaklu solver is not installed") def test_sensitivities(self): - def exact_diff_a(y, a, b): - return np.array([ - [y[0]**2 + 2 * a], - [y[0]] - ]) + return np.array([[y[0] ** 2 + 2 * a], [y[0]]]) def exact_diff_b(y, a, b): return np.array([[y[0]], [0]]) - for convert_to_format in ['', 'python', 'casadi', 'jax']: + for convert_to_format in ["", "python", "casadi", "jax"]: model = pybamm.BaseModel() v = pybamm.Variable("v") u = pybamm.Variable("u") a = pybamm.InputParameter("a") b = pybamm.InputParameter("b") - model.rhs = {v: a * v**2 + b * v + a**2} + model.rhs = {v: a * v ** 2 + b * v + a ** 2} model.algebraic = {u: a * v - u} model.initial_conditions = {v: 1, u: a * 1} model.convert_to_format = convert_to_format - solver = pybamm.IDAKLUSolver(root_method='lm') - model.calculate_sensitivities = ['a', 'b'] - solver.set_up(model, inputs={'a': 0, 'b': 0}) + solver = pybamm.IDAKLUSolver(root_method="lm") + model.calculate_sensitivities = ["a", "b"] + solver.set_up(model, inputs={"a": 0, "b": 0}) all_inputs = [] for v_value in [0.1, -0.2, 1.5, 8.4]: for u_value in [0.13, -0.23, 1.3, 13.4]: @@ -357,24 +368,20 @@ def exact_diff_b(y, a, b): for b_value in [0.82, 1.9]: y = np.array([v_value, u_value]) t = 0 - inputs = {'a': a_value, 'b': b_value} + inputs = {"a": a_value, "b": b_value} all_inputs.append((t, y, inputs)) for t, y, inputs in all_inputs: - if model.convert_to_format == 'casadi': + if model.convert_to_format == "casadi": use_inputs = casadi.vertcat(*[x for x in inputs.values()]) else: use_inputs = inputs - sens = model.sensitivities_eval( - t, y, use_inputs - ) + sens = model.sensitivities_eval(t, y, use_inputs) np.testing.assert_allclose( - sens['a'], - exact_diff_a(y, inputs['a'], inputs['b']) + sens["a"], exact_diff_a(y, inputs["a"], inputs["b"]) ) np.testing.assert_allclose( - sens['b'], - exact_diff_b(y, inputs['a'], inputs['b']) + sens["b"], exact_diff_b(y, inputs["a"], inputs["b"]) )