Scaling of sensitivities in PyBaMM against finite difference check

[AlexGuibert]

Hello,

I am trying to use PyBaMM to compute the sensitivities of some functions of interest (e.g., average temperature) with a custom discharge profile. It is not complaining about the implementation and the code is running. However, when I am trying to validate the sensitivities with a finite difference scheme, it seems that the sensitivities from PyBaMM are scaled compared to the one obtained with the finite difference (~ a factor of 3 difference.). For instance here for the average temperature w.r.t. positive electrode porosity:

Is there a reason why this is happening and am I missing something here? Please find below a minimum working example I have used to generate the plot above.

import pybamm as pb
import numpy as np
import matplotlib.pyplot as plt

# Creation of the models
options = {"thermal": "lumped"}                     # electrochemical model coupled to a lumped thermal model
model = pb.lithium_ion.DFN(options=options)         # Doyle Fuller Newman (P2D) model
geometry = pb.battery_geometry(form_factor="cylindrical")
pb.set_logging_level("NOTICE")

# Define flight profile
segment1 = 8 # [A]
segment2 = 4 # [A]

t_1 = 30. # [s]
t_2 = 30. # [s]

kH = 15 # sharpness of the transition for the Heaviside function
def H(t, t_start, t_end, load):
    H_start = 0.5*load + 0.5*load*pb.tanh(kH*(t - t_start))
    H_end = 0.5*load + 0.5*load*pb.tanh(kH*(t - t_end))
    return (H_start - H_end) 

def my_current(t):
    flight = H(t, 0, t_1, segment1) 
    flight += H(t, t_1, t_2, segment2)
    return flight

# Time discretization
t_eval = np.arange(0, t_1 + t_2, 1)

# Definition of the chemistries
chemistry = pb.parameter_sets.Chen2020         
parameters = pb.ParameterValues(chemistry=chemistry)
parameters["Current function [A]"] = my_current
parameters.update({"Positive electrode porosity": "[input]"}) # design variables
parameters.process_geometry(geometry)
parameters.process_model(model)

# Discretization and notice options
var_pts = {"x_n": 30, "x_s": 30, "x_p": 30, "r_n": 10, "r_p": 10}
mesh = pb.Mesh(geometry, model.default_submesh_types, var_pts)
discretization = pb.Discretisation(mesh, model.default_spatial_methods)
discretization.process_model(model)


# Solving the model
solver = pb.CasadiSolver(mode="fast")
solution = solver.solve(
    model, t_eval, inputs={"Positive electrode porosity": 0.335},
    calculate_sensitivities=True
)

# solution vectors
time = solution["Time [s]"].entries
volt = solution["Volume-averaged cell temperature [K]"].entries
dTdp = solution["Volume-averaged cell temperature [K]"].sensitivities["Positive electrode porosity"]

# Solving the "perturbed" model (0.01%)
perturbation = 0.0001
solution_pert = solver.solve(
    model, t_eval, inputs={"Positive electrode porosity": 0.335*(1 - perturbation)},
    calculate_sensitivities=True)
volt_perturb = solution_pert["Volume-averaged cell temperature [K]"].entries
dTdp_pert_1   = (volt - volt_perturb) / perturbation

# Postprocessing
# time vs temperature sensitivities
plt.figure(figsize=(7.0, 4.8))
plt.xlabel(r"Time [s]")
plt.ylabel(r"dTdp")
plt.plot(time, dVdp, linewidth=4)
plt.plot(time, dTdp_pert_1, linewidth=4)
plt.legend(["PyBaMM", "Finite diff. 0.01 %"])
plt.show()

Thanks a lot!

[martinjrobins]

Hi @AlexGuibert, the difference you see here is related to the base value you have chosen for the parameter (i.e. 0.335), the finite difference formula should be dTdt = (volt - volt_perturb) / (0.335 * perturbation)

[AlexGuibert]

Hi, yup I missed it (I shouldn't code on a Friday evening). Thanks a lot!