Obtaining the EIS for a sinusoidal input current

[Akila1993]

Hi All,

I am trying to obtain the EIS for a given AC input current signal as indicated in the attached file. Can anyone suggested how to extract the terminal voltage and current entries from this in order to obtain the EIS.

I tried to obtain the EIS via #pbeis but it seems the results are wrong.

Please provide your insights in this regard.

impoting libraries

import pybamm
import numpy as np
from functools import partial
import matplotlib.pyplot as plt
import pbeis

Define the DFN model

model = pybamm.lithium_ion.DFN()

params = pybamm.ParameterValues("Chen2020")

def current(t, amp=100e-3, freq=1):
return amp * pybamm.sin(2 * np.pi * freq * t)

set up frequency params

freq = 2e3
samples_per_period = 128
num_periods = 3

set up t_eval

dt = 1 / freq / samples_per_period
t_eval = np.array(range(0, samples_per_period*num_periods)) * dt

Update current function with 2 kHz sinusoid

params['Current function [A]'] = partial(current, amp=100e-3, freq=freq)

Set up solvers

solver_fast = pybamm.CasadiSolver(mode='fast')

Get different solutions

sol_fast = pybamm.Simulation(model=model, parameter_values=params).
solve(t_eval=t_eval, solver=solver_fast)

sim_step = pybamm.Simulation(model=model, parameter_values=params, solver=solver_fast)
t, t_max, dt = 0, t_eval[-1], t_eval[1]
while t < t_max:
sol_step = sim_step.step(dt)
t += dt

fig, axs = plt.subplots(2, 1, figsize=(10, 10/1.6), sharex=True,
constrained_layout=True)
for y_key, ax in zip(['Current [A]', 'Terminal voltage [V]'], axs):
for sol, label in zip([sol_fast],
['fast']):
ax.plot(sol['Time [s]'].data, sol[y_key].data, lw=2, label=label)
ax.legend(bbox_to_anchor=(1, 0.5), loc='center left', fontsize=16)
ax.set_ylabel(y_key, fontsize=18)
ax.tick_params(labelsize=16)
ax.grid('on')

axs[1].set_xlabel('Time [s]', fontsize=18)
fig.suptitle('Solver diffs', fontsize=20)
fig.show()

EIS plot

eis_sim = pbeis.EISSimulation(model=model, parameter_values=params)
eis_sim.solve(pbeis.logspace(-4,4,100)) # calculate impedance at log-spaced frequencies
eis_sim.nyquist_plot()

[rtimms]

Here is an example using SPM and the Chen 2020 parameter set (DFN works too). Note:

1. you need to include the option "surface form": "differential" to include capacitance in your model
2. pbeis has only been tested on PyBaMM version 22.10. It was worked on as part of a summer project and I haven't had the bandwidth to continue working on it since then.
3. There is a bug with the CasadiSolver that might give incorrect results when doing high-frequency simulations (see this discussion, so best to use the SicpySolver for now. The IDAKLUSolver should work too, but I haven't tested it.

import pbeis
import pybamm
import numpy as np
import matplotlib.pyplot as plt
import time as timer
from scipy.fft import fft

# Set up
model = pybamm.lithium_ion.SPM(options={"surface form": "differential"}, name="SPM")
parameter_values = pybamm.ParameterValues("Chen2020")
frequencies = np.logspace(-4, 2, 30)

# Time domain
I = 50 * 1e-3
number_of_periods = 20
samples_per_period = 16


def current_function(t):
    return I * pybamm.sin(2 * np.pi * pybamm.InputParameter("Frequency [Hz]") * t)


parameter_values["Current function [A]"] = current_function

start_time = timer.time()

sim = pybamm.Simulation(
    model, parameter_values=parameter_values, solver=pybamm.ScipySolver()
)

impedances_time = []
for frequency in frequencies:
    # Solve
    period = 1 / frequency
    dt = period / samples_per_period
    t_eval = np.array(range(0, 1 + samples_per_period * number_of_periods)) * dt
    sol = sim.solve(t_eval, inputs={"Frequency [Hz]": frequency})
    # Extract final two periods of the solution
    time = sol["Time [s]"].entries[-3 * samples_per_period - 1 :]
    current = sol["Current [A]"].entries[-3 * samples_per_period - 1 :]
    voltage = sol["Terminal voltage [V]"].entries[-3 * samples_per_period - 1 :]
    # FFT
    current_fft = fft(current)
    voltage_fft = fft(voltage)
    # Get index of first harmonic
    idx = np.argmax(np.abs(current_fft))
    impedance = -voltage_fft[idx] / current_fft[idx]
    impedances_time.append(impedance)

end_time = timer.time()
time_elapsed = end_time - start_time
print("Time domain method: ", time_elapsed, "s")

# Frequency domain
methods = ["direct"]
impedances_freqs = []
for method in methods:
    start_time = timer.time()
    eis_sim = pbeis.EISSimulation(model, parameter_values=parameter_values)
    impedances_freq = eis_sim.solve(frequencies, method)
    end_time = timer.time()
    time_elapsed = end_time - start_time
    print(f"Frequency domain ({method}): ", time_elapsed, "s")
    impedances_freqs.append(impedances_freq)

# Compare
_, ax = plt.subplots()
ax = pbeis.nyquist_plot(impedances_time, ax=ax, label="Time", alpha=0.7)
for i, method in enumerate(methods):
    ax = pbeis.nyquist_plot(
        impedances_freqs[i], ax=ax, label=f"Frequency ({method})", alpha=0.7
    )
ax.legend()
plt.suptitle(f"{model.name}")
plt.show()

[brosaplanella]

My guess here is that OKane2022 has many parameters depending on concentration. EIS uses a linearised version of the model so it is only an approximation. One way to check if this is true is that the two solutions should get closer together as the current amplitude decreases.

[Akila1993]

Hi @brosaplanella your suggestion is correct I change current amplitude 1 and it gives me the following result.

It is possible for you to support me for this discussion also please:

[Akila1993]

In this sample code if I need to calculate SOC throughout the simulation for both frequency domain and time domain,

Can I use the following for SOC calculation?
SoC_init = 1;
Q = parameter_values["Nominal cell capacity [A.h]"];
SoC = SoC_init - sol["Discharge Capacity [A.h]"]/Q;

This will work for the Time domain as it stores solutions in sol. In that case, how can I get this from pbeis.EISSimulation.

[brosaplanella]

You mean defining the SoC you want to use or computing what the SoC is for a given simulation?

The actual implementation will be slightly different for each case, but you can update the initial concentrations to match your desired SoC.

[Akila1993]

Hi @brosaplanella Yes, I will need to define the SOC that I want to use. In that case how I can approach it? Could you please advise or direct me to such a discussion?

[Akila1993]

Here is an example using SPM and the Chen 2020 parameter set (DFN works too). Note:

1. you need to include the option "surface form": "differential" to include capacitance in your model
2. pbeis has only been tested on PyBaMM version 22.10. It was worked on as part of a summer project and I haven't had the bandwidth to continue working on it since then.
3. There is a bug with the CasadiSolver that might give incorrect results when doing high-frequency simulations (see this discussion, so best to use the SicpySolver for now. The IDAKLUSolver should work too, but I haven't tested it.

import pbeis
import pybamm
import numpy as np
import matplotlib.pyplot as plt
import time as timer
from scipy.fft import fft

# Set up
model = pybamm.lithium_ion.SPM(options={"surface form": "differential"}, name="SPM")
parameter_values = pybamm.ParameterValues("Chen2020")
frequencies = np.logspace(-4, 2, 30)

# Time domain
I = 50 * 1e-3
number_of_periods = 20
samples_per_period = 16


def current_function(t):
    return I * pybamm.sin(2 * np.pi * pybamm.InputParameter("Frequency [Hz]") * t)


parameter_values["Current function [A]"] = current_function

start_time = timer.time()

sim = pybamm.Simulation(
    model, parameter_values=parameter_values, solver=pybamm.ScipySolver()
)

impedances_time = []
for frequency in frequencies:
    # Solve
    period = 1 / frequency
    dt = period / samples_per_period
    t_eval = np.array(range(0, 1 + samples_per_period * number_of_periods)) * dt
    sol = sim.solve(t_eval, inputs={"Frequency [Hz]": frequency})
    # Extract final two periods of the solution
    time = sol["Time [s]"].entries[-3 * samples_per_period - 1 :]
    current = sol["Current [A]"].entries[-3 * samples_per_period - 1 :]
    voltage = sol["Terminal voltage [V]"].entries[-3 * samples_per_period - 1 :]
    # FFT
    current_fft = fft(current)
    voltage_fft = fft(voltage)
    # Get index of first harmonic
    idx = np.argmax(np.abs(current_fft))
    impedance = -voltage_fft[idx] / current_fft[idx]
    impedances_time.append(impedance)

end_time = timer.time()
time_elapsed = end_time - start_time
print("Time domain method: ", time_elapsed, "s")

# Frequency domain
methods = ["direct"]
impedances_freqs = []
for method in methods:
    start_time = timer.time()
    eis_sim = pbeis.EISSimulation(model, parameter_values=parameter_values)
    impedances_freq = eis_sim.solve(frequencies, method)
    end_time = timer.time()
    time_elapsed = end_time - start_time
    print(f"Frequency domain ({method}): ", time_elapsed, "s")
    impedances_freqs.append(impedances_freq)

# Compare
_, ax = plt.subplots()
ax = pbeis.nyquist_plot(impedances_time, ax=ax, label="Time", alpha=0.7)
for i, method in enumerate(methods):
    ax = pbeis.nyquist_plot(
        impedances_freqs[i], ax=ax, label=f"Frequency ({method})", alpha=0.7
    )
ax.legend()
plt.suptitle(f"{model.name}")
plt.show()

Hi @rtimms if I need to validate this results what would be possibility to check the results from PyBaMM is accurate enough? Could you please guide me.