Merge pull request #2478 from pybamm-team/issue-1143-equivalent-circuit-model

Issue 1143 equivalent circuit model

Author: Scottmar93

.gitignore

@@ -13,6 +13,7 @@
 !CMakeLists.txt
 !pybamm/CITATIONS.txt
 !pybamm/input/**/*.csv
+!tests/unit/test_parameters/*.csv
 !benchmarks/benchmark_images/*.png
 
 # running files

CHANGELOG.md

@@ -2,6 +2,7 @@
 
 ## Features
 
+- Equivalent circuit models ([#2478](https://github.com/pybamm-team/PyBaMM/pull/2478))
 - Added `scale` and `reference` attributes to `Variable` objects, which can be use to make the ODE/DAE solver better conditioned ([#2440](https://github.com/pybamm-team/PyBaMM/pull/2440))
 - SEI reactions can now be asymmetric ([#2425](https://github.com/pybamm-team/PyBaMM/pull/2425))
 - New Idaklu solver options for jacobian type and linear solver, support Sundials v6 ([#2444](https://github.com/pybamm-team/PyBaMM/pull/2444))

docs/source/models/equivalent_circuit/index.rst

@@ -0,0 +1,7 @@
+Equivalent Circuit Models
+=========================
+
+.. toctree::
+  :maxdepth: 1
+
+  thevenin

docs/source/models/equivalent_circuit/thevenin.rst

@@ -0,0 +1,5 @@
+Thevenin Model
+==============
+
+.. autoclass:: pybamm.equivalent_circuit.Thevenin
+    :members:

docs/source/models/index.rst

@@ -14,4 +14,5 @@ to see how these models can be solved, and compared, using PyBaMM.
   base_models/index
   lithium_ion/index
   lead_acid/index
+  equivalent_circuit/index
   submodels/index

docs/source/models/submodels/equivalent_circuit_elements/index.rst

@@ -0,0 +1,10 @@
+Equivalent Circuit Elements
+===========================
+
+.. toctree::
+
+  ocv_element
+  resistor_element
+  rc_element
+  thermal
+  voltage_model

docs/source/models/submodels/equivalent_circuit_elements/ocv_element.rst

@@ -0,0 +1,6 @@
+OCV Element
+===========
+
+.. autoclass:: pybamm.equivalent_circuit_elements.OCVElement
+    :members:
+

docs/source/models/submodels/equivalent_circuit_elements/rc_element.rst

@@ -0,0 +1,6 @@
+RC Element
+==========
+
+.. autoclass:: pybamm.equivalent_circuit_elements.RCElement
+    :members:
+

docs/source/models/submodels/equivalent_circuit_elements/resistor_element.rst

@@ -0,0 +1,6 @@
+Resistor Element
+================
+
+.. autoclass:: pybamm.equivalent_circuit_elements.ResistorElement
+    :members:
+

docs/source/models/submodels/equivalent_circuit_elements/thermal.rst

@@ -0,0 +1,6 @@
+Thermal SubModel
+================
+
+.. autoclass:: pybamm.equivalent_circuit_elements.ThermalSubModel
+    :members:
+

docs/source/models/submodels/equivalent_circuit_elements/voltage_model.rst

@@ -0,0 +1,6 @@
+Voltage Model
+=============
+
+.. autoclass:: pybamm.equivalent_circuit_elements.VoltageModel
+    :members:
+

docs/source/models/submodels/index.rst

@@ -19,3 +19,4 @@ Submodels
     porosity/index
     thermal/index
     transport_efficiency/index
+    equivalent_circuit_elements/index

docs/source/parameters/index.rst

@@ -10,3 +10,4 @@ Parameters
   lithium_ion_parameters
   lead_acid_parameters
   parameter_sets
+  process_parameter_data

docs/source/parameters/process_parameter_data.rst

@@ -0,0 +1,10 @@
+Process Parameter Data
+======================
+
+.. autofunction:: pybamm.parameters.process_1D_data
+
+.. autofunction:: pybamm.parameters.process_2D_data
+
+.. autofunction:: pybamm.parameters.process_2D_data_csv
+
+.. autofunction:: pybamm.parameters.process_3D_data_csv

examples/scripts/run_ecm.py

@@ -0,0 +1,44 @@
+import pybamm
+
+pybamm.set_logging_level("INFO")
+
+# options = {"number of rc elements": 2}
+options = {}
+model = pybamm.equivalent_circuit.Thevenin(options=options)
+
+parameter_values = model.default_parameter_values
+# parameter_values.update(
+#     {
+#         "R2 [Ohm]": 0.3e-3,
+#         "C2 [F]": 1000 / 0.3e-3,
+#         "Element-2 initial overpotential [V]": 0,
+#     },
+#     check_already_exists=False,
+# )
+
+experiment = pybamm.Experiment(
+    [
+        (
+            "Discharge at C/10 for 10 hours or until 3.3 V",
+            "Rest for 1 hour",
+            "Charge at 100 A until 4.1 V (1 second period)",
+            "Hold at 4.1 V until 5 A (1 seconds period)",
+            "Rest for 1 hour",
+        ),
+    ]
+)
+
+sim = pybamm.Simulation(model, experiment=experiment, parameter_values=parameter_values)
+sim.solve()
+sim.plot(
+    output_variables=[
+        "SoC",
+        "Open circuit voltage [V]",
+        "Current [A]",
+        "Cell temperature [degC]",
+        "Entropic change [V/K]",
+        "R0 [Ohm]",
+        "R1 [Ohm]",
+        "C1 [F]",
+    ]
+)

pybamm/__init__.py

@@ -111,6 +111,7 @@
 )
 from .models.full_battery_models import lead_acid
 from .models.full_battery_models import lithium_ion
+from .models.full_battery_models import equivalent_circuit
 
 #
 # Submodel classes
@@ -132,6 +133,7 @@
     thermal,
     transport_efficiency,
     particle_mechanics,
+    equivalent_circuit_elements,
 )
 from .models.submodels.interface import kinetics
 from .models.submodels.interface import sei
@@ -161,6 +163,7 @@
 from .parameters.thermal_parameters import thermal_parameters, ThermalParameters
 from .parameters.lithium_ion_parameters import LithiumIonParameters
 from .parameters.lead_acid_parameters import LeadAcidParameters
+from .parameters.ecm_parameters import EcmParameters
 from .parameters.size_distribution_parameters import *
 from .parameters.parameter_sets import parameter_sets
 from .parameters_cli import add_parameter, remove_parameter, edit_parameter

pybamm/expression_tree/broadcasts.py

@@ -549,7 +549,9 @@ def full_like(symbols, fill_value):
         return array_type(entries, domains=sum_symbol.domains)
 
     except NotImplementedError:
-        if sum_symbol.shape_for_testing == (1, 1):
+        if sum_symbol.shape_for_testing == (1, 1) or sum_symbol.shape_for_testing == (
+            1,
+        ):
             return pybamm.Scalar(fill_value)
         if sum_symbol.evaluates_on_edges("primary"):
             return FullBroadcastToEdges(

pybamm/expression_tree/interpolant.py

@@ -68,12 +68,12 @@ def __init__(
             x1, x2 = x
             if y.ndim != 2:
                 raise ValueError("y should be two-dimensional if len(x)=2")
-            if x1.shape[0] != y.shape[1]:
+            if x1.shape[0] != y.shape[0]:
                 raise ValueError(
                     "len(x1) should equal y=shape[1], "
                     f"but x1.shape={x1.shape} and y.shape={y.shape}"
                 )
-            if x2 is not None and x2.shape[0] != y.shape[0]:
+            if x2 is not None and x2.shape[0] != y.shape[1]:
                 raise ValueError(
                     "len(x2) should equal y=shape[0], "
                     f"but x2.shape={x2.shape} and y.shape={y.shape}"
@@ -153,9 +153,18 @@ def __init__(
                     "interpolator should be 'linear' or 'cubic' if x is two-dimensional"
                 )
             else:
-                interpolating_function = interpolate.interp2d(
-                    x1, x2, y, kind=interpolator
+                if extrapolate:
+                    fill_value = None
+                else:
+                    fill_value = np.nan
+                interpolating_function = interpolate.RegularGridInterpolator(
+                    (x1, x2),
+                    y,
+                    method=interpolator,
+                    bounds_error=False,
+                    fill_value=fill_value,
                 )
+
         elif len(x) == 3:
             self.dimension = 3
 
@@ -241,13 +250,7 @@ def _function_evaluate(self, evaluated_children):
                 children_eval_flat.append(child)
         if self.dimension == 1:
             return self.function(*children_eval_flat).flatten()[:, np.newaxis]
-        elif self.dimension == 2:
-            res = self.function(*children_eval_flat)
-            if res.ndim > 1:
-                return np.diagonal(res)[:, np.newaxis]
-            else:
-                return res[:, np.newaxis]
-        elif self.dimension == 3:
+        elif self.dimension in [2, 3]:
 
             # If the children are scalars, we need to add a dimension
             shapes = []