I tried to populate different initial electrolyte concentration across the cells in the pack. I have an array of initial electrolyte concentration (initial_eleConc) similar to total heat transfer coefficient (htc). The changes I made to utils.py, simulations.py, solver_utils.py are given below.
Changed line 220, and added line 229 in utils.py.
# -*- coding: utf-8 -*-
"""
Created on Thu Sep 23 10:33:13 2021
@author: Tom
"""
from scipy.interpolate import interp1d, interp2d
import numpy as np
import pandas as pd
import os
import liionpack
from skspatial.objects import Plane
from skspatial.objects import Points
ROOT_DIR = os.path.dirname(os.path.abspath(liionpack.__file__))
CIRCUIT_DIR = os.path.join(ROOT_DIR, "circuits")
DATA_DIR = os.path.join(ROOT_DIR, "data")
INIT_FUNCS = os.path.join(ROOT_DIR, "init_funcs")
def interp_current(df):
r"""
Returns an interpolation function for current w.r.t time
Parameters
----------
df : pandas.DataFrame or Dict
Contains data for 'Time' and 'Cells Total Current' from which to
construct an interpolant function
Returns
-------
f : function
interpolant function of total cell current with time.
"""
t = df["Time"]
I = df["Cells Total Current"]
f = interp1d(t, I)
return f
def _z_from_plane(X, Y, plane):
r"""
Given X and Y and a plane provide Z
X - temperature
Y - flow rate
Z - heat transfer coefficient
Parameters
----------
X : float (array)
x-coordinate.
Y : float (array)
z-coordinate.
plane : skspatial.object.Plane
plane returned from read_cfd_data.
Returns
-------
z : float (array)
z-coordinate.
"""
a, b, c = plane.normal
d = plane.point.dot(plane.normal)
z = (d - a * X - b * Y) / c
return z
def _fit_plane(xv, yv, dbatt):
r"""
Private method to fit plane to CFD data
Parameters
----------
xv : ndarray
temperature meshgrid points.
yv : ndarray
flow_rate meshgrid points.
dbatt : ndarray
cfd data for heat transfer coefficient.
Returns
-------
plane : skspatial.object.Plane
htc varies linearly with temperature and flow rate so relationship
describes a plane
"""
nx, ny = xv.shape
pts = []
for i in range(nx):
for j in range(ny):
pts.append([xv[i, j], yv[i, j], dbatt[i, j]])
points = Points(pts)
plane = Plane.best_fit(points, tol=1e-6)
return plane
def read_cfd_data(data_dir=None, filename="cfd_data.xlsx", fit="linear"):
r"""
A very bespoke function to read heat transfer coefficients from an excel
file
Parameters
----------
data_dir : str, optional
Path to data file. The default is None. If unspecified the module
liionpack.DATA_DIR folder will be used
filename : str, optional
The default is 'cfd_data.xlsx'.
fit : str
options are 'linear' (default) and 'interpolated'.
Returns
-------
funcs : list
an interpolant is returned for each cell in the excel file.
"""
if data_dir is None:
data_dir = liionpack.DATA_DIR
fpath = os.path.join(data_dir, filename)
ncells = 32
flow_bps = np.array(pd.read_excel(fpath, sheet_name="massflow_bps", header=None))
temp_bps = np.array(pd.read_excel(fpath, sheet_name="temperature_bps", header=None))
xv, yv = np.meshgrid(temp_bps, flow_bps)
data = np.zeros([len(temp_bps), len(flow_bps), ncells])
fits = []
for i in range(ncells):
data[:, :, i] = np.array(
pd.read_excel(fpath, sheet_name="cell" + str(i + 1), header=None)
)
# funcs.append(interp2d(xv, yv, data[:, :, i], kind='linear'))
if fit == "linear":
fits.append(_fit_plane(xv, yv, data[:, :, i]))
elif fit == "interpolated":
fits.append(interp2d(xv, yv, data[:, :, i], kind="linear"))
return data, xv, yv, fits
def get_linear_htc(planes, T, Q):
r"""
A very bespoke function that is called in the solve process to update the
heat transfer coefficients for every battery - assuming linear relation
between temperature, flow rate and heat transfer coefficient.
Parameters
----------
planes : list
each element of the list is a plane equation describing linear relation
between temperature, flow rate and heat transfer coefficient.
T : float array
The temperature of each battery.
Q : float
The flow rate for the system.
Returns
-------
htc : float
Heat transfer coefficient for each battery.
"""
ncell = len(T)
htc = np.zeros(ncell)
for i in range(ncell):
htc[i] = _z_from_plane(T[i], Q, planes[i])
return htc
def get_interpolated_htc(funcs, T, Q):
r"""
A very bespoke function that is called in the solve process to update the
heat transfer coefficients for every battery
Parameters
----------
funcs : list
each element of the list is an interpolant function.
T : float array
The temperature of each battery.
Q : float
The flow rate for the system.
Returns
-------
htc : float
Heat transfer coefficient for each battery.
"""
ncell = len(T)
htc = np.zeros(ncell)
for i in range(ncell):
htc[i] = funcs[i](T[i], Q)
return htc
def build_inputs_dict(I_batt, htc,initial_eleConc):
r"""
Function to convert inputs and external_variable arrays to list of dicts
As expected by the casadi solver. These are then converted back for mapped
solving but stored individually on each returned solution.
Can probably remove this process later
Parameters
----------
I_batt : float array
The input current for each battery.
htc : float array
the heat transfer coefficient for each battery.
Returns
-------
inputs_dict : list
each element of the list is an inputs dictionary corresponding to each
battery.
"""
inputs_dict = []
for i in range(len(I_batt)):
inputs_dict.append(
{
# 'Volume-averaged cell temperature': T_batt[i],
"Current": I_batt[i],
"Total heat transfer coefficient [W.m-2.K-1]": htc[i],
"Initial concentration in electrolyte [mol.m-3]": initial_eleConc[i],
}
)
return inputs_dict
Added line 73 in simulations.py.
# -*- coding: utf-8 -*-
"""
Created on Wed Sep 22 15:37:51 2021
@author: tom
"""
import pybamm
def _current_function(t):
r"""
Internal function to make current an input parameter
Parameters
----------
t : float
Dummy time parameter.
Returns
-------
TYPE
DESCRIPTION.
"""
return pybamm.InputParameter("Current")
def create_simulation(parameter_values=None, experiment=None, make_inputs=False):
r"""
Create a PyBaMM simulation set up for interation with liionpack
Parameters
----------
parameter_values : pybamm.ParameterValues class
DESCRIPTION. The default is None.
experiment : pybamm.Experiment class
DESCRIPTION. The default is None.
make_inputs : bool, optional
Changes "Current function [A]" and "Total heat transfer coefficient
[W.m-2.K-1]" to be inputs that are controlled by liionpack.
The default is False.
Returns
-------
sim : pybamm.Simulation
A simulation that can be solved individually or passed into the
liionpack solve method
"""
# Create the pybamm model
model = pybamm.lithium_ion.SPMe(
options={
"thermal": "lumped",
}
)
# geometry = model.default_geometry
if parameter_values is None:
# load parameter values and process model and geometry
chemistry = pybamm.parameter_sets.Chen2020
parameter_values = pybamm.ParameterValues(chemistry=chemistry)
# Change the current function to be an input as this is set by the external circuit
if make_inputs:
parameter_values.update(
{
"Current function [A]": _current_function,
}
)
parameter_values.update(
{
"Current": "[input]",
"Total heat transfer coefficient [W.m-2.K-1]": "[input]",
"Initial concentration in electrolyte [mol.m-3]":"[input]",
},
check_already_exists=False,
)
solver = pybamm.CasadiSolver(mode="safe")
sim = pybamm.Simulation(
model=model,
experiment=experiment,
parameter_values=parameter_values,
solver=solver,
)
return sim
if __name__ == "__main__":
sim = create_simulation()
sim.solve([0, 1800])
sim.plot()
Changed lines 86, 122,163,256,266 in solver_utils.py.
# -*- coding: utf-8 -*-
"""
Created on Thu Sep 23 10:44:31 2021
@author: Tom
"""
import casadi
import pybamm
import numpy as np
import time as ticker
import liionpack as lp
def _mapped_step(model, solutions, inputs_dict, integrator, variables, t_eval):
r"""
Internal function to process the model for one timestep in a mapped way.
Mapped versions of the integrator and variables functions should already
have been made.
Parameters
----------
model : pybamm.Model
The built model
solutions : list of pybamm.Solution objects for each battery
Used to get the last state of the system and use as x0 and z0 for the
casadi integrator
inputs_dict : list of inputs_dict objects for each battery
DESCRIPTION.
integrator : mapped casadi.integrator
Produced by _create_casadi_objects
variables : mapped variables evaluator
Produced by _create_casadi_objects
t_eval : float array of times to evaluate
Produced by _create_casadi_objects
Returns
-------
sol : list
solutions that have been stepped forward by one timestep
var_eval : list
evaluated variables for final state of system
"""
len_rhs = model.concatenated_rhs.size
N = len(solutions)
if solutions[0] is None:
# First pass
x0 = casadi.horzcat(*[model.y0[:len_rhs] for i in range(N)])
z0 = casadi.horzcat(*[model.y0[len_rhs:] for i in range(N)])
else:
x0 = casadi.horzcat(*[sol.y[:len_rhs, -1] for sol in solutions])
z0 = casadi.horzcat(*[sol.y[len_rhs:, -1] for sol in solutions])
# t_min = [0.0]*N
t_min = 0.0
inputs = []
for temp in inputs_dict:
inputs.append(casadi.vertcat(*[x for x in temp.values()] + [t_min]))
ninputs = len(temp.values())
inputs = casadi.horzcat(*inputs)
# p = casadi.horzcat(*zip(inputs, external_variables, [t_min]*N))
# inputs_with_tmin = casadi.vertcat(inputs, np.asarray(t_min))
# Call the integrator once, with the grid
timer = pybamm.Timer()
tic = timer.time()
casadi_sol = integrator(x0=x0, z0=z0, p=inputs)
integration_time = timer.time()
nt = len(t_eval)
xf = casadi_sol["xf"]
# zf = casadi_sol["zf"]
sol = []
xend = []
for i in range(N):
start = i * nt
y_sol = xf[:, start:start + nt]
xend.append(y_sol[:, -1])
# Not sure how to index into zf - need an example
sol.append(pybamm.Solution(t_eval, y_sol, model, inputs_dict[i]))
sol[-1].integration_time = integration_time
toc = timer.time()
lp.logger.debug(f"Mapped step completed in {toc - tic}")
xend = casadi.horzcat(*xend)
var_eval = variables(0, xend, 0, inputs[0:ninputs, :])
return sol, var_eval
def _create_casadi_objects(I_init, htc,initial_eleConc, sim, dt, Nspm, nproc, variable_names):
r"""
Internal function to produce the casadi objects in their mapped form for
parallel evaluation
Parameters
----------
I_init : float
initial guess for current of a battery (not used for simulation).
htc : float
initial guess for htc of a battery (not used for simulation).
sim : pybamm.Simulation
A PyBaMM simulation object that contains the model, parameter_values,
solver, solution etc.
dt : float
The time interval for a single timestep. Fixed throughout the simulation
Nspm : int
Number of individual batteries in the pack.
nproc : int
Number of parallel processes to map to.
variable_names : list
Variables to evaluate during solve. Must be a valid key in the
model.variables
Returns
-------
integrator : mapped casadi.integrator
Solves an initial value problem (IVP) coupled to a terminal value
problem with differential equation given as an implicit ODE coupled
to an algebraic equation and a set of quadratures
variables_fn : mapped variables evaluator
evaluates the simulation and output variables. see casadi function
t_eval : float array of times to evaluate
times to evaluate in a single step, starting at zero for each step
"""
inputs = {"Current": I_init, "Total heat transfer coefficient [W.m-2.K-1]": htc,"Initial concentration in electrolyte [mol.m-3]": initial_eleConc}
solver = sim.solver
# solve model for 1 second to initialise the circuit
t_eval = np.linspace(0, 1, 2)
# Initial solution - this builds the model behind the scenes
sim.solve(t_eval, inputs=inputs)
# step model
# Code to create mapped integrator
t_eval = np.linspace(0, dt, 11)
t_eval_ndim = t_eval / sim.model.timescale.evaluate()
inp_and_ext = inputs
# No external variables - Temperature solved as lumped model in pybamm
# External variables could (and should) be used if battery thermal problem
# Includes conduction with any other circuits or neighboring batteries
# inp_and_ext.update(external_variables)
integrator = solver.create_integrator(
sim.built_model, inputs=inp_and_ext, t_eval=t_eval_ndim
)
integrator = integrator.map(Nspm, "thread", nproc)
# Variables function for parallel evaluation
casadi_objs = sim.built_model.export_casadi_objects(variable_names=variable_names)
variables = casadi_objs["variables"]
t, x, z, p = (
casadi_objs["t"],
casadi_objs["x"],
casadi_objs["z"],
casadi_objs["inputs"],
)
variables_stacked = casadi.vertcat(*variables.values())
variables_fn = casadi.Function("variables", [t, x, z, p], [variables_stacked])
variables_fn = variables_fn.map(Nspm, "thread", nproc)
return integrator, variables_fn, t_eval
def solve(
netlist=None,
parameter_values=None,
experiment=None,
I_init=1.0,
htc=None,initial_eleConc=None,
initial_soc=0.5,
nproc=12,
output_variables=None,
):
r"""
Solves a pack simulation
Parameters
----------
netlist : pandas.DataFrame
A netlist of circuit elements with format. desc, node1, node2, value.
Produced by liionpack.read_netlist or liionpack.setup_circuit
parameter_values : pybamm.ParameterValues class
A dictionary of all the model parameters
experiment : pybamm.Experiment class
The experiment to be simulated. experiment.period is used to
determine the length of each timestep.
I_init : float, optional
Initial guess for single battery current [A]. The default is 1.0.
htc : float array, optional
Heat transfer coefficient array of length Nspm. The default is None.
initial_soc : float
The initial state of charge for every battery. The default is 0.5
nproc : int, optional
Number of processes to start in parallel for mapping. The default is 12.
output_variables : list, optional
Variables to evaluate during solve. Must be a valid key in the
model.variables
Raises
------
Exception
DESCRIPTION.
Returns
-------
output : ndarray shape [# variable, # steps, # batteries]
simulation output array
"""
if netlist is None or parameter_values is None or experiment is None:
raise Exception("Please supply a netlist, paramater_values, and experiment")
# Get netlist indices for resistors, voltage sources, current sources
Ri_map = netlist["desc"].str.find("Ri") > -1
V_map = netlist["desc"].str.find("V") > -1
I_map = netlist["desc"].str.find("I") > -1
Nspm = np.sum(V_map)
protocol = lp.generate_protocol_from_experiment(experiment)
dt = experiment.period
Nsteps = len(protocol)
# Solve the circuit to initialise the electrochemical models
V_node, I_batt = lp.solve_circuit(netlist)
sim = lp.create_simulation(parameter_values, make_inputs=True)
lp.update_init_conc(sim, SoC=initial_soc)
v_cut_lower = parameter_values["Lower voltage cut-off [V]"]
v_cut_higher = parameter_values["Upper voltage cut-off [V]"]
# The simulation output variables calculated at each step for each battery
# Must be a 0D variable i.e. battery wide volume average - or X-averaged for 1D model
variable_names = [
"Terminal voltage [V]",
"Measured battery open circuit voltage [V]",
"Local ECM resistance [Ohm]",
]
if output_variables is not None:
for out in output_variables:
if out not in variable_names:
variable_names.append(out)
# variable_names = variable_names + output_variables
Nvar = len(variable_names)
# Storage variables for simulation data
shm_i_app = np.zeros([Nsteps, Nspm], dtype=float)
shm_Ri = np.zeros([Nsteps, Nspm], dtype=float)
output = np.zeros([Nvar, Nsteps, Nspm], dtype=float)
# Initialize currents in battery models
shm_i_app[0, :] = I_batt * -1
time = 0
# step = 0
end_time = dt * Nsteps
step_solutions = [None] * Nspm
V_terminal = []
record_times = []
integrator, variables_fn, t_eval = _create_casadi_objects(
I_init, htc[0],initial_eleConc[0], sim, dt, Nspm, nproc, variable_names
)
sim_start_time = ticker.time()
for step in range(Nsteps):
step_solutions, var_eval = _mapped_step(
sim.built_model,
step_solutions,
lp.build_inputs_dict(shm_i_app[step, :], htc,initial_eleConc),
integrator,
variables_fn,
t_eval,
)
output[:, step, :] = var_eval
time += dt
# Calculate internal resistance and update netlist
temp_v = output[0, step, :]
temp_ocv = output[1, step, :]
temp_Ri = np.abs(output[2, step, :])
shm_Ri[step, :] = temp_Ri
netlist.loc[V_map, ("value")] = temp_ocv
netlist.loc[Ri_map, ("value")] = temp_Ri
netlist.loc[I_map, ("value")] = protocol[step]
# print('Stepping time', np.around(ticker.time()-tic, 2), 's')
if np.any(temp_v < v_cut_lower):
print("Low V limit reached")
break
if np.any(temp_v > v_cut_higher):
print("High V limit reached")
break
# step += 1
if time <= end_time:
record_times.append(time)
V_node, I_batt = lp.solve_circuit(netlist)
V_terminal.append(V_node.max())
if time < end_time:
shm_i_app[step + 1, :] = I_batt[:] * -1
all_output = {}
all_output["Time [s]"] = np.asarray(record_times)
all_output["Pack current [A]"] = np.asarray(protocol[: step + 1])
all_output["Pack terminal voltage [V]"] = np.asarray(V_terminal)
all_output["Cell current [A]"] = shm_i_app[: step + 1, :]
for j in range(Nvar):
all_output[variable_names[j]] = output[j, : step + 1, :]
toc = ticker.time()
pybamm.logger.notice(
"Solve circuit time " + str(np.around(toc - sim_start_time, 3)) + "s"
)
return all_output
The example code I ran is:
import liionpack as lp
import numpy as np
import pybamm
# Generate the netlist
netlist = lp.setup_circuit(Np=16, Ns=2, Rb=1e-4, Rc=1e-2, Ri=5e-2, V=3.2, I=80.0)
output_variables = [
'X-averaged total heating [W.m-3]',
'Volume-averaged cell temperature [K]',
'X-averaged negative particle surface concentration [mol.m-3]',
'X-averaged positive particle surface concentration [mol.m-3]',
'X-averaged electrolyte concentration [mol.m-3]',
]
# Heat transfer coefficients
htc = np.ones(32) * 10
#initial_eleConc=np.ones(32) * 1000.0
initial_eleConc=np.array([1000., 1000., 2000., 2000., 1000., 1000., 1000., 1000., 1000.,
1000., 1000., 1000., 1000., 1000., 1000., 1000., 1000., 1000.,
1000., 1000., 1000., 1000., 1000., 3000., 3000., 3000., 1000.,
1000., 1000., 1000., 1000., 1000.])
# Cycling experiment, using PyBaMM
experiment = pybamm.Experiment(
["Charge at 50 A for 30 minutes", "Rest for 15 minutes", "Discharge at 50 A for 30 minutes", "Rest for 30 minutes"],
period="10 seconds",
)
# PyBaMM parameters
chemistry = pybamm.parameter_sets.Chen2020
parameter_values = pybamm.ParameterValues(chemistry=chemistry)
# Solve pack
output = lp.solve(netlist=netlist,
parameter_values=parameter_values,
experiment=experiment,
output_variables=output_variables,
htc=htc,initial_eleConc=initial_eleConc)
lp.plot_output(output)
However, in the simulation the first value of initial electrolyte concentration, i.e., initial_eleConc[0] is populating across all the cells. Here is the output showing that
:
I think the issue is in line 256 of solver_utils.py, i.e.
integrator, variables_fn, t_eval = _create_casadi_objects(
I_init, htc[0],initial_eleConc[0], sim, dt, Nspm, nproc, variable_names
)
where it is only taking the first element of the array initial_eleConc[0]. I gave initial_eleConc array similar to htc(heat transfer coefficient) array. As htc array is going into lumped model, so every cell can have their own individual heat transfer coefficient. In contrast, the parameters in the base model seems not populating across the cells.
![[Bug]: lp make lcapy circuit](https://avatars.githubusercontent.com/u/102300401?v=4)
![[Bug]: pip install in Windows](https://avatars.githubusercontent.com/u/28443643?v=4)
![[Bug]: can't conda install via environment.yml on mac](https://avatars.githubusercontent.com/u/20817509?v=4)
![[Bug]: Prada2013 not working](https://avatars.githubusercontent.com/u/4000822?v=4)
2.8k Dec 30, 2022
3 Feb 1, 2022
1 Jan 24, 2022
52 Oct 16, 2022
28 Dec 23, 2022
20 Dec 31, 2022
1 Dec 29, 2021
1 Jan 9, 2022
4 Jan 19, 2022
111 Dec 31, 2022
3 Jul 29, 2022
3 Apr 8, 2022
9 Apr 26, 2022
147 Jan 4, 2023
5 Aug 17, 2022
28 Nov 21, 2022
3 Nov 3, 2022
1 Dec 20, 2021
3 Jan 4, 2022
7 Oct 9, 2022