Easy statistical distributions

[TomTranter]

Description

Add a nice way to pass a statistical distribution as an input for varying circuit and battery parameters. For example if you want to create a pack with a single particle model and vary the particle radius according to a statistical distribution we could add a helper function or an example of how to do that properly.

Motivation

This will help to model packs with real world statistical variation in parameters

Possible Implementation

It is already possible to pass an array of values as an input. We could add an example showing how to populate this array from a distribution or could go one step further and hand the selection of values from the distribution to the solver algorithms by passing the distribution object. This might be useful for comparing multiple different variations...

Additional context

No response

[valentinsulzer]

I agree this would be very useful to have. It would be good to implement this in a way that doesn't require duplicate changes to both PyBaMM and liionpack.

A possible non-invasive solution for this would be to implement helper functions to do the following:

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

model = pybamm.lithium_ion.SPM()

param = model.default_parameter_values
param.update(
    {
        "Negative electrode diffusivity [m2.s-1]": pybamm.InputParameter("D_s_n"),
        "Positive electrode diffusivity [m2.s-1]": pybamm.InputParameter("D_s_p"),
    }
)
sim = pybamm.Simulation(model, parameter_values=param)

Vs = []
n = 100
Dsn_set = np.random.normal(loc=-13, scale=0.5, size=n)
Dsp_set = np.random.normal(loc=-13, scale=0.5, size=n)
for i, (D_s_n, D_s_p) in enumerate(zip(Dsn_set, Dsp_set)):
    print(i, D_s_n, D_s_p)
    inputs = {
        "D_s_n": 10**D_s_n,
        "D_s_p": 10**D_s_p,
    }
    solution = sim.solve([0, 3600], inputs=inputs)
    Vs.append(solution["Terminal voltage [V]"])

Vs_2000 = np.concatenate([V(1000).flatten() for V in Vs])
print(Vs_2000)
fig, axes = plt.subplots(2, 2)

axes.flat[0].hist(Dsn_set)
axes.flat[1].hist(Dsp_set)
axes.flat[2].hist(Vs_2000)
axes.flat[0].set_title("Dsn")
axes.flat[1].set_title("Dsp")
axes.flat[2].set_title("V")
plt.show()

[valentinsulzer]

Got thoroughly nerdsniped by this. Here are some cooler plots.

#
# Example showing how to load and solve the DFN
#

import pybamm
import numpy as np
import pandas as pd
from scipy.interpolate import interp1d

model = pybamm.lithium_ion.SPM()

param = model.default_parameter_values
param.search("active material volume fraction")
param.update(
    {
        "Negative electrode active material volume fraction": pybamm.InputParameter(
            "eps_n"
        ),
        "Positive electrode active material volume fraction": pybamm.InputParameter(
            "eps_p"
        ),
    }
)
sim = pybamm.Simulation(model, parameter_values=param)

n = 200
epsn_set = np.random.normal(loc=0.6, scale=0.05, size=n)
epsp_set = np.random.normal(loc=0.5, scale=0.05, size=n)


t_eval = np.linspace(0, 5000, 100)
results = []
summary = []
V_interp = np.linspace(
    param["Upper voltage cut-off [V]"], param["Lower voltage cut-off [V]"], 100
)
for i, (eps_n, eps_p) in enumerate(zip(epsn_set, epsp_set)):
    print(i, eps_n, eps_p)
    inputs = {"eps_n": eps_n, "eps_p": eps_p}
    solution = sim.solve(t_eval, inputs=inputs)
    capacity = solution["Discharge capacity [A.h]"].entries[-1]
    summary.append({**inputs, "Capacity [A.h]": capacity})

    t = solution["Time [s]"].data
    V = solution["Terminal voltage [V]"].data
    t_interp = interp1d(V, t, kind="linear", bounds_error=False)(V_interp)
    for t, V in zip(t_interp, V_interp):
        results.append(
            {
                **inputs,
                "Time [s]": t,
                "Voltage [V]": V,
                "run": i,
                "Capacity [A.h]": capacity,
            }
        )

df = pd.DataFrame(results)
df_summary = pd.DataFrame(summary)

fig, ax = plt.subplots(2, 3, figsize=(10, 6))
for i, var in [(0, "eps_n"), (1, "eps_p")]:
    sns.distplot(
        df[var],
        hist=True,
        kde=True,
        bins=int(180 / 5),
        color="darkblue",
        hist_kws={"edgecolor": "black"},
        kde_kws={"linewidth": 4},
        ax=ax[0][i],
    )

sns.scatterplot(
    data=df_summary, x="eps_n", y="eps_p", hue="Capacity [A.h]", ax=ax[0][2]
)

sns.lineplot(
    data=df, x="Time [s]", y="Voltage [V]", orient="y", errorbar=("sd", 1), ax=ax[1][0]
)
ax[1][0].set_ylabel("Voltage with one standard deviation")
sns.lineplot(
    data=df, x="Time [s]", y="Voltage [V]", orient="y", hue="eps_n", ax=ax[1][1]
)
sns.lineplot(
    data=df, x="Time [s]", y="Voltage [V]", orient="y", hue="eps_p", ax=ax[1][2]
)
fig.tight_layout()
plt.show()

[TomTranter]

Haha that is beautiful but yes you did! Leave something for the hackers to do

…

On Thu, 22 Sep 2022 at 20:01, Valentin Sulzer ***@***.***> wrote: Got thoroughly nerdsniped by this. Here are some cooler plots. ## Example showing how to load and solve the DFN# import pybammimport numpy as npimport pandas as pdfrom scipy.interpolate import interp1d model = pybamm.lithium_ion.SPM() param = model.default_parameter_valuesparam.search("active material volume fraction")param.update( { "Negative electrode active material volume fraction": pybamm.InputParameter( "eps_n" ), "Positive electrode active material volume fraction": pybamm.InputParameter( "eps_p" ), } )sim = pybamm.Simulation(model, parameter_values=param) n = 200epsn_set = np.random.normal(loc=0.6, scale=0.05, size=n)epsp_set = np.random.normal(loc=0.5, scale=0.05, size=n) t_eval = np.linspace(0, 5000, 100)results = []summary = []V_interp = np.linspace( param["Upper voltage cut-off [V]"], param["Lower voltage cut-off [V]"], 100 )for i, (eps_n, eps_p) in enumerate(zip(epsn_set, epsp_set)): print(i, eps_n, eps_p) inputs = {"eps_n": eps_n, "eps_p": eps_p} solution = sim.solve(t_eval, inputs=inputs) capacity = solution["Discharge capacity [A.h]"].entries[-1] summary.append({**inputs, "Capacity [A.h]": capacity}) t = solution["Time [s]"].data V = solution["Terminal voltage [V]"].data t_interp = interp1d(V, t, kind="linear", bounds_error=False)(V_interp) for t, V in zip(t_interp, V_interp): results.append( { **inputs, "Time [s]": t, "Voltage [V]": V, "run": i, "Capacity [A.h]": capacity, } ) df = pd.DataFrame(results)df_summary = pd.DataFrame(summary) fig, ax = plt.subplots(2, 3, figsize=(10, 6))for i, var in [(0, "eps_n"), (1, "eps_p")]: sns.distplot( df[var], hist=True, kde=True, bins=int(180 / 5), color="darkblue", hist_kws={"edgecolor": "black"}, kde_kws={"linewidth": 4}, ax=ax[0][i], ) sns.scatterplot( data=df_summary, x="eps_n", y="eps_p", hue="Capacity [A.h]", ax=ax[0][2] ) sns.lineplot( data=df, x="Time [s]", y="Voltage [V]", orient="y", errorbar=("sd", 1), ax=ax[1][0] )ax[1][0].set_ylabel("Voltage with one standard deviation")sns.lineplot( data=df, x="Time [s]", y="Voltage [V]", orient="y", hue="eps_n", ax=ax[1][1] )sns.lineplot( data=df, x="Time [s]", y="Voltage [V]", orient="y", hue="eps_p", ax=ax[1][2] )fig.tight_layout()plt.show() [image: Figure_1] <https://user-images.githubusercontent.com/20817509/191829935-e9767172-2404-4224-9fc9-27b7187e2402.png> — Reply to this email directly, view it on GitHub <#190 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ABV5YRLZU6OYVPZQ6PXTKMLV7SUHPANCNFSM6AAAAAAQS27DZU> . You are receiving this because you authored the thread.Message ID: ***@***.***>