import numpy as np
from impedance.models.circuits.fitting import rmse
from impedance.models.circuits.elements import circuit_elements, K # noqa
[docs]def linKK(f, Z, c=0.85, max_M=50, fit_type='real', add_cap=False):
""" A method for implementing the Lin-KK test for validating linearity [1]
Parameters
----------
f: np.ndarray
measured frequencies
Z: np.ndarray of complex numbers
measured impedances
c: np.float
cutoff for mu
max_M: int
the maximum number of RC elements
fit_type: str
selects which components of data are fit ('real', 'imag', or
'complex')
add_cap: bool
option to add a serial capacitance that helps validate data with no
low-frequency intercept
Returns
-------
M: int
number of RC elements used
mu: np.float
under- or over-fitting measure
Z_fit: np.ndarray of complex numbers
impedance of fit at input frequencies
resids_real: np.ndarray
real component of the residuals of the fit at input frequencies
resids_imag: np.ndarray
imaginary component of the residuals of the fit at input frequencies
Notes
-----
The lin-KK method from Schönleber et al. [1] is a quick test for checking
the
validity of EIS data. The validity of an impedance spectrum is analyzed by
its reproducibility by a Kramers-Kronig (KK) compliant equivalent circuit.
In particular, the model used in the lin-KK test is an ohmic resistor,
:math:`R_{Ohm}`, and :math:`M` RC elements.
.. math::
\\hat Z = R_{Ohm} + \\sum_{k=1}^{M} \\frac{R_k}{1 + j \\omega \\tau_k}
The :math:`M` time constants, :math:`\\tau_k`, are distributed
logarithmically,
.. math::
\\tau_1 = \\frac{1}{\\omega_{max}} ; \\tau_M = \\frac{1}{\\omega_{min}}
; \\tau_k = 10^{\\log{(\\tau_{min}) + \\frac{k-1}{M-1}\\log{{(
\\frac{\\tau_{max}}{\\tau_{min}}}})}}
and are not fit during the test (only :math:`R_{Ohm}` and :math:`R_{k}`
are free parameters).
In order to prevent under- or over-fitting, Schönleber et al. propose using
the ratio of positive resistor mass to negative resistor mass as a metric
for finding the optimal number of RC elements.
.. math::
\\mu = 1 - \\frac{\\sum_{R_k \\ge 0} |R_k|}{\\sum_{R_k < 0} |R_k|}
The argument :code:`c` defines the cutoff value for :math:`\\mu`. The
algorithm starts at :code:`M = 3` and iterates up to :code:`max_M` until a
:math:`\\mu < c` is reached. The default of 0.85 is simply a heuristic
value based off of the experience of Schönleber et al., but a lower value
may give better results.
If the argument :code:`c` is :code:`None`, then the automatic determination
of RC elements is turned off and the solution is calculated for
:code:`max_M` RC elements. This manual mode should be used with caution as
under- and over-fitting should be avoided.
[1] Schönleber, M. et al. A Method for Improving the Robustness of
linear Kramers-Kronig Validity Tests. Electrochimica Acta 131, 20–27 (2014)
`doi: 10.1016/j.electacta.2014.01.034
<https://doi.org/10.1016/j.electacta.2014.01.034>`_.
"""
if c is not None:
M = 0
mu = 1
while mu > c and M < max_M:
M += 1
ts = get_tc_distribution(f, M)
elements, mu = fit_linKK(f, ts, M, Z, fit_type, add_cap)
if M % 10 == 0:
print(M, mu, rmse(eval_linKK(elements, ts, f), Z))
else:
M = max_M
ts = get_tc_distribution(f, M)
elements, mu = fit_linKK(f, ts, M, Z, fit_type, add_cap)
Z_fit = eval_linKK(elements, ts, f)
resids_real = residuals_linKK(elements, ts, Z, f, residuals='real')
resids_imag = residuals_linKK(elements, ts, Z, f, residuals='imag')
return M, mu, Z_fit, resids_real, resids_imag
[docs]def get_tc_distribution(f, M):
""" Returns the distribution of time constants for the linKK method """
t_max = 1/(2 * np.pi * np.min(f))
t_min = 1/(2 * np.pi * np.max(f))
ts = np.zeros(shape=(M,))
ts[0] = t_min
ts[-1] = t_max
if M > 1:
for k in range(2, M):
ts[k-1] = 10**(np.log10(t_min) +
((k-1)/(M-1))*np.log10(t_max/t_min))
return ts
[docs]def fit_linKK(f, ts, M, Z, fit_type='real', add_cap=False):
""" Fits the linKK model using linear regression
Parameters
----------
f: np.ndarray
measured frequencies
ts: np.ndarray
logarithmically spaced time constants of RC elements
M: int
the number of RC elements
Z: np.ndarray of complex numbers
measured impedances
fit_type: str
selects which components of data are fit ('real', 'imag', or
'complex')
add_cap: bool
option to add a serial capacitance that helps validate data with no
low-frequency intercept
Returns
-------
elements: np.ndarray
values of fit :math:`R_k` in RC elements and series :math:`R_0`,
L, and optionally C.
mu: np.float
under- or over-fitting measure
Notes
-----
Since we have a system of equations, :math:`Ax ~= b`, that's linear wrt
:math:`R_k`, we can fit the model by calculating the pseudo-inverse of A.
:math:`Ax` is our model fit, :math:`\\hat{Z}`, and :math:`b` is the
normalized real or imaginary component of the impedance data,
:math:`Re(Z)/|Z|` or :math:`Im(Z)/|Z|`, respectively.
:math:`\\hat{Z} = R_0 + \\sum^M_{k=1}(R_k / |Z|(1 + j * w * \\tau_k))`.
:math:`x` is an (M+1) :math:`\\times` 1 matrix where the first row
contains :math:`R_0` and subsequent rows contain :math:`R_k` values.
A is an N :math:`\\times` (M+1) matrix, where N is the number of data
points, and M is the number of RC elements.
Examples
--------
Fitting the real part of data, the first column of A contains
values of :math:`\\frac{1}{|Z|}`, the second column contains
:math:`Re(1 / |Z| (1 + j * w * \\tau_1))`, the third contains
:math:`Re(1 / |Z| (1 + j * w * \\tau_2))` and so on. The :math:`R_k` values
within the x matrix are found using :code:`numpy.linalg.pinv` when
fit_type = 'real' or 'imag'. When fit_type = 'complex' the coefficients are
found "manually" using :math:`r = ||A'x - b'||^2 + ||A''x - b'||^2`
according to Eq 14 of Schonleber [1].
[1] Schönleber, M. et al. A Method for Improving the Robustness of
linear Kramers-Kronig Validity Tests. Electrochimica Acta 131, 20–27 (2014)
`doi: 10.1016/j.electacta.2014.01.034
<https://doi.org/10.1016/j.electacta.2014.01.034>`_.
"""
w = 2 * np.pi * f
# Fitting model has M RC elements plus 1 series resistance and 1 series
# inductance
a_re = np.zeros((f.size, M+2))
a_im = np.zeros((f.size, M+2))
if add_cap:
a_re = np.zeros((f.size, M+3))
a_im = np.zeros((f.size, M+3))
# Column for series capacitance. Real part = 0.
a_im[:, -2] = - 1 / (w * np.abs(Z))
# Column for series resistance, R_o in model. Imaginary part = 0.
a_re[:, 0] = 1 / np.abs(Z)
# Column for series inductance to capture inevitable contributions from
# the measurement system. Real part = 0.
a_im[:, -1] = w / np.abs(Z)
# Columns for series RC elements
for i, tau in enumerate(ts):
a_re[:, i+1] = K([1, tau], f).real / np.abs(Z)
a_im[:, i+1] = K([1, tau], f).imag / np.abs(Z)
if fit_type == 'real':
elements = np.linalg.pinv(a_re).dot(Z.real / np.abs(Z))
# After fitting real part, need to use imaginary component of fit to
# find values of series inductance and capacitance
a_im = np.zeros((f.size, 2))
a_im[:, -1] = w / np.abs(Z)
if add_cap:
a_im[:, -2] = -1 / (w * np.abs(Z))
elements[-2] = 1e-18 # nullifies series C without dividing by 0
Z_fit_re = eval_linKK(elements, ts, f)
coefs = np.linalg.pinv(a_im).dot((Z.imag - Z_fit_re.imag) / np. abs(Z))
if add_cap:
elements[-2:] = coefs
else:
elements[-1] = coefs[-1]
elif fit_type == 'imag':
elements = np.linalg.pinv(a_im).dot(Z.imag / np.abs(Z))
# Calculates real part of impedance from fitting to imaginary parts
# without ohmic resistance, i.e. only the real parts of series RC
# elements.
z_re = eval_linKK(elements, ts, f)
# Weighting used in Boukamp et al - "A Linear Kronig-Kramers
# Transform for Immittance Data Validation" 1995, J. Electrochem
# Soc. 142 (6)
ws = 1 / (Z.real**2 + Z.imag**2)
# Finds ohmic resistance for imaginary part fit according to Eq 7
# of Boukamp et al.
elements[0] = np.sum(ws * (Z.real - z_re.real)) / np.sum(ws)
elif fit_type == 'complex':
# x = (A*•A)^-1
# y = A*•b
# Pseudoinsverse, A^+ = (A*•A)^-1•A* and R_k values = A^+•b
x = np.linalg.inv(a_re.T.dot(a_re) + a_im.T.dot(a_im))
y = a_re.T.dot(Z.real / np.abs(Z)) + a_im.T.dot(Z.imag / np.abs(Z))
elements = x.dot(y)
else:
raise ValueError("Invalid choice of fit_type, please choose from "
"\'real\', \'imag\', or \'complex\'")
if add_cap:
mu = calc_mu(elements[1:-2])
else:
mu = calc_mu(elements[1:-1])
return elements, mu
[docs]def eval_linKK(elements, ts, f):
""" Builds a circuit of RC elements to be used in LinKK """
circuit_string = 's([R({},{}),'.format([elements[0]], f.tolist())
for (Rk, tk) in zip(elements[1:], ts):
circuit_string += f'K({[Rk, tk]},{f.tolist()}),'
circuit_string += 'L({},{}),'.format([elements[-1]], f.tolist())
if elements.size == (ts.size + 3):
circuit_string += 'C({},{}),'.format([1/elements[-2]], f.tolist())
circuit_string = circuit_string.strip(',')
circuit_string += '])'
return eval(circuit_string, circuit_elements)
[docs]def residuals_linKK(elements, ts, Z, f, residuals='real'):
""" Calculates the residual between the data and a LinKK fit """
err = (Z - eval_linKK(elements, ts, f))/np.abs(Z)
if residuals == 'real':
return err.real
elif residuals == 'imag':
return err.imag
elif residuals == 'both':
z1d = np.zeros(Z.size*2, dtype=np.float64)
z1d[0:z1d.size:2] = err.real
z1d[1:z1d.size:2] = err.imag
return z1d
[docs]def calc_mu(Rs):
""" Calculates mu for use in LinKK """
neg_sum = sum(abs(x) for x in Rs if x < 0)
pos_sum = sum(abs(x) for x in Rs if x >= 0)
return 1 - neg_sum/pos_sum