Tip

An interactive online version of this notebook is available, which can be accessed via Open this notebook in Google Colab

Alternatively, you may download this notebook and run it offline.

Using submodels in PyBaMM#

In this notebook we show how to modify existing models by swapping out submodels, and how to build your own model from scratch using existing submodels. To see all of the models and submodels available in PyBaMM, please take a look at the documentation here.

Changing a submodel in an exisiting battery model#

PyBaMM is designed to be a flexible modelling package that allows users to easily compare different models and numerical techniques within a common framework. Battery models within PyBaMM are built up using a number of submodels that describe different physics included within the model, such as mass conservation in the electrolyte or charge conservation in the solid. For ease of use, a number of popular battery models are pre-configured in PyBaMM. As an example, we look at the Single Particle Model (for more information see here).

First we import pybamm

[1]:

%pip install pybamm -q    # install PyBaMM if it is not installed
import pybamm

Note: you may need to restart the kernel to use updated packages.

Then we load the SPM

[2]:

model = pybamm.lithium_ion.SPM()

We can look at the submodels that make up the SPM by accessing model.submodels, which is a dictionary of the submodel’s name (i.e. the physics it represents) and the submodel that is selected

[3]:

for name, submodel in model.submodels.items():
    print(name, submodel)

external circuit <pybamm.models.submodels.external_circuit.explicit_control_external_circuit.ExplicitCurrentControl object at 0x1514783d0>
porosity <pybamm.models.submodels.porosity.constant_porosity.Constant object at 0x151478790>
Negative interface utilisation <pybamm.models.submodels.interface.interface_utilisation.full_utilisation.Full object at 0x151478940>
Positive interface utilisation <pybamm.models.submodels.interface.interface_utilisation.full_utilisation.Full object at 0x1514789d0>
negative particle mechanics <pybamm.models.submodels.particle_mechanics.no_mechanics.NoMechanics object at 0x151478a00>
positive particle mechanics <pybamm.models.submodels.particle_mechanics.no_mechanics.NoMechanics object at 0x151478a60>
negative primary active material <pybamm.models.submodels.active_material.constant_active_material.Constant object at 0x1514787f0>
positive primary active material <pybamm.models.submodels.active_material.constant_active_material.Constant object at 0x151478a90>
electrolyte transport efficiency <pybamm.models.submodels.transport_efficiency.bruggeman_transport_efficiency.Bruggeman object at 0x1514788e0>
electrode transport efficiency <pybamm.models.submodels.transport_efficiency.bruggeman_transport_efficiency.Bruggeman object at 0x151478ac0>
transverse convection <pybamm.models.submodels.convection.transverse.no_convection.NoConvection object at 0x151478af0>
through-cell convection <pybamm.models.submodels.convection.through_cell.no_convection.NoConvection object at 0x151478a30>
negative primary open-circuit potential <pybamm.models.submodels.interface.open_circuit_potential.single_ocp.SingleOpenCircuitPotential object at 0x151478bb0>
positive primary open-circuit potential <pybamm.models.submodels.interface.open_circuit_potential.single_ocp.SingleOpenCircuitPotential object at 0x151478b50>
negative interface <pybamm.models.submodels.interface.kinetics.inverse_kinetics.inverse_butler_volmer.InverseButlerVolmer object at 0x151478b20>
negative interface current <pybamm.models.submodels.interface.kinetics.inverse_kinetics.inverse_butler_volmer.CurrentForInverseButlerVolmer object at 0x151478c40>
positive interface <pybamm.models.submodels.interface.kinetics.inverse_kinetics.inverse_butler_volmer.InverseButlerVolmer object at 0x151478c70>
positive interface current <pybamm.models.submodels.interface.kinetics.inverse_kinetics.inverse_butler_volmer.CurrentForInverseButlerVolmer object at 0x151478c10>
negative primary particle <pybamm.models.submodels.particle.fickian_diffusion.FickianDiffusion object at 0x151478ca0>
negative primary total particle concentration <pybamm.models.submodels.particle.total_particle_concentration.TotalConcentration object at 0x151478cd0>
positive primary particle <pybamm.models.submodels.particle.fickian_diffusion.FickianDiffusion object at 0x151478b80>
positive primary total particle concentration <pybamm.models.submodels.particle.total_particle_concentration.TotalConcentration object at 0x151478d30>
negative electrode potential <pybamm.models.submodels.electrode.ohm.leading_ohm.LeadingOrder object at 0x151478d00>
positive electrode potential <pybamm.models.submodels.electrode.ohm.leading_ohm.LeadingOrder object at 0x151478d60>
electrolyte diffusion <pybamm.models.submodels.electrolyte_diffusion.constant_concentration.ConstantConcentration object at 0x151478be0>
leading-order electrolyte conductivity <pybamm.models.submodels.electrolyte_conductivity.leading_order_conductivity.LeadingOrder object at 0x151478d90>
negative surface potential difference <pybamm.models.submodels.electrolyte_conductivity.surface_potential_form.explicit_surface_form_conductivity.Explicit object at 0x151478df0>
positive surface potential difference <pybamm.models.submodels.electrolyte_conductivity.surface_potential_form.explicit_surface_form_conductivity.Explicit object at 0x151478e20>
thermal <pybamm.models.submodels.thermal.isothermal.Isothermal object at 0x151478dc0>
current collector <pybamm.models.submodels.current_collector.homogeneous_current_collector.Uniform object at 0x151478e50>
primary sei <pybamm.models.submodels.interface.sei.no_sei.NoSEI object at 0x151478ee0>
primary sei on cracks <pybamm.models.submodels.interface.sei.no_sei.NoSEI object at 0x151478f40>
lithium plating <pybamm.models.submodels.interface.lithium_plating.no_plating.NoPlating object at 0x151478eb0>
total interface <pybamm.models.submodels.interface.total_interfacial_current.TotalInterfacialCurrent object at 0x151454070>

When you load a model in PyBaMM it builds by default. Building the model sets all of the model variables and sets up any variables which are coupled between different submodels: this is the process which couples the submodels together and allows one submodel to access variables from another. If you would like to swap out a submodel in an existing battery model you need to load it without building it by passing the keyword build=False

[4]:

model = pybamm.lithium_ion.SPM(build=False)

This collects all of the submodels which make up the SPM, but doesn’t build the model. Now you are free to swap out one submodel for another. For instance, you may want to assume that diffusion within the negative particles is infinitely fast, so that the PDE describing diffusion is replaced with an ODE for the uniform particle concentration. To change a submodel you simply update the dictionary entry, in this case to the XAveragedPolynomialProfile submodel

[5]:

model.submodels["negative primary particle"] = pybamm.particle.XAveragedPolynomialProfile(
    model.param, "negative", options={**model.options, "particle": "uniform profile"}
)

where we pass in the model parameters, the electrode (negative or positive) the submodel corresponds to, and the name of the polynomial we want to use. In the example we assume uniform concentration within the particle, corresponding to a zero-order polynomial.

Now if we look at the submodels again we see that the model for the negative particle has been changed

[6]:

for name, submodel in model.submodels.items():
    print(name, submodel)

external circuit <pybamm.models.submodels.external_circuit.explicit_control_external_circuit.ExplicitCurrentControl object at 0x151594940>
porosity <pybamm.models.submodels.porosity.constant_porosity.Constant object at 0x151594d00>
Negative interface utilisation <pybamm.models.submodels.interface.interface_utilisation.full_utilisation.Full object at 0x151594eb0>
Positive interface utilisation <pybamm.models.submodels.interface.interface_utilisation.full_utilisation.Full object at 0x151594f40>
negative particle mechanics <pybamm.models.submodels.particle_mechanics.no_mechanics.NoMechanics object at 0x151594f70>
positive particle mechanics <pybamm.models.submodels.particle_mechanics.no_mechanics.NoMechanics object at 0x151594fd0>
negative primary active material <pybamm.models.submodels.active_material.constant_active_material.Constant object at 0x151594d60>
positive primary active material <pybamm.models.submodels.active_material.constant_active_material.Constant object at 0x1515a2070>
electrolyte transport efficiency <pybamm.models.submodels.transport_efficiency.bruggeman_transport_efficiency.Bruggeman object at 0x151594e50>
electrode transport efficiency <pybamm.models.submodels.transport_efficiency.bruggeman_transport_efficiency.Bruggeman object at 0x151594fa0>
transverse convection <pybamm.models.submodels.convection.transverse.no_convection.NoConvection object at 0x1515a20a0>
through-cell convection <pybamm.models.submodels.convection.through_cell.no_convection.NoConvection object at 0x1515a2040>
negative primary open-circuit potential <pybamm.models.submodels.interface.open_circuit_potential.single_ocp.SingleOpenCircuitPotential object at 0x1515a2160>
positive primary open-circuit potential <pybamm.models.submodels.interface.open_circuit_potential.single_ocp.SingleOpenCircuitPotential object at 0x1515a2100>
negative interface <pybamm.models.submodels.interface.kinetics.inverse_kinetics.inverse_butler_volmer.InverseButlerVolmer object at 0x1515a20d0>
negative interface current <pybamm.models.submodels.interface.kinetics.inverse_kinetics.inverse_butler_volmer.CurrentForInverseButlerVolmer object at 0x1515a21f0>
positive interface <pybamm.models.submodels.interface.kinetics.inverse_kinetics.inverse_butler_volmer.InverseButlerVolmer object at 0x1515a2220>
positive interface current <pybamm.models.submodels.interface.kinetics.inverse_kinetics.inverse_butler_volmer.CurrentForInverseButlerVolmer object at 0x1515a21c0>
negative primary particle <pybamm.models.submodels.particle.x_averaged_polynomial_profile.XAveragedPolynomialProfile object at 0x15154bb50>
negative primary total particle concentration <pybamm.models.submodels.particle.total_particle_concentration.TotalConcentration object at 0x1515a2280>
positive primary particle <pybamm.models.submodels.particle.fickian_diffusion.FickianDiffusion object at 0x1515a2130>
positive primary total particle concentration <pybamm.models.submodels.particle.total_particle_concentration.TotalConcentration object at 0x1515a22e0>
negative electrode potential <pybamm.models.submodels.electrode.ohm.leading_ohm.LeadingOrder object at 0x1515a22b0>
positive electrode potential <pybamm.models.submodels.electrode.ohm.leading_ohm.LeadingOrder object at 0x1515a2310>
electrolyte diffusion <pybamm.models.submodels.electrolyte_diffusion.constant_concentration.ConstantConcentration object at 0x1515a2190>
leading-order electrolyte conductivity <pybamm.models.submodels.electrolyte_conductivity.leading_order_conductivity.LeadingOrder object at 0x1515a2340>
negative surface potential difference <pybamm.models.submodels.electrolyte_conductivity.surface_potential_form.explicit_surface_form_conductivity.Explicit object at 0x1515a23a0>
positive surface potential difference <pybamm.models.submodels.electrolyte_conductivity.surface_potential_form.explicit_surface_form_conductivity.Explicit object at 0x1515a23d0>
thermal <pybamm.models.submodels.thermal.isothermal.Isothermal object at 0x1515a2370>
current collector <pybamm.models.submodels.current_collector.homogeneous_current_collector.Uniform object at 0x1515a2400>
primary sei <pybamm.models.submodels.interface.sei.no_sei.NoSEI object at 0x1515a2490>
primary sei on cracks <pybamm.models.submodels.interface.sei.no_sei.NoSEI object at 0x1515a24f0>
lithium plating <pybamm.models.submodels.interface.lithium_plating.no_plating.NoPlating object at 0x1515a2460>
total interface <pybamm.models.submodels.interface.total_interfacial_current.TotalInterfacialCurrent object at 0x1515a2550>

Building the model also sets up the equations, boundary and initial conditions for the model. For example, if we look at model.rhs before building we see that it is empty

[7]:

model.rhs

[7]:

{}

If we try to use this empty model, PyBaMM will give an error. So, before proceeding we must build the model

[8]:

model.build_model()

Now if we look at model.rhs we see that it contains an entry relating to the concentration in each particle, as expected for the SPM

[9]:

model.rhs

[9]:

{Variable(0x620af1e84efc93fa, Discharge capacity [A.h], children=[], domains={}): Multiplication(0x3098e50eb9cc5275, *, children=['0.0002777777777777778', 'Current function [A]'], domains={}),
 Variable(-0x5e5303cde5e32a1d, Average negative particle concentration [mol.m-3], children=[], domains={'primary': ['current collector']}): MatrixMultiplication(0x26225f38ea92e2a0, @, children=['mass(Average negative particle concentration [mol.m-3])', '-3.109280896985319e-05 * Current function [A] / (Number of electrodes connected in parallel to make a cell * Electrode width [m] * Electrode height [m]) / Negative electrode thickness [m] / x-average(3.0 * Negative electrode active material volume fraction / Negative particle radius [m]) / x-average(Negative particle radius [m])'], domains={'primary': ['current collector']}),
 Variable(0x2e6e9aee084a77f, X-averaged positive particle concentration [mol.m-3], children=[], domains={'primary': ['positive particle'], 'secondary': ['current collector']}): Divergence(-0x6a7c97412b8b9861, div, children=['Positive electrode diffusivity [m2.s-1] * grad(X-averaged positive particle concentration [mol.m-3])'], domains={'primary': ['positive particle'], 'secondary': ['current collector']})}

Now the model can be used in a simulation and solved in the usual way, and we still have access to model defaults such as the default geometry and default spatial methods which are used in the simulation

[10]:

simulation = pybamm.Simulation(model)
simulation.solve([0, 3600])
simulation.plot()

[10]:

<pybamm.plotting.quick_plot.QuickPlot at 0x1515c5a00>

Building a custom model from submodels#

Instead of editing a pre-existing model, you may wish to build your own model from scratch by combining existing submodels of you choice. In this section, we build a Single Particle Model in which the diffusion is assumed infinitely fast in both particles.

To begin, we load a base lithium-ion model. This sets up the basic model structure behind the scenes, and also sets the default parameters to be those corresponding to a lithium-ion battery. Note that the base model does not select any default submodels, so there is no need to pass build=False.

[11]:

model = pybamm.lithium_ion.BaseModel()

Submodels can be added to the model.submodels dictionary in the same way that we changed the submodels earlier.

We use the simplest model for the external circuit, which is the explicit “current control” submodel

[12]:

model.submodels["external circuit"] = pybamm.external_circuit.ExplicitCurrentControl(model.param, model.options)

We want to build a 1D model, so select the Uniform current collector model (if the current collectors are behaving uniformly, then a 1D model is appropriate). We also want the model to be isothermal, so select the thermal model accordingly. Further, we assume that the porosity and active material are constant in space and time.

[13]:

model.submodels["current collector"] = pybamm.current_collector.Uniform(model.param)
model.submodels["thermal"] = pybamm.thermal.isothermal.Isothermal(model.param)
model.submodels["porosity"] = pybamm.porosity.Constant(model.param, model.options)
model.submodels["negative active material"] = pybamm.active_material.Constant(
    model.param, "negative", model.options
)
model.submodels["positive active material"] = pybamm.active_material.Constant(
    model.param, "positive", model.options
)

We assume that the current density varies linearly in the electrodes. This corresponds the the leading-order terms in Ohm’s law in the limit in which the SPM is derived in [3]

[14]:

model.submodels["negative electrode potentials"] = pybamm.electrode.ohm.LeadingOrder(
    model.param, "negative"
)
model.submodels["positive electrode potentials"] = pybamm.electrode.ohm.LeadingOrder(
    model.param, "positive"
)

We assume uniform concentration in both the negative and positive particles. We also have to separately specify a model for the total concentration in each electrode, which is calculated from the concentration in the particles (not a separate ODE)

[22]:

options = {**model.options, "particle": "uniform profile"}
model.submodels["negative primary particle"] = pybamm.particle.XAveragedPolynomialProfile(model.param, "negative", options)
model.submodels["positive primary particle"] = pybamm.particle.XAveragedPolynomialProfile(model.param, "positive", options)

model.submodels["negative total particle concentration"] = pybamm.particle.TotalConcentration(model.param, "negative", options)
model.submodels["positive total particle concentration"] = pybamm.particle.TotalConcentration(model.param, "positive", options)

In the Single Particle Model, the overpotential can be obtained by inverting the Butler-Volmer relation, so we choose the InverseButlerVolmer submodel for the interface, with the “main” lithium-ion reaction (and default lithium ion options). Because of how the current is implemented, we also need to separately specify the CurrentForInverseButlerVolmer submodel. We also need to specify the submodel for open-circuit potential.

[16]:

model.submodels[
    "negative open-circuit potential"
] = pybamm.open_circuit_potential.SingleOpenCircuitPotential(
    model.param, "negative", "lithium-ion main", options=model.options
)
model.submodels[
    "positive open-circuit potential"
] = pybamm.open_circuit_potential.SingleOpenCircuitPotential(
    model.param, "positive", "lithium-ion main", options=model.options
)
model.submodels[
    "negative interface"
] = pybamm.kinetics.InverseButlerVolmer(
    model.param, "negative", "lithium-ion main", options=model.options
)
model.submodels[
    "positive interface"
] = pybamm.kinetics.InverseButlerVolmer(
    model.param, "positive", "lithium-ion main", options=model.options
)
model.submodels[
    "negative interface current"
] = pybamm.kinetics.CurrentForInverseButlerVolmer(
    model.param, "negative", "lithium-ion main"
)
model.submodels[
    "positive interface current"
] = pybamm.kinetics.CurrentForInverseButlerVolmer(
    model.param, "positive", "lithium-ion main"
)
model.submodels["negative interface utilisation"] = pybamm.interface_utilisation.Full(
    model.param, "negative", model.options
)
model.submodels["positive interface utilisation"] = pybamm.interface_utilisation.Full(
    model.param, "positive", model.options
)

We don’t want any particle mechanics, SEI formation or lithium plating in this model

[17]:

model.submodels[
    "Negative particle mechanics"
] = pybamm.particle_mechanics.NoMechanics(model.param, "negative", model.options)
model.submodels[
    "Positive particle mechanics"
] = pybamm.particle_mechanics.NoMechanics(model.param, "positive", model.options)
model.submodels["sei"] = pybamm.sei.NoSEI(model.param, model.options)
model.submodels["sei on cracks"] = pybamm.sei.NoSEI(model.param, model.options, cracks=True)
model.submodels["lithium plating"] = pybamm.lithium_plating.NoPlating(model.param)

Finally, for the electrolyte we assume that diffusion is infinitely fast so that the concentration is uniform, and also use the leading-order model for charge conservation, which leads to a linear variation in ionic current in the electrodes

[18]:

model.submodels["electrolyte diffusion"] = pybamm.electrolyte_diffusion.ConstantConcentration(
    model.param
)
model.submodels["electrolyte conductivity"] = pybamm.electrolyte_conductivity.LeadingOrder(
    model.param
)

Now that we have set all of the submodels we can build the model

[19]:

model.build_model()

We can then use the model in a simulation in the usual way

[20]:

simulation = pybamm.Simulation(model)
simulation.solve([0, 3600])
simulation.plot()

[20]:

<pybamm.plotting.quick_plot.QuickPlot at 0x1518a5df0>

References#

The relevant papers for this notebook are:

[21]:

pybamm.print_citations()

[1] Joel A. E. Andersson, Joris Gillis, Greg Horn, James B. Rawlings, and Moritz Diehl. CasADi – A software framework for nonlinear optimization and optimal control. Mathematical Programming Computation, 11(1):1–36, 2019. doi:10.1007/s12532-018-0139-4.
[2] Charles R. Harris, K. Jarrod Millman, Stéfan J. van der Walt, Ralf Gommers, Pauli Virtanen, David Cournapeau, Eric Wieser, Julian Taylor, Sebastian Berg, Nathaniel J. Smith, and others. Array programming with NumPy. Nature, 585(7825):357–362, 2020. doi:10.1038/s41586-020-2649-2.
[3] Scott G. Marquis, Valentin Sulzer, Robert Timms, Colin P. Please, and S. Jon Chapman. An asymptotic derivation of a single particle model with electrolyte. Journal of The Electrochemical Society, 166(15):A3693–A3706, 2019. doi:10.1149/2.0341915jes.
[4] Venkat R. Subramanian, Vinten D. Diwakar, and Deepak Tapriyal. Efficient macro-micro scale coupled modeling of batteries. Journal of The Electrochemical Society, 152(10):A2002, 2005. doi:10.1149/1.2032427.
[5] Valentin Sulzer, Scott G. Marquis, Robert Timms, Martin Robinson, and S. Jon Chapman. Python Battery Mathematical Modelling (PyBaMM). Journal of Open Research Software, 9(1):14, 2021. doi:10.5334/jors.309.