bmlite

Summary

Battery Analysis and Training Models for Optimization and Degradation Studies (BATMODS) is a Python package with an API for pre-built battery models. The original purpose of the package was to quickly generate synthetic data for machine learning models to train with. However, the models are generally useful for any battery simulation or analysis. BATMODS-lite includes the following:

1. A library and API for pre-built battery models

2. Kinetic/transport properties for common battery materials

Accessing the documentation

Documentation is accessible via Python’s help() function which prints docstrings from the specified package, module, function, class, etc. You can also access the documentation by visiting the website, hosted on Read the Docs. The website includes search functionality and more detailed examples.

Submodules

• bmlite.P2D
• bmlite.SPM
• bmlite.materials
• bmlite.mathutils
• bmlite.mesh
• bmlite.plotutils

Classes

Constants	Physical constants.
Experiment	Experiment builder.
IDAResult	Results container.
IDASolver	SUNDIALS IDA solver.

Functions

templates(model[, file])

Print simulation templates.

Package Contents

class bmlite.Constants

Physical constants.

property F: float

Faraday’s constant [C/kmol].

property R: float

Gas constant [J/kmol/K].

class bmlite.Experiment(**kwargs)

Experiment builder.

A class to define an experimental protocol. Use the add_step() method to add a series of sequential steps. Each step defines a control mode, a constant or time-dependent load profile, a time span, and optional limiting criteria to stop the step early if a specified event/state is detected.

Parameters:

**kwargs (dict, optional) – IDASolver keyword arguments that span all steps.

See also

IDASolver

The solver class, with documentation for most keyword arguments that you might want to adjust.

add_step(mode, value, tspan, limits=None, **kwargs)

Add a step to the experiment.

Parameters:

• mode (str) – Control mode, {‘current_A’, ‘current_C’, ‘voltage_V’, ‘power_W’}.

• value (float | Callable) – Value of boundary contion mode, in the appropriate units. Note that negative and positive values of current and power reference charge and discharge directions, respectively.

• tspan (float or tuple[float, float] or 1D array) – Relative times for recording solution [s]. Providing a float will result in the solver picking time steps to save on its own. A tuple is interpreted as (tmax, dt) where the first element is the max time for the step (in seconds) and the second is the time interval between steps (also seconds). You can also provide any custom array of times at which to save the solution by providing a 1D np.array; however, the first element must be zero and the array must be in a monotonically increasing order, and there must be at least three elements. An array like np.array([0, tmax]) will will result in the solver choosing its own time steps, similar to just providing a float. See notes for more information.

• limits (tuple[str, float], optional) – Stopping criteria for the new step, must be entered in sequential name/value pairs. Allowable names are {‘current_A’, ‘current_C’, ‘voltage_V’, ‘power_W’, ‘capacity_Ah’,’time_s’, ‘time_min’, ‘time_h’}. Multiple limits are allowed by entering consecutive pairs of names and values. Capacity limits track the throughput of the step and are calculated by integrating current over time. The time limits are in reference to total experiment time. Step times are controlled using the ‘tspan’ argument instead of the ‘limits’ input. The default is None. Current and power limits should follow the same sign convention as the mode, i.e., negative values for charge and positive for discharge.

• **kwargs (dict, optional) – IDASolver keyword arguments specific to the new step only.

Raises:

• ValueError – ‘mode’ is invalid.

• ValueError – A ‘limits’ name is invalid.

• TypeError – ‘tspan’ must be type float, tuple, or np.array.

• ValueError – ‘tspan’ tuple must be length 2.

• TypeError – ‘tspan’ tuple values must be type float.

• ValueError – ‘tspan[1]’ must be less than ‘tspan[0]’ when given a tuple.

• ValueError – ‘tspan’ arrays must be one-dimensional.

• ValueError – ‘tspan[0]’ must be zero when given an array.

• ValueError – ‘tspan’ arrays must be monotonically increasing.

See also

IDASolver

The solver class, with documentation for most keyword arguments that you might want to adjust.

Notes

For time-dependent loads, use a Callable for value with a function signature like def load(t: float) -> float, where t is the step’s relative time, in seconds.

When tspan is given as a 2-tuple, like (tmax, dt), the time span is constructed as:

tspan = np.arange(0., tspan[0], tspan[1])

In the case where tmax is not an integer multiple of dt, a final time point is appended to ensure that tspan[-1] == tmax. If this is too restrictive, you can instead provide a custom 1D np.array for the tspan argument. However, the array is checked to make sure the first element is zero and the array is monotonically increasing. If either of these checks fail, a ValueError is raised.

print_steps()

Print a formatted/readable list of steps.

property num_steps: int

Return number of steps.

Returns:

num_steps (int) – Number of steps.

property steps: list[dict]

Return steps list.

Returns:

steps (list[dict]) – List of the step dictionaries.

class bmlite.IDAResult(**kwargs)

Results container.

Inherits from RichResult. The solution class groups output from IDA into an object with the fields:

Parameters:

• message (str) – Human-readable description of the status value.

• success (bool) – True if the solver was successful (status >= 0). False otherwise.

• status (int) – Reason for the algorithm termination. Negative values correspond to errors, and non-negative values to different successful criteria.

• t (ndarray, shape(n,)) – Solution time(s). The dimension depends on the method. Stepwise solutions will only have 1 value whereas solutions across a full ‘tspan’ will have many.

• y (ndarray, shape(n, m)) – State variable values at each solution time. Rows correspond to indices in ‘t’ and columns match indexing from ‘y0’.

• yp (ndarray, shape(n, m)) – State variable time derivate values at each solution time. Row and column indexing matches ‘y’.

• i_events (ndarray, shape(k, num_events) or None) –

  Provides an array for each detected event ‘k’ specifying indices for which event(s) occurred. i_events[k,i] != 0 if ‘events[i]’ occurred at ‘t_events[k]’. The sign of ‘i_events’ indicates the direction of zero-crossing:

  • -1 indicates ‘events[i]’ was decreasing

  • +1 indicates ‘events[i]’ was increasing

  Output for ‘i_events’ will be None when either ‘eventsfn’ was None or if no events occurred during the solve.

• t_events (ndarray, shape(k,) or None) – Times at which events occurred or None if ‘eventsfn’ was None or no events were triggered during the solve.

• y_events (ndarray, shape(k, m) or None) – State variable values at each ‘t_events’ value or None. Rows and columns correspond to ‘t_events’ and ‘y0’ indexing, respectively.

• yp_events (ndarray, shape(k, m) or None) – State variable time derivative values at each ‘t_events’ value or None. Row and column indexing matches ‘y_events’.

• nfev (int) – Number of times that ‘resfn’ was evaluated.

• njev (int) – Number of times the Jacobian was evaluated, ‘jacfn’ or internal finite difference method.

Notes

Terminal events are appended to the end of ‘t’, ‘y’, and ‘yp’. However, if an event was not terminal then it will only appear in *_events outputs and not within the main output arrays.

‘nfev’ and ‘njev’ are cumulative for stepwise solution approaches. The values are reset each time ‘init_step’ is called.

class bmlite.IDASolver(resfn, **options)

SUNDIALS IDA solver.

A class to wrap the implicit differential algebraic (IDA) solver from SUNDIALS [1] [2]. IDA solves both ordinary differential equations (ODEs) and differiential agebraic equations (DAEs).

Parameters:

• resfn (Callable) – Residual function with signature f(t, y, yp, res[, userdata]). See the notes for more information.

• **options (dict, optional) – Keyword arguments to describe the solver options. A full list of names, types, descriptions, and defaults is given below.

• userdata (object or None, optional) – Additional data object to supply to all user-defined callables. Cannot be None (default) if ‘resfn’ takes in 5 arguments.

• calc_initcond ({'y0', 'yp0', None}, optional) – Determines whether or not ‘y0’ or ‘yp0’ are corrected before the first step. Requires ‘calc_init_dt’ if not None (default).

• calc_init_dt (float, optional) – Step size for initial condition correction. Positive for forward, negative for backward integration. Default is 0.01.

• algebraic_idx (array_like[int] or None, optional) – Indices ‘i’ for ‘y[i]’ that are for purely algebraic variables. If None (default) the problem is assumed to be a pure ODE.

• first_step (float, optional) – The initial step size. The default is 0, which uses an estimated value internally determined by SUNDIALS.

• min_step (float, optional) – Minimum allowable step size. The default is 0.

• max_step (float, optional) – Maximum allowable step size. Use 0 (default) for unbounded steps.

• rtol (float, optional) – Relative tolerance. It is recommended to not use values larger than 1e-3 or smaller than 1e-15. The default is 1e-5.

• atol (float or array_like[float], optional) – Absolute tolerance. A scalar will apply to all variables equally, while an array (matching ‘y’ length) sets specific tolerances for eqch variable. The default is 1e-6.

• linsolver ({'dense', 'band', 'sparse', ...}, optional) – Choice of linear solver, defaults to ‘dense’. ‘band’ requires both ‘lband’ and ‘uband’. ‘sparse’ uses SuperLU_MT [3] and requires ‘sparsity’. When using an iterative method (‘gmres’, ‘bicgstab’, ‘tfqmr’) the number of Krylov dimensions is set using ‘krylov_dim’. ‘lapackdense’ and ‘lapackband’ can also be used as alternatives to ‘dense’ and ‘band’. They use OpenBLAS-linked LAPACK [4] routines, but can have noticeable overhead for small (<100) systems.

• lband (int or None, optional) – Lower Jacobian bandwidth. Given a DAE system 0 = F(t, y, yp), the Jacobian is J = dF_i/dy_j + cj*dF_i/dyp_j. Required when ‘linsolver’ is ‘band’. Use zero if no values are below the main diagonal. Defaults to None.

• uband (int or None, optional) – Upper Jacobian bandwidth. Required when ‘linsolver’ is ‘band’. Use zero if no elements are above the main diagonal. Defaults to None.

• sparsity (2D np.array or sparse matrix or None, optional) – Jacobian sparsity pattern. Required when ‘linsolver’ is ‘sparse’. The shape must be (N, N) where N is the size of the system. Zero entries indicate fixed zeros in the Jacobian. If ‘jacfn’ is None, this argument activates a custom Jacobian routine (not part of the original SUNDIALS package). The routine works with all direct linear solvers but may increase step count. Reduce ‘max_step’ to help with this, if needed. Defaults to None.

• nthreads (int or None, optional) – Number of threads to use with the ‘sparse’ linear solver. If None (default), 1 is used. Use -1 to use all available threads.

• krylov_dim (int or None, optional) – Maximum number of Krylov basis vectors for iterative solvers. Will default to 5 if invalid/None when required. Larger values improve convergence but increase memory usage. Only applies to the ‘gmres’, ‘bicgstab’, and ‘tfqmr’ linear solvers.

• max_order (int, optional) – Specifies the maximum order for the linear multistep BDF method. The value must be in the range [1, 5]. The default is 5.

• max_num_steps (int, optional) – The maximum number of steps taken by the solver in each attempt to reach the next output time. The default is 500.

• max_nonlin_iters (int, optional) – Specifies the maximum number of nonlinear solver iterations in one step. The default is 4.

• max_conv_fails (int, optional) – Specifies the max number of nonlinear solver convergence failures in one step. The default is 10.

• constraints_idx (array_like[int] or None, optional) – Specifies indices ‘i’ in the ‘y’ state variable array for which inequality constraints should be applied. Constraint types must be specified in ‘constraints_type’, see below. The default is None.

• constraints_type (array_like[int] or None, optional) – If ‘constraints_idx’ is not None, then this option must include an array of equal length specifying the types of constraints to apply. Values should be in {-2, -1, 1, 2} which apply y[i] < 0, y[i] <= 0, y[i] >= 0, and y[i] > 0, respectively. The default is None.

• eventsfn (Callable or None, optional) –

  Events function with signature g(t, y, yp, events[, userdata]). If None (default), no events are tracked. See the notes for more information. Requires ‘num_events’ be set when not None.

  The function may also have these optional attributes:

  terminal: list[bool, int], optional

  Specifies solver behavior for each event. A boolean stops the solver (True) or just records the event (False). An integer stops the solver after than many occurrences. The default is [True]*num_events.

  direction: list[int], optional

  Determines which event slopes to track: -1 (negative), +1 (positive), or 0 (both). If not provided the default [0]*num_events is used.

  You can assign attributes like eventsfn.terminal = [True] to any function in Python, after it has been defined.

• num_events (int, optional) – Number of events to track. The default is 0.

• jacfn (Callable or None, optional) – Jacobian function like J(t, y, yp, res, cj, JJ[, userdata]). The function should fill the pre-allocated memory for JJ with the values defined by JJ[i,j] = dres_i/dy_j + cj*dres_i/dyp_j. An internal finite difference method is applied when None (default). Note that the template for JJ is determined by the linear solver choice: a 1D array for ‘sparse’ or 2D array for ‘dense’ and ‘band’. See the notes for more information.

• precond (IDAPrecond or None, optional) – Preconditioner functions. Only compatible with iterative linear solvers. Must be an instance of IDAPrecond if not None (default).

• jactimes (IDAJacTimes or None, optional) – Jacobian-vector product functions. Only compatible with iterative linear solvers. Must be an instance of IDAJacTimes when provided. Difference quotient approximations are used with iterative solvers if None (default).

Notes

Return values from all user-defined function (e.g., ‘resfn’, ‘eventsfn’, and ‘jacfn’) are ignored by the solver. Instead the solver directly reads from pre-allocated memory. Output arrays (e.g., ‘res’, ‘events’, and ‘JJ’) from each user-defined callable should be filled within each respective function. When setting values across the entire array/matrix at once, don’t forget to use [:] to fill the existing array rather than overwriting it. For example, using res[:] = F(t, y, yp) is correct whereas res = F(t, y, yp) is not.

When any user-defined function require data outside of their normal arguments, you can supply optional ‘userdata’. When given, ‘userdata’ must appear in ALL function signatures (‘rhsfn’, ‘eventsfn’, ‘jacfn’, ‘precond’, and ‘jactimes’) even if it is not used in all functions. Of course this does not apply to functions that are provided as None. Note that ‘userdata’ only takes up one argument position; however, ‘userdata’ can be any Python object. Therefore, to pass more than one extra argument you should pack all data into a single tuple, dict, dataclass, etc. and pass them all together as ‘userdata’. The data can be unpacked as needed within the functions.

When using user-denfined Jacobians supplied via ‘jacfn’, the allocated memory for JJ depends on the chosen linear solver. For the ‘dense’ and ‘band’ solver options, JJ is a normal 2D array with shape (N, N), where N is the size of the system. This means that the template takes up more memory than is actually necessary for the banded solver. This is a choice of convenience, however, to avoid requiring the user to deal with complex indexing in compressed formats. There is a plan for future releases to allow users to decide between this convenient, but less memory-efficient template and a more memory-efficient 1D template that only tracks entries within the bandwidth. While users currently have access to write to indexes outside the bandwidth in the present implementation for the banded solver, any values written outside the given bandwidth are internally ignored by the solver. The JJ template details for the ‘sparse’ solver option do not follow the same 2D format, and are discussed below.

When users chose to provide their own Jacobian via ‘jacfn’ and specify use of the ‘sparse’ linear solver, the memory allocated for JJ is a 1D array with length equal to the number of non-zero entries in the supplied sparsity pattern. This is a choice made for memory efficiency, but does require users to understand the indexing of the JJ entries. The exposed 1D array is in compressed sparse column (CSC) format, which means that the non-zero entries are ordered first by column and then by row. A minimal example is provided below for a 2x2 Jacobian with three non-zero entries.

import numpy as np
import sksundae as sun

def resfn(t, y, yp, res):
    # Simple DAE residual expressions.
    res[0] = yp[0] + y[0]
    res[1] = y[0] - y[1]

def jacfn(t, y, yp, res, cj, JJ):
    # Analytical Jacobian. JJ is a 1D array of length 3 due to only
    # three non-zero entries. Values are indexed in a CSC format.
    JJ[0] = 1.0 + cj  # dres0/dy0 + cj*dres0/dyp0
    JJ[1] = 1.0       # dres1/dy0 + cj*dres1/dyp0
    JJ[2] = -1.0      # dres1/dy1 + cj*dres1/dyp1

sparsity = np.array([[1, 0], [1, 1]])  # Jacobian sparsity pattern

solver = sun.ida.IDA(
    resfn,
    jacfn=jacfn,
    algebraic_idx=[1],  # y[1] is purely algebraic
    sparsity=sparsity,
    linsolver='sparse',
)

soln = solver.solve([0, 1], [1.0, 1.0], [0.0, 0.0])

References

[1]

A. C. Hindmarsh, P. N. Brown, K. E. Grant, S. L. Lee, R. Serban, D. E. Shumaker, and C. S. Woodward, “SUNDIALS: Suite of Nonlinear and Differential/Algebraic Equation Solvers,” ACM TOMS, 2005, DOI: 10.1145/1089014.1089020

[2]

D. J. Gardner, D. R. Reynolds, C. S. Woodward, C. J. Balos, “Enabling new flexibility in the SUNDIALS suite of nonlinear and differential/algebraic equation solvers,” ACM TOMS, 2022, DOI: 10.1145/3539801

[3]

J. W. Demmel, J. R. Gilbert, and X. S. Li, “An Asynchronous Parallel Supernodal Algorithm for Sparse Gaussian Elimination,” SIMAX, 1999, DOI: 10.1137/S0895479897317685

[4]

E. Anderson, Z. Bai, C. Bischof, S. Blackford, J. Demmel, J. Dongarra, J. Du Croz, A. Greenbaum, S. Hammarling, A. McKenney, D. Sorensen, “LAPACK Users’ Guide,” Society for Industrial and Applied Mathematics, 1999, Philidelphia, PA.

Examples

The following example solves the Robertson problem, which is a classic test problem for programs that solve stiff ODEs. A full description of the problem is provided by MATLAB. While initializing the solver, algebraic_idx=[2] specifies y[2] is purely algebraic, and calc_initcond='yp0' tells the solver to determine the values for ‘yp0’ at ‘tspan[0]’ before starting to integrate. That is why ‘yp0’ can be initialized as an array of zeros even though plugging in ‘y0’ to the residuals expressions actually gives yp0 = [-0.04, 0.04, 0]. The initialization is checked against the correct answer after solving.

import numpy as np
import sksundae as sun
import matplotlib.pyplot as plt

def resfn(t, y, yp, res):
    res[0] = yp[0] + 0.04*y[0] - 1e4*y[1]*y[2]
    res[1] = yp[1] - 0.04*y[0] + 1e4*y[1]*y[2] + 3e7*y[1]**2
    res[2] = y[0] + y[1] + y[2] - 1.0

solver = sun.ida.IDA(resfn, algebraic_idx=[2], calc_initcond='yp0')

tspan = np.hstack([0, 4*np.logspace(-6, 6)])
y0 = np.array([1, 0, 0])
yp0 = np.zeros_like(y0)

soln = solver.solve(tspan, y0, yp0)
assert np.allclose(soln.yp[0], [-0.04, 0.04, 0], rtol=1e-3)

soln.y[:, 1] *= 1e4  # scale y[1] so it is visible in the figure
plt.semilogx(soln.t, soln.y)
plt.show()

bmlite.templates(model, file=None)

Print simulation templates.

If no input, a list of available templates will print. Otherwise, given a file name or index, a template will print.

Parameters:

• model (str) – Model package name, e.g., ‘SPM’. Input is case insensitive. All model names are converted to uppercase within this function.

• file (str | int, optional) – File name or index. The default is None.

Raises:

• AttributeError – ‘model’ is not a valid subpackage.

• FileNotFoundError – ‘model’ has no ‘templates’ directory.