Notebook

Equivalent Circuit Parameter Identification¶

This notebook provides example usage for identifying stationary parameters for a two RC branch Thevenin model. The Thevenin model represents an electrochemical battery through an empirical circuit model capable of capturing the electrical response of the battery. This model can be extended with a thermal submodel, as well as additional parallel resistor-capacitor branches.

Setting up the Environment¶

Before we begin, we need to ensure that we have all the necessary tools. We will install PyBOP from its development branch and upgrade some dependencies:

In [1]:

%pip install --upgrade pip ipywidgets
%pip install pybop -q

Requirement already satisfied: pip in /Users/engs2510/.pyenv/versions/pybop-3.12/lib/python3.12/site-packages (24.0)
Requirement already satisfied: ipywidgets in /Users/engs2510/.pyenv/versions/pybop-3.12/lib/python3.12/site-packages (8.1.2)
Requirement already satisfied: comm>=0.1.3 in /Users/engs2510/.pyenv/versions/pybop-3.12/lib/python3.12/site-packages (from ipywidgets) (0.2.2)
Requirement already satisfied: ipython>=6.1.0 in /Users/engs2510/.pyenv/versions/pybop-3.12/lib/python3.12/site-packages (from ipywidgets) (8.23.0)
Requirement already satisfied: traitlets>=4.3.1 in /Users/engs2510/.pyenv/versions/pybop-3.12/lib/python3.12/site-packages (from ipywidgets) (5.14.2)
Requirement already satisfied: widgetsnbextension~=4.0.10 in /Users/engs2510/.pyenv/versions/pybop-3.12/lib/python3.12/site-packages (from ipywidgets) (4.0.10)
Requirement already satisfied: jupyterlab-widgets~=3.0.10 in /Users/engs2510/.pyenv/versions/pybop-3.12/lib/python3.12/site-packages (from ipywidgets) (3.0.10)
Requirement already satisfied: decorator in /Users/engs2510/.pyenv/versions/pybop-3.12/lib/python3.12/site-packages (from ipython>=6.1.0->ipywidgets) (5.1.1)
Requirement already satisfied: jedi>=0.16 in /Users/engs2510/.pyenv/versions/pybop-3.12/lib/python3.12/site-packages (from ipython>=6.1.0->ipywidgets) (0.19.1)
Requirement already satisfied: matplotlib-inline in /Users/engs2510/.pyenv/versions/pybop-3.12/lib/python3.12/site-packages (from ipython>=6.1.0->ipywidgets) (0.1.6)
Requirement already satisfied: prompt-toolkit<3.1.0,>=3.0.41 in /Users/engs2510/.pyenv/versions/pybop-3.12/lib/python3.12/site-packages (from ipython>=6.1.0->ipywidgets) (3.0.43)
Requirement already satisfied: pygments>=2.4.0 in /Users/engs2510/.pyenv/versions/pybop-3.12/lib/python3.12/site-packages (from ipython>=6.1.0->ipywidgets) (2.17.2)
Requirement already satisfied: stack-data in /Users/engs2510/.pyenv/versions/pybop-3.12/lib/python3.12/site-packages (from ipython>=6.1.0->ipywidgets) (0.6.3)
Requirement already satisfied: pexpect>4.3 in /Users/engs2510/.pyenv/versions/pybop-3.12/lib/python3.12/site-packages (from ipython>=6.1.0->ipywidgets) (4.9.0)
Requirement already satisfied: parso<0.9.0,>=0.8.3 in /Users/engs2510/.pyenv/versions/pybop-3.12/lib/python3.12/site-packages (from jedi>=0.16->ipython>=6.1.0->ipywidgets) (0.8.4)
Requirement already satisfied: ptyprocess>=0.5 in /Users/engs2510/.pyenv/versions/pybop-3.12/lib/python3.12/site-packages (from pexpect>4.3->ipython>=6.1.0->ipywidgets) (0.7.0)
Requirement already satisfied: wcwidth in /Users/engs2510/.pyenv/versions/pybop-3.12/lib/python3.12/site-packages (from prompt-toolkit<3.1.0,>=3.0.41->ipython>=6.1.0->ipywidgets) (0.2.13)
Requirement already satisfied: executing>=1.2.0 in /Users/engs2510/.pyenv/versions/pybop-3.12/lib/python3.12/site-packages (from stack-data->ipython>=6.1.0->ipywidgets) (2.0.1)
Requirement already satisfied: asttokens>=2.1.0 in /Users/engs2510/.pyenv/versions/pybop-3.12/lib/python3.12/site-packages (from stack-data->ipython>=6.1.0->ipywidgets) (2.4.1)
Requirement already satisfied: pure-eval in /Users/engs2510/.pyenv/versions/pybop-3.12/lib/python3.12/site-packages (from stack-data->ipython>=6.1.0->ipywidgets) (0.2.2)
Requirement already satisfied: six>=1.12.0 in /Users/engs2510/.pyenv/versions/pybop-3.12/lib/python3.12/site-packages (from asttokens>=2.1.0->stack-data->ipython>=6.1.0->ipywidgets) (1.16.0)
Note: you may need to restart the kernel to use updated packages.
Note: you may need to restart the kernel to use updated packages.

Importing Libraries¶

With the environment set up, we can now import PyBOP alongside other libraries we will need:

In [2]:

import numpy as np

import pybop

Importing Parameters¶

This can be completed by importing a JSON representation, such as the one in the PyBOP examples. To import via JSON, either download the example file, or create your own and update the path below to reference the corresponding file.

In [3]:

# parameter_set = pybop.ParameterSet(
# json_path="examples/scripts/parameters/initial_ecm_parameters.json"
# )
# parameter_set.import_parameters()

Alternatively, define the initial parameter set with a dictionary. Ensure you have definitions for all R's, C's, and initial overpotentials for any additional RC elements.

In this example, we use the default parameter value for the "Open-circuit voltage [V] as provided by the original PyBaMM class. To update this, provide a function definition that matches this function.

In [4]:

parameter_set = pybop.ParameterSet(
    params_dict={
        "chemistry": "ecm",
        "Initial SoC": 0.5,
        "Initial temperature [K]": 25 + 273.15,
        "Cell capacity [A.h]": 5,
        "Nominal cell capacity [A.h]": 5,
        "Ambient temperature [K]": 25 + 273.15,
        "Current function [A]": 5,
        "Upper voltage cut-off [V]": 4.2,
        "Lower voltage cut-off [V]": 3.0,
        "Cell thermal mass [J/K]": 1000,
        "Cell-jig heat transfer coefficient [W/K]": 10,
        "Jig thermal mass [J/K]": 500,
        "Jig-air heat transfer coefficient [W/K]": 10,
        "Open-circuit voltage [V]": pybop.empirical.Thevenin().default_parameter_values[
            "Open-circuit voltage [V]"
        ],
        "R0 [Ohm]": 0.001,
        "Element-1 initial overpotential [V]": 0,
        "Element-2 initial overpotential [V]": 0,
        "R1 [Ohm]": 0.0002,
        "R2 [Ohm]": 0.0003,
        "C1 [F]": 13000,
        "C2 [F]": 36924,
        "Entropic change [V/K]": 0.0004,
    }
)

Identifying the Parameters¶

Now that the initial parameter set is constructed, we can start the PyBOP fitting process. First, we define the model class with two RC elements.

In [5]:

model = pybop.empirical.Thevenin(
    parameter_set=parameter_set, options={"number of rc elements": 2}
)

In this example, we are going to try to fit all five parameters at once. This isn't recommend for real-life application as identifiablity is challenging to guarantee with this large a parameter space. To do this, we define the pybop.parameters as,

In [6]:

parameters = [
    pybop.Parameter(
        "R0 [Ohm]",
        prior=pybop.Gaussian(0.0002, 0.0001),
        bounds=[1e-4, 1e-2],
    ),
    pybop.Parameter(
        "R1 [Ohm]",
        prior=pybop.Gaussian(0.0001, 0.0001),
        bounds=[1e-5, 1e-2],
    ),
    pybop.Parameter(
        "R2 [Ohm]",
        prior=pybop.Gaussian(0.0001, 0.0001),
        bounds=[1e-5, 1e-2],
    ),
    pybop.Parameter(
        "C1 [F]",
        prior=pybop.Gaussian(10000, 2500),
        bounds=[2500, 5e4],
    ),
    pybop.Parameter(
        "C1 [F]",
        prior=pybop.Gaussian(10000, 2500),
        bounds=[2500, 5e4],
    ),
]

Let's create some synthetic data to identify the parameters. This data is then corrupted with a small amount of Gaussian noise to represent some additional uncertainty in the measured values. We can then form the pybop.Dataset from this data.

In [7]:

sigma = 0.001
t_eval = np.arange(0, 900, 2)
values = model.predict(t_eval=t_eval)
corrupt_values = values["Voltage [V]"].data + np.random.normal(0, sigma, len(t_eval))

# Form dataset
dataset = pybop.Dataset(
    {
        "Time [s]": t_eval,
        "Current function [A]": values["Current [A]"].data,
        "Voltage [V]": corrupt_values,
    }
)

The FittingProblem class provides us with a single class that holds all of the objects we need to evaluate our selected SumSquaredError cost function.

In [8]:

problem = pybop.FittingProblem(model, parameters, dataset)
cost = pybop.SumSquaredError(problem)

The cost function can be interrogated manually via the cost([params]) API. In this example, that would look like the following,

In [9]:

cost([0.001, 0.001, 0.001, 5000, 5000])

Out[9]:

0.024944621550803514

Next, we construct the optimisation class with our algorithm of choice and run it. In this case, we select the CMA-ES method as it provides global optimisation capability. For the sake of reducing the runtime of this example, we limit the maximum iterations to 100; however, feel free to update this value.

In [10]:

optim = pybop.Optimisation(cost, optimiser=pybop.CMAES)
optim.set_max_iterations(300)
x, final_cost = optim.run()
print("Initial parameters:", cost.x0)
print("Estimated parameters:", x)

Initial parameters: [2.29696565e-04 3.53341865e-05 1.63145688e-05 1.07259649e+04
 9.73352990e+03]
Estimated parameters: [1.00367942e-03 3.80907013e-04 1.00242178e-04 1.07259636e+04
 9.73353073e+03]

Plotting and Visualisation¶

PyBOP provides various plotting utilities to visualize the results of the optimisation.

In [11]:

pybop.quick_plot(problem, parameter_values=x, title="Optimised Comparison");

Convergence and Parameter Trajectories¶

To assess the optimisation process, we can plot the convergence of the cost function and the trajectories of the parameters:

In [12]:

pybop.plot_convergence(optim)
pybop.plot_parameters(optim);

Conclusion¶

This notebook illustrates how to perform parameter estimation using CMA-ES in PyBOP, providing insights into the optimisation process through various visualisations.