Battery response sensitivities for an electric aircraft discharge profile

[AlexGuibert]

Hello,
I am trying to use PyBaMM to model the battery response to a flight profile of an electric aircraft. I would like to do it such that I can also compute sensitivities of some quantities of interest (e.g., final temperature) with respect to some battery parameters (e.g., particle radius). Usually, I would use the Experiment module from PyBaMM to model the battery response to flight but it seems that experiments do not allow for sensitivities. So, I was thinking of dividing the flight into simple segments where the load is constant and computing the sensitivities segment-wise. However, to do this I need to be able to initialize each simulation (segment) from the previous one but I am unsure how to do it.

I was thinking of a class such as the following to be called in a main file:

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

class electro_model:
    def run(input):
        # input
        # - 0: solution to start from (from previous run) | None if the first run
        # - 1: current for the segment [A]
        # - 2: duration of the segment [s]

        # output
        # solution for next run which contains also info such as voltage/temperature increase and sensitivities

        # 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")

        # Definition of the chemistries
        chemistry = pb.parameter_sets.Chen2020         
        parameters = pb.ParameterValues(chemistry=chemistry)
        parameters.process_geometry(geometry)
        parameters.process_model(model)

        # Set up sensitivity analysis - Define design variables
        parameters.update({"Ambient temperature [K]": "[input]"})

        # 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)

        # Time discretization
        tstep = 0.1 # [s]
        t_eval = np.linspace(0, input[2], int(input[2]/tstep))

        # Solving the model
        discharge = input[1] # [A]
        solver = pb.CasadiSolver(mode="fast")
        solution = solver.solve(
            model, t_eval, inputs={'Ambient temperature [K]': 298, 'Current function [A]': discharge},
            calculate_sensitivities=True, initial_conditions=input[0]
        )

        # output solution 
        return solution

The main file would then call the function run for each segment such that

from electro import electro_model
import matplotlib.pyplot as plt
import numpy as np

# The flight profile is the following
#       Discharge at 4 A for 60 s
#       Discharge at 2 A for 600 s
#       Discharge at 4 A for 60 s
loads       = [4.0, 4.0, 4.0, 4.0, 4.0]       # in Amps
duration    = [6.0, 6.0, 6.0, 6.0, 6.0]   # in s

# input for first run 
input = []
input.append(None)  # initial solution is None
input.append(loads[0])     # current for first run [A]
input.append(duration[0])    # duration of the first run [s]

# first run
solution   = electro_model.run(input)

# inialize the solution vectors
voltage     = solution["Terminal voltage [V]"].entries
temperature = solution["Volume-averaged cell temperature [K]"].entries
time        = solution["Time [s]"].entries

# run with the remaining segment
for i in range(len(loads)-1):
    input[0] = solution
    input[1] = loads[i+1]
    input[2] = duration[i+1]
    # run
    solution = electro_model.run(input)
    # append the results
    voltage     = np.append(voltage, solution["Terminal voltage [V]"].entries)
    temperature = np.append(temperature, solution["Volume-averaged cell temperature [K]"].entries)
    time        = np.append(time, solution["Time [s]"].entries)

I am unsure if I am passing properly the previous solution to initialize a model. Also, is it possible within PyBaMM? What would be the best way to enable segment-wise sensitivities?

Thank you! Any help or advice would be greatly appreciated.

[TomTranter]

Hi @AlexGuibert yeah this seems a bit tricky to do as sensitivities are calculated by the solver and there isn't a way to control them through a simulation which is what you would use to run an experiment and also start a simulation from an old solution. I think there are two options: Run a simulation and experiment then manually calculate sensitivities yourself or change the Current function to be a function instead of a constant and run that through a solver. Maybe @martinjrobins could advise on whether this assessment is the true - he knows more about sensitivities and might have another suggestion

[TomTranter]

There's an example here on using an interpolant as Current function [A]
https://docs.pybamm.org/en/stable/source/examples/notebooks/getting_started/tutorial-4-setting-parameter-values.html

[AlexGuibert]

Hi @TomTranter,
Thanks a lot for the quick response! It solved my problem, I created a single function for the current as you suggested with all the segments, and relaxed the discontinuity by using smooth Heaviside functions with pybamm.tanh() .

I am now trying to validate the sensitivities given by PyBaMM against a simple finite difference scheme and it seems that there is some scaling issue between the finite difference results and PyBaMM (~ a factor of 3 difference) as you can see in the picture. Any idea where it might be coming from?

Here I am computing the sensitivity of the average temperature with respect to the porosity and the percentages in the legend indicate the perturbation size with respect to the nominal value. I observe a similar trend for the voltage instead of the temperature.

Thanks!

[TomTranter]

Can you open a new discussion and post the code?

[AlexGuibert]

I have opened a new discussion regarding this topic, here is the link to it #3814

Thanks!

[davidhowey]

FYI PyBOP should be able to do this (eventually), although it's still a work in progress. cc @BradyPlanden see https://github.com/pybop-team/PyBOP for the repo

[AlexGuibert]

Thanks for sharing, looking forward to it!