{
"cells": [
{
"cell_type": "markdown",
"id": "0",
"metadata": {},
"source": [
"# Empirical models with scipy trust-constr:\n",
"## Identifying equivalent circuit parameters using bounds, nonlinear constraints, and gradient-based optimisers\n",
"\n",
"Here, we provide a short example of how to identify equivalent-circuit parameters using the [scipy trust-constr](https://docs.scipy.org/doc/scipy-1.14.0/reference/generated/scipy.optimize.minimize.html) 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.\n",
"\n",
"### Importing libraries\n",
"\n",
"If you don't already have PyBOP installed, check out the [installation guide](https://pybop-docs.readthedocs.io/en/latest/installation.html) first.\n",
"\n",
"We begin by importing the necessary libraries."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "1",
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import scipy.optimize\n",
"\n",
"import pybop\n",
"\n",
"pybop.plot.PlotlyManager().pio.renderers.default = \"notebook_connected\""
]
},
{
"cell_type": "markdown",
"id": "2",
"metadata": {},
"source": [
"### Initialising the model and parameters\n",
"\n",
"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.\n",
"\n",
"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```.\n",
"\n",
"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."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "3",
"metadata": {},
"outputs": [],
"source": [
"parameter_set = {\n",
" \"chemistry\": \"ecm\",\n",
" \"Initial SoC\": 0.5,\n",
" \"Initial temperature [K]\": 25 + 273.15,\n",
" \"Cell capacity [A.h]\": 5,\n",
" \"Nominal cell capacity [A.h]\": 5,\n",
" \"Ambient temperature [K]\": 25 + 273.15,\n",
" \"Current function [A]\": 5,\n",
" \"Upper voltage cut-off [V]\": 4.2,\n",
" \"Lower voltage cut-off [V]\": 3.0,\n",
" \"Cell thermal mass [J/K]\": 1000,\n",
" \"Cell-jig heat transfer coefficient [W/K]\": 10,\n",
" \"Jig thermal mass [J/K]\": 500,\n",
" \"Jig-air heat transfer coefficient [W/K]\": 10,\n",
" \"Open-circuit voltage [V]\": pybop.empirical.Thevenin().default_parameter_values[\n",
" \"Open-circuit voltage [V]\"\n",
" ],\n",
" \"R0 [Ohm]\": 0.01,\n",
" \"Element-1 initial overpotential [V]\": 0,\n",
" \"Element-2 initial overpotential [V]\": 0,\n",
" \"R1 [Ohm]\": 0.005,\n",
" \"R2 [Ohm]\": 0.0003,\n",
" \"C1 [F]\": 10000,\n",
" \"C2 [F]\": 5000,\n",
" \"Entropic change [V/K]\": 0.0004,\n",
"}\n",
"\n",
"model = pybop.empirical.Thevenin(\n",
" parameter_set=parameter_set, options={\"number of rc elements\": 2}\n",
")"
]
},
{
"cell_type": "markdown",
"id": "4",
"metadata": {},
"source": [
"### Generating synthetic data\n",
"\n",
"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."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "5",
"metadata": {},
"outputs": [],
"source": [
"sigma = 0.001\n",
"t_eval = np.linspace(0, 500, 500)\n",
"values = model.predict(t_eval=t_eval)\n",
"corrupt_values = values[\"Voltage [V]\"].data + np.random.normal(0, sigma, len(t_eval))\n",
"\n",
"# Form dataset\n",
"dataset = pybop.Dataset(\n",
" {\n",
" \"Time [s]\": t_eval,\n",
" \"Current function [A]\": values[\"Current [A]\"].data,\n",
" \"Voltage [V]\": corrupt_values,\n",
" }\n",
")"
]
},
{
"cell_type": "markdown",
"id": "6",
"metadata": {},
"source": [
"### Identifying parameters\n",
"\n",
"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."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "7",
"metadata": {},
"outputs": [],
"source": [
"parameters = pybop.Parameters(\n",
" pybop.Parameter(\n",
" \"R0 [Ohm]\",\n",
" prior=pybop.Gaussian(2e-3, 1e-4),\n",
" bounds=[1e-4, 1e-1],\n",
" ),\n",
" pybop.Parameter(\n",
" \"R1 [Ohm]\",\n",
" prior=pybop.Gaussian(1e-3, 1e-4),\n",
" bounds=[1e-5, 1e-2],\n",
" ),\n",
" pybop.Parameter(\n",
" \"C1 [F]\",\n",
" prior=pybop.Gaussian(5000, 300),\n",
" bounds=[2500, 5e4],\n",
" ),\n",
")"
]
},
{
"cell_type": "markdown",
"id": "8",
"metadata": {},
"source": [
"We can now set up a problem and cost for identifying these parameters from our synthetic dataset."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "9",
"metadata": {},
"outputs": [],
"source": [
"problem = pybop.FittingProblem(model, parameters, dataset)\n",
"cost = pybop.SumSquaredError(problem)"
]
},
{
"cell_type": "markdown",
"id": "10",
"metadata": {},
"source": [
"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. \n",
"\n",
"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$.\n",
"\n",
"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."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "11",
"metadata": {},
"outputs": [],
"source": [
"tau_constraint = scipy.optimize.NonlinearConstraint(lambda x: x[1] * x[2], 5, 750)"
]
},
{
"cell_type": "markdown",
"id": "12",
"metadata": {},
"source": [
"Let's dig into what's going on here.\n",
"\n",
"The nonlinear constraint will receive a parameter vector from PyBOP. We know from our parameters list that parameter 0 will be $R_0$, \n",
"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).\n",
"\n",
"Now all we need to do is set up an optimiser to use this."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "13",
"metadata": {},
"outputs": [],
"source": [
"optim = pybop.SciPyMinimize(\n",
" cost,\n",
" method=\"trust-constr\",\n",
" constraints=tau_constraint,\n",
" max_iterations=100,\n",
")"
]
},
{
"cell_type": "markdown",
"id": "14",
"metadata": {},
"source": [
"Note that ```COBYLA```, ```COBYQA```, and ```SLSQP``` can also be used as the method; these are all different Scipy optimisers with constraint capabilities.\n",
"\n",
"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."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "15",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"OptimisationResult:\n",
" Initial parameters: [2.01470532e-03 1.05171055e-03 4.57653422e+03]\n",
" Optimised parameters: [8.41410880e-03 6.45916310e-03 4.57841684e+03]\n",
" Final cost: 0.0013884834785099927\n",
" Optimisation time: 8.70180082321167 seconds\n",
" Number of iterations: 100\n",
" SciPy result available: Yes\n"
]
},
{
"data": {
"text/html": [
" \n",
" "
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/html": [
""
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/html": [
""
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/html": [
""
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/html": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"results = optim.run()\n",
"\n",
"# Plot the time series\n",
"pybop.plot.dataset(dataset)\n",
"\n",
"# Plot the timeseries output\n",
"pybop.plot.quick(problem, problem_inputs=results.x, title=\"Optimised Comparison\")\n",
"\n",
"# Plot convergence\n",
"pybop.plot.convergence(optim)\n",
"# Plot the parameter traces\n",
"pybop.plot.parameters(optim);"
]
},
{
"cell_type": "markdown",
"id": "16",
"metadata": {},
"source": [
"### Conclusion\n",
"\n",
"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.\n",
"\n",
"A difficulty here is that fitting RC parameters is a [sloppy problem](https://sethna.lassp.cornell.edu/Sloppy/FittingExponentials.html) -- 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.\n",
"\n",
"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...\n",
"\n",
"+ If good initial guesses are available for the RC timescales, use a gradient-based method in the way we've done here, as it'll typically give the most accurate fit;\n",
"+ Otherwise, use a gradient-free method to estimate timescales, then maybe consider moving back to a gradient-based method!"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.12.2"
}
},
"nbformat": 4,
"nbformat_minor": 5
}