Merge pull request #1612 from pybamm-team/issue-1611-CCCV

Issue 1611 cccv

docs/source/models/submodels/external_circuit/function_control_external_circuit.rst

@@ -8,4 +8,7 @@ Function control external circuit
     :members:
 
 .. autoclass:: pybamm.external_circuit.PowerFunctionControl
+    :members:
+
+.. autoclass:: pybamm.external_circuit.CCCVFunctionControl
     :members:

pybamm/experiments/experiment.py

@@ -55,7 +55,11 @@ class Experiment:
         faster at simulating individual steps but has higher set-up overhead
     drive_cycles : dict
         Dictionary of drive cycles to use for this experiment.
-
+    cccv_handling : str, optional
+        How to handle CCCV. If "two-step" (default), then the experiment is run in
+        two steps (CC then CV). If "ode", then the experiment is run in a single step
+        using an ODE for current: see
+        :class:`pybamm.external_circuit.CCCVFunctionControl` for details.
     """
 
     def __init__(
@@ -66,7 +70,11 @@ def __init__(
         termination=None,
         use_simulation_setup_type="new",
         drive_cycles={},
+        cccv_handling="two-step",
     ):
+        if cccv_handling not in ["two-step", "ode"]:
+            raise ValueError("cccv_handling should be either 'two-step' or 'ode'")
+        self.cccv_handling = cccv_handling
 
         self.period = self.convert_time_to_seconds(period.split())
         operating_conditions_cycles = []
@@ -75,9 +83,28 @@ def __init__(
             if (isinstance(cycle, tuple) or isinstance(cycle, str)) and all(
                 [isinstance(cond, str) for cond in cycle]
             ):
-                operating_conditions_cycles.append(
-                    cycle if isinstance(cycle, tuple) else (cycle,)
-                )
+                if isinstance(cycle, str):
+                    processed_cycle = (cycle,)
+                else:
+                    processed_cycle = []
+                    idx = 0
+                    finished = False
+                    while not finished:
+                        step = cycle[idx]
+                        if idx < len(cycle) - 1:
+                            next_step = cycle[idx + 1]
+                        else:
+                            next_step = None
+                            finished = True
+                        if self.is_cccv(step, next_step):
+                            processed_cycle.append(step + " then " + next_step)
+                            idx += 2
+                        else:
+                            processed_cycle.append(step)
+                            idx += 1
+                        if idx >= len(cycle):
+                            finished = True
+                operating_conditions_cycles.append(tuple(processed_cycle))
             else:
                 try:
                     # Condition is not a string
@@ -153,7 +180,20 @@ def read_string(self, cond, drive_cycles):
             must be numbers, 'C' denotes the unit of the external circuit (can be A for
             current, C for C-rate, V for voltage or W for power), and 'hours' denotes
             the unit of time (can be second(s), minute(s) or hour(s))
+        drive_cycles: dict
+            A map specifying the drive cycles
         """
+        if " then " in cond:
+            # If the string contains " then ", then this is a two-step CCCV experiment
+            # and we need to split it into two strings
+            cond_CC, cond_CV = cond.split(" then ")
+            op_CC, _ = self.read_string(cond_CC, drive_cycles)
+            op_CV, event_CV = self.read_string(cond_CV, drive_cycles)
+            return {
+                "electric": op_CC["electric"] + op_CV["electric"],
+                "time": op_CV["time"],
+                "period": op_CV["period"],
+            }, event_CV
         # Read period
         if " period)" in cond:
             cond, time_period = cond.split("(")
@@ -165,20 +205,12 @@ def read_string(self, cond, drive_cycles):
         if "Run" in cond:
             cond_list = cond.split()
             if "at" in cond:
-                raise ValueError(
-                    """Instruction must be
-                    For example: {}""".format(
-                        examples
-                    )
-                )
+                raise ValueError(f"Instruction must be of the form: {examples}")
             dc_types = ["(A)", "(V)", "(W)"]
             if all(x not in cond for x in dc_types):
                 raise ValueError(
-                    """Type of drive cycle must be
-                    specified using '(A)', '(V)' or '(W)'.
-                    For example: {}""".format(
-                        examples
-                    )
+                    "Type of drive cycle must be specified using '(A)', '(V)' or '(W)'."
+                    f" For example: {examples}"
                 )
             # Check for Events
             elif "for" in cond:
@@ -208,7 +240,7 @@ def read_string(self, cond, drive_cycles):
                 time = drive_cycles[cond_list[1]][:, 0][-1]
                 period = np.min(np.diff(drive_cycles[cond_list[1]][:, 0]))
                 events = None
-        elif "Run" not in cond:
+        else:
             if "for" in cond and "or until" in cond:
                 # e.g. for 3 hours or until 4.2 V
                 cond_list = cond.split()
@@ -238,7 +270,8 @@ def read_string(self, cond, drive_cycles):
                         examples
                     )
                 )
-        return electric + (time,) + (period,), events
+
+        return {"electric": electric, "time": time, "period": period}, events
 
     def extend_drive_cycle(self, drive_cycle, end_time):
         "Extends the drive cycle to enable for event"
@@ -404,3 +437,25 @@ def read_termination(self, termination):
                     "e.g. '80% capacity' or '4 Ah capacity'"
                 )
         return termination_dict
+
+    def is_cccv(self, step, next_step):
+        """
+        Check whether a step and the next step indicate a CCCV charge
+        """
+        if self.cccv_handling == "two-step" or next_step is None:
+            return False
+        # e.g. step="Charge at 2.0 A until 4.2V"
+        # next_step="Hold at 4.2V until C/50"
+        if (
+            step.startswith("Charge")
+            and "until" in step
+            and "V" in step
+            and "Hold at " in next_step
+            and "V until" in next_step
+        ):
+            _, events = self.read_string(step, None)
+            next_op, _ = self.read_string(next_step, None)
+            # Check that the event conditions are the same as the hold conditions
+            if events == next_op["electric"]:
+                return True
+        return False

pybamm/models/full_battery_models/base_battery_model.py

@@ -164,7 +164,7 @@ def __init__(self, extra_options):
             "lithium plating": ["none", "reversible", "irreversible"],
             "lithium plating porosity change": ["true", "false"],
             "loss of active material": ["none", "stress-driven", "reaction-driven"],
-            "operating mode": ["current", "voltage", "power"],
+            "operating mode": ["current", "voltage", "power", "CCCV"],
             "particle": [
                 "Fickian diffusion",
                 "fast diffusion",
@@ -410,7 +410,7 @@ def default_parameter_values(self):
     def default_geometry(self):
         return pybamm.battery_geometry(
             options=self.options,
-            current_collector_dimension=self.options["dimensionality"]
+            current_collector_dimension=self.options["dimensionality"],
         )
 
     @property
@@ -440,12 +440,8 @@ def default_submesh_types(self):
             "positive electrode": pybamm.MeshGenerator(pybamm.Uniform1DSubMesh),
             "negative particle": pybamm.MeshGenerator(pybamm.Uniform1DSubMesh),
             "positive particle": pybamm.MeshGenerator(pybamm.Uniform1DSubMesh),
-            "negative particle size": pybamm.MeshGenerator(
-                pybamm.Uniform1DSubMesh
-            ),
-            "positive particle size": pybamm.MeshGenerator(
-                pybamm.Uniform1DSubMesh
-            ),
+            "negative particle size": pybamm.MeshGenerator(pybamm.Uniform1DSubMesh),
+            "positive particle size": pybamm.MeshGenerator(pybamm.Uniform1DSubMesh),
         }
         if self.options["dimensionality"] == 0:
             base_submeshes["current collector"] = pybamm.MeshGenerator(pybamm.SubMesh0D)
@@ -731,7 +727,7 @@ def build_model(self):
         pybamm.logger.info("Finish building {}".format(self.name))
 
     def new_empty_copy(self):
-        """ See :meth:`pybamm.BaseModel.new_empty_copy()` """
+        """See :meth:`pybamm.BaseModel.new_empty_copy()`"""
         new_model = self.__class__(name=self.name, options=self.options)
         new_model.use_jacobian = self.use_jacobian
         new_model.convert_to_format = self.convert_to_format
@@ -756,6 +752,10 @@ def set_external_circuit_submodel(self):
             self.submodels[
                 "external circuit"
             ] = pybamm.external_circuit.PowerFunctionControl(self.param)
+        elif self.options["operating mode"] == "CCCV":
+            self.submodels[
+                "external circuit"
+            ] = pybamm.external_circuit.CCCVFunctionControl(self.param)
         elif callable(self.options["operating mode"]):
             self.submodels[
                 "external circuit"

pybamm/models/submodels/external_circuit/__init__.py

@@ -4,6 +4,7 @@
     FunctionControl,
     VoltageFunctionControl,
     PowerFunctionControl,
+    CCCVFunctionControl,
     LeadingOrderFunctionControl,
     LeadingOrderVoltageFunctionControl,
     LeadingOrderPowerFunctionControl,

pybamm/models/submodels/external_circuit/current_control_external_circuit.py

@@ -1,6 +1,7 @@
 #
 # External circuit with current control
 #
+import pybamm
 from .base_external_circuit import BaseModel, LeadingOrderBaseModel
 
 
@@ -17,6 +18,7 @@ def get_fundamental_variables(self):
         I = self.param.dimensional_current_with_time
 
         variables = {
+            "Current density variable": pybamm.Scalar(1, name="i_cell"),
             "Total current density": i_cell,
             "Total current density [A.m-2]": i_cell_dim,
             "Current [A]": I,

pybamm/models/submodels/external_circuit/function_control_external_circuit.py

@@ -6,22 +6,41 @@
 
 
 class FunctionControl(BaseModel):
-    """External circuit with an arbitrary function."""
+    """
+    External circuit with an arbitrary function, implemented as a control on the current
+    either via an algebraic equation, or a differential equation.
+
+    Parameters
+    ----------
+    param : parameter class
+        The parameters to use for this submodel
+    external_circuit_function : callable
+        The function that controls the current
+    control : str, optional
+        The type of control to use. Must be one of 'algebraic' (default)
+        or 'differential'.
+    """
 
-    def __init__(self, param, external_circuit_function):
+    def __init__(self, param, external_circuit_function, control="algebraic"):
         super().__init__(param)
         self.external_circuit_function = external_circuit_function
+        self.control = control
 
     def get_fundamental_variables(self):
         param = self.param
         # Current is a variable
-        i_cell = pybamm.Variable("Total current density")
+        i_var = pybamm.Variable("Current density variable")
+        if self.control == "algebraic":
+            i_cell = i_var
+        elif self.control == "differential":
+            i_cell = pybamm.maximum(i_var, param.current_with_time)
 
         # Update derived variables
-        I = i_cell * pybamm.AbsoluteValue(param.I_typ)
+        I = i_cell * abs(param.I_typ)
         i_cell_dim = I / (param.n_electrodes_parallel * param.A_cc)
 
         variables = {
+            "Current density variable": i_var,
             "Total current density": i_cell,
             "Total current density [A.m-2]": i_cell_dim,
             "Current [A]": I,
@@ -36,15 +55,25 @@ def get_fundamental_variables(self):
     def set_initial_conditions(self, variables):
         super().set_initial_conditions(variables)
         # Initial condition as a guess for consistent initial conditions
-        i_cell = variables["Total current density"]
+        i_cell = variables["Current density variable"]
         self.initial_conditions[i_cell] = self.param.current_with_time
 
+    def set_rhs(self, variables):
+        super().set_rhs(variables)
+        # External circuit submodels are always equations on the current
+        # The external circuit function should provide an update law for the current
+        # based on current/voltage/power/etc.
+        if self.control == "differential":
+            i_cell = variables["Current density variable"]
+            self.rhs[i_cell] = self.external_circuit_function(variables)
+
     def set_algebraic(self, variables):
         # External circuit submodels are always equations on the current
         # The external circuit function should fix either the current, or the voltage,
         # or a combination (e.g. I*V for power control)
-        i_cell = variables["Total current density"]
-        self.algebraic[i_cell] = self.external_circuit_function(variables)
+        if self.control == "algebraic":
+            i_cell = variables["Current density variable"]
+            self.algebraic[i_cell] = self.external_circuit_function(variables)
 
 
 class VoltageFunctionControl(FunctionControl):
@@ -53,7 +82,7 @@ class VoltageFunctionControl(FunctionControl):
     """
 
     def __init__(self, param):
-        super().__init__(param, self.constant_voltage)
+        super().__init__(param, self.constant_voltage, control="algebraic")
 
     def constant_voltage(self, variables):
         V = variables["Terminal voltage [V]"]
@@ -66,7 +95,7 @@ class PowerFunctionControl(FunctionControl):
     """External circuit with power control."""
 
     def __init__(self, param):
-        super().__init__(param, self.constant_power)
+        super().__init__(param, self.constant_power, control="algebraic")
 
     def constant_power(self, variables):
         I = variables["Current [A]"]
@@ -76,11 +105,29 @@ def constant_power(self, variables):
         )
 
 
+class CCCVFunctionControl(FunctionControl):
+    """External circuit with constant-current constant-voltage control."""
+
+    def __init__(self, param):
+        super().__init__(param, self.cccv, control="differential")
+
+    def cccv(self, variables):
+        # Multiply by the time scale so that the votage overshoot only lasts a few
+        # seconds
+        K_aw = 1 * self.param.timescale  # anti-windup
+        K_V = 1 * self.param.timescale
+        i_var = variables["Current density variable"]
+        i_cell = variables["Total current density"]
+        V = variables["Terminal voltage [V]"]
+        V_CCCV = pybamm.Parameter("Voltage function [V]")
+        return -K_aw * (i_var - i_cell) + K_V * (V - V_CCCV)
+
+
 class LeadingOrderFunctionControl(FunctionControl, LeadingOrderBaseModel):
     """External circuit with an arbitrary function, at leading order."""
 
-    def __init__(self, param, external_circuit_class):
-        super().__init__(param, external_circuit_class)
+    def __init__(self, param, external_circuit_class, control="algebraic"):
+        super().__init__(param, external_circuit_class, control=control)
 
     def _get_current_variable(self):
         return pybamm.Variable("Leading-order total current density")
@@ -93,7 +140,7 @@ class LeadingOrderVoltageFunctionControl(LeadingOrderFunctionControl):
     """
 
     def __init__(self, param):
-        super().__init__(param, self.constant_voltage)
+        super().__init__(param, self.constant_voltage, control="algebraic")
 
     def constant_voltage(self, variables):
         V = variables["Terminal voltage [V]"]
@@ -106,7 +153,7 @@ class LeadingOrderPowerFunctionControl(LeadingOrderFunctionControl):
     """External circuit with power control, at leading order."""
 
     def __init__(self, param):
-        super().__init__(param, self.constant_power)
+        super().__init__(param, self.constant_power, control="algebraic")
 
     def constant_power(self, variables):
         I = variables["Current [A]"]