# Hybrid stochastic-deterministic algorithms for the interpretation of Electrochemical Impedance Spectroscopy spectra of Proton Exchange Membrane Fuel Cells

P. Wu, S. Touhami, W. Aït-Idir, H. Yelloz, C. Marty, F. Micoud, J. Dillet, O.

Lottin, Julia Mainka

► To cite this version:

P. Wu, S. Touhami, W. Aït-Idir, H. Yelloz, C. Marty, et al., Hybrid stochastic-deterministic algorithms for the interpretation of Electrochemical Impedance Spectroscopy spectra of Proton Exchange Membrane Fuel Cells. Electrochimica Acta, 2025, 518, pp.145673. 10.1016/j.electacta.2025.145673. hal-04930033

HAL Id: hal-04930033

https://hal.univ-lorraine.fr/hal-04930033v1

Submitted on 5 Feb 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Hybrid stochastic-deterministic algorithms for the interpretation of Electrochemical Impedance Spectroscopy spectra of Proton Exchange Membrane Fuel Cells

P. Wu<sup>a ID</sup>, S. Touhami<sup>a</sup>, W. Aït-Idir<sup>a</sup>, H. Yelloz<sup>a</sup>, C. Marty<sup>b</sup>, F. Micoud<sup>b</sup>, J. Dillet<sup>a ID</sup>, O. Lottin<sup>a ID</sup>, J. Mainka<sup>a ID</sup>

## ARTICLE INFO

### Keywords:

Proton exchange membrane fuel cell (PEMFC)
Electrochemical impedance spectroscopy (eis)
Equivalent electrical circuit (eec) modeling
Global Optimization
Hybrid Algorithms

## ABSTRACT

This work presents a case study on the use of three hybrid algorithms combining a stochastic part -Genetic Algorithms (GA), Particle Swarm Optimization (PSO), or Simulated Annealing (SA)- with a deterministic Nelder-Mead (NM) algorithm for the estimation of the equivalent electrical circuit (EEC) parameters for the interpretation of Proton Exchange Membrane Fuel Cell (PEMFC) impedance data. These hybrid methods were evaluated on mathematical test functions as well as for the interpretation of simulated and experimental PEMFC impedance spectra using EEC of different complexity.

The three stochastic/deterministic methods were compared in terms of stability, efficiency, ability to explore multiple solutions, and computing resources. The results showed that all hybrid methods were able to improve the interpretation of experimental EIS data by identifying satisfying and physically meaningful solutions, with low least-square residuals and by reducing the sensitivity to initial conditions while accelerating convergence. All methods allowed an improvement compared to the use of one single type of algorithm alone -deterministic and stochastic.

## 1. Introduction

Electrochemical Impedance Spectroscopy (EIS) is an efficient tool for characterizing electrochemical systems during operation, with only minor perturbations [1–3]. It is based on the application of small amplitude sinusoidal variations of the electrical voltage (potentiostatic mode) or current (galvanostatic mode) of varying frequency f around a given operation point so that the response of the system can be considered as linear (Fig. 1, left). The response is recorded in terms of current or voltage to determine the transfer function, i.e. the impedance:

$$Z(\omega) = \frac{\tilde{u}(t)}{i(t)} = Z_{re}(\omega) + iZ_{im}(\omega) \quad (1)$$

With $i^2 = -1$ being the imaginary unit number and $\omega = 2\pi f$ the angular frequency (in rad/s).

The impedance is a complex function of the frequency at which the perturbation is applied, as well as of the system operating point. It enables the separate analysis of the different phenomena occurring in the system at the material scale in the frequency domain according to the shape, or other typical characteristics (inflexion points, radii, etc.) of the impedance spectra as shown in Fig. 1, right, in the case of the Nyquist

plot of an ideal Polymer Electrolyte Membrane Fuel Cell (PEMFC) spectrum.

One way to model the impedance is under the form of equivalent electrical circuits (EEC) composed of -more or less simple- electrical components (such as resistors, capacitances, inductances, etc.) aiming at representing the charge and mass transport phenomena in the individual components of the system. To interpret the impedance of a PEMFC, the most extensively used EEC are **Randles-type circuits** [4–6] and **Transmission Line Models** (TLM) [5,6]. TLM are a generalization of the Randles circuit [1], taking into account the electrical potential variation through the electrode thickness: they must be used with thick electrodes. To represent mass transport limitations at low frequencies, the most common way is to add a finite **Warburg element** in the EEC [7]. This element can be placed in series with the charge transfer resistance (in the case of the Randles circuit) or in series with the **Cathode Catalyst Layer** (CCL) circuit, depending on the most limiting layer for oxygen diffusion: the CCL or the **Gas Diffusion Layer** (GDL), respectively [1].

The identification of appropriate EEC parameters consists in solving a **Complex Non-linear Least-Square** (CNLS) problem. This requires to use optimization algorithms [8,9] to minimize the square of the total

<footnote><sup>a</sup> LEMTA, CNRS, Université de Lorraine, Nancy, France</footnote>

<footnote><sup>b</sup> Univ. Grenoble Alpes, CEA, LITEN, F-38054, Grenoble, France</footnote>
<figure><img src="image_3.png" /><figcaption>Fig. 1. Principle of EIS measurement and ideal PEMFC spectrum in the Nyquist plot with several representative parameters.</figcaption></figure>

<figure><img src="image_4.png" /></figure>

difference between the experimental data and the EEC impedance, represented by an error function defined herein as:

$$err = \frac{1}{n} \sum \frac{|Z_{\text{theo}} - Z_{\text{exp}}|^2}{|Z_{\text{exp}}|^2} \qquad (2)$$

Where $Z_{\text{theo}}$ is the expression of the impedance model, a function of multiple parameters: $Z_{\text{theo}} = f(r, \omega)$, with $r = (r_1, r_2, r_3...)$ the vector of EEC parameters. $n$ stands for the number of experimental points and $Z_{\text{exp}}(\omega)$ is the measured impedance of the system.

There exists a large number of CNLS solving algorithms, which can generally be classified in the following categories [10–12]:

* Local algorithms, converging only towards the nearest (local) minimum from a given starting point and global algorithms, able to search for the global minimum in each domain.

* Deterministic algorithms, always following the same path and converging toward identical results in given conditions and stochastic algorithms, with routines that generate and use random numbers so that they necessarily take different paths and may lead to different solutions when run several times in identical conditions. Often, stochastic algorithms are also global algorithms because in practice, the global search ability of an algorithm relies on stochastic routines.

Deterministic algorithms are generally efficient in terms of convergence and repeatability. This is why they are widely used for the identification of EEC parameters [13–16]. However, their result may be very sensitive to the values used to initiate the algorithm: when the initial values are too far from the optimum solution, the algorithm may converge to a local minimum. This becomes critical when the problem dimension (the number of parameters) is high, so that the (empirical) search for appropriate initial values may become difficult. Therefore, deterministic algorithms are difficult to run reliably in an automated manner, for instance to monitor the ageing of a system or to characterize new materials.

Stochastic algorithms are much better suited to global optimization problems than deterministic ones: since they generate random parameter values at the beginning of each execution, they generally do not need an initial vector but a parameter range instead. Among the most common stochastic algorithms, simulated annealing algorithms is an exception that still require an initial vector, but the influence of the starting point is lower than with deterministic methods since the algorithm will explore the possible solutions within a predefined search space [13,17–21]. However, stochastic algorithms may return different results when run multiple times, even more so when the local minima lead to similar values of the error function: this can be considered as a weakness in some case or conversely, as a way to put forward various

possible solutions. However, convergence could be slow when approaching a minimum, so that converging to the minimum of the error function requires high computational costs.

In this work, three of the most common stochastic algorithms, i.e. Genetic Algorithms (GA), Particle Swarm Optimization (PSO) algorithms and Simulated Annealing (SA) algorithms, are used in combination with a deterministic Nelder-Mead algorithm for the identification of the impedance parameters of PEMFC EEC. All those algorithms were taken from MATLAB standard libraries. The output of the stochastic algorithm is used as the input of the deterministic algorithm (Fig. 2) -run as a second step- to benefit from the advantages of both. The algorithms as well as the EEC and the associated model parameters are introduced in section 2. These hybrid methods were evaluated on mathematical test functions (section 3.1) and for the interpretation of simulated and experimental impedance data obtained on PEMFC (section 3.2), which are modeled by Randles-type circuits or by TLM, depending on the shape of the impedance spectra at high frequencies. In our case, such EEC include at least 5, and sometimes up to 11 parameters.

In all cases, it turned out that these hybrid methods are globally efficient enough in terms of precision and robustness, and to explore multiple possible solutions. The different algorithm combinations are compared with the final objective to select the best concept for interpreting PEMFC impedance data.

## 2. Algorithms

### 2.1. Concept of hybrid algorithms

Stochastic and deterministic algorithms each have their advantages and weaknesses. Deterministic algorithms are generally fast and stable, but they are very sensitive to the starting point, all the more so when there are multiple possible solutions, i.e. local minima. Conversely, stochastic algorithms are powerful tools to identify a global optimum, but they are much slower and less precise<fnref n="1" />: they may converge to different solutions for multiple executions with the same input conditions. These issues become more significant when the dimension of the problem (the number of model parameters) is high.

The concept of hybrid optimization algorithms -combining multiple optimization techniques or algorithms-, is known for several decades. It is difficult to pinpoint the exact first proposal of a hybrid optimization algorithm, as different researchers have independently developed and combined various optimization methods over time. A remarkable work was proposed by Moscato in 2000 [22] who introduced a “memetic

<footnote n="1">¹ In the context of optimization, the term "precision" refers to the degree of accuracy or closeness to the true optimum that an optimization algorithm can achieve [48].</footnote>
<figure><img src="image_5.png" /><figcaption>Fig. 2. Outline of the principle of hybrid stochastic-determinist algorithms.</figcaption></figure>

algorithm” (MA) which combines a genetic algorithm (GA) with deterministic local refinement strategies (GA-MA).

Hybrid optimization algorithms have been proposed and applied to a considerable range of systems. Zhang and his team [23] tested a PSO-LM combination to fit the multi-peaks Brillouin scattering spectrum and results showed that the recognition effect of this method is obviously better than that of traditional algorithms and it can improve the accuracy of information analysis. Wang and al. [10] presented a hybrid GA-Levenberg-Marquardt (LM) method for cell suspension measurement using EIS: the effectiveness of the GA-LM method was demonstrated through results on simulations and experimental data of yeast suspensions. A. Kabouche [24] used a Genetic Algorithm (GA) and Simulated Annealing (SA) coupled with Nelder-Mead (NM) method for the parameter estimation of a liquid-liquid phase equilibrium model. T. Zhou et al. [12] did a case study by combining GA with a deterministic mixed-integer nonlinear programming (MINLP) algorithm for integrated solvent and process design. Many other studies [12,22,24,25] also showed that hybrid algorithms have a better performance in global optimization compared to stochastic [26,27] and deterministic algorithms [1,16,28–30] alone. This work proposes to go one step further to reveal the strength of hybrid method in practical scenarios and discusses different aspects of algorithm performances within the context of EEC parameter identification of fuel cells.

In our previous works, the characterization of PEMFC impedance spectra was made using deterministic NM algorithms [1,16,28–30]. They yield satisfying results in the case of EEC with a low number of parameters only when the starting point corresponds to approximate values of the parameters, which complicates the study of new materials (i.e. when the initial values are difficult to assess). When a deeper analysis of PEMFC impedance spectra and thus more elaborated models with a higher number of parameters (>10) is required, NM algorithms alone become instable. This is why, in this work several alternative hybrid methods are proposed. As the main objective is to compare the performance of those hybrid methods to deterministic algorithms, EEC and experimental data of some of our previous studies are used [28,29, 31–33].

## 2.2. Optimization algorithms

Zero order algorithms do not require the calculation of the function derivative: this is the case of all algorithms used in this work.

### 2.2.1. Deterministic Nelder-Mead algorithm

The Nelder-Mead (NM) algorithm [34] is based on the evolution of a polytope of $n + 1$ points surrounding a $n$-dimensional vector of parameters $x$, called simplex. The simplex is modified at each iteration until its vertices are close enough to each other, i.e. when $x$ approaches the nearest local minimum.

Initially, the algorithm creates a simplex around the initial guess $x_0$ by increasing each component $x_0(i)$ by a certain percentage (5 % in MATLAB). The simplex is modified at each iteration according to a set of four operations: shrinkage, contraction, reflection, and expansion. The sequence and order of these operations depends on the evolution of the

function value (Fig. 3).

### 2.2.2. Stochastic algorithms

Genetic Algorithms (GA) are global stochastic algorithms inspired by Charles Darwin's theory of natural evolution [35]. A group of individuals (i.e. a set of parameter values) called **population** is generated randomly at the beginning of each execution: in the earlier version of the algorithms, each one was associated with a binary code representing the “chromosome” of the individual, and a fitness value calculated using the function defined by the user (the error function in our case, Eq. (2)). The fitness value shows how well each individual is adapted to its environment. At each iteration, new individuals are generated from the information contained in the chromosomes of the individuals of the previous generation. This is achieved using three to four elementary operators:

* **Selection:** individuals with better fitness value have higher chance to survive, i.e. to be selected to generate the next population. The probability of being selected varies between (0,1], the best individual being always selected. A ratio called generation gap determines how many individuals will be eliminated at each generation, which is sometimes considered as a fourth operator: **elitism**.

* **Crossover:** crossover is the main method for generating new individuals. Selected individuals can mate with each other, exchanging a part of their chromosomes to create a new individual added to the population. The crossover probability is denoted $P_c$.

* **Mutation:** mutation ensures the ability of searching the global optimum. To maintain a relative stability of the population, the mutation probability $P_m$ should be much smaller than $P_c$. It consists for instance of swapping some binary code of the chromosome between 0 and 1 after the crossover phase [36,37].

These steps are repeated at each generation, i.e. iteration (Fig. 4). Hence, the population evolves to improve the individual's fitness score from one generation to the other. The algorithm stops when the generation number $N_{gene}$ exceeds the maximum number $Max_{gene}$ or when the average change of best individual's fitness value $\Delta$ is less than the user-defined tolerance $\epsilon$ over a certain number of generations which is called the maximum stall generation. Genetic algorithms were originally

<figure><img src="image_6.png" /><figcaption>Fig. 3. Structure of the simplex used in Nelder-Mead algorithms. The simplex is updated at each iteration to approach the nearest local minimum.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Fig. 4. Flowchart of GA, individuals are selected and updated in each generation.</figcaption></figure>

designed to solve combinatorial problems and have gradually been applied to continuous problems using real instead of binary numbers [38].

**Particle Swarm Optimization (PS or PSO)** algorithms are also based on the evolution of a population. They were first proposed by J. Kennedy and R. Eberhart [39] and were later modified by E. Mezura-Montes [40] and M-E. Pedersen [41]. This family of algorithm is inspired by the behavior of a flock of birds searching for food or escaping a predator. Socio-biologists believe that birds can share information between all members, which helps them to find the best ways.

Like GA, PSO algorithms are based on a population of particles (equivalent to individuals of GA, i.e. a set of parameter values) randomly created at the beginning of the procedure and then modified to find the best location (equivalent to the fitness value of GA, i.e. the value of the error function). Besides the location, each particle is also associated with an initial velocity, determining the change that is allowed for the value of a parameter between 2 iterations. The algorithm iteratively updates each particle's location and velocity according to 3 criteria:

* the current location and velocity of the particle: current flying.

* the best location the particle has ever found: flying memory.

* and the best location the entire population has ever found: the swarm influence.

The algorithm is executed as follows:

* **Initialization**: PS optimization algorithms create randomly n particles in the problem area with their positions and velocities. Then, they evaluate the objective function (here the error function) value for each particle and record their current position $p(i)$ and velocity $v(i)$, which will be replaced by the position with the best function value $i$ has ever found in the following procedures. With $b$ the individual with the best function value among all particles: $b(i) = \min(f(p(i)))$ with $i = 1, ..., n$.

* **Iteration steps**: the velocity and the position of the particles are updated as:

$$p_{t+1}(i) = p_t(i) + v_{t+1}(i) \times \Delta t \quad (3)$$

With $\Delta t = 1$ by definition and

$$v_{t+1}(i) = Wv_t(i) + y_1u_1(p-x) + y_2u_2(g-x) \quad (4)$$

$v_{t+1}(i)$ is a weighted sum of the previous velocity $v_t(i)$, the difference between the best location $p$ the particle $i$ has ever found and its current location $x$ (i.e. $p-x$), and the difference between the best location the entire population has found $g$ and the current location of particle $i$ (i.e. $g-x$).

$W$, $y_1$ and $y_2$ are respectively 3 parameters called the inertia weight, the self-adjustment weight, and the social adjustment weight, they determine how much weight should be given between refining the search result of the particle itself and recognizing the search result of the swarm. These parameters control the tradeoff between **exploration** and **exploitation** (see in 2.4). $u_1$ and $u_2$ are random numbers between 0 and 1. In MATLAB, $W$ is higher at the beginning of the routine and decreases through the iterations to improve exploitation.

By repeating these steps (Fig. 5), the particle with the best position $b(i)$ is updated at each iteration, the algorithm stops when $b$ does not change anymore (change smaller than the tolerance several times in a row) or the max number of iterations is reached. The associated function implanted in MATLAB is “particleswarm”.

**Simulated Annealing (SA)** algorithms were invented by S. Kirkpatrick, C.D. Gelatt and M.P. Vecchi [42] and are inspired from metallurgy: annealing consists in heating a material above its temperature of recrystallization, maintaining this temperature for a given time, and then cooling it down. As the temperature decrease is relatively slow, the particles are likely to find a better position to minimize their Gibbs free energy.

Considering the target function (function “err” in Eq. (1)) as the Gibbs energy, the algorithm systematically lowers the “temperature” $T$ of the system during the execution:

* **Initialization**: unlike other stochastic methods, SA requires a starting vector instead of search intervals as input to start the execution. However, it is much less sensitive to the starting point than a deterministic algorithm.

* At the beginning, the algorithm creates a random trial point (solution) and chooses a “distance” of the trial point from the current point (the starting point for the first iteration) based on the current temperature. The higher the temperature, the larger the exploration domain around the current location.

* The algorithm then compares the new point with the current point and can accept the new point when it is not better than the current point: The acceptance probability depends on how worse the new point is and on the current temperature. **The higher the temperature, or the better the new point, the more likely the new point is to be accepted.** The acceptance probability $p$ should not be >0.5.

As the temperature lowers, the SA algorithm reanneals when it has accepted a certain number of points by raising again the temperature in each dimension (Fig. 6).

Like other stochastics algorithms, the SA algorithm stops when the average change is less than the tolerance $\epsilon$ during a certain number of iterations or if the maximum number of iterations is reached. The associated function implanted in MATLAB is “simulannealbnd”.

It is worth mentioning that like the NM algorithm, GA, PS, SA are also heuristic algorithms: they use problem-specific rules and strategies to look for the solutions of complex problems in a reasonably efficient manner. Of course, it can never be guaranteed that the true global solution [43] is found, for both stochastic and deterministic algorithms. Deterministic and stochastic optimization methods can be complementary when combining efficient deterministic guided search strategies with random exploration of stochastic methods to enhance the search process and improve the quality of the solution.

Many other algorithms exist for CNLS problems, with various versions. The number of optimization methods and their implementation is
<figure><img src="image_8.png" /><figcaption>Fig. 5. Flowchart of PSO algorithm.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 6. Flowchart of the SA algorithm.</figcaption></figure>
constantly increasing with comparative studies evaluating their performance [14,44–47]. Such benchmarking studies can be of great value in helping end-users for choosing the most suitable optimization method for their problems. The present work is rather to be considered as a case study or a proof of concept [48] since it is more oriented toward a specific problem.

The performance of hybrid methods is mostly determined by the stochastic part(s). Therefore, although there exists few research on benchmarking of hybrid methods, the results from benchmarking of the stochastic part(s) can offer a first indication of the overall performance.

PSO algorithms are relatively recent. Like GA, they are population-based and depend on information sharing. Results with mathematical functions that are used for general benchmarking [49] showed that PSO algorithms have a larger differential in computational efficiency<fnref n="2" /> than GA when used to solve unconstrained nonlinear problems with continuous variables. The differential efficiency of PSO algorithms is however lower when applied to constrained nonlinear problems with discontinuous or discrete variables.

However, on average, PSO algorithms and GA yield the same efficiency for most problems. SA being based on one “individual” instead of a whole population, this algorithm is thus expected to be less efficient in terms of precision, but should converge much faster, as confirmed by a study of F. Jia and D. Lichti [45]. On the other hand, PSO algorithms exhibit certain drawbacks, including the risk of premature convergence [50,51]. Ongoing research is actively exploring further developments and enhancements to address these limitations.

## 2.3. Input parameters

GA, PS and SA algorithms have some common or equivalent control parameters.

* **MaxIteration**, standing for the maximum number of iterations. Note that it is called “max generation” in the case of GA, considering a generation as an iteration.

* **FunctionTolerance** and **MaxStallIteration**: the algorithm stops if the average relative change in the target function value over MaxStallIteration iterations/generations is less than or equal to FunctionTolerance. These two parameters and the MaxIteration constitute the algorithm stopping condition. Note that FunctionTolerance should be small enough compared to the target function value (i.e. the residual). On the other hand, a very low value for FunctionTolerance and a very high value for MaxStallIteration will require more computing resources without improving the algorithm's performance. In our case, setting FunctionTolerance to $10^{-12}$ (being $10^{-7}$ times the residual) and MaxStallIteration to 50 gave satisfying results for all algorithms and studied cases.

* **SwarmSize/PopulationSize** represents the number of individuals/particles of the population/swarm for PS algorithms and GA, respectively. Increasing this number can improve the algorithm performances, up to a certain value when a critical density of individuals is reached (saturation). This parameter is not applicable to SA algorithms since it is not population-based.

* The **Bound limits** (upper and lower bound) are also important parameters that can influence the algorithm's performance. Conservative bound limits can reduce the search space and thus the convergence time. On the other hand, in many real-world problems, certain parameters have physical or practical limits and violating these limits may render the solution unrealistic [52–54]. In this work, the limits of the EEC parameters were set according to a

<table><caption>**Table 1**
Fixed input parameters for all optimization algorithms and test sets.</caption><thead><tr><td>MaxIter</td><td>Function Tolerance</td><td>MaxStall Iteration</td><td>Swarm/Population Size</td></tr></thead><tbody><tr><td>5000</td><td>10<sup>-12</sup></td><td>50</td><td>5000</td></tr></tbody></table>

graphical estimation starting from the impedance spectra in the Nyquist plot (Fig. 1).

Table 1 gives the values of the parameters listed above. They were identical for all the test sets considered in the following, so that the algorithm's performances can be compared on similar bases.

In MATLAB, GA and PS algorithms generate automatically a random **initial population**, while the SA algorithm needs user-specified initial points. For a better comparison between the three algorithms, a code generating random initial values within the bound limits was used for SA.

In addition, some other parameters are specific to each algorithm: for instance, crossover and mutation for the GA, the inertia, self-adjustment and social adjustment weights for PS algorithm, and the acceptance probability for SA algorithm. Their modification requires extreme precaution and was considered out of the scope of this work, so that default values were used.

## 2.4. Evaluation criteria

The performance of optimization algorithms is generally evaluated according to three criteria: efficiency, reliability, and quality of the output [48].

The **efficiency** of an optimization algorithm refers to the computational effort required to obtain a solution. In mathematical programming, there are two main indicators of efficiency: the number of fundamental evaluations and the execution time. A third and less common indicator is the amount of memory used [48]. Among them, the **number of fundamental evaluations (nEval)** can be used as a standard unit of time and is by definition hardware independent [48]. This is the main reason why it is chosen here for the assessment of algorithm's efficiency.

The **reliability** or robustness of an optimization algorithm is defined as the ability to perform well over a wide range of problems [20,21]. The most common criteria to evaluate the reliability is the success rate [55]. For stochastic algorithms, the result of different executions can vary when there are several local minima of the target function, all the more so when these local minima correspond to function values that are close to each other. It is therefore important to run the algorithm several times to include as many solutions as possible. Finding the best solution(s) and exploring as many local minima as possible results from a tradeoff between exploitation (finding the best solution) and exploration (finding a maximum number of solutions).

It must be noted that these 2 objectives are not necessarily contradictory. Two parameters are used for their evaluation:

* **Exploitation**: when the analytical global minimum is known, the **success rate** is used to evaluate the reliability of the algorithm. Otherwise, the **occurrence** of each solution (number of times they appear in a number of runs) is used.

* **Exploration**: i.e. the **number of different solutions** found. The higher the number of distinct solutions found, the more efficient the algorithm is, in terms of exploration.

The **quality** of an algorithm can be evaluated by the value of the target function at the output, i.e. the **residual** (Eq. (2)), which introduces a normalized version of the average relative error, allowing the comparison between different test sets. Different solutions can be distinguished according to the values of the residual and those of each model parameter.

<footnote n="2">² The differential in computational efficiency [48] in optimization indicates how the efficiency of the optimization algorithm changes when the input changes (algorithm parameter, target function...).</footnote>
<figure><img src="image_10.png" /><figcaption>Fig. 7. Performance criteria of algorithms and the corresponding evaluation parameters.</figcaption></figure>

Fig. 7 gives an illustration of the performance criteria and their corresponding evaluation parameters.

The improvement of the performances of a hybrid approach compared to purely deterministic algorithms depends mainly on the stochastic part. Therefore, the three hybrid methods presented in the following use all one identical deterministic algorithm, i.e. NM, as second step, resulting in: GA-NM, PS-NM and SA-NM.

## 3. Results & discussion

The main objective of this work is to validate the relevance of hybrid algorithms for least-square fitting of PEMFC EI spectra in different operating conditions, and for various MEA materials and flow field geometry. Before that, a first assessment of their performance was done using several usual test functions (§ 3.1) and simulated PEMFC EIS data (§ 3.2). The research data repository [56] contains all of the data and source codes used to generate the results of this article.

### 3.1. Results with test functions

A preliminary study with test functions enables to estimate the performance of optimization algorithms in controlled conditions widely discussed in the literature, by solving well-defined mathematical problems with known solutions.

Several reference test functions [50] of different complexity and characteristics have been used. They are listed in Table 2. For simplicity, the graphics of these functions are omitted.

The dimension of the function was set to $d = 10$ in all cases. Their analytical global minima are known, so that the success rates and the average number of fundamental evaluations can be estimated over 100 executions.

It can be seen in Table 3 that the performance of the three hybrid algorithms depends on the considered problem. However, the GA-NM algorithm outperforms PS-NM and SA-NM in most cases in terms of success rate, while it takes more fundamental function evaluations to converge. The exact value of the evaluation criteria of each method depends on many factors including the algorithmic parameters such as the size of search space, the dimension of the target functions, etc.

### Evolution of residuals between stochastic and deterministic steps

Fig. 8 shows the evolution of the residuals with the number of iterations performed for each stochastic algorithm alone (top) and when hybridized with deterministic Nelder-Mead method (bottom), in the case of the Dixon-Price function which is considered the most coherent with experimental results among test functions (see in conclusion). Note that the last values in the graphs at the top corresponds to the first values

<caption>Table 2</caption>

Mathematical test functions: name, expression and analytical global minimum, all functions are in 10D version.

<table><thead><tr><th>Name</th><th>Expression</th><th>Search area &amp; Global minimum</th><th>Specificity</th></tr></thead><tbody><tr><td>Ackey function<br>[57]</td><td>f(x) = - aexp(-<br>b√(1/d ∑<sub>i=1</sub><sup>d</sup> x<sub>i</sub><sup>2</sup>)) -<br>exp(1/d ∑<sub>i=1</sub><sup>d</sup> cos(cxi)) +<br>a + exp(1)<br>a = 20, b = 0.2, c = 2π</td><td>[-10;10] for<br>all x<sub>i</sub><br>f(x*) = 0, at<br>x* = (0,...,0)</td><td>a nearly flat<br>outer region,<br>and a large<br>hole at the<br>center</td></tr><tr><td>Dixon-Price<br>function<br>[58]</td><td>f(x) = (x<sub>1</sub> - 1)<sup>2</sup> +<br>∑<sub>i=2</sub><sup>d</sup> i(2x<sub>i</sub><sup>2</sup> - x<sub>i-1</sub>)<sup>2</sup> +<br>a*sin*<sup>b</sup>(x - x*)</td><td>[-10;10] for<br>all x<sub>i</sub><br>f(x*) = 0,<br>at x<sub>i</sub> = 2<sup>2i-2</sup>/2,<br>for i = 1,...,d</td><td>a valley<br>shaped<br>function, with<br>a sinusoidal<br>noise</td></tr><tr><td>Griewank<br>function<br>[51]</td><td>f(x) = ∑<sub>i=1</sub><sup>d</sup> x<sub>i</sub><sup>2</sup> / 4000 -<br>Π<sub>i=1</sub><sup>d</sup> cos(x<sub>i</sub>/√i) + 1</td><td>[-10;10] for<br>all x<sub>i</sub><br>f(x*) = 0, at<br>x* = (0,...,0)</td><td>many<br>widespread<br>local minima<br>regularly<br>distributed</td></tr><tr><td>Michalewicz<br>function<br>[59]</td><td>f(x) = -<br>∑<sub>i=1</sub><sup>d</sup> sin(xi)sin<sup>2m</sup>(bx<sub>i</sub><sup>2</sup>/π)</td><td>[-5;5] for all x<sub>i</sub><br>at d = 10,<br>f(x*) = -9.66015</td><td>Steep ridges<br>and drops</td></tr><tr><td>Rastrigin<br>function<br>[60]</td><td>f(x) = 10d + ∑<sub>i=1</sub><sup>d</sup> [x<sub>i</sub><sup>2</sup> -<br>10cos(2πx<sub>i</sub>)]</td><td>[-5;5] for all x<sub>i</sub><br>f(x*) = 0, at<br>x* = (0,...,0)</td><td>Many<br>regularly<br>distributed<br>local minima<br>local minima<br>with different<br>amplitudes</td></tr><tr><td>Schwefel<br>function<br>[61]</td><td>f(x) = 418.9829d -<br>∑<sub>i=1</sub><sup>d</sup> x<sub>i</sub>sin(√|x<sub>i</sub>|)</td><td>[0;π] for all x<sub>i</sub><br>f(x*) = 0,<br>at x* = (420.9687,<br>...420.9687)</td><td>A hypercube<br>shape</td></tr><tr><td>Sum of<br>different<br>power<br>(sumpow<br>functions<br>[50]</td><td>f(x) = ∑<sub>i=1</sub><sup>d</sup> |x<sub>i</sub>|<sup>i+1</sup></td><td>[-1;1] for all x<sub>i</sub><br>f(x*) = 0, at<br>x* = (0,...,0).</td><td></td></tr><tr><td>Trid function<br>[62]</td><td>f(x) = ∑<sub>i=1</sub><sup>d</sup> (x<sub>i</sub> - 1)<sup>2</sup> -<br>∑<sub>i=2</sub><sup>d</sup> x<sub>i</sub>x<sub>i-1</sub></td><td>[-100;100] for<br>all x<sub>i</sub><br>f(x*) = -d(d + 4)(d -<br>1)/6, at x<sub>i</sub> =<br>i(d + 1 - i),<br>for all i = 1, 2,<br>...d</td><td>bowl shaped<br>with single<br>minimum</td></tr></tbody></table>
<table><caption>**Table 3**
Results of the hybrid algorithms with the test functions of Table 2: success rate and fundamental evaluations.</caption><thead><tr><th>Function</th><th>GA-NM Success rate</th><th>GA-NM average fundamental evaluation</th><th>PS-NM Success rate</th><th>PS-NM average fundamental evaluation</th><th>SA-NM Success rate</th><th>SA-NM average fundamental evaluation</th></tr></thead><tbody><tr><td>Ackley</td><td>100 %</td><td>897,250</td><td>100 %</td><td>652,671</td><td>27.00 %</td><td>8741</td></tr><tr><td>Dixon-Price</td><td>22.00 %</td><td>718,250</td><td>70 %</td><td>467,955</td><td>5 %</td><td>7469</td></tr><tr><td>Griewank</td><td>60 %</td><td>411,450</td><td>32 %</td><td>410,719</td><td>30 %</td><td>6858</td></tr><tr><td>Michealwicz</td><td>100 %</td><td>482,375</td><td>5 %</td><td>578,323</td><td>0 %</td><td></td></tr><tr><td>Rasrigin</td><td>100 %</td><td>571,210</td><td>100 %</td><td>815,430</td><td>0 %</td><td></td></tr><tr><td>Schweel</td><td>100 %</td><td>518,244</td><td>100 %</td><td>505,712</td><td>20 %</td><td>6912</td></tr><tr><td>Sumpow</td><td>100 %</td><td>411,273</td><td>100 %</td><td>375,741</td><td>100 %</td><td>10,890</td></tr><tr><td>Trid</td><td>100 %</td><td>625,939</td><td>100 %</td><td>482,714</td><td>100 %</td><td>5953</td></tr></tbody></table>

<figure><img src="image_11.png" /><figcaption>**Fig. 8.** Example of evolution of residuals as functions of the number of iterations performed for Dixon-Price function. Top: stochastic algorithms alone (from left to right: GA, PS and SA), bottom: hybrid stochastic/deterministic algorithms (from left to right: GA-NM, PS-NM and SA-NM).</figcaption></figure>

in the graphs at the bottom. With the same function tolerance for all algorithms, a significant improvement of the residuals is obtained by combining stochastic to deterministic algorithms: stochastic algorithms converge slowly due to their random nature, especially in the later stages of optimization when precision is required. Hybrid algorithms take advantage of the deterministic steps to accelerate convergence as the solution approaches an optimum.

Note that the behavior of the algorithm when applied to real-world problems could differ from that observed with test functions, that may not be sufficiently complex and diverse compared to real-world applications. Nevertheless, test functions are useful to assess the efficiency and robustness of optimization algorithms in tackling various challenges commonly encountered in real-world optimization scenarios [48,61,63, 64], such as smoothness, noise, multimodality, correlation between parameters, etc. This point is discussed in the following sections.

## 3.2. Results from simulated & experimental PEMFC impedance data

To validate the efficiency of the hybrid algorithms as well as compare their performances for the interpretation of PEMFC impedance spectra, they were tested on a variety of simulated and experimental data using equivalent electrical circuits (EEC) of different complexity.

To be representative and showing the versatility of the numerical methods, previously published experimental data sets were used [1, 29–31]. Those data sets were measured on PEMFC with various **Membrane-Electrode Assemblies** (MEA), flow field geometry, and operating conditions. Depending on the shape of the spectra, different EEC were tested with 6 (low difficulty) to 11 (high difficulty) parameters. The following sections present the results obtained with two representative and distinct application examples:

* **Case A** [29]: a 30 cm² (10 mm × 30 mm) cell with 5 straight parallel gas channels using a commercial MEA based on Nafion® XL 100 membrane, Pt catalyst supported on high surface area carbon (HSAC, 800 m² g⁻¹) with loadings of 0.1 mgPt·cm⁻² at the anode and 0.3 mgPt·cm⁻² at the cathode. Sigracet 28 BCE gas diffusion layers were used compressed to 175 µm using PTFE gaskets. The cell was operated at 80 °C and fed with H₂ and air in counter-flow with stoichiometries set to 1.5 and 2, respectively and humidified to 50 % RH at the cell inlet. The outlet (absolute) pressure was set to 1.5 bar on both sides. The MEA was subjected to 240 h of Accelerated Stress Test (AST) that significantly damaged the cathode and the anode with noticeable impact on the EI spectra. Impedance data were measured in galvanostatic mode at 0.5 A·cm⁻² with a perturbation amplitude
limited to 50 mA (peak to peak), with frequencies ranging from 20 mHz to 10 kHz, and 10 points per decade.

* **Case B** [33]: a 25 cm² (50 mm × 50 mm) cell with 5 serpentine parallel channels using a MEA made by decal transfer by CEA-LITEN on 820.15 GORE membranes. The catalyst layers were made of commercial Tanaka carbon supported catalyst (Vulcan XC72, 50 wt % Pt, TEC10V50E) with a Pt loading of 0.1 mgpt.cm⁻² for both electrodes. Sigracet 22 BB gas diffusion layers were used compressed by 20 % using PTFE gaskets. The fuel cell temperature was maintained at 80 °C and it was fed with hydrogen and air in counter-flow with an inlet relative humidity of 50 % for both gases. The stoichiometry ratio was set to 4 and 1.5 for air and hydrogen, respectively. The impedance data were measured in galvanostatic mode at 1 A. cm⁻² with a perturbation amplitude of 10 % and for frequencies ranging from 10 mHz to 10 kHz. More information can be found in [31].

### 3.2.1. Equivalent electrical circuits (EEC)

The shape of the impedance spectra depended sensitively if they were measured with cell A or cell B giving a first hint of the most appropriate model. Different EEC were thus used to interpret the data measured with each of both cells, but all consist in series and parallel combinations of several unit elements, as described in the following.

The major contributions to the impedance of a PEMFC stem from the reaction kinetics of the Oxygen Reduction Reaction (ORR), oxygen transport limitation at the cathode side, ion resistance (mainly) through the electrolyte membrane and electronic resistance through the other components and interfaces in series. These contributions are usually described using a Randles EEC [1,4], consisting in:

* **A high frequency resistance R<sub>hf</sub>** accounting for the membrane and other ohmic resistances of and between the MEA components, connected in series with the circuit representing the CCL.

* **A Warburg element Z<sub>w</sub>** accounting for oxygen diffusion limitations, in series with a **charge transfer resistance R<sub>ct</sub>**, both in parallel with a **double-layer capacitance C<sub>dl</sub>**, the two latter representing the reaction kinetics of the ORR.

Introduced by E.Warburg [7], the finite Warburg element (bounded diffusion impedance) is based on simple assumptions: Fickian binary diffusion, surface reaction at the active interface (with the membrane) with a constant flux (driven by the cell current) and a constant gas concentration at the passive interface, i.e. the interface with the gas diffusion medium or channel. Its impedance can be expressed as [2]:

$$Z_{W}(\omega) = \frac{R_{d}}{\sqrt{i\omega\tau}} \tanh\sqrt{i\omega\tau} \quad (5)$$

With τ the characteristic diffusion time and R<sub>d</sub> the oxygen diffusion resistance. In the standard/default Randles EEC, the Warburg element is set within the ORR circuit (Fig. 9, left), implying that the CCL is at the

origin of oxygen diffusion limitations [1,65]. The circuit can also be modified with the Warburg in series with the ORR circuit (Fig. 9, right), which is more consistent when/if the GDL is the most limiting layer for oxygen diffusion [1,29,66].

The impedance expression of the default Randles circuit (Randles-CCL, Fig. 9, left) is [4]:

$$Z_{\text{Randles-CCL}}(\omega) = R_{\text{hf}} + \left( \frac{1}{R_{\text{ct}} + Z_{w}(\omega)} + i\omega C_{\text{dl}} \right)^{-1} \quad (6)$$

And the one of the modified Randles EEC with the Warburg in series (Randles-GDL, Fig. 9, right) is [4]:

$$Z_{\text{Randles-GDL}}(\omega) = Z_{\tau}(\omega) + Z_{\text{GDL}}(\omega) = R_{\text{hf}} + \left( \frac{1}{R_{\text{ct}}} + i\omega C_{\text{dl}} \right)^{-1} + Z_{w}(\omega) \quad (7)$$

With Z<sub>τ</sub> (ω) = R<sub>hf</sub> + (1/R<sub>ct</sub> + iωC<sub>dl</sub>)<sup>-1</sup> the reduced Randles circuit, and Z<sub>w</sub>(ω) the Warburg impedance (Eq. (7)).

An alternative to the Warburg impedance for oxygen diffusion impedance in the GDL is the Kulikovsky impedance Z<sub>K</sub> [67]. This expression considers that the double-layer capacitance of the porous electrode influences the faradaic current density and thus the oxygen flux through the adjacent GDL. It can be written as a function of the finite Warburg element Z<sub>W</sub>(ω) and the impedance parameters of the ORR, i.e. R<sub>ct</sub> and C<sub>dl</sub>:

$$Z_{\text{K}}(\omega) = \frac{Z_{\text{W}}(\omega)}{1 + i\omega C_{\text{dl}} R_{\text{ct}}} \quad (8)$$

Independently of the position and expression of the diffusion impedance, the EIS spectra associated with Randles-type circuits look like a combination of two semicircles, associated respectively to the mass transfer phenomena at low frequencies and the cathode electrochemical reaction kinetics at intermediate frequencies. The intercept with the real axis at high frequencies corresponds to the value of R<sub>hf</sub> (Fig. 10, top).

Although adapted for most commercial MEA, the Randles EEC does not account for the volumetric effects linked to the electrode thickness since it is based on the assumption of surface reaction: this means that it is suitable only for thin electrodes [68,69]. This results in an 90° intercept of the EIS spectra with the horizontal axis in the Nyquist plot at high frequency (Fig. 10, top).

When the volumetric nature of the electrodes (i.e. its thickness) is no longer negligible, a 45° line in the high-frequency region of the EIS spectra (Fig. 10, bottom) is visible. It is generally attributed to ion transfer through the CCL. To describe such EIS spectra, Transmission Line Models (TLM) are used. They were first introduced by de Levie in 1963 [68] to describe the electrochemical behavior of porous thick electrodes. TLM assume that the electrode consists of an assembly of elementary horizontal cylindrical pores of identical diameter and properties, with a length equal to the electrode thickness. TLM EEC are composed of an infinite assembly of unit circuits connected in series, i.e. a transmission line. The circuit introduced by de Levie was used to model a blocked electrode (zero net current) like those of

<figure><img src="image_12.png" /><figcaption>Fig. 9. Randles EEC in default configuration, Randles-CCL (left) and in GDL configuration, Randles-GDL (right).</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Fig. 10. Nyquist diagram of impedance data measured on PEMFC with commercial (top) and hot-pressed MEA using Gas Diffusion Electrodes GDE (bottom), and corresponding fitting curves obtained with a Randles-type (top) and a TLM-like EEC (bottom) [1].</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 11. TLM-type EEC of a PEMFC. The electrode thickness is represented by a succession of unit Randles circuits connected by electric resistances, the latter attributed here to ion conduction through the ionomer in the electrode pores [1].</figcaption></figure>

supercapacitors. An example of a TLM-like EEC suitable for FC electrodes, such as those of a PEMFC is shown in Fig. 11 [1]. Each elementary circuit consists of a Randles EEC interconnected via ionic resistances Rion accounting for ion transport through the electrolyte phase in the catalyst layer. Note that the resistance of the electron conducting phase (carbon) in PEMFC electrodes are usually considered as negligible compared that of the ion conducting phase (ionomer) [1].

$$Z_{\text{cell}} = R_{\text{HF}} + Z_{\text{CCL}}(\omega) = R_{\text{HF}} + Z_{\text{TLM}}(\omega) \quad (9)$$

With:

$$Z_{\text{TLM}}(\omega) = \frac{\sqrt{R_{\text{ion}}}}{\sqrt{i\omega C_{\text{dl}} + \left(\frac{1}{(R_{\text{ct}} + Z_{\text{W}}(\omega))}\right)}} \coth\left(\sqrt{i\omega R_{\text{ion}} C_{\text{dl}} + \left(\frac{R_{\text{ion}}}{(R_{\text{ct}} + Z_{\text{W}}(\omega))}\right)}\right) \quad (10)$$

When Rion is sufficiently small compared to Z(Rct||Cdl), meaning ion transport being negligible compared to the reaction kinetics limitations,

the expression of TLM impedance ZTLM(ω) in Eq. (10) tends to the expression of the Randles circuit in Eq. (6) [1]. In this case the 45° slope at high frequencies tends to a 90° slope.

Like in the case of Randles EEC, the diffusion impedance can also be in series with the TLM if the GDL is the limiting layer for oxygen diffusion:

$$Z_{\text{cell}} = R_{\text{HF}} + Z_{\text{CCL}}(\omega) + Z_{\text{GDL}}(\omega) = R_{\text{HF}} + Z_{\text{TLM}}(\omega) + Z_{\text{GDL}}(\omega) \quad (11)$$

$$Z_{\text{TLM}}(\omega) = \frac{\sqrt{R_{\text{ion}}}}{\sqrt{i\omega C_{\text{dl}} + \left(\frac{1}{R_{\text{ct}}}\right)}} \coth\left(\sqrt{i\omega R_{\text{ion}} C_{\text{dl}} + \left(\frac{R_{\text{ion}}}{R_{\text{ct}}}\right)}\right) \quad (12)$$

In this case, ZGDL(ω) can be expressed either by a Warburg impedance ZW (Eq. (5)) or by a Kulikovsky impedance ZK (Eq. (8)).

Note that it is also possible to include two oxygen diffusion impedances in the EEC: one for the CCL and one for the GDL. However, there would be strong correlations between their characteristic parameters. These questions are beyond the scope of this work, like the discussion on
the relative importance of these diffusion impedances. Hence, for simplicity, only one oxygen diffusion impedance is used in the following.

In some other cases, the anode has also to be considered for modeling properly the PEMFC impedance [29]. The impedance of hydrogen diffusion being always negligible [70], the anode is modelled through a simple parallel connection of a charge transfer resistance $R_{\text{cta}}$ and a double-layer capacitance $C_{\text{dla}}$ representing the kinetics of the Hydrogen Oxidation Reaction (HOR):

$$Z_{\text{anode}}(\omega) = \left( \frac{1}{R_{\text{cta}}} + i\omega C_{\text{dla}} \right)^{-1} \quad (13)$$

Finally, an inductance $Z_L(\omega) = i\omega L$ is sometimes needed at the highest frequencies as inductive phenomena of the electric connections can alter the shape of the spectra [71].

### 3.2.2. Case A: PEMFC with aged anode-Randles vs. TLM-type circuits with simulated & experimental data

#### Performance tests with simulated data

Performance validation and evaluation of algorithms can be carried out by solving already known problems. First, the 3 hybrid methods (GA-NM, PS-NM and SA-NM) were tested on simulated EIS spectra generated from the parameter values in Table 4. Those values were identified from experimental PEMFC data using standard Nelder-Mead algorithm with a Randles-Anode model (noted as RA in the following parts, Fig. 12) with 7 parameters [1,29]. The corresponding impedance expression is:

$$Z_{\text{cell}} = Z_{\text{Randles}} + Z_{\text{anode}} \quad (14)$$

With $Z_{\text{Randles}}$ and $Z_{\text{anode}}$ given by Eqs. (5) and (13), respectively. Eq. (4) can thus be considered as a test function for EIS. The spectra were simulated with various noise levels following a gaussian distribution with a standard deviation up to about 50 %. This enables to evaluate the stability of the algorithms against noise and the sensitivity of each parameter.

The bounds of $R_{ct}$ and $R_d$ were set to ± 50 % of the values estimated visually, i.e. the diameter of the corresponding semicircle in the Nyquist plot (Fig. 1). The upper bound of $R_{ct,a}$ is set to lower bound value of $R_{ct}$, i. e. 0.1 Ω.cm². This choice is made since the high frequency loop associated with the HOR the anode appears relatively small in the spectrum compared to that at intermediate frequencies associated to ORR kinetics at the cathode (Fig. 13 and Fig. 16). The search intervals of the other parameters were set to [0,1], keeping in mind that most of the actual values are much lower than one: large intervals slow down the convergence but may help to highlight the advantages and differences between the optimization algorithms.

**Table 6:** Best solution with minimum residual and occurrences of 100 runs for each hybrid algorithm. shows the best solution with minimum residuals and its occurrence number when repeating each hybrid algorithm 100 times, with a noise level varying between 5 % and 50 %. For clarity only the best solution (minimum residuals) is shown.

The three hybrid algorithms return identical parameter values and show identical fit quality (i.e. residuals) up to 50 % noise. The estimated values are almost identical to the input parameters of the simulated spectra, with a tiny difference (lower than 2 %) as the noise increases. In terms of reliability, GA-NM and PS-NM algorithms converge more often to the correct solution than SA-NM (about 70 % vs. 30 % occurrences). Whereas in terms of efficiency, SA-NM needs globally less function

<caption>**Table 4**</caption>

Parameter values of the RA EEC in Fig. 12 identified from experimental data of a PEMFC subjected to 240 h of AST (Figs. 13 and 16) using standard NM algorithm [29]. These values are used for the simulated spectra of case A.

<table><thead><tr><th>R<sub>hf</sub><br>(Ω.cm<sup>2</sup>)</th><th>R<sub>ct</sub><br>(Ω.cm<sup>2</sup>)</th><th>C<sub>dla</sub><br>(F/cm<sup>2</sup>)</th><th>R<sub>d</sub><br>(Ω.cm<sup>2</sup>)</th><th>τ(s)</th><th>R<sub>ct,a</sub><br>(Ω*cm<sup>2</sup>)</th><th>C<sub>dla,a</sub><br>(F/cm<sup>2</sup>)</th></tr></thead><tbody><tr><td>0.0728</td><td>0.250</td><td>0.0179</td><td>0.325</td><td>0.114</td><td>0.0555</td><td>0.0051</td></tr></tbody></table>

evaluations than GA-NM and PS-NM, these later behaving similarly: about 10³ for SA-NM vs. 10⁷ to 10⁸ fundamental evaluations for the two others (not shown in Table 6: Best solution with minimum residual and occurrences of 100 runs for each hybrid algorithm. for readability reasons). This is consistent with the observations on test functions of section 3.1 and those reported in the literature [45].

#### Performance tests with experimental data

The estimation is then repeated starting from the experimental dataset [29] -with identical bound limits as shown in Table 5- allowing to account for real noise distribution and measurement errors.

The 3 hybrid algorithms converged toward the 2 solutions shown in Table 7 sorted by ascending order of residuals, i.e. quality. The experimental and fitted spectra are shown in Fig. 13.

The first solution is identical to that found in our previous work using NM algorithm alone (Table 4), while the second solution has 2 parameters reaching the bound limits, which makes it invalid. As shown in Table 8 and in the pie chart in Fig. 14, the highest success rates for the solution of best quality (lowest residuals) are obtained with GA-NM and PS-NM algorithms (about 70 % and 80 %, respectively), while it is lower with SA-NM (≈ 18 %). This confirms the higher reliability of the two former algorithms compared to the latter. The lowest number of fundamental evaluations is observed with SA-NM algorithm (nEval in Table 8) which confirms again the higher efficiency of this method compared to the two others. This agrees with our observations on simulated data (Table 6).

Table 8 also shows that the stochastic step (GA/PS/SA) takes significantly more fundamental evaluations than the deterministic step (NM), confirming that NM can significantly accelerate the convergence. Note that in combination with SA, NM needs about twice as many evaluations than with GA or PS, which may be explained by the lower reliability of this stochastic algorithm.

#### Towards more complex EEC models

As shown in Fig. 16, using a Randles circuit with an anode in series (RA model) does not make it possible to fit satisfyingly experimental spectrum at high frequencies. In our previous work [29], the NM algorithm used for data interpolation did not allow the use of a more complex EEC.

Thanks to their robustness, the hybrid algorithms enable to fit the data using more complex EEC that are numerically and -most probably- also physically more appropriate. A variety of combinations were tested. It turned out that for the present case, the circuit described in Fig. 15 was the most suited for the data of case A: it consists of a TLM-like circuit for the cathode accounting only for the ORR kinetics (without diffusion impedance), a high frequency resistance and a Warburg element in series representing diffusion through the GDL (Equation (Z<sub>cell</sub> = R<sub>HF</sub> + Z(ω) + ZGDL(ω) = R<sub>HF</sub> + Ztlm(ω) + ZGDLω (11)). The circuit is completed by an anode (13) and an inductance. This TLM-WarburgGDL-Anode-Inductance model is called tWAL -for short- in the following. It contains two more parameters than the RA circuit, R<sub>ion</sub> and L, accounting for ion conduction through the cathode thickness and inductive phenomena of the connecting wires, respectively. These parameters impact the shape of the spectrum at high frequencies. The oxygen diffusion impedance was put in series (in the GDL) for a better consistency with our previous results [1,31]. The position of the diffusion impedance impacts mostly the values of the diffusion impedance parameters, rather than the shape of the spectrum.

The interval limits were set to [0;1] for both, R<sub>ion</sub> and L. 100 repetitions of the hybrid algorithms resulted in 6 main solutions, found with at least 1 of the 3 algorithms (Table 9). However, only the first solution is numerically and physically coherent: lowest residual value, consistent fitting curves, as well as parameter values within the bound limits.

**Table 10** shows the values of the evaluation criteria (reliability/occurrence and efficiency/nEval) of the solutions identified for case A for each hybrid algorithms using the TLM-like model. In addition, Fig. 17 shows the pie charts of the solutions found with each method. All 3 methods return the same set of solutions, keeping in mind that these
<figure><img src="image_15.png" /><figcaption>Fig. 12. Randles-Anode (RA) circuit [27].</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 13. Experimental and fitted impedance spectra of the different numerical solutions for case A, shown in Table 7 obtained with the RA model in Fig. 12, experimental data taken from [27].</figcaption></figure>

<caption>**Table 5**</caption>

Lower (lb) and upper bound (ub) limits of the EEC parameters for case A: around ±50 % of the approximative values of $R_{ct}$ and $R_d$ estimated visually from the Nyquist plot, a maximum value for $R_{ct,a}$ corresponding to the lb value of $R_{ct}$, and [0,1] (large intervals) for the others.

<table><thead><tr><th></th><th>R<sub>hf</sub><br>(Ω.cm<sup>2</sup>)</th><th>R<sub>ct</sub><br>(Ω.cm<sup>2</sup>)</th><th>C<sub>dl</sub><br>(F/cm<sup>2</sup>)</th><th>R<sub>d</sub><br>(Ω.cm<sup>2</sup>)</th><th>τ (s)</th><th>R<sub>ct,a</sub><br>(Ω<sup>2</sup>cm<sup>2</sup>)</th><th>C<sub>dl,a</sub><br>(F/cm<sup>2</sup>)</th></tr></thead><tbody><tr><td>lb</td><td>0</td><td>0.1</td><td>0</td><td>0.1</td><td>0</td><td>0</td><td>0</td></tr><tr><td>ub</td><td>1</td><td>0.3</td><td>1</td><td>0.4</td><td>1</td><td>0.1</td><td>1</td></tr></tbody></table>

solutions cannot be compared or tested on absolute basis without any physical interpretation.

Again, the three algorithms were able to find the best solution in terms of fitting quality (lowest residual), this time with relatively low and comparable reliability (occurrences): 9 % for GA-NM and 8 % with PS-NM and SA-NM. Note that this time SA-NM shows comparable performances to the two other methods in terms of reliability, and once again the lowest number of fundamental evaluations, i.e. highest efficiency.

### 3.2.3. Case B: PEMFC with low Pt-cathode loading-TLM-type circuits with different oxygen diffusion impedances applied to experimental data

Here, the objective was to identify the most appropriate circuit to fit

satisfyingly the entire EI spectrum of case B (Fig. 18). In a previous work [31] based on NM algorithm alone, the data had to be fitted using a Randles-type circuit, since NM alone was not able to identify a solution with the more complex TLM-like models. As shown on Fig. 18, using Randles-type circuit did not allow to satisfyingly fit the impedance data at high frequencies. This was acceptable as the work focused on oxygen diffusion impedance influencing the spectra at low frequencies. It is now possible to fit the data with more complex TLM-type circuits thanks to the robustness of the hybrid algorithms.

The data of case B were fitted using two TLM-like circuits: the same as for case A, i.e. the tWAL model shown in Fig. 15, and a second, called tKAL in the following in which the Warburg diffusion impedance (Eq. (5)) is replaced by a Kulikovsky impedance (Eq. (8)), since this

<figure><img src="image_17.png" /><figcaption>Fig. 14. Pie charts for the 3 hybrid-algorithms of the occurrence of the two solutions for case A given in Table 7 identified with the RA model.</figcaption></figure>
<table><caption>**Table 6**
Best solution with minimum residual and occurrences of 100 runs for each hybrid algorithm.</caption><thead><tr><th>Noise</th><th>R<sub>hf</sub> (Ω. cm<sup>2</sup>)</th><th>R<sub>ct</sub> (Ω. cm<sup>2</sup>)</th><th>C<sub>dl</sub> (F/cm<sup>2</sup>)</th><th>R<sub>d</sub> (Ω. cm<sup>2</sup>)</th><th>τ (s)</th><th>R<sub>ct.a</sub> (Ω. cm<sup>2</sup>)</th><th>C<sub>dl.a</sub> (F/cm<sup>2</sup>)</th><th>residuals</th><th>Exploit Occurrence GA-NM</th><th>Exploit Occurrence PS-NM</th><th>Exploit Occurrence SA-NM</th></tr></thead><tbody><tr><td>5 %</td><td>0.0728</td><td>0.251</td><td>0.018</td><td>0.324</td><td>0.114</td><td>0.0557</td><td>5.10 × 10<sup>-3</sup></td><td>1.50 × 10<sup>-4</sup></td><td>72</td><td>77</td><td>31</td></tr><tr><td>10 %</td><td>0.0727</td><td>0.249</td><td>0.0178</td><td>0.326</td><td>0.113</td><td>0.0553</td><td>5.13 × 10<sup>-3</sup></td><td>6.78 × 10<sup>-4</sup></td><td>67</td><td>69</td><td>24</td></tr><tr><td>15 %</td><td>0.0726</td><td>0.248</td><td>0.0181</td><td>0.325</td><td>0.115</td><td>0.0576</td><td>5.10 × 10<sup>-3</sup></td><td>1.66 × 10<sup>-3</sup></td><td>75</td><td>70</td><td>27</td></tr><tr><td>20 %</td><td>0.0727</td><td>0.25</td><td>0.0176</td><td>0.326</td><td>0.114</td><td>0.0543</td><td>4.97 × 10<sup>-3</sup></td><td>2.63 × 10<sup>-3</sup></td><td>74</td><td>73</td><td>28</td></tr><tr><td>30 %</td><td>0.0727</td><td>0.25</td><td>0.0176</td><td>0.326</td><td>0.114</td><td>0.0543</td><td>4.97 × 10<sup>-3</sup></td><td>2.63 × 10<sup>-3</sup></td><td>73</td><td>73</td><td>28</td></tr><tr><td>50 %</td><td>0.0723</td><td>0.254</td><td>0.0178</td><td>0.323</td><td>0.113</td><td>0.0544</td><td>5.01 × 10<sup>-3</sup></td><td>5.65 × 10<sup>-3</sup></td><td>75</td><td>76</td><td>27</td></tr></tbody></table>

<table><caption>**Table 7**
Solutions found for case A starting from experimental data (Fig. 13) with corresponding residuals [29]. Red = parameter value reached its bound limit.</caption><thead><tr><th>Solution</th><th>R<sub>hf</sub> (Ω.cm<sup>2</sup>)</th><th>R<sub>ct</sub> (Ω.cm<sup>2</sup>)</th><th>C<sub>dl</sub> (F/cm<sup>2</sup>)</th><th>R<sub>d</sub> (Ω.cm<sup>2</sup>)</th><th>τ (s)</th><th>R<sub>ct.a</sub> (Ω.cm<sup>2</sup>)</th><th>C<sub>dl.a</sub> (F/cm<sup>2</sup>)</th><th>residual</th></tr></thead><tbody><tr><td>1</td><td>0.0728</td><td>0.25</td><td>0.0179</td><td>0.325</td><td>0.114</td><td>0.056</td><td>0.0051</td><td>0.0017</td></tr><tr><td>2</td><td>0.0866</td><td>0.125</td><td>0.0099</td><td>0.4</td><td>0.111</td><td>0.1</td><td>0.0744</td><td>0.0121</td></tr></tbody></table>

<table><caption>**Table 8**
Evaluation criteria (reliability/occurrence and efficiency/nEval) of each hybrid method for the solutions of case A, shown in Table 7 using the RA model (Fig. 12).</caption><thead><tr><th>Minimum</th><th>GA-NM<br>Occurrence</th><th>Efficiency<br>(nEval_GA)</th><th>Efficiency<br>(nEval_NM)</th><th>PS-NM<br>Occurrence</th><th>Efficiency<br>(nEval_PS)</th><th>Efficiency<br>(nEval_NM)</th><th>SA-NM<br>Occurrence</th><th>Efficiency<br>(nEval_SA)</th><th>Efficiency<br>(nEval_NM)</th></tr></thead><tbody><tr><td>1 (best)</td><td>71</td><td>1 639 663</td><td>609</td><td>80</td><td>400 316</td><td>488</td><td>18</td><td>5 170</td><td>1 361</td></tr><tr><td>6</td><td>20</td><td>1 942 203</td><td>572</td><td>20</td><td>343 571</td><td>495</td><td>82</td><td>4 923</td><td>2 021</td></tr></tbody></table>

component is physically more adapted than a Warburg impedance to describe oxygen diffusion impedance of the GDL [24,31]. However, even if the fundamental EEC parameters and their number are identical

for both diffusion impedances, the Kulikovsky impedance depends on $R_{ct}$ and $C_{dl}$, which increases the correlation between the model parameters, making their estimation more difficult. Therefore, more restricted

<figure><img src="image_18.png" /><figcaption>**Fig. 15.** TLM-like model with Warburg element in the GDL, anode and inductance (tWAL model) used for fitting the impedance data of case A, Fig. 16.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>**Fig. 16.** Nyquist plot of experimental impedance data of case A and fitting curves of the best solution of RA and tWAL models, with zoom to high frequencies.</figcaption></figure>

**Table 9**

Main solutions for case A when fitting the experimental data with the TLM-like circuit in Fig. 15.

<table><thead><tr><th>Minimum</th><th>R<sub>hf</sub><br>(Ω.cm<sup>2</sup>)</th><th>R<sub>ct</sub><br>(Ω.cm<sup>2</sup>)</th><th>C<sub>dl</sub><br>(F/cm<sup>2</sup>)</th><th>R<sub>ion</sub><br>(Ω.cm<sup>2</sup>)</th><th>R<sub>ds</sub><br>(Ω.cm<sup>2</sup>)</th><th>τ<sub>s</sub><br>(s)</th><th>R<sub>ct,a</sub><br>(Ω*cm<sup>2</sup>)</th><th>C<sub>dl,a</sub><br>(F/cm<sup>2</sup>)</th><th>L<br>(H*cm<sup>2</sup>)</th><th>residual</th></tr></thead><tbody><tr><td>1</td><td>0.054</td><td>0.315</td><td>0.170</td><td>0.157</td><td>0.269</td><td>0.0164</td><td>0.0144</td><td>9.05×10<sup>-3</sup></td><td>5.09×10<sup>-7</sup></td><td>8.46×10<sup>-4</sup></td></tr><tr><td>2</td><td>0.0525</td><td>0.197</td><td>0.0323</td><td>0.170</td><td>0.4</td><td>0.132</td><td>9.33×10<sup>-3</sup></td><td>0.0106</td><td>5.27×10<sup>-7</sup></td><td>8.98×10<sup>-4</sup></td></tr><tr><td>3</td><td>0.0571</td><td>0.193</td><td>0.0305</td><td>0.191</td><td>0.342</td><td>0.127</td><td>0.0515</td><td>1</td><td>3.60×10<sup>-7</sup></td><td>1.30×10<sup>-3</sup></td></tr><tr><td>4</td><td>0.0572</td><td>0.181</td><td>0.032</td><td>0.181</td><td>0.4</td><td>0.13</td><td>0.0187</td><td>1</td><td>3.58×10<sup>-7</sup></td><td>1.40×10<sup>-3</sup></td></tr><tr><td>5</td><td>0.0571</td><td>0.1</td><td>0.118</td><td>0.725</td><td>0.368</td><td>0.136</td><td>4.95×10<sup>-3</sup></td><td>9.91×10<sup>-3</sup></td><td>4.91×10<sup>-7</sup></td><td>1.50×10<sup>-3</sup></td></tr><tr><td>6</td><td>0.0565</td><td>0.1</td><td>0.109</td><td>0.823</td><td>0.251</td><td>0.136</td><td>0.1</td><td>0.535</td><td>3.63×10<sup>-7</sup></td><td>1.70×10<sup>-3</sup></td></tr></tbody></table>

**Table 10**

Evaluation criteria (reliability/occurrence and efficiency/nEval) of each hybrid method for the solutions of case A, shown in Table 9 using the tWAL model (Fig. 15).

<table><thead><tr><th>Minimum</th><th>GA-NM<br>Occurrence</th><th>Efficiency<br>(nEval GA)</th><th>Efficiency<br>(nEval NM)</th><th>PS-NM<br>Occurrence</th><th>Efficiency<br>(nEval PS)</th><th>Efficiency<br>(nEval NM)</th><th>SA-NM<br>Occurrence</th><th>Efficiency<br>(nEval SA)</th><th>Efficiency<br>(nEval NM)</th></tr></thead><tbody><tr><td>1 (best)</td><td>9</td><td>4 311 803</td><td>3 385</td><td>8</td><td>579 375</td><td>2 091</td><td>8</td><td>5 116</td><td>5 472</td></tr><tr><td>2</td><td>23</td><td>4 136 717</td><td>5 816</td><td>0</td><td>-</td><td>-</td><td>7</td><td>5 075</td><td>6 787</td></tr><tr><td>3</td><td>4</td><td>1 850 150</td><td>4 369</td><td>5</td><td>2 301 000</td><td>3 918</td><td>15</td><td>5 077</td><td>5 067</td></tr><tr><td>4</td><td>16</td><td>6 443 875</td><td>2 756</td><td>22</td><td>2 291 818</td><td>2 354</td><td>8</td><td>5 109</td><td>6 162</td></tr><tr><td>5</td><td>17</td><td>1 803 544</td><td>3 401</td><td>25</td><td>496 000</td><td>1 734</td><td>27</td><td>5 071</td><td>5 038</td></tr><tr><td>6</td><td>20</td><td>3 606 937</td><td>7 980</td><td>11</td><td>5 140</td><td>5 845</td><td>24</td><td>378 333</td><td>5 858</td></tr></tbody></table>

bound limits were used for the diffusion impedance parameters $R_{ds}$ (Table 11) that was set to ±50 % of the approximative diameter of the low frequency loop as shown in Fig. 1. $R_{hf}$ and $\tau_s$ were set to the values given in [31]. For comparison purpose, the data were also fitted with a Randles-type circuit equivalent to the tKAL model, a Kulikovsky impedance in series with the OKK kinetics circuit (without diffusion impedance), an anode circuit and a high frequency resistance without inductance, called rKA which is an EEC used in [31]. Same bound limits are used for the parameters of the Randles- and TLM-type circuits.

Fig. 18 shows the experimental spectrum and fitting curves obtained with the 3 models, and Table 12 shows the identified parameters. Numerically better solutions are found with the TLM models compared to the Randles-type circuit. The fact that the value of $R_{ion}$ is of the same order than $R_{ct}$ confirms once more the importance to account for ion

conduction through the cathode volume influencing the spectra at high frequencies. It can be observed that the fitting curves of the two TLM-type circuits are visually identical and yield close residuals (2.10 $10^{-6}$ for tWAL and 2.14 $10^{-6}$ for tKAL) for the computationally and physically best solution. The tKAL model is however physically more reasonable for representing gas diffusion through the GDL, which is in agreement with the former conclusions obtained with the Randles-like model [31]. That is why the following discussion focuses on the results obtained with the tKAL.

Table 13 depicts the main solutions found with all the 3 hybrid algorithms for case B using the tKAL model. The particularity of this case is that the first solution, which has the smallest residual value is not correct: the value of $R_{ct}$ has reached the lower bound. With a residual value only 1 % larger, solution 2 is the only one for which all parameter values
<figure><img src="image_20.png" /><figcaption>Fig. 17. Pie charts of the occurrence percentage of the six solutions for case A (Table 9) identified with the tWAL model with the 3 hybrid-algorithms.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Fig. 18. EIS spectra of case B, experimental data and fitting curves of best solutions of each model (rKA, tWAL and tKAL), experimental data taken from [29].</figcaption></figure>

<caption>**Table 11**</caption>

Lower (lb) and upper bound (ub) limits of the EEC parameters for case B: ±50 % of the approximative values of $R_{ct}$ and $R_{ds}$ estimated visually from the Nyquist plot, a maximum value for $R_{ct,a}$ corresponding to the lb value of $R_{ct}$, a maximum value of 0.05 for $\tau_s$, and [0,1] (large intervals) for the others.

<table><caption>**Table 12**</caption><thead><tr><th></th><th>R<sub>hf</sub> (Ω.cm<sup>2</sup>)</th><th>R<sub>ct</sub> (Ω.cm<sup>2</sup>)</th><th>C<sub>dl</sub> (F/cm<sup>2</sup>)</th><th>R<sub>ion</sub> (Ω.cm<sup>2</sup>)</th><th>R<sub>ds</sub> (Ω.cm<sup>2</sup>)</th><th>τ<sub>s</sub> (s)</th><th>R<sub>ctu</sub> (Ω.cm<sup>2</sup>)</th><th>C<sub>dla</sub> (F/cm<sup>2</sup>)</th><th>L (H.cm<sup>2</sup>)</th></tr></thead><tbody><tr><td>Lb2</td><td>0.05</td><td>0.0875</td><td>0</td><td>0</td><td>0.0325</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Ub2</td><td>0.15</td><td>0.2625</td><td>1</td><td>1</td><td>0.13</td><td>0.05</td><td>0.1</td><td>1</td><td>1</td></tr></tbody></table>

Parameter values identified from experimental data of case B (Fig. 18) with different EEC models: rKA, tWAL and tKAL.

<table><thead><tr><th>Solution</th><th>R<sub>hf</sub> (Ω.cm<sup>2</sup>)</th><th>R<sub>ct</sub> (Ω.cm<sup>2</sup>)</th><th>C<sub>dl</sub> (F/cm<sup>2</sup>)</th><th>R<sub>ion</sub> (Ω.cm<sup>2</sup>)</th><th>R<sub>ds</sub> (Ω.cm<sup>2</sup>)</th><th>τ<sub>s</sub> (s)</th><th>R<sub>ctu</sub> (Ω.cm<sup>2</sup>)</th><th>C<sub>dla</sub> (F/cm<sup>2</sup>)</th><th>L (H.cm<sup>2</sup>)</th><th>Residual</th></tr></thead><tbody><tr><td>rKA</td><td>0.111</td><td>0.0857</td><td>0.0034</td><td></td><td>0.0757</td><td>1.70 × 10<sup>-3</sup></td><td>0.0668</td><td>0.238</td><td></td><td>1,24 × 10<sup>-3</sup></td></tr><tr><td>tWAL</td><td>0.0688</td><td>0.1</td><td>0.0070</td><td>0.11</td><td>0.0696</td><td>7.73 × 10<sup>-4</sup></td><td>0.0664</td><td>0.24</td><td>1,45 × 10<sup>-6</sup></td><td>3,99 × 10<sup>-4</sup></td></tr><tr><td>tKAL</td><td>0.0633</td><td>0.0879</td><td>0.0036</td><td>0.188</td><td>0.0665</td><td>1.74 × 10<sup>-3</sup></td><td>0.0664</td><td>0.237</td><td>1.62 × 10<sup>-6</sup></td><td>4.04 × 10<sup>-4</sup></td></tr></tbody></table>

stay within the bound limits and can thus be considered as the most consistent solution.

Concerning algorithmic performance, similar tendency can be observed as for the previous cases (Fig. 19, Table 14): PS-NM found the best solution with the highest number of occurrences, while GA-NM identified a greater variety of solutions overall. Meanwhile, SA-NM demonstrated the lowest number of fundamental evaluations. Furthermore, PS-NM algorithm appears capable of eliminating the physically meaningless solution 1.

## 4. Conclusion & perspectives

Three hybrid stochastic-deterministic optimization algorithms, i.e. GA-NM, PS-NM and SA-NM, have been tested and compared to classical deterministic Nelder-Mead method, and the usual non-hybridized stochastic algorithms for interpreting PEMFC EIS data. To cover a wide range of scenarios, >10 MEA bench setups and 100 datasets were tested, incorporating various MEA materials (membranes, electrodes, GDLs)

and operating conditions (temperature, relative humidity, current density, etc.). The results across these diverse cases were found to be consistent and comparable.

These results lead to the following main conclusions:

* The hybrid approach combines the **global optimization capabilities** of stochastic algorithms with the **local refinement efficiency** of the deterministic Nelder-Mead algorithm. This reduces sensitivity to initial conditions and accelerates convergence, offering a balance between exploration and exploitation.

* With **suitable parameter search intervals**, the 3 methods were able to achieve the main objective, i.e. find at least one physically and numerically good solution over 100 executions.

* The **setting of the parameter search intervals** can be done starting from graphical interpretations of the EI spectra to assess the orders of magnitude of the main parameters. This is much easier and practical than choosing an appropriate starting point when using NM or any deterministic algorithm alone.
<table><caption>**Table 13**
Main solutions found with the 3 hybrid algorithms for case B using the tKAL model.</caption><thead><tr><th>Solution</th><th>R<sub>hf</sub><br>(Ω.cm<sup>2</sup>)</th><th>R<sub>ct</sub><br>(Ω.cm<sup>2</sup>)</th><th>C<sub>dl</sub><br>(F/cm<sup>2</sup>)</th><th>R<sub>ton</sub><br>(Ω.cm<sup>2</sup>)</th><th>R<sub>ds</sub><br>(Ω.cm<sup>2</sup>)</th><th>τ<sub>s</sub><br>(s)</th><th>R<sub>ct,a</sub><br>(Ω.cm<sup>2</sup>)</th><th>C<sub>dl,a</sub><br>(F/cm<sup>2</sup>)</th><th>L<br>(H.cm<sup>2</sup>)</th><th>Residual</th></tr></thead><tbody><tr><td>1</td><td>0.0684</td><td>0.0875</td><td>7.15×10<sup>-3</sup></td><td>0.3</td><td>0.034</td><td>3.34×10<sup>-4</sup></td><td>0.0667</td><td>0.0234</td><td>1.54×10<sup>-6</sup></td><td>2.12×10<sup>-6</sup></td></tr><tr><td>2</td><td>0.0705</td><td>0.127</td><td>4.24×10<sup>-3</sup></td><td>0.127</td><td>0.0778</td><td>0.0396</td><td>0.0241</td><td>6.52×10<sup>-3</sup></td><td>1.55×10<sup>-6</sup></td><td>2.14×10<sup>-6</sup></td></tr><tr><td>3</td><td>0.0687</td><td>0.0875</td><td>5.40×10<sup>-3</sup></td><td>0.202</td><td>0.0782</td><td>0.0395</td><td>0.0462</td><td>0.0138</td><td>1.57×10<sup>-6</sup></td><td>2.15×10<sup>-6</sup></td></tr><tr><td>4</td><td>0.0592</td><td>0.135</td><td>3.98×10<sup>-3</sup></td><td>0.231</td><td>0.0594</td><td>0.0391</td><td>0.016</td><td>1</td><td>1.76×10<sup>-6</sup></td><td>2.16×10<sup>-6</sup></td></tr><tr><td>5</td><td>0.0554</td><td>0.0875</td><td>8.71×10<sup>-3</sup></td><td>0.511</td><td>0.0325</td><td>0.0396</td><td>0.0369</td><td>0.45</td><td>1.98×10<sup>-6</sup></td><td>2.60×10<sup>-6</sup></td></tr><tr><td>6</td><td>0.15</td><td>0.261</td><td>3.71×10<sup>-6</sup></td><td>1</td><td>0.13</td><td>7.04×10<sup>-7</sup></td><td>0.095</td><td>2.63×10<sup>-7</sup></td><td>1</td><td>1.64×10<sup>-2</sup></td></tr></tbody></table>

<figure><img src="image_22.png" /><figcaption>**Fig. 19.** Pie charts of the occurrence of the six solutions for case B (Table 14) identified with the tKAL model with the 3 hybrid-algorithms.</figcaption></figure>

<table><caption>**Table 14**
Occurrences and number of function evaluations of the 3 hybrid algorithms for each solution of case B obtained with the tKAL model.</caption><thead><tr><th>Minimum</th><th>GA-NM<br>Recurrence</th><th>Efficiency<br>(nEval_GA)</th><th>Efficiency<br>(nEval_NM)</th><th>PS-NM<br>recurrence</th><th>Efficiency<br>(nEval_PS)</th><th>Efficiency<br>(nEval_NM)</th><th>SA-NM<br>Recurrence</th><th>Efficiency<br>(nEval_SA)</th><th>Efficiency<br>(nEval_NM)</th></tr></thead><tbody><tr><td>1 (invalid)</td><td>1</td><td>3 154 500</td><td>5 632</td><td>0</td><td></td><td></td><td>7</td><td>5 001</td><td>7 612</td></tr><tr><td>2 (best)</td><td>16</td><td>8 621 453</td><td>12 065</td><td>22</td><td>1 322 500</td><td>9 267</td><td>13</td><td>5 267</td><td>10 302</td></tr><tr><td>3</td><td>9</td><td>4 106 611</td><td>7 800</td><td>30</td><td>2444 666</td><td>10 772</td><td>14</td><td>5 305</td><td>13 831</td></tr><tr><td>4</td><td>11</td><td>2 170 836</td><td>8 305</td><td>11</td><td>1 446 363</td><td>9 104</td><td>20</td><td>4 951</td><td>6 445</td></tr><tr><td>5</td><td>10</td><td>306 711</td><td>3 768</td><td>31</td><td>426 129</td><td>4 368</td><td>4</td><td>5 332</td><td>7 932</td></tr><tr><td>6</td><td>36</td><td>247 500</td><td>905</td><td>0</td><td></td><td></td><td>19</td><td>4 798</td><td>1 285</td></tr></tbody></table>

* Hybrid algorithms can help identify several unexpected possible solutions, allowing researchers to explore potential physical interpretations and select the most relevant models, this can be particularly interesting when characterizing new materials and systems.

The three hybrid methods can therefore help to improve the interpretation of experimental EIS data and compare different EEC to identify the most physically relevant model. However, it has to be noted:

* For complex models or situations, e.g. for new MEA compositions or PEMFC configurations, or otherwise said for unknown order of magnitude of parameters, the best approach consists in using PS-NM or GA-NM as they enable thorough explorations of possible solutions.

* When dealing with a relatively simple optimization problems or when it is possible to impose stricter boundary limits, i.e. for common MEA compositions and PEMFC configurations for which the parameter orders of magnitude are known, employing SA-NM may result in lower computational time.

Finally, it's worth pointing out that, although the use of hybrid algorithms leads to a better estimation of impedance parameters from a numerical quality point of view, the question of the physical representativeness of equivalent electrical circuits has not been addressed in this
work. This aspect remains an important challenge to be tackled in the future, for example by using other characterization methods -such as Distribution of Relaxation Times (DRT)- to determine the main physical phenomena at the origin of the electrical response of the devices analyzed.

### CRediT authorship contribution statement

P. Wu: Writing – original draft, Visualization, Investigation, Formal analysis, Data curation. S. Touhami: Investigation, Data curation. W. Aït-Idir: Data curation. H. Yelloz: Investigation, Data curation. C. Marty: Data curation. F. Micoud: Project administration, Investigation, Formal analysis, Conceptualization. J. Dillet: Investigation, Data curation. O. Lottin: Writing – review & editing, Validation, Supervision, Project administration, Investigation, Formal analysis, Conceptualization. J. Mainka: Writing – review & editing, Visualization, Validation, Supervision, Project administration, Investigation, Formal analysis, Conceptualization.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Acknowledgments

This work has received funding from the European Union's EIT Raw Materials project n° 19247, ALPE: "Advanced Low-Platinum hierarchical Electrocatalysts for low-T fuel cells", as well as from the region Grand Est, the European Fund for Economic and Regional Development FEDER, the French Environmental and Energy Management Agency ADEME and the French National Research Agency ANR-project LOCALI (ANR-17-CE05-0003).

The authors also thank Kévin Mozet from Métro'NRJ – LEMTA (Université de Lorraine-CNRS – https://lemta.univ-lorraine.fr/).

### Data availability

Published on Recherche Data Gouv: DOI: https://doi.org/10.57745/VRFDKM.

### References

[1] S. Touhami, J. Mainka, J. Dillet, S.A.H. Taleb, O. Lottin, Transmission Line Impedance Models Considering Oxygen Transport Limitations in Polymer Electrolyte Membrane Fuel Cells, J. Electrochem. Soc. 166 (15) (Nov. 2019) F1209, https://doi.org/10.1149/2.0891915jes.

[2] J. Mainka, Local impedance in H2/air Proton Exchange Membrane Fuel Cells (PEMFC): theoretical and experimental investigations, HAL CCSD (2011).

[3] T.E. Springer, T.A. Zawodzinski, M.S. Wilson, S. Gottesfeld, Characterization of Polymer Electrolyte Fuel Cells Using AC Impedance Spectroscopy, J. Electrochem. Soc. 143 (2) (Feb. 1996) 587, https://doi.org/10.1149/1.1836485.

[4] J.E.B. Randles, K.W. Somerton, Kinetics of rapid electrode reactions. Part 3.— Electron exchange reactions, Trans. Faraday Soc. 48 (0) (Jan. 1952) 937–950, https://doi.org/10.1039/TF9524800937.

[5] P.M. Gomadam, J.W. Weidner, Analysis of electrochemical impedance spectroscopy in proton exchange membrane fuel cells, Int. J. Energy Res. 29 (12) (2005) 1133–1151, https://doi.org/10.1002/er.1144.

[6] M.A. Rubio, A. Urquia, R. Kuhn, S. Dormido, Electrochemical parameter estimation in operating proton exchange membrane fuel cells, J. Power. Sources. 183 (1) (Aug. 2008) 118–125, https://doi.org/10.1016/j.jpowsour.2008.05.011.

[7] E. Warburg, Ueber das Verhalten sogenannter unpolarisirbarer Elektroden gegen Wechselstrom, Ann. Phys. 303 (3) (1899) 493–499, https://doi.org/10.1002/andp.18993030302.

[8] B.A. Boukamp, A Nonlinear Least Squares Fit procedure for analysis of immittance data of electrochemical systems, Solid. State Ion. 20 (1) (Feb. 1986) 31–44, https://doi.org/10.1016/0167-2738(86)90031-7.

[9] J.R. Macdonald, J. Schoonman, A.P. Lehnen, Applicability and power of complex nonlinear least squares for the analysis of impedance and admittance data, J. Electroanal. Chem. Interfacial. Electrochem. 131 (Jan. 1982) 77–95, https://doi.org/10.1016/0022-0728(82)87062-9.

[10] L. Wang, et al., A hybrid Genetic Algorithm and Levenberg-Marquardt (GA-LM) method for cell suspension measurement with electrical impedance spectroscopy, Review of scientific instruments 91 (12) (2020) 124104, https://doi.org/10.1063/5.0029491, 124104.

[11] J. Zhou, H.-J. Kang, A hybrid least-squares genetic algorithm-based algorithm for simultaneous identification of geometric and compliance errors in industrial robots, Advances in mechanical engineering 7 (6) (2015) 168781401559028, https://doi.org/10.1177/1687814015590289.

[12] T. Zhou, Y. Zhou, K. Sundmacher, A hybrid stochastic-deterministic optimization approach for integrated solvent and process design, Chem. Eng. Sci. 159 (Feb. 2017) 207–216, https://doi.org/10.1016/j.ces.2016.03.011.

[13] V. Pizarro-Carmona, S. Castano-Solís, M. Cortés-Carmona, J. Fraire-Ardanuy, D. Jimenez-Bermejo, GA-based approach to optimize an equivalent electric circuit model of a Li-ion battery-pack, Expert. Syst. Appl. 172 (Jun. 2021) 114647, https://doi.org/10.1016/j.eswa.2021.114647.

[14] C.G. Moles, P. Mendes, J.R. Banga, Parameter Estimation in Biochemical Pathways: a Comparison of Global Optimization Methods, Genome Res. 13 (11) (Nov. 2003) 2467–2474, https://doi.org/10.1101/gr.1262503.

[15] Q. Yang, J. Xu, B. Cao, X. Li, A simplified fractional order impedance model and parameter identification method for lithium-ion batteries, PLoS. One 12 (2) (2017) e0172424, https://doi.org/10.1371/journal.pone.0172424 e0172424.

[16] J. Mainka, G. Maranzana, J. Dillet, S. Didierjean, O. Lottin, On the estimation of high frequency parameters of Proton Exchange Membrane Fuel Cells via Electrochemical Impedance Spectroscopy, J. Power. Sources. 253 (May 2014) 381–391, https://doi.org/10.1016/j.jpowsour.2013.12.064.

[17] B.T. Habte, F. Jiang, Effect of microstructure morphology on Li-ion battery graphite anode performance: electrochemical impedance spectroscopy modeling and analysis, Solid. State Ion. 314 (2018) 81–91, https://doi.org/10.1016/j.ssi.2017.11.024.

[18] E.S. Castillo-Romero, R. Jaimes-López, M. Miranda-Hernández, Estimation of the electrochemical kinetic parameters of redox couples, on activated porous electrodes, based on impedance spectroscopy and particle swarm optimization, J. Energy Storage 55 (2022) 105574, https://doi.org/10.1016/j.est.2022.105574.

[19] B. Wang, S.E. Li, H. Peng, Z. Liu, Fractional-order modeling and parameter identification for lithium-ion batteries, J. Power. Sources. 293 (2015) 151–161, https://doi.org/10.1016/j.jpowsour.2015.05.059.

[20] S. Sharifi-Asl, M.L. Taylor, Z. Lu, G.R. Engelhardt, B. Kursten, D.D. Macdonald, Modeling of the electrochemical impedance spectroscopic behavior of passive iron using a genetic algorithm approach, Electrochim. Acta 102 (Jul. 2013) 161–173, https://doi.org/10.1016/j.electacta.2013.03.143.

[21] A.P. Bhondekar, R. Kaur, R. Kumar, R. Vig, P. Kapur, A novel approach using Dynamic Social Impact Theory for optimization of impedance-Tongue (iTongue), Chemometrics and intelligent laboratory systems 109 (1) (2011) 65–76, https://doi.org/10.1016/j.chemolab.2011.08.002.

[22] P. Moscato, "On Evolution, Search, Optimization, Genetic Algorithms and Martial Arts Towards Memetic Algorithms".

[23] Y. Zhang, W. Jia, X. Fu, D. Li, C. Yu, [A Multi-Peak Brillouin Scattering Spectrum Feature Extraction Method Based on Multi-Criteria Decision-Making and Particle Swarm Optimization-Levenberg Marquardt Hybrid Optimization Algorithm], Guang. Pu Xue Yu Guang. Pu Fen. Xi. 36 (7) (Jul. 2016) 2183–2188.

[24] A. Kabouche, A. Boultif, A. Abidi, N. Gherraf, Interaction parameter estimation in liquid–liquid phase equilibrium modeling using stochastic and hybrid algorithms, Fluid. Phase Equilib. 336 (Dec. 2012) 113–121, https://doi.org/10.1016/j.fluid.2012.09.002.

[25] M. Kowski, M.J. Palys, Genetic Algorithm-Based Improvement to Equivalent Circuit Analysis of Experimental Data in Electrochemical Impedance Spectroscopy, ChemElectroChem. 8 (15) (2021) 2956–2967, https://doi.org/10.1002/cec2.202100778.

[26] D. Yu, X. Li, S.S. Araya, S.L. Sahlin, V. Liso, Data-driven diagnosis of high temperature PEM fuel cells based on the electrochemical impedance spectroscopy: robustness improvement and evaluation, J. Energy Chem. 96 (Sep. 2024) 544–558, https://doi.org/10.1016/j.jechem.2024.05.014.

[27] J. Lv, Z. Yu, G. Sun, J. Liu, Deep learning-based fault diagnosis and Electrochemical Impedance Spectroscopy frequency selection method for Proton Exchange Membrane Fuel Cell, J. Power. Sources. 591 (Jan. 2024) 233815, https://doi.org/10.1016/j.jpowsour.2023.233815.

[28] J. Mainka, W. Gao, N. He, J. Dillet, O. Lottin, A General Equivalent Electrical Circuit Model for the characterization of MXene/graphene oxide hybrid-fiber supercapacitors by electrochemical impedance spectroscopy – Impact of fiber length, Electrochim. Acta 404 (Feb. 2022) 139740, https://doi.org/10.1016/j.electacta.2021.139740.

[29] S. Touhami, L. Dubau, J. Mainka, J. Dillet, M. Chatenet, O. Lottin, Anode aging in polymer electrolyte membrane fuel Cells I: anode monitoring by ElectroChemical impedance spectroscopy, J. Power. Sources. 481 (2021) 228908, https://doi.org/10.1016/j.jpowsour.2020.228908.

[30] M. Daoudi, et al., Impact of sulfonated poly(ether ether ketone) membranes pretreatments on their physicochemical properties and fuel cell performances, Journal of Polymer Science (Jan. 2023), https://doi.org/10.1002/pol.20220649 pol.20220649.

[31] W. Aït-Idir, et al., Oxygen diffusion impedance in proton exchange membrane fuel cells – insights into electrochemical impedance spectra and equivalent electrical circuit modeling, Electrochim. Acta 472 (Dec. 2023) 143430, https://doi.org/10.1016/j.electacta.2023.143430.

[32] M. Daoudi et al., "Impact of sulfonated poly(ether ether ketone) membranes pretreatments on their physicochemical properties and fuel cell performances," *Journal of Polymer Science*, vol. n/a, no. n/a, doi: 10.1002/pol.20220649.
[33] R. Sgarbi, et al., Does the platinum-loading in proton-exchange membrane fuel cell cathodes influence the durability of the membrane-electrode assembly? Ind. Chem. Mater. 1 (4) (2023) 501–515, https://doi.org/10.1039/D3IM00059A.

[34] J.C. Lagarias, J.A. Reeds, M.H. Wright, P.E. Wright, Convergence Properties of the Nelder-Mead Simplex Method in Low Dimensions, SIAM journal on optimization 9 (1) (1998) 112–147, https://doi.org/10.1137/S1052623496303470.

[35] C. Reeves, Genetic Algorithms, in: Handbook of Metaheuristics, 146, 2010, pp. 109–139, https://doi.org/10.1007/978-1-4419-1665-5_5.

[36] Goldberg David Edward, Genetic Algorithms in search, optimization, and Machine Learning /David E. Goldberg, Addison-Wesley publishing company, 1989. Reading (Mass.) [etc.].

[37] A.R. CONN, N.I.M. GOULD, P.L. TOINT, A globally convergent augmented Lagrangian algorithm for optimization with general constraints and simple bounds, SIAM. J. Numer. Anal. 28 (2) (1991) 545–572, https://doi.org/10.1137/0728030.

[38] L. Momayez, P. Dupont, G. Delacourt, O. Lottin, H. Peerhossaini, Genetic algorithm based correlations for heat transfer calculation on concave surfaces, Appl. Therm. Eng. 29 (17) (Dec. 2009) 3476–3481, https://doi.org/10.1016/j.applthermaleng.2009.05.025.

[39] J. Kennedy, R. Eberhart, Particle swarm optimization, in: Proceedings of ICNN'95-International Conference on Neural Networks, IEEE, Perth, WA, Australia, 1995, pp. 1942–1948, https://doi.org/10.1109/ICNN.1995.488968.

[40] E. Mezura-Montes, C.A. Coello Coello, Constraint-handling in nature-inspired numerical optimization: past, present and future, Swarm. Evol. Comput. 1 (4) (Dec. 2011) 173–194, https://doi.org/10.1016/j.swevo.2011.10.001.

[41] M. Pedersen, Good Parameters for Particle Swarm Optimization, 2010, https://doi.org/10.13140/RG.2.2.23771.49445.

[42] S. Kirkpatrick, J. Gelatt, M.P. Vecchi, Optimization by Simulated Annealing, Science (American Association for the Advancement of Science) 220 (4598) (1983) 671–680, https://doi.org/10.1126/science.220.4598.671.

[43] É.D. Taillard, Design of Heuristic Algorithms For Hard Optimization : With Python Codes for the Travelling Salesman Problem, Springer Nature, 2023. Accessed: Jun. 06, 2023[Online]. Available: https://directory.doabooks.org/handle/20.500.12854/93984.

[44] A. Kerr, K. Mullen, A comparison of genetic algorithms and simulated annealing in maximizing the thermal conductance of harmonic lattices, Comput. Mater. Sci. 157 (2019) 31–36, https://doi.org/10.1016/j.commatsci.2018.10.007. Feb.

[45] F. Jia, D. Lichti, A COMPARISON OF SIMULATED ANNEALING, GENETIC ALGORITHM AND PARTICLE SWARM OPTIMIZATION IN OPTIMAL FIRST-ORDER DESIGN OF INDOOR TLS NETWORKS, ISPRS Ann. Photogramm. Remote Sens. Spatial Inf. Sci. (Sep. 2017) 75–82, https://doi.org/10.5194/isprs-annals-IV-2-W4-75-2017. IV-2/W4.

[46] M. Ghazanfari, S. Alizadeh, M. Fathian, D.E. Koulouriotis, Comparing simulated annealing and genetic algorithm in learning FCM, Appl. Math. Comput. 192 (1) (Sep. 2007) 56–68, https://doi.org/10.1016/j.amc.2007.02.144.

[47] J.J. Moré, B.S. Garbow, K.E. Hillstrom, Testing Unconstrained Optimization Software, ACM Trans. Math. Softw. 7 (1) (Mar. 1981) 17–41, https://doi.org/10.1145/355934.355936.

[48] V. Beiranvand, W. Hare, Y. Lucet, Best Practices for Comparing Optimization Algorithms, Optim. Eng. 18 (Dec. 2017), https://doi.org/10.1007/s11081-017-9366-1.

[49] R. Hassan, B. Cohanim, and O. de Weck, "A Comparison of Particle Swarm Optimization and the Genetic Algorithm," Apr. 2005. doi: 10.2514/6.2005-1897.

[50] "Virtual Library of Simulation Experiments: test Functions and Datasets." Accessed: Nov. 21, 2023. [Online]. Available: http://www.sfu.ca/~ssurjano/.

[51] A.O. Griewank, Generalized descent for global optimization, J. Optim. Theory. Appl. 34 (1) (May 1981) 11–39, https://doi.org/10.1007/BF00933356.

[52] J.A. Bilmes, "A Gentle Tutorial of the EM Algorithm and its Application to Parameter Estimation for Gaussian Mixture and Hidden Markov Models".

[53] J.S. Rustagi, Optimization Techniques in Statistics, Academic Press, 2014. Accessed: Jan. 04, 2024. [Online]. Available: https://www.perlego.com/book/1896141/optimization-techniques-in-statistics-pdf.

[54] J.S. Kowalik, Practical Optimization, SIAM Review 25 (2) (1983) 283–284, https://doi.org/10.1137/1025065.

[55] A. Törn and A. Zilinskas, Eds., "Global Optimization," vol. 350, 1989, doi: 10.1007/978-1-540-50871-6.

[56] P. Wu, S. Touhami, W. Aït-Idir, J. Dillet, O. Lottin, J. Mainka, Hybrid stochastic-deterministic algorithms for the interpretation of Electrochemical Impedance Spectroscopy spectra of Proton Exchange Membrane Fuel Cells: data for replication, Recherche Data Gouv (2024), https://doi.org/10.57745/VRFDKM. May 06.

[57] D.H. Ackley, A connectionist machine for genetic hillclimbing. The Kluwer International Series in Engineering and Computer Science. Boston, [Massachusetts], Kluwer Academic Publishers, 1987, https://doi.org/10.1007/978-1-4613-1997-9.

[58] L.C.W. Dixon, R.C. Price, Truncated Newton method for sparse unconstrained optimization using automatic differentiation, J. Optim. Theory. Appl. 60 (2) (Feb. 1989) 261–275, https://doi.org/10.1007/BF00940007.

[59] *Michalewicz Zbigniew, Genetic Algorithms + Data Structures = Evolution Programs /ZbigniewMichalewicz, 3rd Revised and Extended edition*, Springer, Berlin Heidelberg, 1996.

[60] H. Mühlenbein, M. Schomisch, J. Born, The parallel genetic algorithm as function optimizer, Parallel. Comput. 17 (6) (Sep. 1991) 619–632, https://doi.org/10.1016/S0167-8191(05)80052-3.

[61] M. Jamil, X.-S. Yang, A Literature Survey of Benchmark Functions For Global Optimization Problems, arXiv.org (2013), https://doi.org/10.48550/arxiv.1308.4008.

[62] E.P. Adorio, "MVF-Multivariate Test Functions Library in C for Unconstrained Global Optimization".

[63] A. Fischbach, M. Zaefferer, J. Stork, M. Echtenbruck, and T. Bartz-Beielstein, "From Real World Data to Test Functions," Jan. 2016, p. 159.

[64] M. Zaefferer, A. Fischbach, B. Naujoks, T. Bartz-Beielstein, Simulation-based test functions for optimization algorithms, in: Proceedings of the Genetic and Evolutionary Computation Conference, Berlin Germany, ACM, Jul. 2017, pp. 905–912, https://doi.org/10.1145/3071178.3071190.

[65] A.Z. Weber, A. Kusoglu, Unexplained transport resistances for low-loaded fuel-cell catalyst layers, J. Mater. Chem. A 2 (41) (2014) 17207–17211, https://doi.org/10.1039/C4TA02952F.

[66] S. Cruz-Manzo, P. Greenwood, Study of oxygen diffusion in the cathode catalyst layer and gas diffusion layer for polymer electrolyte fuel cells with EIS, Journal of Electroanalytical Chemistry 892 (Jul. 2021) 115270, https://doi.org/10.1016/j.jelechem.2021.115270.

[67] A. Kulikovsky, Why impedance of the gas diffusion layer in a PEM fuel cell differs from the Warburg finite-length impedance? Electrochem. commun. 84 (2017) 28–31, https://doi.org/10.1016/j.elecom.2017.09.014.

[68] R. de Levie, On porous electrodes in electrolyte solutions—IV, Electrochim. Acta 9 (9) (Sep. 1964) 1231–1245, https://doi.org/10.1016/0013-4686(64)85015-5.

[69] R. de Levie, On porous electrodes in electrolyte solutions: I. Capacitance effects, Electrochim. Acta 8 (10) (Oct. 1963) 751–780, https://doi.org/10.1016/0013-4686(63)80042-0.

[70] M. Ciureanu, H. Wang, Electrochemical Impedance Study of Electrode-Membrane Assemblies in PEM Fuel Cells: I. Electro-oxidation of H 2 and H 2 /CO Mixtures on Pt-Based Gas-Diffusion Electrodes, J. Electrochem. Soc. 146 (11) (Nov. 1999) 4031, https://doi.org/10.1149/1.1392588.

[71] A.R. Thompson, A.V. Call, D.J. Cumming, D.C. Sinclair, R.H. Elder, Electrochemical Impedance Spectroscopy Data from Solid Oxide Cells Undergoing Co-Electrolysis: the Influence of Rig Inductance, ECS. Trans. 68 (1) (Jun. 2015) 3417, https://doi.org/10.1149/06801.3417ecst.