Issues when working on parameter estimation using calculate_sensitivities method.

[Mrzhang-hub]

Recently I am working on parameter estimation using calculate_sensitivities method. But there are some issues which I have encountered.

1. Because pybamm has not provided calculate_sensitivities for experiment. Then I have to use the following method, which is not convenient. Is there a better way?

#experiment data
data = pd.read_csv('LGM50_789_1C_25degC.csv')
t = data['Time [s]']
y = data['Voltage [V]']
#design the same condition:["Discharge at 1C until 2.5 V (5 seconds period)", "Rest for 2 hours"]
Q = param['Nominal cell capacity [A.h]']
current = 0 * t
current[0:3597] = Q
param["Current function [A]"] = pybamm.Interpolant(np.array(t), np.array(current), pybamm.t)

2. There is a bug after running:

import numpy as np
import os
import pybamm
import pbparam
import pandas as pd
os.chdir(os.path.join(
    pbparam.__path__[0],
    'input',
    'data'
))
data = pd.read_csv('LGM50_789_1C_25degC.csv')
t = data['Time [s]']
y = data['Voltage [V]']
model = pybamm.lithium_ion.SPMe(
    options={
        "thermal": "lumped",
        "dimensionality": 0,
        "cell geometry": "arbitrary",
        "electrolyte conductivity": "integrated",
    },
    name="TSPMe",
)
param = pybamm.ParameterValues("Chen2020")
param.update({"Negative electrode diffusivity [m2.s-1]": '[input]'}, check_already_exists=False)
param.update({"Positive electrode diffusivity [m2.s-1]": '[input]'}, check_already_exists=False)
Q = param['Nominal cell capacity [A.h]']
current = 0 * t
current[0:3597] = Q
param["Current function [A]"] = pybamm.Interpolant(
    np.array(t), np.array(current), pybamm.t
)
sim = pybamm.Simulation(model, parameter_values=param, solver=pybamm.CasadiSolver())
# initial paramters
start_params = [3.54e-14, 1.01e-13]
sim.solve(inputs={"Negative electrode diffusivity [m2.s-1]": start_params[0],
                "Positive electrode diffusivity [m2.s-1]": start_params[1]},
                calculate_sensitivities=True)

The bug is:
RuntimeError: .../casadi/core/interpolant.cpp:67: Assertion "is_increasing(g)" failed: Gridpoints must be strictly increasing

[brosaplanella]

1. You could define your current as a function (using a heaviside function). However, this is different from the experiment as in your example you discharge for a given time, while in the experiment you discharged until a certain condition was met. I think he latter is what actually causes trouble with the sensitivities, and I can't think of a simple way of hacking it.

2. Hard to tell what's going wrong. Can you test if setting calculate_sensitivities=False works? Then we will know that the problem is with the sensitivities. Also, can you share the full error log?

[Mrzhang-hub]

Sir, thanks for your reply.

1. I have no idea what the heaviside function is and how it is used(Is there any tutorial?). And I have another question, how to set t_eval to when lower voltage cut-off is reached (because experiment can not be used in calculating sensitivity)? If I set t_eval to a fixed value (such as np.arange(0,3600,100), I could not guarantee it will reach the lower voltage cut-off at the end of t_eval.
2. The problem has solved. The bug is caused by the simulation time which is not linearly increasing.

[brosaplanella]

1. You can define your current with a function:

def my_current_function(t):
   return I0 * (t >= t0)

where I0 and t0 are variables with the parameters to match the step function you want.

2. You mean the bug is fixed for calculate_sensitivities=True or that setting it to False runs?

[Mrzhang-hub]

1. It is exactly what I have defined my current function. But is there a way to define a current function, which can realize the same effect with this experimental condition("Discharge at 1C until 2.5V","Rest for 2 hours") ?
2. Yes, the bug is fixed for calculate_sensitivities=True. But there is still a bug when I calculate the parameter sensitivity of terminal voltage towards negative electrode diffusivity([Bug]: ys too large when computing sensitivities pybamm-team/PyBaMM#3130). However, it works normally when calculating terminal voltage towards some other parameters such as negative electrode porosity.

[brosaplanella]

1. No, the only way of setting a stopping condition based on voltage is through an experiment.
2. Not sure where the issue is. I can't run it on my end as the code you share requires additional files, could you please simplify it to a minimum working example (i.e. strip all functionality that is not needed to reproduce the error)?

[Mrzhang-hub]

1. But when using experiment, calculate_sensitivities could not be used. It is so contradictory...
2. the minimum working example:

import pybamm
import numpy as np
model = pybamm.lithium_ion.SPMe()
param = pybamm.ParameterValues('Chen2020')
param.update({'Negative electrode diffusivity [m2.s-1]':'[input]'})
sim = pybamm.Simulation(model,parameter_values=param,solver=pybamm.CasadiSolver())
t_eval = np.arange(0,3600,100)
sol = sim.solve(t_eval,inputs={'Negative electrode diffusivity [m2.s-1]':3.3e-14},calculate_sensitivities=True)

[brosaplanella]

1. Yes, sensitivities cannot be computed with an experiment as discussed in create calculate_sensitivities for experiment pybamm-team/PyBaMM#3133 (basically, experiments combine multiple models so it is not possible to compute a single sensitivity).

2. This seems to be a PyBaMM issue, I have now update the MWE in [Bug]: ys too large when computing sensitivities pybamm-team/PyBaMM#3130, we can continue the discussion there.

[Mrzhang-hub]

Thanks for you patient reply. I am working on gradient-based parameter optimization using sensitivity calculation. When It is done, I would like to upload it to this repository(maybe someone will be interested).

[brosaplanella]

Sounds good! It is amongst the features we want to implement, but not very high on our priority list at the moment. Note that issue #15 is about sensitivities, so let's have future discussions there.