Add GroupedSPMe model

README.md

@@ -92,10 +92,10 @@ The table below lists the currently supported [models](https://github.com/pybop-
 | Single Particle Model with Electrolyte (SPMe) | Particle Swarm Optimization (PSO)                        | Root Mean Squared Error (RMSE) <tr></tr> |
 | Doyle-Fuller-Newman (DFN)                     | Exponential Natural Evolution Strategy (xNES)            | Minkowski <tr></tr>                      |
 | Many Particle Model (MPM)                     | Separable Natural Evolution Strategy (sNES)              | Sum of Power <tr></tr>                   |
-| Multi-Species Multi-Reactants (MSMR)          | Adaptive Moment Estimation with Weight Decay (AdamW)     | Gaussian Log Likelihood <tr></tr>        |
+| Multi-Species Multi-Reaction (MSMR)           | Adaptive Moment Estimation with Weight Decay (AdamW)     | Gaussian Log Likelihood <tr></tr>        |
 | Weppner-Huggins                               | Improved Resilient Backpropagation (iRProp-)             | Log Posterior <tr></tr>                  |
 | Equivalent Circuit Models (ECM)               | SciPy Minimize & Differential Evolution                  | Unscented Kalman Filter (UKF) <tr></tr>  |
-|                                               | Cuckoo Search                                            | Gravimetric Energy Density <tr></tr>     |
+| Grouped-parameter SPMe (GroupedSPMe)          | Cuckoo Search                                            | Gravimetric Energy Density <tr></tr>     |
 |                                               | Gradient Descent                                         | Volumetric Energy Density<tr></tr>       |
 |                                               | Nelder-Mead                                              | <tr></tr>                                |
 

examples/scripts/comparison_examples/grouped_SPMe.py

@@ -0,0 +1,120 @@
+import matplotlib.pyplot as plt
+import numpy as np
+
+import pybop
+from pybop.models.lithium_ion.basic_SPMe import convert_physical_to_grouped_parameters
+
+# Prepare figure
+layout_options = dict(
+    xaxis_title="Time / s",
+    yaxis_title="Voltage / V",
+)
+plot_dict = pybop.plot.StandardPlot(layout_options=layout_options)
+
+# Unpack parameter values from Chen2020
+parameter_set = pybop.ParameterSet.pybamm("Chen2020")
+ce0 = parameter_set["Initial concentration in electrolyte [mol.m-3]"]
+T = parameter_set["Ambient temperature [K]"]
+parameter_set["Electrolyte diffusivity [m2.s-1]"] = parameter_set[
+    "Electrolyte diffusivity [m2.s-1]"
+](ce0, T)
+parameter_set["Electrolyte conductivity [S.m-1]"] = parameter_set[
+    "Electrolyte conductivity [S.m-1]"
+](ce0, T)
+
+# Make the conductivities artificially large
+parameter_set["Electrolyte conductivity [S.m-1]"] = 1e16  # 0.9487
+parameter_set["Negative electrode conductivity [S.m-1]"] = 1e16
+parameter_set["Positive electrode conductivity [S.m-1]"] = 1e16
+
+# Define a test protocol
+initial_state = {"Initial SoC": 0.9}
+experiment = pybop.Experiment(
+    [
+        "Discharge at 1C until 2.5 V (5 seconds period)",
+        "Rest for 30 minutes (5 seconds period)",
+        # "Charge at 2C until 4.1 V (5 seconds period)",
+        # "Rest for 30 minutes (5 seconds period)",
+    ],
+)
+
+# Run an example SPMe simulation
+model_options = {"surface form": "differential", "contact resistance": "true"}
+time_domain_SPMe = pybop.lithium_ion.SPMe(
+    parameter_set=parameter_set,
+    options=model_options,
+)
+simulation = time_domain_SPMe.predict(
+    initial_state=initial_state, experiment=experiment
+)
+dataset = pybop.Dataset(
+    {
+        "Time [s]": simulation["Time [s]"].data,
+        "Current function [A]": simulation["Current [A]"].data,
+        "Voltage [V]": simulation["Voltage [V]"].data,
+    }
+)
+plot_dict.add_traces(dataset["Time [s]"], dataset["Voltage [V]"])
+
+# Test model in the time domain
+grouped_parameter_set = convert_physical_to_grouped_parameters(parameter_set)
+time_domain_grouped = pybop.lithium_ion.GroupedSPMe(
+    parameter_set=grouped_parameter_set,
+    options=model_options,
+    build=True,
+)
+time_domain_grouped.set_initial_state(initial_state)
+time_domain_grouped.set_current_function(dataset)
+simulation = time_domain_grouped.predict(t_eval=dataset["Time [s]"])
+dataset = pybop.Dataset(
+    {
+        "Time [s]": simulation["Time [s]"].data,
+        "Current function [A]": simulation["Current [A]"].data,
+        "Voltage [V]": simulation["Voltage [V]"].data,
+    }
+)
+plot_dict.add_traces(dataset["Time [s]"], dataset["Voltage [V]"], line_dash="dash")
+plot_dict()
+
+# Set up figure
+fig, ax = plt.subplots()
+ax.grid()
+
+# Compare models in the frequency domain
+freq_domain_SPMe = pybop.lithium_ion.SPMe(
+    parameter_set=parameter_set, options=model_options, eis=True
+)
+freq_domain_grouped = pybop.lithium_ion.GroupedSPMe(
+    parameter_set=grouped_parameter_set,
+    options=model_options,
+    eis=True,
+    build=True,
+)
+
+for i, model in enumerate([freq_domain_SPMe, freq_domain_grouped]):
+    NSOC = 11
+    Nfreq = 60
+    fmin = 4e-4
+    fmax = 1e3
+    SOCs = np.linspace(0, 1, NSOC)
+    frequencies = np.logspace(np.log10(fmin), np.log10(fmax), Nfreq)
+
+    impedances = 1j * np.zeros((Nfreq, NSOC))
+    for ii, SOC in enumerate(SOCs):
+        model.set_initial_state({"Initial SoC": SOC})
+        simulation = model.simulateEIS(inputs=None, f_eval=frequencies)
+        impedances[:, ii] = simulation["Impedance"]
+
+        if i == 0:
+            ax.plot(np.real(impedances[:, ii]), -np.imag(impedances[:, ii]), "b")
+        if i == 1:
+            ax.plot(np.real(impedances[:, ii]), -np.imag(impedances[:, ii]), "r--")
+
+# Show figure
+ax.set(xlabel=r"$Z_r(\omega)$ [$\Omega$]", ylabel=r"$-Z_j(\omega)$ [$\Omega$]")
+ax.set_aspect("equal", "box")
+plt.show()
+fig.savefig("Nyquist.png")
+
+# mdic = {"Z": impedances, "f": frequencies, "SOC": SOCs}
+# savemat("Simulated data SPMe/Z_SPMe_SOC_Pybop_chen2020.mat", mdic)

pybop/models/base_model.py

@@ -241,30 +241,7 @@ def set_initial_state(self, initial_state: dict, inputs: Optional[Inputs] = None
         """
         self.clear()
 
-        initial_state = self.convert_to_pybamm_initial_state(initial_state)
-
-        if isinstance(self.pybamm_model, pybamm.equivalent_circuit.Thevenin):
-            initial_state = self.get_initial_state(initial_state, inputs=inputs)
-            self._unprocessed_parameter_set.update({"Initial SoC": initial_state})
-
-        else:
-            if not self.pybamm_model._built:  # noqa: SLF001
-                self.pybamm_model.build_model()
-
-            # Temporary construction of attributes for PyBaMM
-            self._model = self.pybamm_model
-            self._unprocessed_parameter_values = self._unprocessed_parameter_set
-
-            # Set initial state via PyBaMM's Simulation class
-            pybamm.Simulation.set_initial_soc(self, initial_state, inputs=inputs)
-
-            # Update the default parameter set for consistency
-            self._unprocessed_parameter_set = self._parameter_values
-
-            # Clear the pybamm objects
-            del self._model
-            del self._unprocessed_parameter_values
-            del self._parameter_values
+        self._set_initial_state(initial_state=initial_state, inputs=inputs)
 
         # Use a copy of the updated default parameter set
         self._parameter_set = self._unprocessed_parameter_set.copy()
@@ -331,11 +308,11 @@ def set_up_for_eis(self, model):
             model.param, None, model.options, control="algebraic"
         ).get_fundamental_variables()
 
-        # Perform the replacement
-        symbol_replacement_map = {
-            model.variables[name]: variable
-            for name, variable in external_circuit_variables.items()
-        }
+        # Define the variables to replace
+        symbol_replacement_map = {}
+        for name, variable in external_circuit_variables.items():
+            if name in model.variables.keys():
+                symbol_replacement_map[model.variables[name]] = variable
 
         # Don't replace initial conditions, as these should not contain
         # Variable objects

pybop/models/empirical/base_ecm.py

@@ -1,3 +1,5 @@
+from typing import Optional
+
 import numpy as np
 import pybamm
 
@@ -121,6 +123,22 @@ def _check_params(
             return self.param_checker(inputs, allow_infeasible_solutions)
         return True
 
+    def _set_initial_state(self, initial_state: dict, inputs: Optional[Inputs] = None):
+        """
+        Set the initial state of charge or concentrations for the battery model.
+
+        Parameters
+        ----------
+        initial_state : dict
+            A valid initial state, e.g. the initial state of charge or open-circuit voltage.
+        inputs : Inputs
+            The input parameters to be used when building the model.
+        """
+        initial_state = self.convert_to_pybamm_initial_state(initial_state)
+
+        initial_state = self.get_initial_state(initial_state, inputs=inputs)
+        self._unprocessed_parameter_set.update({"Initial SoC": initial_state})
+
     def get_initial_state(
         self,
         initial_value,

pybop/models/lithium_ion/__init__.py

@@ -2,4 +2,4 @@
 # Import lithium ion based models
 #
 from .base_echem import EChemBaseModel
-from .echem import SPM, SPMe, DFN, MPM, MSMR, WeppnerHuggins
+from .echem import SPM, SPMe, DFN, MPM, MSMR, WeppnerHuggins, GroupedSPMe

pybop/models/lithium_ion/base_echem.py

@@ -2,7 +2,7 @@
 import warnings
 from typing import Optional
 
-from pybamm import LithiumIonParameters
+from pybamm import LithiumIonParameters, Simulation
 from pybamm import lithium_ion as pybamm_lithium_ion
 
 from pybop.models.base_model import BaseModel, Inputs
@@ -154,6 +154,37 @@ def _check_params(
 
         return True
 
+    def _set_initial_state(self, initial_state: dict, inputs: Optional[Inputs] = None):
+        """
+        Set the initial state of charge or concentrations for the battery model.
+
+        Parameters
+        ----------
+        initial_state : dict
+            A valid initial state, e.g. the initial state of charge or open-circuit voltage.
+        inputs : Inputs
+            The input parameters to be used when building the model.
+        """
+        initial_state = self.convert_to_pybamm_initial_state(initial_state)
+
+        if not self.pybamm_model._built:  # noqa: SLF001
+            self.pybamm_model.build_model()
+
+        # Temporary construction of attributes for PyBaMM
+        self._model = self.pybamm_model
+        self._unprocessed_parameter_values = self._unprocessed_parameter_set
+
+        # Set initial state via PyBaMM's Simulation class
+        Simulation.set_initial_soc(self, initial_state, inputs=inputs)
+
+        # Update the default parameter set for consistency
+        self._unprocessed_parameter_set = self._parameter_values
+
+        # Clear the pybamm objects
+        del self._model
+        del self._unprocessed_parameter_values
+        del self._parameter_values
+
     def cell_volume(self, parameter_set: Optional[ParameterSet] = None):
         """
         Calculate the total cell volume in m3.