Here, we provide a short example of how to identify equivalent-circuit parameters using the scipy trust-constr method -- a trust-region based optimiser that can handle parameter bounds, and both linear and nonlinear constraints. As shown here, these turn out to be useful tools for fitting empirical battery models.
If you don't already have PyBOP installed, check out the installation guide first.
We begin by importing the necessary libraries.
import numpy as np
import scipy.optimize
import pybop
PyBOP needs to know some model parameters in order to run. Where these are unknown, and to be fitted, any sensible initial guess can be provided.
Model parameters can either be defined using a PyBOP ParameterSet object, or by importing from a JSON file. For clarity, we will define the initial parameters from a ParameterSet.
Here, we'll fit a Thevenin model with two RC pairs. The initial parameter set must contain all the parameters that the PyBaMM model will use, and therefore it will need definitions for each circuit component, and each RC pair overpotential.
parameter_set = {
"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.01,
"Element-1 initial overpotential [V]": 0,
"Element-2 initial overpotential [V]": 0,
"R1 [Ohm]": 0.005,
"R2 [Ohm]": 0.0003,
"C1 [F]": 10000,
"C2 [F]": 5000,
"Entropic change [V/K]": 0.0004,
}
model = pybop.empirical.Thevenin(
parameter_set=parameter_set, options={"number of rc elements": 2}
)
Ordinarily, we would want to fit models to experimental data. For simplicity we'll use our model to generate some synthetic data, by simulating a discharge and corrupting the results with random noise. We will then use this synthetic data to form a PyBOP Dataset, from which we can demonstrate the model fitting procedure.
sigma = 0.001
t_eval = np.linspace(0, 500, 500)
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,
}
)
Now for the fun part! We begin by defining the parameters that we want to fit. In this case, it's the series resistance $R_0$, and the RC parameters $R_1$ and $C_1$. Since we don't want to make things too easy, we've set up our parameters to have quite different values to those used in simulating the data.
parameters = pybop.Parameters(
pybop.Parameter(
"R0 [Ohm]",
prior=pybop.Gaussian(2e-3, 1e-4),
bounds=[1e-4, 1e-1],
),
pybop.Parameter(
"R1 [Ohm]",
prior=pybop.Gaussian(1e-3, 1e-4),
bounds=[1e-5, 1e-2],
),
pybop.Parameter(
"C1 [F]",
prior=pybop.Gaussian(5000, 300),
bounds=[2500, 5e4],
),
)
We can now set up a problem and cost for identifying these parameters from our synthetic dataset.
problem = pybop.FittingProblem(model, parameters, dataset)
cost = pybop.SumSquaredError(problem)
The voltage over any given RC branch will evolve over a timescale $\tau = R \times C$. The data we work with will place some limits on which parameters can realistically be identified. We can't expect to fit fast timescales with low sample-rate data; similarly, we can't get good parameters for a long timescale if we only have short amounts of data.
To illustrate, imagine trying to fit a timescale of $\tau=0.1$ s from data with a sample rate of 1 Hz. Virtually all the interesting dynamics will have happened over the course of a single sample, so there's not enough information to work from for parameter identification. Similarly, it would be unreasonable to fit a timescale of $\tau=1000$ s from only one minute of data; across the range of data that we have, the dynamics will have barely changed, giving us very little to work from for fitting $R$ and $C$.
In general therefore, we need to be careful to make sure our parameter values are sensible. To make sure that the optimiser doesn't propose excessively long or short timescales, we can use nonlinear constraints. The scipy.optimize.NonlinearConstraint function handles this for us.
tau_constraint = scipy.optimize.NonlinearConstraint(lambda x: x[1] * x[2], 5, 750)
Let's dig into what's going on here.
The nonlinear constraint will receive a parameter vector from PyBOP. We know from our parameters list that parameter 0 will be $R_0$,
parameter 1 will be $R_1$, and parameter 2 will be $C_1$. We can therefore compute the RC timescale $R_1 \times C_1$ using lambda x: x[1] * x[2]. Next, we tell scipy that we want this to be constrained to the range $5\leq R_1 \times C_1 \leq 750$. There's a lot of flexibility in our choice of lower and upper bounds, however these are reasonable numbers given our sampling rate (1 Hz), and data time-span (500 s).
Now all we need to do is set up an optimiser to use this.
optim = pybop.SciPyMinimize(
cost,
method="trust-constr",
constraints=tau_constraint,
max_iterations=100,
)
Note that COBYLA, COBYQA, and SLSQP can also be used as the method; these are all different Scipy optimisers with constraint capabilities.
Finally, we run the parameteriser, and plot some results. Don't worry if the solver sometimes terminates early -- this happens when the model receives a bad set of parameters, but PyBOP will catch and handle these cases automatically.
x, final_cost = optim.run()
print("Estimated parameters:", x)
# Plot the time series
pybop.plot_dataset(dataset)
# Plot the timeseries output
pybop.quick_plot(problem, problem_inputs=x, title="Optimised Comparison")
# Plot convergence
pybop.plot_convergence(optim)
# Plot the parameter traces
pybop.plot_parameters(optim);
--------------------------------------------------------------------------- KeyboardInterrupt Traceback (most recent call last) Cell In[9], line 1 ----> 1 x, final_cost = optim.run() 2 print("Estimated parameters:", x) 4 # Plot the time series File ~/Documents/Git/PyBOP/pybop/optimisers/base_optimiser.py:177, in BaseOptimiser.run(self) 166 def run(self): 167 """ 168 Run the optimisation and return the optimised parameters and final cost. 169 (...) 175 The final cost associated with the best parameters. 176 """ --> 177 self.result = self._run() 179 # Store the optimised parameters 180 x = self.result.x File ~/Documents/Git/PyBOP/pybop/optimisers/scipy_optimisers.py:71, in _run(self) 64 def _run(self): 65 """ 66 Internal method to run the optimization using a PyBOP optimiser. 67 68 Returns 69 ------- 70 result : pybop.Result ---> 71 The result of the optimisation including the optimised parameter values and cost. 72 """ 73 result = self._run_optimiser() 75 try: File ~/Documents/Git/PyBOP/pybop/optimisers/scipy_optimisers.py:229, in _run_optimiser(self) 222 if np.isinf(self._cost0): 223 raise ValueError( 224 "The initial parameter values return an infinite cost." 225 ) 227 return minimize( 228 self.cost_wrapper, --> 229 self.x0, 230 bounds=self._scipy_bounds, 231 callback=callback, 232 **self._options, 233 ) File _pydevd_bundle/pydevd_cython_darwin_312_64.pyx:1187, in _pydevd_bundle.pydevd_cython_darwin_312_64.SafeCallWrapper.__call__() File _pydevd_bundle/pydevd_cython_darwin_312_64.pyx:627, in _pydevd_bundle.pydevd_cython_darwin_312_64.PyDBFrame.trace_dispatch() File _pydevd_bundle/pydevd_cython_darwin_312_64.pyx:937, in _pydevd_bundle.pydevd_cython_darwin_312_64.PyDBFrame.trace_dispatch() File _pydevd_bundle/pydevd_cython_darwin_312_64.pyx:928, in _pydevd_bundle.pydevd_cython_darwin_312_64.PyDBFrame.trace_dispatch() File _pydevd_bundle/pydevd_cython_darwin_312_64.pyx:585, in _pydevd_bundle.pydevd_cython_darwin_312_64.PyDBFrame.do_wait_suspend() File /Applications/PyCharm.app/Contents/plugins/python/helpers/pydev/pydevd.py:1201, in PyDB.do_wait_suspend(self, thread, frame, event, arg, send_suspend_message, is_unhandled_exception) 1198 from_this_thread.append(frame_id) 1200 with self._threads_suspended_single_notification.notify_thread_suspended(thread_id, stop_reason): -> 1201 self._do_wait_suspend(thread, frame, event, arg, suspend_type, from_this_thread) File /Applications/PyCharm.app/Contents/plugins/python/helpers/pydev/pydevd.py:1216, in PyDB._do_wait_suspend(self, thread, frame, event, arg, suspend_type, from_this_thread) 1213 self._call_mpl_hook() 1215 self.process_internal_commands() -> 1216 time.sleep(0.01) 1218 self.cancel_async_evaluation(get_current_thread_id(thread), str(id(frame))) 1220 # process any stepping instructions KeyboardInterrupt:
Here, we have considered parameter estimation for an equivalent circuit model. We have introduced some issues of practical identifiability, whereby particularly fast or slow timescales can't be parameterised very well. To work around this challenge, nonlinear constraints have been implemented.
A difficulty here is that fitting RC parameters is a sloppy problem -- there's a range of timescales that will all produce reasonable looking results, and an optimiser will struggle to choose any one timescale over another. Gradient-free methods such as pybop.CMAES may do a better job of fitting timescales. Nevertheless, the resistance values can be fitted very accurately with gradient-based methods, whereas gradient-free methods don't always come up with such a good solution.
The ideal approach would be to fit resistances with gradient-based or ordinary least-squares methods, and timescales using gradient-free optimisers. In the absence of this approach, general recommendations are as follows...