Problem with implementing EIS

[alirahdarian]

Hello,
I am using the following code for generating EIS data. The problem is that I do not have the semicircle in my plot. Can someone please let me know what is the problem of my code?

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

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

#EIS parameters
I = 1 * 1e-3
number_of_periods = 20
samples_per_period = 100

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

parameter_values["Current function [A]"] = current_function
#parameter_values["Contact resistance [Ohm]"] = 0.023
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[-5 * samples_per_period - 1 :]
#print('time=',sol["Time [s]"])
#print('stable time=',time)
current = sol["Current [A]"].entries[-5 * samples_per_period - 1 :]
voltage = sol["Terminal voltage [V]"].entries[-5 * 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)

#Nyquist Plot
print('impedances:', impedances_time)
data = np.array(impedances_time)
print('data after array=', data)
plt.plot(data.real, -data.imag,linestyle='None' , marker='o')
plt.xlabel(r"$Z_\mathrm{Re}$ [Ohm]")
plt.ylabel(r"$-Z_\mathrm{Im}$ [Ohm]")
plt.show()

sim.plot(output_variables=["Current [A]",
"Terminal voltage [V]", 'Battery open-circuit voltage [V]'])

[brosaplanella]

I believe the artifact at high frequencies is due to numerical accuracy. If you increase the current amplitude it goes away (if you wanted to keep small amplitudes you would have to tighten the solver tolerances to improve accuracy). However, I still do not get the semicircle. I am not too familiar with EIS for the SPM model and not sure if the semicircle should be there or is it due to physics not included in SPM.

There is also this side package for EIS: https://github.com/rish31415/pybamm-eis. It is based on a quite old PyBaMM version so not sure if it still works, but worth a try. It uses a better method for EIS.