The Probabilistic Deconvolution of the Distribution of Relaxation Times with Finite Gaussian Processes

Electrochemical impedance spectroscopy (EIS) is a tool widely used to study the properties of electrochemical systems. The distribution of relaxation times (DRT) has emerged as one of the main methods for the analysis of EIS spectra. Gaussian processes can be used to regress EIS data, quantify uncertainty, and deconvolve the DRT, but current implementations do not constrain the DRT to be positive and can only use the imaginary part of EIS spectra. Herein, we overcome both issues by using a finite Gaussian process approximation to develop a new framework called the finite Gaussian process distribution of relaxation times (fGP-DRT). The analysis on artificial EIS data shows that the fGP-DRT method consistently recovers exact DRT from noise-corrupted EIS spectra while accurately regressing experimental data. Furthermore, the fGP-DRT framework is used as a machine learning tool to provide probabilistic estimates of the impedance at unmeasured frequencies. The method is further validated against experimental data from fuel cells and batteries. In short, this work develops a novel probabilistic approach for the analysis of EIS data based on Gaussian process, opening a new stream of research for the deconvolution of DRT.

Deconvolving (log ) from EIS measurements is challenging because an ill-posed inverse problem needs to be solved [1,[39][40][41][42][43][44][45][46][47][48]. Several approaches, including maximum entropy, genetic algorithms, ridge regression, evolutionary algorithms, elastic net, Ridge and Lasso regression, and Bayesian methods have been used [4,41,43,[45][46][47][49][50][51][52][53]. Among them, probabilistic methods have emerged as the method of choice for DRT deconvolution due to their ability to quantify uncertainty. In that context, our research group recently developed the Gaussian process distribution of relaxation times (GP-DRT) method, which assumes the DRT to be a Gaussian process (GP) [1]. The GP-DRT was shown to effectively deconvolve the DRT, fitting both synthetic and real EIS data. One particularly appealing trait of the GP-DRT method is that its hyperparameters can be obtained by maximizing the evidence, ensuring maximal consistency with experiments. The GP-DRT method has, however, two limitations. First, the (log ) estimated using the GP-DRT is unconstrained and, therefore, can be negative (which is unphysical). Second, the current implementation assumes that only the imaginary part of the impedance spectrum can be used [54]. To overcome these limitations, we developed a new model based on a finite approximation of GPs [55]. Herein, this method will be called the finite Gaussian process distribution of relaxation times (fGP-DRT). The fGP-DRT was set up so that it can recover both the real and imaginary parts of the impedance under the constraint that the DRT is non negative, i.e., (log ) ≥ 0. The fGP-DRT method was successfully validated against artificial and real EIS data. Furthermore, like the GP-DRT, the fGP-DRT was used as a machine learning toolbox, allowing the prediction of EIS spectra at unmeasured frequencies. This work extends the interpretation of EIS data through the lens of probabilistic modeling, and will likely inspire more research works that leverage GPs for DRT deconvolution.

Gaussian Process
A Gaussian process (GP) can be loosely understood as an infinite collection of random variables such that any finite subset of these random variables has a joint Gaussian distribution [1,56,57].
After defining = ( re im ) and = ( re 2 im ), we assume that the experimental impedance is given by where = ( 0 , ∞ , ( ) ⊤ ) ⊤ , ∼ ( , 2 ), is the standard deviation of the random error, and is the 2 × 2 identity matrix. We now suppose that where and are the standard deviations of 0 and ∞ , respectively ( 0 ∼ (0, 2 ) and ∞ ∼ (0, 2 )). It follows from (10) and (11) that the joint distribution of and is the following multivariate Gaussian: Using (13), we can compute the distribution of conditioned to the experimental data exp [1,56]: where

Positivity Constraint
As outlined in the introduction, the current GP-DRT model does not impose any constraint on [1].
However, starting from (14), it is straightforward to impose a constraint on | by taking that | is a truncated multinormal distribution bound by zero from below, i.e., To sample from the truncated multinormal distribution (14), Hamiltonian Monte Carlo sampling can be used [58]. We generated 10,000 samples and discarded the first 1,000 as burn-in. The impedance can then be obtained from the sampled by matrix multiplication, i.e., = . In the figures in the remainder of the article, the means of DRT and impedance are shown as a black line, while the 3 × σ credible intervals are shown as gray regions.

Results
We first validated the consistency of the fGP-DRT model against artificial experiments generated by corrupting exact EIS spectra characterized by known DRTs. Then, the model was tested against real EIS data of fuel cells and batteries.

Artificial experiments
All artificial EIS data were generated, unless otherwise specified, using log-spaced frequencies in the 10 −4 − 10 4 Hz range with 10 points per decade. We first considered a single ZARC model, for which we analyzed the influence of the experimental error, the number of collocation points (N as defined in equation (3)), and truncated spectra. After that, we studied models with overlapping (2×ZARC) and discontinuous timescales (piecewise constant and fractal elements).

Single ZARC Model
A single ZARC circuit is an ohmic resistor ( ∞ ) in series with a parallel circuit consisting of a resistor ( ct ) and a constant phase element (CPE), with characteristic time 0 and parameter .
The impedance, exact ( ), and its corresponding DRT are given by [1,4] and the parameters used in the model are reported in Table 1. For the synthetic experiments, we generated data using (12) where the error standard deviation is given by Nyquist plots of the exact and noise-corrupted impedance spectra are shown in Figure 2  We also considered a single ZARC plus an inductor (values of the parameters in Table 1), whose Nyquist and Bode plots for both unconstrained and constrained DRT are reported in Figure S1 of the SI. We observe that the unconstrained fGP-DRT oscillates about the exact DRT with widening error bands for > 10 3 s ( Figure S1(b)). Instead, the constrained fGP-DRT model exhibits a better recovery of (log ). The real and imaginary components of the impedance obtained using the constrained fGP-DRT match well with the exact impedance spectrum.  Figure S2. The matrix is nearly diagonal for ℓ = 10 −2 , indicating that a value of (log ) at a particular log is only influenced by adjacent log timescales. In contrast, we note in Figure S2 a strong dependency across a broader ranges of timescales as ℓ increases. The results reported in Figure S3, which shows the recovered real and imaginary components of the impedance, the estimated DRT and the values of the entries in the matrix for small (ℓ = 10 −5 ) and moderately high (ℓ = 2.5 and 5.0) length scales, also support these conclusions and highlight the strong influence of the kernel and, in particular, the correlation log-timescale ℓ on the recovered DRT. In the remainder of the article, the hyperparameters are obtained by optimizing the experimental evidence, that is the probability of the experimental data given the model.

Influence of the Experimental Errors
We now illustrate the influence of the experimental noise on the predictions produced by the fGP- Using equation (10) artificial experiments (1,000 for each exp ).

Data Truncation
We tested how well the constrained fGP-DRT model deconvolves the DRT when the EIS data is For further validation, we considered the ECM described in [61] with the data being truncated below 10 −3 Hz. The specific parameters of this ECM are reported in Table S2. Our model successfully recovers the real and imaginary components of the impedance. Additionally, the DRT is estimated in a satisfactory manner together with its error bands ( Figure S5). Hz ( exp = 0.5 Ω) relative to the data presented in Figure 2 with the corresponding recovered DRT in (d), (e), and (f), and real and imaginary parts of the impedance in (g), (h), and (i).

Error of Estimation and Discrepancies
We also explored the influence of the number of collocation points, N in (3) where ̂ is the value of the mean DRT obtained by sampling (16) using a synthetic EIS spectrum.
We generated 100 artificial experiments, 10 each for N = 20, 40, 60, …, 200, and computed the normalized squared residuals 2 for each. Figure 7(a) shows the boxplot of 2 as a function of N.
As N increases, the descrepancies significantly reduce, evidenced by the smaller box for higher values of N, leading to an improved DRT recovery (Figure 7(b)). Specifically, the curve of the recovered DRT becomes smoother as N increases ( Figure S6). A representative recovered DRT for N = 200 ( 2 = 8.25 × 10 −5 ) is shown in Figure 7

2×ZARC Model
In many electrochemical systems, numerous physical processes take place concurrently.
Therefore, we investigated the capability of the fGP-DRT model to capture overlapping features by studying a 2×ZARC model. The DRT and impedance response are given by [1,5,49]  where the specific parameters used are reported in Table S1. We set the value 1 = 0.1 s and investigated two values for 2 , namely 2 = 1 and 10 s. Figure 8 shows the DRT, real, and imaginary parts of the impedance. For separated timescales with 1 = 0.1 s and 2 = 10 s, the Nyquist plot consists in two semicircles (Figure 8(a)). The DRT recovered using the fGP-DRT method correctly shows two distinct peaks (Figure 8(c)). Similarly, for 1 = 0.1 s and 2 = 1.0 s, the Nyquist plot of the impedance resembles a single semicircle with close peaks in the DRT as shown in Figure 8(b) and (d), respectively. The impedance is well regressed as shown in Figure   8

Piecewise Constant and Fractal Models
We tested how the fGP-DRT model handles discontinuities in the DRT. To that end, the piecewise constant (PWC) [36,63]  Corresponding DRT recovered using the fGP-DRT method under the non-negativity constraint (c) and (d), and real and imaginary parts of the impedance (e) and (f).

Real Experiments
Having tested the consistency of the fGP-DRT model with synthetic experiments, we evaluated its performance against real data from three symmetric solid oxide fuel cells and two batteries.

Ba0.95La0.05FeO3--based Symmetric Solid Oxide Cells
EIS data was obtained from symmetric cells with Ba0.95La0.05FeO3- (BLF) [64] as the electrode material and samarium-doped ceria as the electrolyte. The tests were performed in the frequency range from 100 mHz to 20 kHz with five points per decade using a VSP potentiostat (BioLogic).
The symmetric cell was tested at 500 and 550 °C in an atmosphere consisting of N2 and O2 with an oxygen partial pressure pO2 = 60%. The experimental data for each temperature was regressed against a 2×ZARC (the best-fitting parameters are reported in Table S3). The regressed ECM and deconvolved fGP-DRT model match closely, as shown in Figure 10 Additionally, we considered the EIS of a symmetric cell with samarium-doped ceria as the electrolyte and the electrode consisting of a composite containing BLF and Ag2O (molar ratio 20:1) [65]. EIS data was collected at 550 and 700 °C in a mixture of N2 and O2 with pO2 = 60%.
We found that the fGP-DRT is consistent with the ECM model ( Figure S7), except for a small deviation for > 10 s ( Figure S7(d)).

Sr0.9Ce0.1Fe0.8Ni0.2O3--based Symmetric Protonic Ceramic Cells
We analyzed the EIS data from a symmetric cell with Sr0.9Ce0.1Fe0.8Ni0.2O3- (SCFN) as the electrode material and BaZr0.1Ce0.7Y0.2-x O3- as the electrolyte. The cell was tested at 500 °C in two atmospheres, namely 97% Air-3% H 2 O and 94% Air-6% H 2 O [67]. The data was obtained for frequencies ranging from 0.01 Hz to 200 kHz and a 10 mV amplitude of the input voltage. The measured EIS data was regressed against a 2×ZARC ECM whose parameters are reported in Table   S4. The Nyquist plots of the SCFN with the deconvolved fGP-DRT model and the regressed ECM are presented in Figure 11

Composite Polymer Electrolytes
We analyzed batteries with Li-metal as the anode, a composite polymer electrolyte, and LiFePO4 as the cathode [68,69]. The EIS spectrum was regressed using the fGP-DRT model and a 3×ZARC ECM (the obtained ECM parameters are reported in Table S5). The Nyquist plot, recovered DRT with the regressed ECM, and real and imaginary parts of the impedance are displayed in the left panel of Figure S8. We observe that the DRT from the fGP-DRT and the ECM are consistent and that the fGP-DRT matches closely the experimental impedance.

Solid-like Dual-salt Polymer Electrolytes
We used experimental data obtained from a battery with the same anode and cathode as used in composite polymer electrolyte (see Section 3.2.2.1 above), and a solid-like dual-salt polymer electrolyte [67][68][69][70]. EIS experiments were carried out for frequencies between 1 Hz and 7 MHz.
We repeated the same procedure as described Section 3.2.2.1, the results are displayed in the right panel of Figure S8. As above, the DRT obtained using the ECM and fGP-DRT method are consistent and the experimental impedance is closely regressed by the fGP-DRT method.

Conclusion
In this work, we developed a novel fGP-DRT model based on a finite GP approximation that accurately deconvolves the DRT from EIS data. The new model inherits several traits of the previously developed GP-DRT since: 1) it assumes that the DRT is a GP; 2) analysis of synthetic and real experiments shows that this method is consistent; 3) its hyperparameters can be selected by maximizing the experimental evidence; and 4) the fGP-DRT is robust against experimental noise and data truncation. In addition, the fGP-DRT model outperforms the GP-DRT as it can: 1) be constrained to produce only non-negative DRTs; and 2) use both the real and imaginary parts of the impedance. In short, this article develops a new method for the probabilistic analysis of EIS spectra, opening up new research avenues that leverage finite GPs for DRT deconvolution.