Error with inputs

[amirsh213]

I am trying to follow pybamm youtube instructions on parameterization through optimization but I ran into an error about the inputs. For current, I am using input experiment and interpolation. The script is long and I paste the part that raises the error.

teval= time_exp[0:320]
current_exp= current_exp[0:320]
Voltage_exp= Voltage_exp[0:320]

timescale = parameter_values.evaluate(model.timescale)
current_interpolant = pybamm.Interpolant(teval, -current_exp, timescale * pybamm.t)

param.update({"Current function [A]":current_interpolant})
param.update({'Initial concentration in negative electrode [mol.m-3]': '[input]',
             'Initial concentration in positive electrode [mol.m-3]': '[input]',
             'Negative electrode reference diffusivity [m2.s-1]': '[input]',
             'Positive electrode reference diffusivity [m2.s-1]': '[input]',
             'Negative electrode reference constant [m2.5.mol-0.5.s-1]':'[input]',
             'Positive electrode reference constant [m2.5.mol-0.5.s-1]':'[input]'})

def sumsq(parameters):
    print('solving for "parameters"=', parameters)
    simulation= solver.solve(model, teval, inputs={'Initial concentration in negative electrode [mol.m-3]': parameters[0],
             'Initial concentration in positive electrode [mol.m-3]': parameters[1],
             'Negative electrode reference diffusivity [m2.s-1]': parameters[2],
             'Positive electrode reference diffusivity [m2.s-1]': parameters[3],
             'Negative electrode reference constant [m2.5.mol-0.5.s-1]': parameters[4],
             'Positive electrode reference constant [m2.5.mol-0.5.s-1]': parameters[5]}
                            )["Terminal voltage [V]"](teval)
    
    return np.sum((simulation-Voltage_exp)**2)

print ("starting parameter is=", x0)
res= scipy.optimize.minimize(sumsq, x0, bounds=bounds)

Something that I don't get is that if I use the whole arrays of teval, current_exp, and voltage_exp it does not raise any error although it does not converge. When I use only the first 320 elements it raises the error below.

error

starting parameter is= [800, 35000.0, 1e-13, 1e-13, 1e-11, 1e-11]
solving for "parameters"= [8.0e+02 3.5e+04 1.0e-13 1.0e-13 1.0e-11 1.0e-11]
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
~\AppData\Local\Temp/ipykernel_14904/1874768652.py in <module>
     31 x0=[800, 3.5e4, 1e-13, 1e-13, 1e-11, 1e-11]
     32 print ("starting parameter is=", x0)
---> 33 res= scipy.optimize.minimize(sumsq, x0, bounds=bounds)
     34 print("recovered parameter is=", res.x)
     35 

~\anaconda3\lib\site-packages\scipy\optimize\_minimize.py in minimize(fun, x0, args, method, jac, hess, hessp, bounds, constraints, tol, callback, options)
    621                                   **options)
    622     elif meth == 'l-bfgs-b':
--> 623         return _minimize_lbfgsb(fun, x0, args, jac, bounds,
    624                                 callback=callback, **options)
    625     elif meth == 'tnc':

~\anaconda3\lib\site-packages\scipy\optimize\lbfgsb.py in _minimize_lbfgsb(fun, x0, args, jac, bounds, disp, maxcor, ftol, gtol, eps, maxfun, maxiter, iprint, callback, maxls, finite_diff_rel_step, **unknown_options)
    304             iprint = disp
    305 
--> 306     sf = _prepare_scalar_function(fun, x0, jac=jac, args=args, epsilon=eps,
    307                                   bounds=new_bounds,
    308                                   finite_diff_rel_step=finite_diff_rel_step)

~\anaconda3\lib\site-packages\scipy\optimize\optimize.py in _prepare_scalar_function(fun, x0, jac, args, bounds, epsilon, finite_diff_rel_step, hess)
    259     # ScalarFunction caches. Reuse of fun(x) during grad
    260     # calculation reduces overall function evaluations.
--> 261     sf = ScalarFunction(fun, x0, args, grad, hess,
    262                         finite_diff_rel_step, bounds, epsilon=epsilon)
    263 

~\anaconda3\lib\site-packages\scipy\optimize\_differentiable_functions.py in __init__(self, fun, x0, args, grad, hess, finite_diff_rel_step, finite_diff_bounds, epsilon)
    138 
    139         self._update_fun_impl = update_fun
--> 140         self._update_fun()
    141 
    142         # Gradient evaluation

~\anaconda3\lib\site-packages\scipy\optimize\_differentiable_functions.py in _update_fun(self)
    231     def _update_fun(self):
    232         if not self.f_updated:
--> 233             self._update_fun_impl()
    234             self.f_updated = True
    235 

~\anaconda3\lib\site-packages\scipy\optimize\_differentiable_functions.py in update_fun()
    135 
    136         def update_fun():
--> 137             self.f = fun_wrapped(self.x)
    138 
    139         self._update_fun_impl = update_fun

~\anaconda3\lib\site-packages\scipy\optimize\_differentiable_functions.py in fun_wrapped(x)
    132             # Overwriting results in undefined behaviour because
    133             # fun(self.x) will change self.x, with the two no longer linked.
--> 134             return fun(np.copy(x), *args)
    135 
    136         def update_fun():

~\AppData\Local\Temp/ipykernel_14904/1874768652.py in sumsq(parameters)
      6 def sumsq(parameters):
      7     print('solving for "parameters"=', parameters)
----> 8     simulation= solver.solve(model, teval, inputs={'Initial concentration in negative electrode [mol.m-3]': parameters[0],
      9              'Initial concentration in positive electrode [mol.m-3]': parameters[1],
     10              'Negative electrode reference diffusivity [m2.s-1]': parameters[2],

~\anaconda3\lib\site-packages\pybamm\solvers\base_solver.py in solve(self, model, t_eval, external_variables, inputs, initial_conditions, nproc, calculate_sensitivities)
   1021             ninputs = len(ext_and_inputs_list)
   1022             if ninputs == 1:
-> 1023                 new_solution = self._integrate(
   1024                     model,
   1025                     t_eval_dimensionless[start_index:end_index],

~\anaconda3\lib\site-packages\pybamm\solvers\casadi_solver.py in _integrate(self, model, t_eval, inputs_dict)
    180             # Check if the sign of an event changes, if so find an accurate
    181             # termination point and exit
--> 182             solution = self._solve_for_event(solution, init_event_signs)
    183             return solution
    184         elif self.mode in ["safe", "safe without grid"]:

~\anaconda3\lib\site-packages\pybamm\solvers\casadi_solver.py in _solve_for_event(self, coarse_solution, init_event_signs)
    412         if event_idx_lower is None:
    413             # Flag "final time" for termination
--> 414             self.check_interpolant_extrapolation(model, coarse_solution)
    415             coarse_solution.termination = "final time"
    416             return coarse_solution

~\anaconda3\lib\site-packages\pybamm\solvers\casadi_solver.py in check_interpolant_extrapolation(self, model, solution)
    490                             == pybamm.EventType.INTERPOLANT_EXTRAPOLATION
    491                             and (
--> 492                                 event.expression.evaluate(
    493                                     solution.t[-1],
    494                                     solution.y[:, -1].full(),

~\anaconda3\lib\site-packages\pybamm\expression_tree\functions.py in evaluate(self, t, y, y_dot, inputs, known_evals)
    180             return known_evals[self.id], known_evals
    181         else:
--> 182             evaluated_children = [
    183                 child.evaluate(t, y, y_dot, inputs) for child in self.children
    184             ]

~\anaconda3\lib\site-packages\pybamm\expression_tree\functions.py in <listcomp>(.0)
    181         else:
    182             evaluated_children = [
--> 183                 child.evaluate(t, y, y_dot, inputs) for child in self.children
    184             ]
    185             return self._function_evaluate(evaluated_children)

~\anaconda3\lib\site-packages\pybamm\expression_tree\binary_operators.py in evaluate(self, t, y, y_dot, inputs, known_evals)
    151                 return value, known_evals
    152         else:
--> 153             left = self.left.evaluate(t, y, y_dot, inputs)
    154             right = self.right.evaluate(t, y, y_dot, inputs)
    155             return self._binary_evaluate(left, right)

~\anaconda3\lib\site-packages\pybamm\expression_tree\binary_operators.py in evaluate(self, t, y, y_dot, inputs, known_evals)
    151                 return value, known_evals
    152         else:
--> 153             left = self.left.evaluate(t, y, y_dot, inputs)
    154             right = self.right.evaluate(t, y, y_dot, inputs)
    155             return self._binary_evaluate(left, right)

~\anaconda3\lib\site-packages\pybamm\expression_tree\symbol.py in evaluate(self, t, y, y_dot, inputs, known_evals)
    748             return known_evals[self.id], known_evals
    749         else:
--> 750             return self._base_evaluate(t, y, y_dot, inputs)
    751 
    752     def evaluate_for_shape(self):

~\anaconda3\lib\site-packages\pybamm\expression_tree\input_parameter.py in _base_evaluate(self, t, y, y_dot, inputs)
     68             if inputs == "shape test":
     69                 return self.evaluate_for_shape()
---> 70             raise TypeError("inputs should be a dictionary")
     71         try:
     72             input_eval = inputs[self.name]

TypeError: inputs should be a dictionary

[valentinsulzer]

Please post a complete script we can run to diagnose this error.

[amirsh213]

Below is the complete script. I defined my own parameters using Chen2020.

# Importing experiment

data_cur= np.loadtxt(r'C:\Users\plati\Anaconda\Python\Input data\time_current_1C.txt')
time_exp= data_cur[:,0]
current_exp= 0.098*data_cur[:,1]

data_vol= np.loadtxt(r'C:\Users\plati\Anaconda\Python\Input data\time_voltage_1C.txt')
time_exp= data_vol[:,0]
Voltage_exp= data_vol[:,1] 

teval= time_exp[0:320]
current_exp= current_exp[0:320]
Voltage_exp= Voltage_exp[0:320]

# Building model to solve by solver class

model= pybamm.lithium_ion.SPMe()
geometry = model.default_geometry

param = pybamm.ParameterValues(chemistry=chemistry)
timescale = parameter_values.evaluate(model.timescale)
current_interpolant = pybamm.Interpolant(teval, -current_exp, timescale * pybamm.t)

param.update({"Current function [A]":current_interpolant})
param.update({'Initial concentration in negative electrode [mol.m-3]': '[input]',
             'Initial concentration in positive electrode [mol.m-3]': '[input]',
             'Negative electrode reference diffusivity [m2.s-1]': '[input]',
             'Positive electrode reference diffusivity [m2.s-1]': '[input]',
             'Negative electrode reference constant [m2.5.mol-0.5.s-1]':'[input]',
             'Positive electrode reference constant [m2.5.mol-0.5.s-1]':'[input]'})

param.process_geometry(geometry)
param.process_model(model)
var_pts= model.default_var_pts
mesh = pybamm.Mesh(geometry, model.default_submesh_types, model.default_var_pts)
disc = pybamm.Discretisation(mesh, model.default_spatial_methods)
disc.process_model(model);

solver= pybamm.CasadiSolver(mode='fast', rtol=1e-3, atol=1e-6)

# Finding parameters via optimization

def sumsq(parameters):
    print('solving for "parameters"=', parameters)
    simulation= solver.solve(model, teval, inputs={'Initial concentration in negative electrode [mol.m-3]': parameters[0],
             'Initial concentration in positive electrode [mol.m-3]': parameters[1],
             'Negative electrode reference diffusivity [m2.s-1]': parameters[2],
             'Positive electrode reference diffusivity [m2.s-1]': parameters[3],
             'Negative electrode reference constant [m2.5.mol-0.5.s-1]': parameters[4],
             'Positive electrode reference constant [m2.5.mol-0.5.s-1]': parameters[5]}
                            )["Terminal voltage [V]"](teval)
    
    return np.sum((simulation-Voltage_exp)**2)
    
bounds= ((0,1000), (3e4, 5.2e4), (1e-14, 1e-12), (1e-14, 1e-12), (1e-12, 1e-10), (1e-12, 1e-10)) # scipy.minimize

x0=[800, 3.5e4, 1e-13, 1e-13, 1e-11, 1e-11] 
print ("starting parameter is=", x0)
res= scipy.optimize.minimize(sumsq, x0, bounds= bounds)
print("recovered parameter is=", res.x)


[time_current_1C.txt](https://github.com/pybamm-team/PyBaMM/files/8219094/time_current_1C.txt)
[time_voltage_1C.txt](https://github.com/pybamm-team/PyBaMM/files/8219095/time_voltage_1C.txt)

[valentinsulzer]

This code doesn't run. Please check that the code you provide actually runs. Also, it's not really a MWE - I'm sure there are parts of this code you could remove and still get the bug

[amirsh213]

I edited the script as below and now it runs. It seems that again it has a problem with data interpolation of cathode OCP (attached). If I use the "Kim NCA OCP" function it does not raise any error but it does not like either my input data or Kim's data.

# Importing experiment

data_cur= np.loadtxt('Input data/time_current_1C.txt')
time_exp= data_cur[:,0]
current_exp= 0.098*data_cur[:,1]

data_vol= np.loadtxt('Input data/time_voltage_1C.txt')
time_exp= data_vol[:,0]
Voltage_exp= data_vol[:,1] 

teval= time_exp[0:320]
current_exp= current_exp[0:320]
Voltage_exp= Voltage_exp[0:320]

# Building model to solve by solver class

chemistry=pybamm.parameter_sets.Chen2020 # starting from an existing chemistry
parameter_values = pybamm.ParameterValues(chemistry=chemistry)
param = pybamm.ParameterValues(chemistry=chemistry)

model= pybamm.lithium_ion.SPMe()
geometry = model.default_geometry
timescale = parameter_values.evaluate(model.timescale)
current_interpolant = pybamm.Interpolant(teval, -current_exp, timescale * pybamm.t)

param.update({"Current function [A]":current_interpolant})
param.update({'Initial concentration in negative electrode [mol.m-3]': '[input]',
             'Initial concentration in positive electrode [mol.m-3]': '[input]'})

param.process_geometry(geometry)
param.process_model(model)
var_pts= model.default_var_pts
mesh = pybamm.Mesh(geometry, model.default_submesh_types, model.default_var_pts)
disc = pybamm.Discretisation(mesh, model.default_spatial_methods)
disc.process_model(model);

solver= pybamm.CasadiSolver(mode='fast', rtol=1e-3, atol=1e-6)

# Finding parameters via optimization

def sumsq(parameters):
    print('solving for "parameters"=', parameters)
    simulation= solver.solve(model, teval, inputs={'Initial concentration in negative electrode [mol.m-3]': parameters[0],
             'Initial concentration in positive electrode [mol.m-3]': parameters[1]}
                            )["Terminal voltage [V]"](teval)
    
    return np.sum((simulation-Voltage_exp)**2)
    
bounds= ((0,1000), (3e4, 5.2e4)) # scipy.minimize

x0=[800, 3.5e4] 
print ("starting parameter is=", x0)
res= scipy.optimize.minimize(sumsq, x0, bounds= bounds)
print("recovered parameter is=", res.x)

time_current_1C.txt
time_voltage_1C.txt
nca_ocp.csv

[amirsh213]

I noticed that you are fixing the bug in #1970 which is related to OCP data interpolation. I wanted to add something about the drive cycle (current interpolant). Even if I use a function for the positive OCP, the SPMe model does not simulate it (images attached in #1970). It seems that the current interpolant does not read the current properly in this case. Is this issue related to the other as well?