Repository: pybamm-team/PyBaMM

Issue:
Different constant-power approach
Instead of an extra algebraic equation, constant power could be achieved by specifying the current to be `P_applied / voltage`, where `voltage` is the boundary value of the potential.
This only works for the DFN, where potential is a variable. For SPM(e), potential is defined in terms of the current, so it would be circular, and we still need to include the algebraic equation.

Discussion:

Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 429af9cd6..a932e1c68 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -3,6 +3,7 @@
 ## Features
 
 -   Isothermal models now compute heat source terms (but the temperature remains constant). The models now also account for current collector heating when `dimensionality=0` [#1929](https://github.com/pybamm-team/PyBaMM/pull/1929))
+-   Added new models for power control and resistance control ([#1917](https://github.com/pybamm-team/PyBaMM/pull/1917))
 -   Initial concentrations can now be provided as a function of `r` as well as `x` ([#1866](https://github.com/pybamm-team/PyBaMM/pull/1866))
 
 ## Bug fixes
diff --git a/docs/source/models/submodels/external_circuit/current_control_external_circuit.rst b/docs/source/models/submodels/external_circuit/current_control_external_circuit.rst
deleted file mode 100644
index 6f3f0a001..000000000
--- a/docs/source/models/submodels/external_circuit/current_control_external_circuit.rst
+++ /dev/null
@@ -1,6 +0,0 @@
-Current control external circuit
-================================
-
-.. autoclass:: pybamm.external_circuit.CurrentControl
-    :members:
-
diff --git a/docs/source/models/submodels/external_circuit/explicit_control_external_circuit.rst b/docs/source/models/submodels/external_circuit/explicit_control_external_circuit.rst
new file mode 100644
index 000000000..b24ff1d51
--- /dev/null
+++ b/docs/source/models/submodels/external_circuit/explicit_control_external_circuit.rst
@@ -0,0 +1,13 @@
+Explicit control external circuit
+=================================
+
+Current is explicitly specified, either by a function or in terms of other variables.
+
+.. autoclass:: pybamm.external_circuit.ExplicitCurrentControl
+    :members:
+
+.. autoclass:: pybamm.external_circuit.ExplicitPowerControl
+    :members:
+
+.. autoclass:: pybamm.external_circuit.ExplicitResistanceControl
+    :members:
diff --git a/docs/source/models/submodels/external_circuit/function_control_external_circuit.rst b/docs/source/models/submodels/external_circuit/function_control_external_circuit.rst
index e1876dfc5..e97760066 100644
--- a/docs/source/models/submodels/external_circuit/function_control_external_circuit.rst
+++ b/docs/source/models/submodels/external_circuit/function_control_external_circuit.rst
@@ -10,5 +10,8 @@ Function control external circuit
 .. autoclass:: pybamm.external_circuit.PowerFunctionControl
     :members:
 
+.. autoclass:: pybamm.external_circuit.ResistanceFunctionControl
+    :members:
+
 .. autoclass:: pybamm.external_circuit.CCCVFunctionControl
     :members:
\ No newline at end of file
diff --git a/docs/source/models/submodels/external_circuit/index.rst b/docs/source/models/submodels/external_circuit/index.rst
index 600c7718b..21790e466 100644
--- a/docs/source/models/submodels/external_circuit/index.rst
+++ b/docs/source/models/submodels/external_circuit/index.rst
@@ -11,5 +11,5 @@ variable to be constant.
 .. toctree::
   :maxdepth: 1
 
-  current_control_external_circuit
+  explicit_control_external_circuit
   function_control_external_circuit
diff --git a/examples/notebooks/models/using-submodels.ipynb b/examples/notebooks/models/using-submodels.ipynb
index 4deb7bb72..655801716 100644
--- a/examples/notebooks/models/using-submodels.ipynb
+++ b/examples/notebooks/models/using-submodels.ipynb
@@ -335,7 +335,7 @@
    "source": [
     "Submodels can be added to the `model.submodels` dictionary in the same way that we changed the submodels earlier. \n",
     "\n",
-    "We use the simplest model for the external circuit, which is the \"current control\" submodel"
+    "We use the simplest model for the external circuit, which is the explicit \"current control\" submodel"
    ]
   },
   {
@@ -344,7 +344,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "model.submodels[\"external circuit\"] = pybamm.external_circuit.CurrentControl(model.param)"
+    "model.submodels[\"external circuit\"] = pybamm.external_circuit.ExplicitCurrentControl(model.param)"
    ]
   },
   {
diff --git a/examples/scripts/custom_model.py b/examples/scripts/custom_model.py
index f920ddd31..6121b3fb5 100644
--- a/examples/scripts/custom_model.py
+++ b/examples/scripts/custom_model.py
@@ -11,7 +11,7 @@ pybamm.set_logging_level("INFO")
 model = pybamm.lithium_ion.BaseModel(name="my li-ion model")
 
 # set choice of submodels
-model.submodels["external circuit"] = pybamm.external_circuit.CurrentControl(
+model.submodels["external circuit"] = pybamm.external_circuit.ExplicitCurrentControl(
     model.param
 )
 model.submodels["current collector"] = pybamm.current_collector.Uniform(model.param)
diff --git a/pybamm/models/full_battery_models/base_battery_model.py b/pybamm/models/full_battery_models/base_battery_model.py
index 95e1c1291..80e2aa0b4 100644
--- a/pybamm/models/full_battery_models/base_battery_model.py
+++ b/pybamm/models/full_battery_models/base_battery_model.py
@@ -65,12 +65,21 @@ class BatteryModelOptions(pybamm.FuzzyDict):
                 A 2-tuple can be provided for different behaviour in negative and
                 positive electrodes.
             * "operating mode" : str
-                Sets the operating mode for the model. Can be "current" (default),
-                "voltage" or "power". Alternatively, the operating mode can be
-                controlled with an arbitrary function by passing the function directly
-                as the option. In this case the function must define the residual of
-                an algebraic equation. The applied current will be solved for such
-                that the algebraic constraint is satisfied.
+                Sets the operating mode for the model. This determines how the current
+                is set. Can be:
+
+                - "current" (default) : the current is explicity supplied
+                - "voltage"/"power"/"resistance" : solve an algebraic equation for \
+                    current such that voltage/power/resistance is correct
+                - "differential power"/"differential resistance" : solve a \
+                    differential equation for the power or resistance
+                - "explicit power"/"explicit resistance" : current is defined in terms \
+                    of the voltage such that power/resistance is correct
+                - "CCCV": a special implementation of the common constant-current \
+                    constant-voltage charging protocol, via an ODE for the current
+                - callable : if a callable is given as this option, the function \
+                    defines the residual of an algebraic equation. The applied current \
+                    will be solved for such that the algebraic constraint is satisfied.
             * "particle" : str
                 Sets the submodel to use to describe behaviour within the particle.
                 Can be "Fickian diffusion" (default), "uniform profile",
@@ -186,7 +195,17 @@ class BatteryModelOptions(pybamm.FuzzyDict):
                 "reaction-driven",
                 "stress and reaction-driven",
             ],
-            "operating mode": ["current", "voltage", "power", "CCCV"],
+            "operating mode": [
+                "current",
+                "voltage",
+                "power",
+                "differential power",
+                "explicit power",
+                "resistance",
+                "differential resistance",
+                "explicit resistance",
+                "CCCV",
+            ],
             "particle": [
                 "Fickian diffusion",
                 "fast diffusion",
@@ -897,27 +916,36 @@ class BaseBatteryModel(pybamm.BaseModel):
         e.g. (not necessarily constant-) current, voltage, etc
         """
         if self.options["operating mode"] == "current":
-            self.submodels["external circuit"] = pybamm.external_circuit.CurrentControl(
-                self.param
-            )
+            model = pybamm.external_circuit.ExplicitCurrentControl(self.param)
         elif self.options["operating mode"] == "voltage":
-            self.submodels[
-                "external circuit"
-            ] = pybamm.external_circuit.VoltageFunctionControl(self.param)
+            model = pybamm.external_circuit.VoltageFunctionControl(self.param)
         elif self.options["operating mode"] == "power":
-            self.submodels[
-                "external circuit"
-            ] = pybamm.external_circuit.PowerFunctionControl(self.param)
+            model = pybamm.external_circuit.PowerFunctionControl(
+                self.param, "algebraic"
+            )
+        elif self.options["operating mode"] == "differential power":
+            model = pybamm.external_circuit.PowerFunctionControl(
+                self.param, "differential without max"
+            )
+        elif self.options["operating mode"] == "explicit power":
+            model = pybamm.external_circuit.ExplicitPowerControl(self.param)
+        elif self.options["operating mode"] == "resistance":
+            model = pybamm.external_circuit.ResistanceFunctionControl(
+                self.param, "algebraic"
+            )
+        elif self.options["operating mode"] == "differential resistance":
+            model = pybamm.external_circuit.ResistanceFunctionControl(
+                self.param, "differential without max"
+            )
+        elif self.options["operating mode"] == "explicit resistance":
+            model = pybamm.external_circuit.ExplicitResistanceControl(self.param)
         elif self.options["operating mode"] == "CCCV":
-            self.submodels[
-                "external circuit"
-            ] = pybamm.external_circuit.CCCVFunctionControl(self.param)
+            model = pybamm.external_circuit.CCCVFunctionControl(self.param)
         elif callable(self.options["operating mode"]):
-            self.submodels[
-                "external circuit"
-            ] = pybamm.external_circuit.FunctionControl(
+            model = pybamm.external_circuit.FunctionControl(
                 self.param, self.options["operating mode"]
             )
+        self.submodels["external circuit"] = model
 
     def set_tortuosity_submodels(self):
         self.submodels["electrolyte tortuosity"] = pybamm.tortuosity.Bruggeman(
@@ -1139,10 +1167,11 @@ class BaseBatteryModel(pybamm.BaseModel):
         # Hack to avoid division by zero if i_cc is exactly zero
         # If i_cc is zero, i_cc_not_zero becomes 1. But multiplying by sign(i_cc) makes
         # the local resistance 'zero' (really, it's not defined when i_cc is zero)
-        i_cc_not_zero = ((i_cc > 0) + (i_cc < 0)) * i_cc + (i_cc >= 0) * (i_cc <= 0)
-        i_cc_dim_not_zero = ((i_cc_dim > 0) + (i_cc_dim < 0)) * i_cc_dim + (
-            i_cc_dim >= 0
-        ) * (i_cc_dim <= 0)
+        def x_not_zero(x):
+            return ((x > 0) + (x < 0)) * x + (x >= 0) * (x <= 0)
+
+        i_cc_not_zero = x_not_zero(i_cc)
+        i_cc_dim_not_zero = x_not_zero(i_cc_dim)
 
         self.variables.update(
             {
@@ -1193,9 +1222,16 @@ class BaseBatteryModel(pybamm.BaseModel):
             )
         )
 
-        # Power
+        # Power and resistance
         I_dim = self.variables["Current [A]"]
-        self.variables.update({"Terminal power [W]": I_dim * V_dim})
+        I_dim_not_zero = x_not_zero(I_dim)
+        self.variables.update(
+            {
+                "Terminal power [W]": I_dim * V_dim,
+                "Power [W]": I_dim * V_dim,
+                "Resistance [Ohm]": pybamm.sign(I_dim) * V_dim / I_dim_not_zero,
+            }
+        )
 
     def set_degradation_variables(self):
         """
diff --git a/pybamm/models/full_battery_models/lead_acid/loqs.py b/pybamm/models/full_battery_models/lead_acid/loqs.py
index 886d6effb..602fe4d7f 100644
--- a/pybamm/models/full_battery_models/lead_acid/loqs.py
+++ b/pybamm/models/full_battery_models/lead_acid/loqs.py
@@ -64,7 +64,7 @@ class LOQS(BaseModel):
         if self.options["operating mode"] == "current":
             self.submodels[
                 "leading order external circuit"
-            ] = pybamm.external_circuit.LeadingOrderCurrentControl(self.param)
+            ] = pybamm.external_circuit.LeadingOrderExplicitCurrentControl(self.param)
         elif self.options["operating mode"] == "voltage":
             self.submodels[
                 "leading order external circuit"
diff --git a/pybamm/models/submodels/current_collector/homogeneous_current_collector.py b/pybamm/models/submodels/current_collector/homogeneous_current_collector.py
index 4d65f29cd..a8b726793 100644
--- a/pybamm/models/submodels/current_collector/homogeneous_current_collector.py
+++ b/pybamm/models/submodels/current_collector/homogeneous_current_collector.py
@@ -21,6 +21,11 @@ class Uniform(BaseModel):
     def __init__(self, param):
         super().__init__(param)
 
+    def get_fundamental_variables(self):
+        phi_s_cn = pybamm.PrimaryBroadcast(0, "current collector")
+        variables = self._get_standard_negative_potential_variables(phi_s_cn)
+        return variables
+
     def get_coupled_variables(self, variables):
 
         # TODO: grad not implemented for 2D yet
@@ -28,10 +33,8 @@ class Uniform(BaseModel):
         i_boundary_cc = pybamm.PrimaryBroadcast(
             variables["Total current density"], "current collector"
         )
-        phi_s_cn = pybamm.PrimaryBroadcast(0, "current collector")
 
-        variables = self._get_standard_negative_potential_variables(phi_s_cn)
-        variables.update(self._get_standard_current_variables(i_cc, i_boundary_cc))
+        variables = self._get_standard_current_variables(i_cc, i_boundary_cc)
 
         # Hack to get the leading-order current collector current density
         # Note that this should be different from the actual (composite) current
diff --git a/pybamm/models/submodels/electrolyte_diffusion/full_diffusion.py b/pybamm/models/submodels/electrolyte_diffusion/full_diffusion.py
index a7252f4a1..d3c530c3b 100644
--- a/pybamm/models/submodels/electrolyte_diffusion/full_diffusion.py
+++ b/pybamm/models/submodels/electrolyte_diffusion/full_diffusion.py
@@ -113,13 +113,15 @@ class Full(BaseElectrolyteDiffusion):
             T = variables["Cell temperature"]
             tor = variables["Electrolyte tortuosity"]
             i_boundary_cc = variables["Current collector current density"]
-            dce_dx = (
-                -(1 - param.t_plus(c_e, T))
+            lbc = (
+                pybamm.boundary_value(
+                    -(1 - param.t_plus(c_e, T))
+                    / (tor * param.gamma_e * param.D_e(c_e, T)),
+                    "left",
+                )
                 * i_boundary_cc
                 * param.C_e
-                / (tor * param.gamma_e * param.D_e(c_e, T))
             )
-            lbc = pybamm.boundary_value(dce_dx, "left")
         else:
             # left bc at anode/current collector interface
             lbc = pybamm.Scalar(0)
diff --git a/pybamm/models/submodels/external_circuit/__init__.py b/pybamm/models/submodels/external_circuit/__init__.py
index 6cc5b1422..e84a8c0e8 100644
--- a/pybamm/models/submodels/external_circuit/__init__.py
+++ b/pybamm/models/submodels/external_circuit/__init__.py
@@ -1,9 +1,15 @@
 from .base_external_circuit import BaseModel, LeadingOrderBaseModel
-from .current_control_external_circuit import CurrentControl, LeadingOrderCurrentControl
+from .explicit_control_external_circuit import (
+    ExplicitCurrentControl,
+    ExplicitPowerControl,
+    ExplicitResistanceControl,
+    LeadingOrderExplicitCurrentControl,
+)
 from .function_control_external_circuit import (
     FunctionControl,
     VoltageFunctionControl,
     PowerFunctionControl,
+    ResistanceFunctionControl,
     CCCVFunctionControl,
     LeadingOrderFunctionControl,
     LeadingOrderVoltageFunctionControl,
diff --git a/pybamm/models/submodels/external_circuit/current_control_external_circuit.py b/pybamm/models/submodels/external_circuit/current_control_external_circuit.py
deleted file mode 100644
index c51bab2fb..000000000
--- a/pybamm/models/submodels/external_circuit/current_control_external_circuit.py
+++ /dev/null
@@ -1,38 +0,0 @@
-#
-# External circuit with current control
-#
-import pybamm
-from .base_external_circuit import BaseModel, LeadingOrderBaseModel
-
-
-class CurrentControl(BaseModel):
-    """External circuit with current control."""
-
-    def __init__(self, param):
-        super().__init__(param)
-
-    def get_fundamental_variables(self):
-        # Current is given as a function of time
-        i_cell = self.param.current_with_time
-        i_cell_dim = self.param.dimensional_current_density_with_time
-        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,
-            "C-rate": I / self.param.Q,
-        }
-
-        # Add discharge capacity variable
-        variables.update(super().get_fundamental_variables())
-
-        return variables
-
-
-class LeadingOrderCurrentControl(CurrentControl, LeadingOrderBaseModel):
-    """External circuit with current control, for leading order models."""
-
-    def __init__(self, param):
-        super().__init__(param)
diff --git a/pybamm/models/submodels/external_circuit/explicit_control_external_circuit.py b/pybamm/models/submodels/external_circuit/explicit_control_external_circuit.py
new file mode 100644
index 000000000..ccf4a0e3f
--- /dev/null
+++ b/pybamm/models/submodels/external_circuit/explicit_control_external_circuit.py
@@ -0,0 +1,98 @@
+#
+# External circuit with explicit equations for control
+#
+import pybamm
+from .base_external_circuit import BaseModel, LeadingOrderBaseModel
+
+
+class ExplicitCurrentControl(BaseModel):
+    """External circuit with current control."""
+
+    def __init__(self, param):
+        super().__init__(param)
+
+    def get_fundamental_variables(self):
+        # Current is given as a function of time
+        i_cell = self.param.current_with_time
+        i_cell_dim = self.param.dimensional_current_density_with_time
+        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,
+            "C-rate": I / self.param.Q,
+        }
+
+        # Add discharge capacity variable
+        variables.update(super().get_fundamental_variables())
+
+        return variables
+
+
+class ExplicitPowerControl(BaseModel):
+    """External circuit with current set explicitly to hit target power."""
+
+    def __init__(self, param):
+        super().__init__(param)
+
+    def get_coupled_variables(self, variables):
+        param = self.param
+
+        # Current is given as applied power divided by voltage
+        V = variables["Terminal voltage [V]"]
+        P = pybamm.FunctionParameter(
+            "Power function [W]", {"Time [s]": pybamm.t * self.param.timescale}
+        )
+        I = P / V
+
+        # Update derived variables
+        i_cell = I / abs(param.I_typ)
+        i_cell_dim = I / (param.n_electrodes_parallel * param.A_cc)
+
+        variables = {
+            "Total current density": i_cell,
+            "Total current density [A.m-2]": i_cell_dim,
+            "Current [A]": I,
+            "C-rate": I / param.Q,
+        }
+
+        return variables
+
+
+class ExplicitResistanceControl(BaseModel):
+    """External circuit with current set explicitly to hit target resistance."""
+
+    def __init__(self, param):
+        super().__init__(param)
+
+    def get_coupled_variables(self, variables):
+        param = self.param
+
+        # Current is given as applied voltage divided by resistance
+        V = variables["Terminal voltage [V]"]
+        R = pybamm.FunctionParameter(
+            "Resistance function [Ohm]", {"Time [s]": pybamm.t * self.param.timescale}
+        )
+        I = V / R
+
+        # Update derived variables
+        i_cell = I / abs(param.I_typ)
+        i_cell_dim = I / (param.n_electrodes_parallel * param.A_cc)
+
+        variables = {
+            "Total current density": i_cell,
+            "Total current density [A.m-2]": i_cell_dim,
+            "Current [A]": I,
+            "C-rate": I / param.Q,
+        }
+
+        return variables
+
+
+class LeadingOrderExplicitCurrentControl(ExplicitCurrentControl, LeadingOrderBaseModel):
+    """External circuit with current control, for leading order models."""
+
+    def __init__(self, param):
+        super().__init__(param)
diff --git a/pybamm/models/submodels/external_circuit/function_control_external_circuit.py b/pybamm/models/submodels/external_circuit/function_control_external_circuit.py
index 92cfdb989..567801603 100644
--- a/pybamm/models/submodels/external_circuit/function_control_external_circuit.py
+++ b/pybamm/models/submodels/external_circuit/function_control_external_circuit.py
@@ -30,9 +30,9 @@ class FunctionControl(BaseModel):
         param = self.param
         # Current is a variable
         i_var = pybamm.Variable("Current density variable")
-        if self.control == "algebraic":
+        if self.control in ["algebraic", "differential without max"]:
             i_cell = i_var
-        elif self.control == "differential":
+        elif self.control == "differential with max":
             i_cell = pybamm.maximum(i_var, param.current_with_time)
 
         # Update derived variables
@@ -63,7 +63,7 @@ class FunctionControl(BaseModel):
         # 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":
+        if "differential" in self.control:
             i_cell = variables["Current density variable"]
             self.rhs[i_cell] = self.external_circuit_function(variables)
 
@@ -94,15 +94,43 @@ class VoltageFunctionControl(FunctionControl):
 class PowerFunctionControl(FunctionControl):
     """External circuit with power control."""
 
-    def __init__(self, param):
-        super().__init__(param, self.constant_power, control="algebraic")
+    def __init__(self, param, control):
+        super().__init__(param, self.constant_power, control=control)
 
     def constant_power(self, variables):
         I = variables["Current [A]"]
         V = variables["Terminal voltage [V]"]
-        return I * V - pybamm.FunctionParameter(
+        P = V * I
+        P_applied = pybamm.FunctionParameter(
             "Power function [W]", {"Time [s]": pybamm.t * self.param.timescale}
         )
+        if self.control == "algebraic":
+            return P - P_applied
+        else:
+            # Multiply by the time scale so that the overshoot only lasts a few seconds
+            K_P = 0.01 * self.param.timescale
+            return -K_P * (P - P_applied)
+
+
+class ResistanceFunctionControl(FunctionControl):
+    """External circuit with resistance control."""
+
+    def __init__(self, param, control):
+        super().__init__(param, self.constant_resistance, control=control)
+
+    def constant_resistance(self, variables):
+        I = variables["Current [A]"]
+        V = variables["Terminal voltage [V]"]
+        R = V / I
+        R_applied = pybamm.FunctionParameter(
+            "Resistance function [Ohm]", {"Time [s]": pybamm.t * self.param.timescale}
+        )
+        if self.control == "algebraic":
+            return R - R_applied
+        else:
+            # Multiply by the time scale so that the overshoot only lasts a few seconds
+            K_R = 0.01 * self.param.timescale
+            return -K_R * (R - R_applied)
 
 
 class CCCVFunctionControl(FunctionControl):
@@ -119,7 +147,7 @@ class CCCVFunctionControl(FunctionControl):
     """
 
     def __init__(self, param):
-        super().__init__(param, self.cccv, control="differential")
+        super().__init__(param, self.cccv, control="differential with max")
         pybamm.citations.register("Mohtat2021")
 
     def cccv(self, variables):
diff --git a/pybamm/simulation.py b/pybamm/simulation.py
index 90fb70d38..2f69ec20a 100644
--- a/pybamm/simulation.py
+++ b/pybamm/simulation.py
@@ -390,7 +390,7 @@ class Simulation:
                     # check which kind of external circuit model we need (differential
                     # or algebraic)
                     if op_inputs["CCCV switch"] == 1:
-                        control = "differential"
+                        control = "differential with max"
                     else:
                         control = "algebraic"
                     # create the FunctionControl submodel and extract variables
@@ -478,7 +478,7 @@ class Simulation:
                         new_model.algebraic[
                             i_cell
                         ] = pybamm.external_circuit.PowerFunctionControl(
-                            new_model.param
+                            new_model.param, control="algebraic"
                         ).constant_power(
                             new_model.variables
                         )


Test patch:
diff --git a/tests/integration/test_models/test_submodels/test_external_circuit/test_function_control.py b/tests/integration/test_models/test_submodels/test_external_circuit/test_function_control.py
index 9454d6c73..61b6eff9e 100644
--- a/tests/integration/test_models/test_submodels/test_external_circuit/test_function_control.py
+++ b/tests/integration/test_models/test_submodels/test_external_circuit/test_function_control.py
@@ -101,42 +101,72 @@ class TestFunctionControl(unittest.TestCase):
             return I * V - 4
 
         # load models
+        # use DFN since only DFN allows "explicit power"
         models = [
-            pybamm.lithium_ion.SPM({"operating mode": "power"}),
-            pybamm.lithium_ion.SPM({"operating mode": constant_power}),
+            pybamm.lithium_ion.DFN({"operating mode": "power"}),
+            pybamm.lithium_ion.DFN({"operating mode": "explicit power"}),
+            pybamm.lithium_ion.DFN({"operating mode": constant_power}),
         ]
 
-        # load parameter values and process models and geometry
-        params = [model.default_parameter_values for model in models]
-
-        # First model: 4W charge
-        params[0].update({"Power function [W]": 4}, check_already_exists=False)
-
         # set parameters and discretise models
+        solutions = [None] * len(models)
+        t_eval = np.linspace(0, 3600, 100)
         for i, model in enumerate(models):
             # create geometry
             geometry = model.default_geometry
-            params[i].process_model(model)
-            params[i].process_geometry(geometry)
+            param = model.default_parameter_values
+            param.update({"Power function [W]": 4}, check_already_exists=False)
+            param.process_model(model)
+            param.process_geometry(geometry)
             mesh = pybamm.Mesh(
                 geometry, model.default_submesh_types, model.default_var_pts
             )
             disc = pybamm.Discretisation(mesh, model.default_spatial_methods)
             disc.process_model(model)
+            solutions[i] = model.default_solver.solve(model, t_eval)
 
-        # solve model
+        for var in ["Terminal voltage [V]", "Current [A]"]:
+            for sol in solutions[1:]:
+                np.testing.assert_array_almost_equal(
+                    solutions[0][var].data, sol[var].data
+                )
+
+    def test_constant_resistance(self):
+        def constant_resistance(variables):
+            I = variables["Current [A]"]
+            V = variables["Terminal voltage [V]"]
+            return V / I - 2
+
+        # load models
+        # use DFN since only DFN allows "explicit resistance"
+        models = [
+            pybamm.lithium_ion.DFN({"operating mode": "resistance"}),
+            pybamm.lithium_ion.DFN({"operating mode": "explicit resistance"}),
+            pybamm.lithium_ion.DFN({"operating mode": constant_resistance}),
+        ]
+
+        # set parameters and discretise models
         solutions = [None] * len(models)
         t_eval = np.linspace(0, 3600, 100)
         for i, model in enumerate(models):
+            # create geometry
+            geometry = model.default_geometry
+            param = model.default_parameter_values
+            param.update({"Resistance function [Ohm]": 2}, check_already_exists=False)
+            param.process_model(model)
+            param.process_geometry(geometry)
+            mesh = pybamm.Mesh(
+                geometry, model.default_submesh_types, model.default_var_pts
+            )
+            disc = pybamm.Discretisation(mesh, model.default_spatial_methods)
+            disc.process_model(model)
             solutions[i] = model.default_solver.solve(model, t_eval)
 
-        V0 = solutions[0]["Terminal voltage [V]"].entries
-        V1 = solutions[1]["Terminal voltage [V]"].entries
-        np.testing.assert_array_equal(V0, V1)
-
-        I0 = solutions[0]["Current [A]"].entries
-        I1 = solutions[1]["Current [A]"].entries
-        np.testing.assert_array_equal(I0, I1)
+        for var in ["Terminal voltage [V]", "Current [A]"]:
+            for sol in solutions[1:]:
+                np.testing.assert_array_almost_equal(
+                    solutions[0][var].data, sol[var].data
+                )
 
     def test_cccv(self):
         # load models
diff --git a/tests/unit/test_models/test_full_battery_models/test_base_battery_model.py b/tests/unit/test_models/test_full_battery_models/test_base_battery_model.py
index 73fe28b9d..a6cd53ac5 100644
--- a/tests/unit/test_models/test_full_battery_models/test_base_battery_model.py
+++ b/tests/unit/test_models/test_full_battery_models/test_base_battery_model.py
@@ -25,7 +25,7 @@ PRINT_OPTIONS_OUTPUT = """\
 'lithium plating': 'none' (possible: ['none', 'reversible', 'irreversible'])
 'lithium plating porosity change': 'false' (possible: ['false', 'true'])
 'loss of active material': 'stress-driven' (possible: ['none', 'stress-driven', 'reaction-driven', 'stress and reaction-driven'])
-'operating mode': 'current' (possible: ['current', 'voltage', 'power', 'CCCV'])
+'operating mode': 'current' (possible: ['current', 'voltage', 'power', 'differential power', 'explicit power', 'resistance', 'differential resistance', 'explicit resistance', 'CCCV'])
 'particle': 'Fickian diffusion' (possible: ['Fickian diffusion', 'fast diffusion', 'uniform profile', 'quadratic profile', 'quartic profile'])
 'particle mechanics': 'swelling only' (possible: ['none', 'swelling only', 'swelling and cracking'])
 'particle shape': 'spherical' (possible: ['spherical', 'no particles'])
diff --git a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_tests.py b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_tests.py
index 346ae9e8f..2aa77ce11 100644
--- a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_tests.py
+++ b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_tests.py
@@ -220,6 +220,18 @@ class BaseUnitTestLithiumIon:
         options = {"operating mode": "power"}
         self.check_well_posedness(options)
 
+    def test_well_posed_external_circuit_differential_power(self):
+        options = {"operating mode": "differential power"}
+        self.check_well_posedness(options)
+
+    def test_well_posed_external_circuit_resistance(self):
+        options = {"operating mode": "resistance"}
+        self.check_well_posedness(options)
+
+    def test_well_posed_external_circuit_differential_resistance(self):
+        options = {"operating mode": "differential resistance"}
+        self.check_well_posedness(options)
+
     def test_well_posed_external_circuit_cccv(self):
         options = {"operating mode": "CCCV"}
         self.check_well_posedness(options)
diff --git a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py
index 58a9150c8..c4620d687 100644
--- a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py
+++ b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py
@@ -23,6 +23,14 @@ class TestDFN(BaseUnitTestLithiumIon, unittest.TestCase):
         options = {"particle size": "distribution", "particle": "uniform profile"}
         self.check_well_posedness(options)
 
+    def test_well_posed_external_circuit_explicit_power(self):
+        options = {"operating mode": "explicit power"}
+        self.check_well_posedness(options)
+
+    def test_well_posed_external_circuit_explicit_resistance(self):
+        options = {"operating mode": "explicit resistance"}
+        self.check_well_posedness(options)
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")