In [ ]:
# TODO: move into a file and import
from numpy import array, append
def read_impedance_spectrum(filename):
    '''Read the electrochemical impedance spectrum from a comma-separated file.

    Single line of headers, with resistance in ohm, reactance in ohm, followed 
    by frequency in hertz.
    
    Z(f) = R + jX
    
    Parameters
    ----------
    filename : str
        Comma-separated values (CSV) file that stores the spectrum.

    Returns
    -------
    dict
        The complex impedance and frequency as numpy array objects.
    '''
    frequency = array([], dtype=float)
    impedance = array([], dtype=complex)
    with open(filename) as fin:
        lines = fin.readlines()
        # lines[0] contains the column headers
        for line in lines[1:]:
            (R, X, f) = line.rstrip('\n').split(',')
            frequency = append(frequency, float(f))
            # note that the data has opposite sign for the reactance
            impedance = append(impedance, float(R) - float(X)*1j)
    return {
        'frequency': frequency,
        'impedance': impedance
    }

In [ ]:
experimental_data = read_impedance_spectrum('Maxwell electrodes with TEABF4 electrolyte at open circuit coin cell.csv')

In [ ]:
from pycap import PropertyTree, EnergyStorageDevice
from pycap import ElectrochemicalImpedanceSpectroscopy, NyquistPlot
ptree = PropertyTree()
ptree.put_string('type', 'SeriesRC')
ptree.put_double('series_resistance', 3.3)
ptree.put_double('capacitance', 0.26)
device = EnergyStorageDevice(ptree)
ptree = PropertyTree()
ptree.put_double('dc_voltage', 0)
ptree.put_string('harmonics', '1')
ptree.put_string('amplitudes', '5e-3')
ptree.put_string('phases', '0')
ptree.put_int('steps_per_cycle', 512)
ptree.put_int('cycles', 2)
ptree.put_int('ignore_cycles', 1)
# NOTE: the ``frequency`` argument allows us to compute the same data that in the experiment
eis = ElectrochemicalImpedanceSpectroscopy(ptree, frequencies=experimental_data['frequency'])
eis.run(device)

numerical_data = eis._data

In [ ]:
%matplotlib inline
nyquist = NyquistPlot()
# declare a dummy class that mimics Cap EIS Experiment to call from NyquistPlot
class Dummy:
    def __init__(self, data):
        self._data = data
dummy = Dummy(experimental_data)
nyquist.update(dummy)
nyquist.update(eis)

$\int_0^\infty \frac{|Z - Z_{ref}|}{|Z_{ref}|} d\log(f)$


In [ ]:
# one possible way of computing this
from numpy import trapz, absolute, log
f = experimental_data['frequency']
Z_exp = experimental_data['impedance']
Z_num = numerical_data['impedance']
# minus sign because f goes from upper to lower bound of the range
print(-trapz(absolute(Z_num - Z_exp), log(f)))

In [ ]:
from scipy.optimize import minimize
from numpy import trapz, absolute, log
from pycap import PropertyTree, EnergyStorageDevice
from pycap import ElectrochemicalImpedanceSpectroscopy
def fun(x):
    print(x)

    experimental_data = read_impedance_spectrum('Maxwell electrodes with TEABF4 electrolyte at open circuit coin cell.csv')
    
    ptree = PropertyTree()
    ptree.put_string('type', 'SeriesRC')
    ptree.put_double('series_resistance', x[0])
    ptree.put_double('capacitance', x[1])
    device = EnergyStorageDevice(ptree)

    ptree = PropertyTree()
    ptree.put_double('dc_voltage', 0)
    ptree.put_string('harmonics', '1')
    ptree.put_string('amplitudes', '5e-3')
    ptree.put_string('phases', '0')
    ptree.put_int('steps_per_cycle', 512)
    ptree.put_int('cycles', 2)
    ptree.put_int('ignore_cycles', 1)
    eis = ElectrochemicalImpedanceSpectroscopy(ptree, frequencies=experimental_data['frequency'])
    eis.run(device)

    numerical_data = eis._data

    f = experimental_data['frequency']
    Z_exp = experimental_data['impedance']
    Z_num = numerical_data['impedance']

    return -trapz(absolute(Z_num - Z_exp), log(f))
x0 = [3.3, 0.26]
minimize(fun, x0, bounds=[(0, None), (0, None)])