pyDRTtools

To perform Distribution of Relaxation Times (DRT) analysis on a single Electrochemical Impedance Spectroscopy (EIS) spectrum, follow these steps:

1. Load Data: Load your frequency and impedance data (Real and Imaginary parts) into NumPy arrays.
2. Define Timescales: Create a vector of timescales (tau_vec) using np.logspace over your desired range.
3. Configure Discretization: Choose an rbf_type (e.g., 'Gaussian', 'Piecewise Linear', 'C0 Matern', 'C2 Matern', 'C4 Matern', 'C6 Matern', or 'Inverse Quadratic') and a shape_control method (e.g., 'FWHM Coefficient').
4. Compute Matrices:
   • Use basics.compute_epsilon to find the epsilon parameter.
   • Use basics.assemble_A_re and basics.assemble_A_im to build the discretization matrices for the real and imaginary parts. Note that you must manually augment these matrices to include Ohmic resistance ($R_{\infty}$) and inductance ($L_0$) contributions.
   • Use basics.assemble_M_2 to create the second-order differentiation matrix for ridge regression.
5. Regularization: Find the optimal regularization parameter lambda_value using basics.optimal_lambda with a cross-validation method like 'GCV'.
6. Deconvolution: Use basics.quad_format_combined to prepare the quadratic programming problem, then solve it using cvxopt.solvers.qp to recover the DRT ($\gamma$).

# 1. Setup parameters
rbf_type = 'Gaussian'
shape_control = 'FWHM Coefficient'
coeff = 0.5

# 2. Compute epsilon and discretization matrices
epsilon = basics.compute_epsilon(freq_vec, coeff, rbf_type, shape_control)
A_re = basics.assemble_A_re(freq_vec, tau_vec, epsilon, rbf_type)
# Augment A_re for R_inf and L_0
A_re = np.hstack((np.ones((N_freqs, 1)), np.zeros((N_freqs, 1)), A_re))

A_im = basics.assemble_A_im(freq_vec, tau_vec, epsilon, rbf_type)
# Augment A_im for R_inf and L_0
A_im = np.hstack((np.zeros((N_freqs, 1)), 2*np.pi*freq_vec.reshape((N_freqs, 1)), A_im))

A = np.vstack((A_re, A_im))

# 3. Differentiation matrix
M2 = np.zeros((N_taus+2, N_taus+2))
M2[2:,2:] = basics.assemble_M_2(tau_vec, epsilon, rbf_type)

# 4. Optimal lambda and solve
lambda_value = basics.optimal_lambda(A_re, A_im, np.real(Z_exp), np.imag(Z_exp), M2, -3, 'GCV')
H_combined, c_combined = basics.quad_format_combined(A_re, A_im, np.real(Z_exp), np.imag(Z_exp), M2, lambda_value)

# Solve using cvxopt
from cvxopt import matrix, solvers
# (Assuming G and h are defined for constraints)
sol = solvers.qp(matrix(H_combined), matrix(c_combined), G, h)
x = np.array(sol['x']).flatten()

# Extract results
R_inf_DRT, L_0_DRT = x[0:2]
gamma_DRT = x[2:]

You can generate a batch of EIS datasets by iterating through a list of varying parameters (e.g., varying charge-transfer resistance $R_{ct}$) while keeping other model parameters constant. This is useful for creating datasets for machine learning or parameter sensitivity studies.

The 2xZARC model calculates impedance as: $Z_{total} = R_{\infty} + Z_{ZARC,1} + Z_{ZARC,2}$.

R_inf = 10
R_ct_1_list = [5, 15, 25, 35, 45, 55]
phi_1 = 0.8
tau_1 = 0.1
R_ct_2 = 50
phi_2 = 0.8
tau_2 = 0.0001
sigma_n_exp = 0.5

for n in range(len(R_ct_1_list)):
    Z_exact_1 = R_ct_1_list[n]/(1.+(1j*2.*np.pi*freq_vec*tau_1)**phi_1)
    Z_exact_2 = R_ct_2/(1.+(1j*2.*np.pi*freq_vec*tau_2)**phi_2)
    Z_exact = R_inf + Z_exact_1 + Z_exact_2
    Z_exp = Z_exact + sigma_n_exp*(np.random.normal(0, 1, N_freqs)+1j*np.random.normal(0, 1, N_freqs))
    
    df = pd.DataFrame.from_dict({'Freq': freq_vec, 'Real': np.real(Z_exp), 'Imag': np.imag(Z_exp)})
    df.to_csv(f'2ZARC{n}.csv')

To generate a single Electrochemical Impedance Spectroscopy (EIS) spectrum, you can model the impedance using a ZARC (constant phase element based) model. The total impedance is calculated as the sum of an ohmic resistance ($R_{\infty}$) and one or more ZARC elements.

To simulate realistic experimental data, you can add white noise to the exact impedance using a standard deviation parameter (sigma_n_exp).

Steps:

1. Define model parameters ($R_{\infty}$, $R_{ct}$, $\phi$, $\tau$).
2. Create a log-equispaced frequency vector.
3. Calculate the exact impedance using the formula: $Z_{ZARC} = \frac{R_{ct}}{1 + (j \cdot 2\pi f \tau)^{\phi}}$.
4. Add complex white noise to simulate experimental error.

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd

# 1. Model parameters
R_inf = 10
R_ct_1 = 50
phi_1 = 0.8
tau_1 = 0.1

# 2. Frequencies
log_freq_min, log_freq_max = -2., 6.
N_freqs = 10*int(log_freq_max-log_freq_min)+1
freq_vec = np.logspace(log_freq_min, log_freq_max, num=N_freqs, endpoint=True)

# 3. Exact Impedance (1xZARC)
Z_exact_1 = R_ct_1/(1.+(1j*2.*np.pi*freq_vec*tau_1)**phi_1)
Z_exact = R_inf + Z_exact_1

# 4. Add noise
np.random.seed(12129)
sigma_n_exp = 0.5
Z_exp = Z_exact + sigma_n_exp*(np.random.normal(0, 1, N_freqs)+1j*np.random.normal(0, 1, N_freqs))

# Plotting
plt.plot(np.real(Z_exp), -np.imag(Z_exp), 'o', color='red')