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.
Before we begin, we need to ensure that we have all the necessary tools. We will install PyBOP and upgrade dependencies:
%pip install --upgrade pip ipywidgets -q
%pip install pybop -q
/Users/engs2510/Documents/Git/PyBOP/.nox/notebooks-overwrite/bin/python3: No module named pip Note: you may need to restart the kernel to use updated packages. /Users/engs2510/Documents/Git/PyBOP/.nox/notebooks-overwrite/bin/python3: No module named pip Note: you may need to restart the kernel to use updated packages.
With the environment set up, we can now import PyBOP alongside other libraries we will need:
import numpy as np
import pybop
pybop.PlotlyManager().pio.renderers.default = "notebook_connected"
# 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.
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,
}
)
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.
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,
parameters = pybop.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(
"C2 [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.
sigma = 0.001
t_eval = np.arange(0, 900, 3)
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.
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,
cost([0.001, 0.001, 0.001, 5000, 5000])
0.017319071604517344
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.
optim = pybop.CMAES(cost, max_iterations=300)
x, final_cost = optim.run()
print("Initial parameters:", optim.x0)
print("Estimated parameters:", x)
Initial parameters: [2.38024513e-04 1.33033378e-04 2.63603390e-05 1.15341178e+04 1.55762679e+04] Estimated parameters: [7.72266444e-04 5.85341699e-04 1.25813056e-04 1.15341189e+04 1.55762670e+04]
PyBOP provides various plotting utilities to visualize the results of the optimisation.
pybop.quick_plot(problem, problem_inputs=x, title="Optimised Comparison");
To assess the optimisation process, we can plot the convergence of the cost function and the trajectories of the parameters:
pybop.plot_convergence(optim)
pybop.plot_parameters(optim);
This notebook illustrates how to perform parameter estimation using CMA-ES in PyBOP, providing insights into the optimisation process through various visualisations.