Analytical–Computational Integration of Equivalent Circuit Modeling, Hybrid Optimization, and Statistical Validation for Electrochemical Impedance Spectroscopy

1. Introduction

Electrochemical Impedance Spectroscopy (EIS) is a frequency-domain technique that has become indispensable for the quantitative characterization of electrochemical interfaces and devices [1]. Unlike classical DC methods, EIS applies a small-amplitude sinusoidal perturbation and measures the resulting current response over several orders of magnitude in frequency, thereby resolving the interplay between resistive and reactive contributions at different characteristic timescales [2]. Operating under small-signal, linear, time-invariant conditions, EIS preserves the DC operating point while resolving concurrent faradaic (charge-transfer) and mass-transport phenomena across multiple decades of frequency via separation of characteristic time constants [3].

The complex impedance spectrum, $Z\left(\omega \right)$, encapsulates the system’s dynamic response and can be interpreted by mapping it to an equivalent circuit model in which each element corresponds to a specific physicochemical process [4,5].Classical models, such as the Randles circuit, represent charge-transfer resistance ($R_{ct}$) in parallel with the double-layer capacitance ($C_{dl}$), in series with the uncompensated solution resistance ($R_s$) [6]. More advanced configurations incorporate non-ideal capacitive behavior via constant-phase elements (CPE) or account explicitly for semi-infinite diffusion through the Warburg impedance ($Z_W$), thereby extending the range of electrochemical systems that can be faithfully described [7].

Despite these advances, recent literature still exhibits several unresolved limitations that hinder the reproducibility and reliability of EIS analysis: (i) the absence of a unified framework that integrates analytical derivation, robust parameter estimation, and statistical validation within a single workflow [8]; (ii) limited evaluation of model performance under controlled noise conditions, despite its direct relevance to experimental uncertainty [9]; (iii) scarce incorporation of Kramers–Kronig compliance testing as a systematic physical validity check  [10]; and (iv) insufficient use of quantitative model selection criteria such as Akaike or Bayesian information criteria for objective discrimination among competing equivalent circuits. These gaps often lead to parameter estimates with low reproducibility and model choices driven by heuristic or subjective preference rather than rigorous, quantitative evidence [11].

In this context, we propose an analytical–computational framework that unifies exact derivations of impedance functions for seven representative equivalent circuits with a hybrid fitting pipeline combining global and local optimization. The approach incorporates strict physical parameter bounds, analytical Jacobians for improved numerical stability, bootstrap-based uncertainty quantification, and automated validation against the Kramers–Kronig relations to ensure physical consistency. By systematically evaluating all models under controlled noise conditions, we identify their statistical and physical domains of applicability, thereby providing a reproducible methodology for rigorous EIS model selection.

1.1. Equivalent Circuit Foundations and Phasor Analysis

A rigorous interpretation of Electrochemical Impedance Spectroscopy (EIS) data demands an explicit understanding of how resistive, capacitive, and pseudo-capacitive elements interact within an equivalent circuit [12]. In this context, phasor analysis provides a compact mathematical framework that maps time-domain integro-differential relations into algebraic expressions in the frequency domain, enabling closed-form derivations of complex impedance functions [13].

1.1.1. Parallel Impedance Combination

For elements in parallel, the potential drop is identical across all branches, while the total current is additive [14]. This relation is elegantly expressed through the concept of admittance $Y=1/Z$  [15], yielding the equivalent impedance:

$$Z_{\mathrm{eq}}={\left(\sum _{k=1}^N{\displaystyle \frac{1}{Z_k}}\right)}^{-1}$$

(1)

This formulation is central to archetypal EIS models, such as the Randles circuit, where the charge-transfer resistance ($R_{ct}$) and the double-layer capacitance ($C_{dl}$) or a Constant Phase Element (CPE) operate in parallel [16].

1.2. Charge-Transfer Resistance (${\mathbf{R}}_{\mathbf{ct}}$)

The charge-transfer resistance quantifies the kinetic barrier for electron transfer across the electrode–electrolyte interface. Unlike the solution resistance $R_s$, $R_{ct}$ reflects interfacial phenomena and is modulated by electrode microstructure, activation energy, reactant concentration, and applied overpotential

Linearization of the Butler–Volmer equation under small-signal perturbations ($\eta \ll 25$ mV at 298 K) [17] leads to:

$$R_{ct}={\displaystyle \frac{RT}{nFA{j}_0}}$$

(2)

establishing $R_{ct}$ as a frequency-independent kinetic parameter [18]. In a Nyquist diagram, it is extracted from the mid-frequency semicircle diameter [19]:

$$R_{ct}=Z^{\prime }(\omega \to 0)-Z^{\prime }(\omega \to \infty )$$

(3)

These principles form the analytical backbone for the closed-form derivations presented herein, spanning from the classical Randles topology to advanced non-ideal and diffusion-impedance configurations, and directly supporting the quantitative model selection and validation framework developed in this work.

2. Equivalent Circuit Modeling and Analytical Simulation

The interpretation of Electrochemical Impedance Spectroscopy (EIS) data fundamentally relies on fitting equivalent circuit models, which represent physicochemical phenomena at the electrode–electrolyte interface using ideal and non-ideal electrical elements [20]. The primary components include the solution resistance $R_s$, associated with the ionic conductivity of the electrolyte, separators, and electrical contacts; the charge-transfer resistance $R_{ct}$, which quantifies the interfacial kinetics of the redox reaction; and the ideal double-layer capacitance $C_{dl}$  [21].

The constant-phase element (CPE) is parameterized by Q and n, modeling a non-ideal capacitive behavior with impedance $Z_{\mathrm{CPE}}={\left[Q{\left(j\omega \right)}^n\right]}^{-1}$. To describe mass–transport phenomena, the Warburg impedance $Z_W$ for semi-infinite diffusion is employed, defined as:

$$Z_W\left(\omega \right)={\displaystyle \frac{A_W(1-j)}{\sqrt{2\omega }}}.$$

(4)

Regarding impedance combination rules, in a series configuration the total impedance is $Z_{\mathrm{total}}={\sum }_i{Z}_i$, whereas in a parallel configuration it is obtained as $Z_{\mathrm{total}}={\left({\sum }_i{Y}_i\right)}^{-1}$, where $Y_i=1/Z_i$ denotes the admittance of the i-th element.

2.1. Series RLC (R–L–C)

The series RLC topology, comprising an ohmic resistor (R), an inductor (L), and a capacitor (C) connected in series, represents the canonical form of a lumped-element system with purely discrete resistive, inductive, and capacitive contributions. Although not directly descriptive of interfacial electrochemical processes, it constitutes a benchmark configuration for validating the analytical–computational methodology under controlled noise, owing to its closed-form solution and absence of distributed or diffusive elements [21,22].

The complex impedance is expressed as:

$$Z_{\mathrm{total}}\left(\omega \right)=R+j\omega L-{\displaystyle \frac{j}{\omega C}},$$

(5)

where $\omega =2\pi f$ denotes the angular frequency. The real and imaginary parts are:

$$\begin{array}{cc}\hfill Z^{\prime }\left(\omega \right)& =R,\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill Z^{\prime \prime }\left(\omega \right)& =\omega L-{\displaystyle \frac{1}{\omega C}}.\hfill \end{array}$$

(7)

In the Nyquist representation, the spectrum collapses into a vertical line intersecting the real axis at R, with a positive imaginary component ($Z^{\prime \prime }>0$) in the inductive regime ($\omega \gg {\omega }_0$) and a negative imaginary component ($Z^{\prime \prime }<0$) in the capacitive regime ($\omega \ll {\omega }_0$). The resonant frequency,

$${\omega }_0={\displaystyle \frac{1}{\sqrt{LC}}},$$

(8)

marks the transition point where $Z^{\prime \prime }\left({\omega }_0\right)=0$ and the phase angle passes through $0^{\circ }$.

In the Bode magnitude plot, $\left|Z\right|$ exhibits a well-defined minimum at ${\omega }_0$, corresponding to the cancellation of inductive and capacitive reactances. The Bode phase transitions monotonically from $+{90}^{\circ }$ (purely inductive response) at high frequency to $-{90}^{\circ }$ (purely capacitive response) at low frequency, with the ohmic contribution R governing the real-axis intercept in all regimes.

Relevance to this work: This idealized model provides a stringent test of the proposed fitting and validation framework, enabling direct assessment of numerical stability, parameter identifiability, and noise sensitivity without the confounding effects of dispersion or mass transport.

2.2. Classic Randles ($R_s$–$[R_{ct}\Vert C_{dl}]$)

This circuit consists of $R_s$ in series with the parallel combination $R_{ct}\Vert C_{dl}$, representing charge-transfer resistance coupled with an ideal double-layer capacitance. The impedance of the parallel branch is:

$$Z_p\left(\omega \right)={\displaystyle \frac{R_{ct}}{1+j\omega C_{dl}{R}_{ct}}}.$$

(9)

Separating the real and imaginary components yields:

$$\begin{array}{cc}\hfill Z^{\prime }\left(\omega \right)& =R_s+{\displaystyle \frac{R_{ct}}{1+{\left(\omega C_{dl}{R}_{ct}\right)}^2}}\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill Z^{\prime \prime }\left(\omega \right)& =-{\displaystyle \frac{\omega C_{dl}{R}_{ct}^2}{1+{\left(\omega C_{dl}{R}_{ct}\right)}^2}}\hfill \end{array}$$

(11)

Expected behavior: In the Nyquist plot, a perfect semicircle of diameter $R_{ct}$ is obtained; in the Bode phase plot, a minimum occurs at ${\omega }_{\max }=1/\left(R_{ct}{C}_{dl}\right)$.

2.3. Randles + CPE ($R_s$–$[R_{ct}\Vert \mathrm{CPE}]$)

Here, $C_{dl}$ is replaced by a constant-phase element (CPE):

$$Z_{\mathrm{CPE}}\left(\omega \right)={\displaystyle \frac{1}{Q{\left(j\omega \right)}^n}}.$$

(12)

The parallel branch impedance becomes:

$$Z_p\left(\omega \right)={\displaystyle \frac{R_{ct}}{1+R_{ct}Q{\left(j\omega \right)}^n}}.$$

(13)

Effect: The Nyquist semicircle becomes depressed, with the degree of depression proportional to $(1-n)$, indicating surface heterogeneity [23,24].

2.4. Randles + Warburg [$R_s+(R_{ct}\Vert C_{dl})+Z_W$]

The Randles + Warburg configuration extends the classical Randles circuit by incorporating a semi-infinite diffusion element ($Z_W$) in series with the faradaic branch $R_{ct}\Vert C_{dl}$. This modification enables the simultaneous representation of charge-transfer processes and diffusion limitations that emerge at low frequencies in electrochemical systems governed by mass transport.

The total impedance is expressed as:

$$Z_{\mathrm{total}}\left(\omega \right)=R_s+\left[{\displaystyle \frac{1}{\frac{1}{R_{ct}}+j\omega C_{dl}}}\right]+Z_W\left(\omega \right),$$

(14)

where $Z_W\left(\omega \right)$ denotes the semi-infinite diffusion Warburg impedance:

$$Z_W\left(\omega \right)={\displaystyle \frac{A_W(1-j)}{\sqrt{2\omega }}},$$

(15)

with $A_W$ being the Warburg coefficient and $\omega$ the angular frequency.

Physical interpretation:

• $R_s$: ohmic resistance associated with the electrolyte and electrical contacts.
• $R_{ct}$: charge-transfer resistance at the electrode–electrolyte interface.
• $C_{dl}$: double-layer capacitance.
• $Z_W$: diffusive contribution describing the semi-infinite transport of reactive species.

In the Nyquist plot, this model produces a semicircle at intermediate frequencies (dominated by $R_{ct}\Vert C_{dl}$), followed by a diffusion tail with an approximate slope of ${45}^{\circ }$ at low frequencies. In the Bode plot, the impedance magnitude exhibits a slope of $-0.5$ in the logarithmic scale at low frequencies, while the phase approaches $-{45}^{\circ }$, which constitutes the characteristic signature of a semi-infinite diffusion [25].

2.5. $(R_{ct}+Z_W)\Vert C_{dl}$

In this model, the charge-transfer resistance $R_{ct}$ and the semi-infinite Warburg impedance $Z_W\left(\omega \right)$[ ] (see Equation (4)) are connected in series before entering into parallel with the double-layer capacitance $C_{dl}$:

$$\begin{array}{cc}\hfill Z_{\mathrm{serie}}\left(\omega \right)& =R_{ct}+Z_W\left(\omega \right)\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill Z_{\mathrm{total}}\left(\omega \right)& =R_s+{\left[{\displaystyle \frac{1}{Z_{\mathrm{serie}}\left(\omega \right)}}+j\omega C_{dl}\right]}^{-1}\hfill \end{array}$$

(17)

This topology is particularly relevant for systems in which mass transport of electroactive species is directly coupled to charge-transfer kinetics, such as in porous electrodes or materials with restricted ionic transport.

2.6. $(R_{ct}+Z_W)\Vert \mathrm{CPE}$

This model is a variant of Model , in which the double-layer capacitance $C_{dl}$ is replaced by a constant-phase element (CPE) [16] characterized by parameters Q and n. The semi-infinite Warburg impedance $Z_W\left(\omega \right)$ is defined as in Equation (4):

$$\begin{array}{cc}\hfill Z_{\mathrm{total}}\left(\omega \right)& =R_s+{\left[{\displaystyle \frac{1}{R_{ct}+Z_W\left(\omega \right)}}+Q{\left(j\omega \right)}^n\right]}^{-1}\hfill \end{array}$$

(18)

Application: Suitable for describing rough or heterogeneous interfaces, as well as systems exhibiting a distribution of relaxation times where the diffusion of electroactive species is coupled to a non-ideal capacitive response.

2.7. Complete Randles (CPE + Warburg) [$R_s$–$[R_{ct}\Vert (\mathrm{CPE}+Z_W)]$]

This is the most general topology considered in the present methodology, in which the constant-phase element (CPE) and the semi-infinite Warburg impedance  [16] $Z_W\left(\omega \right)$ (see Equation (4)) are connected in series, and this branch is placed in parallel with the charge-transfer resistance $R_{ct}$, all in series with the solution resistance $R_s$:

$$\begin{array}{cc}\hfill Z_{\mathrm{ser}}\left(\omega \right)& ={\displaystyle \frac{1}{Q{\left(j\omega \right)}^n}}+Z_W\left(\omega \right)\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill Z_{\mathrm{p}}\left(\omega \right)& ={\displaystyle \frac{R_{ct}{Z}_{\mathrm{ser}}\left(\omega \right)}{R_{ct}+Z_{\mathrm{ser}}\left(\omega \right)}}\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill Z_{\mathrm{total}}\left(\omega \right)& =R_s+Z_{\mathrm{p}}\left(\omega \right)\hfill \end{array}$$

(21)

Application: This model simultaneously accounts for (i) non-ideal capacitive behavior due to surface heterogeneity or microscopic roughness, (ii) diffusive effects of electroactive species under semi-infinite transport conditions, and (iii) finite charge-transfer kinetics. It is particularly relevant for complex electrochemical systems such as composite electrodes, porous electrodes coated with non-uniform catalysts, and energy storage devices with coupled ionic and electronic transport.

Comparative Summary

Practical relevance: Models 1–2 are suited for systems dominated by charge-transfer processes with ideal or mildly non-ideal capacitive behavior (e.g., smooth metal electrodes, thin-film catalysts). Models 3–4 are appropriate for electrodes where diffusion plays a significant role, such as porous electrodes or systems with semi-infinite mass transport. Models 5–6 address scenarios with coupled charge-transfer and diffusion under non-ideal double-layer conditions, typical of heterogeneous, rough, or composite electrodes. Model 7 captures the most complex interplay of non-ideal capacitance, diffusion, and finite kinetics, enabling the most general representation of real electrochemical interfaces [26].

3. Methodology

A robust, fully modular computational platform was developed in Python to simulate, fit, and validate electrochemical impedance spectroscopy (EIS) models with complete traceability and rigorous statistical control. The environment automates the entire simulate–fit–evaluate pipeline for seven equivalent circuit topologies, integrating hybrid optimization, non-parametric bootstrapping, physicochemical consistency checks, and comprehensive visualization.

3.1. 1. Analytical Model Definition and Jacobians

Each circuit is implemented as a complex impedance function $Z\left(\omega \right)$ in NumPy, together with its analytical Jacobian $\partial Z/\partial p_i$. The seven models considered are:

(1)

Series RLC: $Z=R+j\omega L+1/\left(j\omega C\right)$

(2)

Classical Randles: $R_s+(R_{ct}\Vert C_{dl})$

(3)

Randles + CPE: $R_s+\left(R_{ct}\Vert {\left[Q{\left(j\omega \right)}^n\right]}^{-1}\right)$

(4)

Randles + Warburg: $R_s+(R_{ct}\Vert C_{dl})+\sigma /\sqrt{\omega }(1-j)$

(5)

$(R_{ct}+Z_W)\Vert C_{dl}$

(6)

$(R_{ct}+Z_W)\Vert \mathrm{CPE}$

(7)

Full Randles: $R_s+\left[R_{ct}\Vert \left({\displaystyle \frac{1}{Q{\left(j\omega \right)}^n}}+{\displaystyle \frac{\sigma }{\sqrt{\omega }}}(1-j)\right)\right]$

Models are registered in a global MODELS dictionary exposing the impedance function, analytical Jacobian, parameter names, physical bounds, and initial relative scaling.

3.2. 2. Synthetic Data Generation with Scaled Gaussian Noise

For each set of true parameters ${\mathbf{p}}_{\mathrm{true}}$, a synthetic complex spectrum is generated and perturbed with additive Gaussian noise proportional to the impedance magnitude:

$$Z_{\exp }\left(\omega \right)=Z({\mathbf{p}}_{\mathrm{true}},\omega )+{\epsilon }_i^{\prime }+j{\epsilon }_i^{\prime \prime },{\epsilon }_i^{\prime },{\epsilon }_i^{\prime \prime }\sim \mathcal{N}\left(0,{\sigma }_i\right),{\sigma }_i={\displaystyle \frac{\alpha }{100}}\left|Z\left({\omega }_i\right)\right|$$

with controlled noise levels $\alpha \in \{2.5,5.0\}\%$.

3.3. 3. Multi–Stage Hybrid Fitting with Internal Scaling

The fitting routine fit_impedance_data(...) executes a robust two–stage optimization pipeline:

(a)

Global search:${\chi }^2\left(\mathbf{p}\right)$ is minimized via differential evolution (Differential Evolution) over dimensionless scaled parameters, enforcing internal physical bounds.

(b)

Local refinement:Levenberg–Marquardt (least_squares) is applied with analytical Jacobians, strict tolerances (xtol=ftol=gtol=1e-8), and convergence verification.

(c)

Automatic scaling: all parameters are normalized by an adaptive vector $\mathbf{s}$ derived from ${\mathbf{p}}_{\mathrm{true}}$ or from physical limits, improving numerical conditioning.

Outputs include optimal parameters, complex residuals, evaluated Jacobian, and cost metrics.

3.4. 4. Uncertainty Estimation via Parallelized Bootstrapping

A non–parametric bootstrap scheme with $N=100$ resamples is executed in parallel (multiprocessing). Each resample is fitted from a neutral point in the scaled parameter space via least_squares:

$${\left\{{\widehat{\mathbf{p}}}^{\left(b\right)}\right\}}_{b=1}^N\Rightarrow {\mathrm{CI}}_{95\%}=\left[{\mathrm{percentile}}_{2.5},{\mathrm{percentile}}_{97.5}\right]$$

3.5. 5. Numerical Evaluation and Fit Diagnostics

Each converged fit is assessed using:

Absolute error: RMSE over Z (both real and imaginary parts). Goodness–of–fit: ${\chi }^2$, $R^2$, and adjusted $R^2$. Model complexity penalty: Akaike (AIC) and Bayesian (BIC) information criteria. Numerical conditioning: condition number of the Jacobian. Parameter correlations: correlation matrix from the pseudo–inverse of $\left(J^{\mathsf{T}}J\right)$, scaled by residual variance

3.6. 6. Systematic Exploration of Models and Noise Levels

For each model–noise combination $\{\mathrm{model}\in [1,7],\alpha \in \{2.5,5.0\left\}\right\}$, 50 independent fits are executed. The fit with the lowest RMSE is retained as the representative optimum. Final results are exported in both CSV and JSON formats for reproducibility and traceability.

3.7. 7. Physical Validation and Automated Visualization

For each optimum fit, the system automatically:

(a)

Generates Nyquist and Bode plots (magnitude and phase) plus residual plots

(b)

Performs Kramers–Kronig validation via the Hilbert transform of $Re\left(Z\right)$

(c)

Exports graphical data in CSV and all associated parameters and metrics in JSON

4. Results

The series RLC model was evaluated at $2.5\%$ and $5.0\%$ noise levels, with results presented as paired panels (left: $2.5\%$; right: $5.0\%$) in Figure 1, Figure 2, Figure 3, Figure 4 and Figure 5.

In the Nyquist diagram (Figure 1), both noise levels reproduce the model’s characteristic arc, with the $5.0\%$ case exhibiting greater scatter in the high-frequency region. Nevertheless, the fitted curves remain in close agreement with the simulated spectra.

The Bode magnitude (Figure 2) shows excellent fidelity across several frequency decades, with only minor deviations at $5.0\%$ noise in the high-frequency limit. The Bode phase (Figure 3) reproduces the expected phase profile without significant distortion from noise.

The residuals (Figure 4) at $2.5\%$ remain centered and nearly homoscedastic. At $5.0\%$, variance increases proportionally to the noise level, yet no systematic bias emerges, indicating the fit’s robustness.

Kramers–Kronig validation via the Hilbert transform (Figure 5) confirms the causality, linearity, and stationarity of the fitted spectra at both noise levels, with no significant violations.

The quantitative summary in Table 1 shows $R^2$ and adjusted $R^2$ values exceeding 0.997 in all cases, with fitted parameters lying within narrow 95% confidence intervals around nominal values. Increasing the noise from $2.5\%$ to $5.0\%$ causes modest increases in RMSE, ${\chi }^2$, AIC, and BIC, as well as slight variations in R, L, and C, without compromising the model’s validity.

Overall, the series RLC model demonstrates high stability and predictive reliability under moderate noise, preserving both statistical robustness and physical consistency.

4.1. Model: Classic Randles $R_s+(R_{ct}\Vert C_{dl})$

The classic Randles model was evaluated under controlled noise levels of $2.5\%$ and $5.0\%$. In both cases, the characteristic semicircle appears slightly depressed due to discrete sampling, while the fit reproduces the experimental trace with high fidelity over several decades of frequency. Figure 6, Figure 7, Figure 8, Figure 9 and Figure 10 present paired panels (left: $2.5\%$ ; right: $5.0\%$) for the Nyquist diagram, Bode plots (magnitude and phase), complex residuals, and the Kramers–Kronig validation.

In the Nyquist plots (Figure 6), a well-defined semicircle is observed, with a diameter corresponding to $R_{ct}$. Increasing noise produces only a larger scatter in the high-frequency region, without altering the global morphology. In the Bode plots (Figure 7 and Figure 8), the phase minimum near ${\omega }_{max}\approx 1/\left(R_{ct}{C}_{dl}\right)$ is accurately captured for both noise levels. The residuals (Figure 9) remain centered and trend-free; at $5.0\%$ noise, the variance increases without introducing systematic bias. The Kramers–Kronig check (Figure 10) confirms causality and linearity across the full frequency domain.

Quantitatively (Table 2), RMSE increases moderately from 3.2200 at $2.5\%$ noise to 6.5700 at $5.0\%$, while $R^2\ge 0.971$ and $R_{\mathrm{adj}}^2\ge 0.971$ are maintained in all cases. The estimated parameters fall within the 95% confidence intervals around their nominal values, confirming robust identifiability even under moderate noise.

Overall, these results support the robustness of the classic Randles model for systems dominated by charge-transfer processes and an ideal double-layer capacitance, preserving predictive accuracy and physical consistency under controlled experimental perturbations.

4.2. Model: Randles + CPE $R_s+\left(R_{ct}\Vert {\displaystyle \frac{1}{Q{\left(j\omega \right)}^n}}\right)$

Increasing the noise level from $2.5\%$ to $5.0\%$produces a controlled growth in error (RMSE: 3.3100 → 6.5200) and a moderate decrease in $R^2$ (0.9903 → 0.9634), while preserving spectral agreement (Table 3). The Nyquist diagram retains the characteristic depressed semicircle associated with the CPE (Figure 11), whereas in the Bode plots (Figure 12 and Figure 13) the fit accurately reproduces the magnitude across several decades and precisely locates the phase minimum.

Complex residuals (Figure 14) remain centered and drift-free for both noise levels; the increase in variance at $5.0\%$ is consistent with the added perturbation. The Kramers–Kronig validation (Figure 15) confirms causality and linearity, with no evidence of nonlinearity or nonstationarity.

Regarding parameters, $R_{ct}$ remains stable and within its 95% confidence intervals (Table 3), while Q and n exhibit the typical collinearity of non-ideal elements (Cond# $\sim {10}^{17}$–${10}^{18}$). For physical interpretation, the effective capacitance is computed using the Brug correction:

$$C_{\mathrm{eff}}={\left[Q{(R_s^{-1}+R_{ct}^{-1})}^{(1-n)}\right]}^{1/n},$$

yielding $C_{\mathrm{eff}}\sim 2$–$3\mu \mathrm{F}$ for $Q\approx 10\mu \mathrm{F}{\mathrm{s}}^{n-1}$, $n\approx 0.90$, $R_s\approx 10\Omega$, and $R_{ct}\approx 100\Omega$. This corresponds to $\tau \sim 2$–$3\times {10}^{-4}\mathrm{s}$ and $f_c\sim {10}^2$–${10}^3\mathrm{Hz}$, in excellent agreement with the experimentally observed phase minimum.

When compared at the same noise level, the classic Randles model is more parsimonious at $2.5\%$ (lower AIC in Table 2; $\Delta$AIC $\approx +8$ relative to CPE), whereas at $5.0\%$ Randles+CPE offers a significant advantage (AIC lower in Table 3; $\Delta$AIC $\approx -19$), capturing non-idealities without overfitting and preserving white residuals.

In summary, Randles+CPE reproduces with high fidelity the depressed arc and fractional phase shift typical of non-ideal interfaces, outperforming the classic Randles model under conditions of pronounced non-ideality or higher noise, while preserving $R_{ct}$ stability under rigorous statistical and physico–mathematical validation (Figure 11, Figure 12, Figure 13, Figure 14 and Figure 15; Table 3).

4.3. Model: Randles + Warburg $R_s+(R_{ct}\Vert C_{dl})+\sigma {\omega }^{-1/2}(1-j)$

The inclusion of the Warburg element introduces a 4.5000 × 10^1∘ tail in the Nyquist plot (Figure 16) and a $-1/2$ slope in the logarithmic magnitude of the Bode plot, accompanied by a phase tending toward $-{45}^{\circ }$ at low frequencies (Figure 17 and Figure 18). Both features are consistently observed at $2.5\%$ and $5.0\%$ noise levels, supporting the presence of a mass-transport–limiting process.

Increasing the noise level raises the RMSE from 3.1900 to 6.3800 and reduces $R^2$ from 0.9930 to 0.9725 (Table 4), while maintaining close agreement in both Nyquist and Bode domains (Figure 16, Figure 17 and Figure 18). Residuals remain centered and trend-free (Figure 19), and Kramers–Kronig validation confirms causality and linearity (Figure 20), ruling out artifacts from nonlinearity or nonstationarity.

$R_s$, $R_{ct}$, and $C_{dl}$ are estimated stably with narrow 95% confidence intervals at both noise levels (Table 4). In contrast, the effective diffusivity ${\sigma }_W$ shows greater uncertainty, especially at $5.0\%$, where the 95% CI reaches the lower bound of the fit (1.0000 × 10⁻⁶–9.4900 × 10⁻⁴). This indicates that, within the available frequency window, the diffusive component is morphologically detectable (Figure 16 and Figure 18) but its magnitude is partially ill-conditioned. Extending the frequency sweep to lower values and/or using a finite-length Warburg variant (e.g., with ${\tau }_D$) would help constrain ${\sigma }_W$.

The phase trend toward $-{45}^{\circ }$ (Figure 18) and the tail in the Nyquist plot (Figure 16) are consistent with semi-infinite diffusion of electroactive species to/from the interface. The stability of $R_{ct}$ under noise confirms that the faradaic kinetics are well identified, with the main variation residing in the mass-transport component (Table 4).

Compared to the classic Randles model, adding the Warburg element improves model likelihood at both noise levels (AIC: 349 vs. 362 at $2.5\%$ ; 348 vs. 353 at $5.0\%$ ; cf. Table 4 vs. Table 2). This indicates that the diffusive signature adds genuine information without overfitting, consistent with white residuals (Figure 19) and the KK validation (Figure 20).

4.4. Model: $(R_{ct}+Z_W)\Vert C_{dl}$

The series combination $R_{ct}+Z_W$, placed in parallel with $C_{dl}$, reproduces an arc with a diffusive tail in the Nyquist plot (Figure 21). In the Bode magnitude representation, the fit accurately follows the simulated data over several frequency decades (Figure 22), while the phase exhibits the expected drift toward the diffusion-controlled regime (approaching $-{45}^{\circ }$ at low frequencies; Figure 23). Increasing noise from $2.5\%$ to $5.0\%$ barely affects the morphological agreement, with discrepancies mainly concentrated at high frequencies, where the capacitive pathway of $C_{dl}$ dominates.

The RMSE increases from 3.3300 to 6.4700, and $R^2$ decreases from 0.9923 to 0.9717 (Table 5), while maintaining close agreement with the reference spectrum. The residuals (real and imaginary parts) remain centered without any trend (Figure 24), and the Kramers–Kronig (Hilbert transform) validation confirms causality and linearity of the simulated spectrum at both noise levels (Figure 25).

The parameters $R_s$, $R_{ct}$, and $C_{dl}$ are stable and have narrow 95% confidence intervals (Table 5). The diffusive amplitude ${\sigma }_W$ is the least constrained: at $2.5\%$ noise it exhibits a wide 95% CI $[1.0000\times {10}^{-6},3.7600\times {10}^{-3}]$, and at $5.0\%$ it reaches the lower bound of the fit $[1.0000\times {10}^{-6},4.6900\times {10}^{-3}]$. This indicates that, within the available frequency window, the diffusive signature is visible (Figure 21 and Figure 23), but its amplitude is partially ill-conditioned due to the capacitive shunt. The condition number is moderate ($\sim 8.0000\times {10}^5$), significantly better than in CPE-based models, supporting the overall numerical stability of the fit.

The parallel arrangement with $C_{dl}$ shunts the faradaic branch at high frequencies (almost purely capacitive trajectories), while at low frequencies $R_{ct}+Z_W$ dominates, producing the 4.5000 × 10^1∘ diffusive tail (Figure 21 and Figure 23). The constancy of $R_{ct}$ under noise suggests that the interfacial kinetics are well-identified, with the main uncertainty arising from the effective mass transport term encapsulated in ${\sigma }_W$.

At $2.5\%$ noise, this model is less statistically supported than Randles+Warburg and Randles+CPE (AIC = 382 vs. 349 and 370; $\Delta$AIC = +33 and +12, respectively), and even less than the classical Randles (362; $\Delta$AIC = +20). At $5.0\%$ , its performance approaches that of the classical Randles (354 vs. 353; $\Delta$AIC = +1), but remains below Randles+Warburg (348; $\Delta$AIC = +6) and particularly Randles+CPE (334; $\Delta$AIC = +20). Overall, the statistical evidence favors alternative descriptions when the available frequency window does not sufficiently constrain ${\sigma }_W$.

The $(R_{ct}+Z_W)\Vert C_{dl}$ topology captures the expected kinetic–diffusive interaction and produces physically consistent fits (Figure 21, Figure 22, Figure 23, Figure 24 and Figure 25); however, under the current measurement window its informational advantage is limited compared to Randles+Warburg or Randles+CPE, due to the poor constraint on ${\sigma }_W$.

4.5. Model: $(R_{ct}+Z_W)\Vert \mathrm{CPE}$

This topology combines semi-infinite diffusion ($Z_W$) with a non-ideal capacitive response (CPE). In the Nyquist plot, a depressed arc appears along with a diffusive extension at low frequencies (Figure 26); in the Bode magnitude plot, the fit reproduces the spectrum across several decades (Figure 27) and captures the fractional phase minimum associated with $n<1$ (Figure 28).

Increasing the noise from $2.5\%$ to $5.0\%$ raises the RMSE from 3.29 to 6.44 and reduces $R^2$ from 0.9904 to 0.9645 (Table 6), while preserving morphological agreement in both representations. The residuals remain centered and trend-free (Figure 29); the Kramers–Kronig test confirms causality and linearity of the spectrum (Figure 30).

The parameters $R_s$ and $R_{ct}$ are stable with narrow 95% confidence intervals at both noise levels (Table 6). The exponent $n\approx 0.90$ indicates a distribution of relaxation times consistent with interfacial heterogeneity/roughness; $Q\sim 9$–$10\mu \mathrm{F}{\mathrm{s}}^{n-1}$ corresponds to an effective capacitance of $\sim 2$–$3\mu \mathrm{F}$ (using Brug’s correction, as in the CPE subsection). The diffusive amplitude ${\sigma }_W$ is the least constrained at $2.5\%$ noise (95% CI reaches the lower bound), but is better defined at $5.0\%$. The condition number is moderate ($\mathrm{Cond}\#\sim 2.2\times {10}^6$), much lower than for models with CPE only, suggesting reduced collinearity between ($Q,n$) and the diffusive branch when $Z_W$ is present.

The CPE accounts for interfacial non-idealities (distribution of time constants), while $Z_W$ captures mass-transport control at low frequencies. The phase transition toward $-{45}^{\circ }$ (Figure 28) and the Nyquist tail (Figure 26) highlight the diffusive contribution; the stability of $R_{ct}$ supports a well-identified faradaic kinetics.

At $2.5\%$ noise, this model markedly outperforms the classical Randles (AIC 352 vs. 362; $\Delta$AIC = $-10$) and Randles+CPE (352 vs. 370; $-18$), and is close to Randles+Warburg (352 vs. 349; $+3$). At $5.0\%$, it becomes the best-performing alternative among those considered (AIC 333), surpassing Randles+CPE (334), Randles+Warburg (348), and the classical Randles (353) (Table 6 and related tables). This suggests that when interfacial non-ideality and a diffusive component are present simultaneously, the $(R_{ct}+Z_W)\Vert \mathrm{CPE}$ architecture offers the best balance between fit quality and parsimony.

The combination of CPE and Warburg in parallel with $R_{ct}$ provides a coherent and statistically preferred representation of systems exhibiting concurrent interfacial non-ideality and diffusion, with white residuals and consistent KK validation (Figure 26, Figure 27, Figure 28, Figure 29 and Figure 30; Table 6).

4.6. Model: Full Randles

The model reliably reproduces the canonical features: a depressed semicircle at intermediate frequencies and a diffusive tail at low frequencies (Nyquist, Figure 31); magnitude decay around the characteristic frequency and a fractional phase minimum due to the CPE (Figure 32 and Figure 33). The Kramers–Kronig validation shows overall agreement between the spectrum and the Hilbert reconstruction (Figure 35), with bounded deviations at the band edges attributable to the spectral window width.

When increasing noise from $2.5\%$ to $5.0\%$, the RMSE rises from 3.2700 to 6.6200 and $R^2$ decreases from 0.9905 to 0.9625 (Table 7), but the fit morphology remains consistent in both representations (Figure 31, Figure 32 and Figure 33). Residuals remain centered and free from systematic trends; only their variance increases, in line with the noise level (Figure 34).

Figure 34. Residuals (real and imaginary parts) for the full Randles model.

Figure 34. Residuals (real and imaginary parts) for the full Randles model.

Figure 35. Kramers–Kronig (Hilbert transform) validation for the full Randles model.

Figure 35. Kramers–Kronig (Hilbert transform) validation for the full Randles model.

Ohmic and kinetic parameters remain stable (narrow 95% CIs in $R_s$ and $R_{ct}$; Table 7). In contrast, Q, n and, particularly, the effective diffusivity $\sigma$ exhibit marked collinearity (Cond# $\sim {10}^{11}$), with $\sigma$ reaching the upper bound of the interval at $2.5\%$ . This indicates that, with the available frequency range, CPE–Warburg separation is not fully excited: the CPE absorbs part of the diffusive curvature and vice versa. A practical way to convey the capacitive response is to report, along with Q and n, the effective capacitance$C_{\mathrm{eff}}$ (Brug correction), which here is on the order of a few $\mu \mathrm{F}$, consistent with the position of the phase minimum (Figure 33).

The exponent $n\approx 0.90$ reflects a distribution of times associated with interfacial roughness/heterogeneity, while the $-{45}^{\circ }$ tail and magnitude slope corroborate a mass transport component (Warburg). The stability of $R_{ct}\approx 100\Omega$ at both noise levels suggests a well-identified faradaic kinetics, decoupled from the uncertainties in Q, n and $\sigma$ (Table 7).

At $2.5\%$ , the AIC of the full Randles (350) is comparable to Randles+Warburg (349) and better than Randles+CPE (370), but does not provide a material advantage over the option without CPE or without Warburg (Figure 31, Figure 32 and Figure 33; Table 7). At $5.0\%$ the model is penalized for complexity (AIC 356; BIC 376) compared to more parsimonious alternatives, especially $(R_{ct}+Z_W)\Vert \mathrm{CPE}$ and Randles+CPE, which achieve lower AIC values with equal or better fidelity. Consequently, if σ is poorly bounded and the AIC improvement is marginal, the full version is not statistically justified.

The full Randles simultaneously captures interfacial non-ideality and diffusion, but with the present spectral range, the gain over partial models is limited and the identifiability of $\sigma$ is fragile. For data with mild non-ideality and moderate diffusion, $(R_{ct}+Z_W)\Vert \mathrm{CPE}$ or Randles+Warburg offer a superior trade-off between physical realism and parsimony; the full model should be reserved for spectra with clear and quantifiable evidence of both contributions (Figure 31, Figure 32, Figure 33, Figure 34 and Figure 35; Table 7).

4.7. Global Model Comparison

Table 8 summarizes the key fitting metrics for all seven evaluated models at both noise levels. At $2.5\%$ , differences in RMSE among the more complex models (Randles full, Randles+Warburg, $(R_{ct}+Z_W)\Vert \mathrm{CPE}$) are marginal, with statistical improvements in AIC and BIC insufficient to justify overparameterization when $\sigma$ or n exhibit broad confidence intervals or strong collinearity. At $5.0\%$ noise, the more parsimonious models (Randles+CPE and $(R_{ct}+Z_W)\Vert \mathrm{CPE}$) yield lower AIC/BIC values with negligible $R^2$ loss, reflecting greater robustness against signal degradation.

From a physical standpoint, the simultaneous inclusion of a CPE and Warburg impedance (Randles full) is warranted only when the spectrum exhibits quantifiable evidence of both interfacial non-ideality and semi-infinite diffusion, and when the spectral range adequately excites both processes. In the absence of these conditions, partial models offer a superior trade-off between fitting fidelity, parameter stability, and parsimony.

5. Conclusions

This work establishes a fully integrated and analytically rigorous framework for Electrochemical Impedance Spectroscopy (EIS) analysis, combining exact impedance formulations, closed-form Jacobians, hybrid global–local optimization, and exhaustive statistical validation. Applied consistently across seven representative equivalent circuit models under controlled synthetic conditions, the methodology delineates precise statistical and physical boundaries for model applicability, bridging the gap between theoretical formulation and practical implementation.

Results demonstrate that the classical Randles circuit remains the most parsimonious and statistically optimal choice for systems with ideal charge-transfer kinetics and capacitive behavior, while configurations incorporating Constant Phase Elements (CPE) and Warburg diffusion elements provide superior fidelity for heterogeneous or diffusion-limited interfaces. The robustness of $R_{ct}$ estimation across noise levels, alongside the ability to quantify non-ideal capacitive parameters and diffusive amplitudes, underscores the stability and reproducibility of the approach.

Critically, automated Kramers–Kronig compliance testing ensures that all fitted spectra satisfy the fundamental principles of linearity and causality, eliminating a recurrent source of interpretative bias in EIS literature. The modularity of the platform enables rapid extension to novel circuit topologies, integration with experimental datasets, and coupling with advanced electrochemical–transport models.

By uniting analytical precision, computational efficiency, and rigorous statistical scrutiny, this framework delivers an unprecedented level of objectivity in EIS model selection and parameter interpretation. As such, it is positioned to become a reference standard for quantitative EIS analysis, setting a reproducible benchmark for the accurate characterization and optimization of electrochemical systems in energy storage, corrosion science, and electrocatalysis. .