Obtaining EIS after running experiment with OKane2022 degradation models

[Akila1993]

Hello,

I am trying to obtain the EIS after running 10 cycles with OKane2022 degradation models by using the code below. But It gives me only one wrong impedance value. See the figure also below. Would someone can correct if my approach is wrong.

Time domain initials

frequencies = np.logspace(-2, 4, 10)
I = 1 * 1e-3

number_of_periods = 20
samples_per_period = 40

def current_function(t):

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

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

start_time = timer.time()

cycle_number = 10
exp = pybamm.Experiment(
["Hold at 4.2 V until C/100",
"Rest for 4 hours",
"Discharge at 0.1C until 2.5 V", # initial capacity check
"Charge at 0.3C until 4.2 V",
"Hold at 4.2 V until C/100",]
+ [("Discharge at 1C until 2.5 V", # ageing cycles
"Charge at 0.3C until 4.2 V",
"Hold at 4.2 V until C/100",)] * cycle_number
+ ["Discharge at 0.1C until 2.5 V"], # final capacity check
period="5 minutes"
)

sim = pybamm.Simulation(model, parameter_values=updated_parameters,experiment=exp,solver=pybamm.CasadiSolver(mode='safe without grid'),var_pts=var_pts)

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

#EIS Nyquist plot
_, ax = plt.subplots()
ax = pbeis.nyquist_plot(impedances_time, ax=ax, label="Time", alpha=0.7)

ax.legend()

plt.suptitle(f"Time domain_Chen2020")
plt.show()