A typical PEMFC includes two electrodes (cathode and anode) and a proton exchange membrane (PEM), while electrodes are mainly composed of a gas diffusion layer and a platinum-based alloy catalyst layer. Hydrogen is decomposed into hydrogen ions H⁺ and electrons e⁻ through catalytic reaction after passing through the anode catalytic layer. Subsequently, oxygen provided by the cathode catalytic layer integrates with electrons and hydrogen ions to generate water and heat. Figure 1 describes the structure of the reaction mechanism diagram, a single PEMFC, and a PEMFC stack of PEMFC. The electrochemical reaction mechanism in PEMFC is expressed as follows:

At anode, H 2 → 2 H + + 2 e −   . (1)

At cathode, 2 H + + 2 e − + 1 2 O 2 → H 2 O . (2)

Overall chemical reaction, 2 H 2 + O 2 → 2 H 2 O + heat   + Electricity . (3)

The net output voltage of the electrochemical model is affected by three kinds of polarization, namely, activation, ohmic, and concentration polarization, while the voltage characteristic function of electrochemical can be written as (Amphlett et al., 1995) V c = E nernst − V act − V ohmic − V con , (4) where   E nernst   indicates the open-circuit voltage (V);   V act   means the activation voltage drop (V);   V ohmic   denotes the ohmic voltage drop (V); and   V con   represents the concentration voltage drop (V).

In addition,   E nernst   calculated from the Nernst equation can be determined by (Ariza et al., 2018) E nernst = Δ G 2 F + Δ S 2 F ( T k − T ref ) + R T 2 F [ ln ( P H 2 ) + 1 2 ( P O 2 ) ] , (5) where Δ G is the increment in Gibbs free energy (J/mol); F indicates Faraday’s constant (96,487 C/mol); Δ S means the increment of entropy (J/mol); R denotes the universal gas constant [8.314 J/(K mol)];   T k   and  T ref   represent the operation ambient and reference temperature (K), respectively; P H 2   and   P O 2   stand for the partial pressures of hydrogen and oxygen, respectively (atm), which are formulated as (Ali et al., 2017) P H 2 = 0.5 × R H a × P H 2 O sat × [ ( R H a × P H 2 O sat P a × exp ( 1.635 ( i cell A ) T k 1.334 ) ) − 1 − 1 ]   , (6) P O 2 = R H c × P H 2 O sat × [ ( R H a × P H 2 O sat P c × exp ( 4.192 ( i cell A ) T k 1.334 ) ) − 1 − 1 ]   , (7) where RH _a and RH _c, respectively, denote the relative humidity of vapor at anode and cathode (atm); P _a and P _c are the inlet pressures of anode and cathode (atm), respectively; and   P H 2 O sat   represents the saturation pressure of water vapor (atm), which is defined as (Ali et al., 2017) T c = T k − 273.15 , (8) log 10 ( P H 2 O sat ) = 2.95 × 10 − 2 × T c − 9.19 × 10 − 5 × T c 2 + 1.44 × 10 − 7 × T c 3 − 2.18 . (9)

Moreover, the activation voltage   V act   demonstrates the slowness of the reaction occurring on the electrode surface, which can be given as follows (Ariza et al., 2018): V act = ε 1 + ε 2 T k + ε 3 T k l n ( C o 2 ) + ε 4 T k l n ( i cell ) , (10) where   ε i   (i = 1,2,3,4) indicates the semi-empirical coefficients and   C o 2   is the oxygen concentration (mol/ cm 3 ), which is determined by (Ali et al., 2017) C O 2 = P O 2 5.08 × 10 6 × e ( − 498 T k ) . (11)

Ohmic voltage drop   V ohmic is mathematically expressed as follows (Ali et al., 2017): V ohmic = i cell ( R m + R c ) , (12) where R _c is the equivalent contact resistance (Ω) and R _m represents the equivalent resistance provided to PEM conduction (Ω), which can be defined by R m = ρ m ( l A ) , (13) where l means the membrane thickness ( μ m ), A is the effective electrode area (cm²), and   ρ m   stands for the specific membrane resistance ( Ω Δ cm) described by (Amphlett et al., 1995) ρ m = 181.6 × [ 1 + 0.03 × ( i cell A ) + 0.062 × ( T k 303 ) 2 ( i cell A ) 2.5 ] [ λ − 0.634 − 3 × ( i cell A ) ] exp [ 4.18 × ( T k − 303 T k ) ] , (14) where   λ   stands for an empirical parameter.

Furthermore, concentration voltage loss   V con   influenced by the concentration can be formulated as (Ariza et al., 2018) V con = − b ln ( ln J A × J max ) , (15) where b means the parametric coefficient (V), J   is the actual current density, and J max   defines the upper bound current density ( A / c m 2 ).

On the whole, after referring to the mentioned Eqs 4–15 for PEMFC, there are seven crucial parameters demand to be identified: ε 1   ,   ε 2   ,   ε 3   ,   ε 4   ,   b ,   λ , and   R c .

To precisely build a mathematical model of PEMFC based on the aforementioned unknown parameters, it is crucial to utilize an objective function to reliably evaluate parameter identification. In addition, the objection function based on the root mean square error (RMSE) can commendably mirror the deviation between the actual and estimated values.

Consequently, RMSE is accounted as the objective function given by R M S E   ( x ) = 1 N ∑ i = 1 N [ V actual ( i ) − V estimate ( i ) ] 2 , (16) where N indicates the total quantity of actual datasets, V _actual is the actual voltage, and V _estimate denotes the estimated voltage.

Additionally, the restraints of crucial parameters can be written as s . t . { ε i , min ≤ ε i ≤ ε i , max λ min ≤ λ ≤ λ max R c , min ≤ R c ≤ R c , max b min ≤ b ≤ b max         , ∀ i ∈ { 1,2,3,4 } . (17)

The BES algorithm (Atlam and Dndar, 2021) is an innovative meta-heuristic algorithm enlightened by biological activities in nature. Individual search strategies can be established through their own mobile or global experience in the group random optimization process to optimize complex problems in the real world. In particular, the self-search and exploration of individuals in the search space are performed by simulating the predation behavior of bald eagles (e.g., prey on salmon living in a specific area). The bald eagles have a specific predation strategy, which can help them consume the least energy while maximizing the probability of predation success.

The bald eagle’s predation strategy has three major aspects: selecting the appropriate search domain, searching in the selected space, and the best chance to swoop on the prey.

In the select stage, the selected search space is related to the last movement of the bald eagle, which can be expressed as P new , i = P best + α × r ( P mean + P i ) , (18) where P best means optimal spatial search position, α is parameters affecting location update, 1.5 ≤ α ≤ 2, r denotes a random number between 0 and 1, P mean indicates average distribution position of the bald eagle after the previous search.

After determining the search space, the bald eagle will search the space that flies spirally to search for prey and find the best-accelerated dive position. The change of the search position is expressed in polar coordinates as follows (Salim et al., 2015): P new , i = P i + y ( i ) × ( P i − P i + 1 ) + x ( i ) × ( P i − P mean ) , (19) x ( i ) = x r ( i ) max ( | x r | ) ,                   y ( i ) = y r ( i ) max ( | y r | ) , (20) x r ( i ) = r ( i ) × sin [ θ ( i ) ] ,                   y r ( i ) = r ( i ) × cos [ θ ( i ) ] , (21) θ ( i ) = a × π × r a n d ,                   r ( i ) = θ ( i ) + R × r a n d , (22) where x ( i ) and y ( i ) are values between 0 and 1 for determining the positions of the bald eagle on the corresponding polar axis, respectively; r ( i ) and θ ( i ) denote the polar diameter and polar angle of the spiral flight in polar coordinates, respectively; a indicates a parameter that takes a value from 5 to 10; R represents a parameter that controls the search period, between 0.2 and 2; and rand is a random number in [0,1].

In the swooping stage, the bald eagle swoops to the prey at the optimal position, while other bald eagles in the population also move to the best position and dive down to attack the prey. The position updates during the swooping are expressed in polar coordinates as (Oliva et al., 2014) P new , i = r a n d × P best + x 1 ( i ) × ( P i − c 1 × P mean ) + y 1 ( i ) × ( P i − c 2 × P best ) , (23) x 1 ( i ) = x r ( i ) max ( | x r | ) ,           y 1 ( i ) = y r ( i ) max ( | y r | ) , (24) x r ( i ) = r ( i ) × sinh [ θ ( i ) ] ,            y r ( i ) = r ( i ) × cos h [ θ ( i ) ] , (25) θ ( i ) = a × π × r a n d ,           r ( i ) = θ ( i ) , (26) where c 1 and c 2 , respectively, stand for moving velocities toward the optimal and central positions and are both between 1 and 2.

On the basis of the BES algorithm, Figure 2 illustrates the flowchart of PEMFC overall parameter identification. Firstly, it is vital to determine which parameters in the steady-state electrochemical model are identified. After that, the output voltage and current data selected as the actual PEMFC will be regarded as the inputs of BES, while the data are transformed into the objective function according to Eq. 16. Furthermore, BES is employed to identify the parameters in accordance with the built model. Finally, the optimal parameter identification results of PEMFC are output after multiple iterations.

When the experimental condition is set as T _k = 333.15K, RH _a = 75%, and RH _c = 75%, the simulation values of seven unknown parameters are generated by BES and GA algorithms. The satisfactory results of parameter identification and minimum RMSE of the PEMFC model are illustrated in Table 2. Particularly, RMSE acquired by BES is 94.11% lower than those acquired by GA, respectively. Hence, the advisable performance of BES is significantly greater than that of GA, which can be attributed to its consideration of accuracy, reliability, and high efficiency.

What is more, Figure 3A depicts the average RMSE acquired by the two algorithms under 10 independent runs. It is transparent that the accuracy of BES is significantly higher than that of GA, which fully demonstrates BES’s superior performance. Given that the average RMSEs of BES and GA algorithms hardly consist in the same order of magnitude, and the value of GA is about 12 times higher than that of BES, the BES’s effect of the actual data approximation is much better than GA.

Figure 3B shows a boxplot of BES and GA, demonstrating the distribution range and upper/lower bounds of simulation explicitly and comprehensively. It can be seen from the chart that the error fluctuation interval and average RMSE of BES are far less than that of GA. Thus, the BES algorithm has accurate searching ability in PEMFC parameter identification and significant global searching ability.

Meanwhile, convergence graphs of the two algorithms are depicted in Figure 3C, while BES attains a stable optimal solution rapidly in approximately 60 iterations based on a global search.

As a single individual optimization and group cooperation mechanism, BES performs excellent local exploitation and global exploration, by which the accuracy and efficiency of parameter identification can be drastically improved.

Figure 3D describes the output V-I fitting curve on the basis of global optimal parameters identification by BES, upon which one can readily discover that the model curve achieved from BES approximates the experimental data to a great extent. It is undeniable that BES has superb performance in PEMFC parameter identification due to superior optimization ability and effectiveness.

When the simulation is in the condition of T _k = 313.15K, RH _a = 75%, and RH _c = 75%, the statistical results of seven unknown parameters and RMSE based on BES and GA are shown in Table 3. In particular, RMSE obtained by BES is significantly smaller than GA, while BES is 96.27% smaller than the RMSE of GA. Upon them, it is obvious that BES can considerably improve the accuracy of PEMFC model parameter identification.

Figure 4A presents the average RMSE obtained by BES and GA, illustrating that BES can acquire a lower average RMSE under 10 independent runs. Thus, BES has better accuracy and stability in the unknown parameter identification of PEMFC.

In addition, the RMSE boxplot obtained by BES and GA is depicted in Figure 4B, upon which the distribution range and upper/lower bound of BES are lower than those of GA. One can easily observe that the balance between global exploration and local exploitation of BES is better than that of GA.

In the meantime, convergence curves of BES and GA are shown in Figure 4C, where BES has about 55 iterations to achieve convergence, while GA needs 200 iterations to achieve convergence stability. Besides, RMSE after BES convergence is smaller than that of GA, while it can be perceived that BES identifies the unknown parameters of the PEMFC model with more accuracy and efficiency.

Figure 4D demonstrates the output V-I fitting curve obtained by the optimal results of BES, where data via BES are highly fitted with experiment data. The result efficiently reflects that BES has a superior ability for PEMFC parameter identification.

Table 4 illustrates the best PEMFC electrochemical model parameters and minimum RMSE under the condition of T _k = 353.15K, RH _a = 50%, and RH _c = 50% through BES and GA algorithm. In light of the result, BES represents the optimal performance owing that the minimum RMSE value acquired by BES is 46.13% lower than GA, which can accomplish parameter identification tasks in PEMFC with higher quality and accuracy.

The average RMSE obtained by BES and GA under 10 independent runs is shown in Figure 5A, while the value of the former is about one-quarter of the latter in the case of multiple operations, which is a lot smaller. Hence, BES outperforms GA in parameter identification and shows good stability and global optimization ability.

As is depicted in Figure 5B, it can be observed from the boxplot of BES and GA that the fluctuation range and average value of RMSE optimized by BES are lower than those of GA with 10 independent runs, which, from the perspective of either its stability or accuracy, effectively confirms the outstanding performance of BES in identifying the accuracy of PEMFC parameters. Moreover, Figure 5C presents the convergence curve of the two algorithms. BES converges at about 90 iterations, while GA does not reach a stable value until 480 iterations. After iterating 270 to 480 times, GA falls into a locally optimal solution and results in low efficiency and reduced accuracy. In general, BES is superior in convergence stability and precision for its exceptional global optimization ability.

In the end, Figure 5D describes convergence graphs of BES algorithms, which indicates that, in the case of parameter identification of BES, the experimental data and model fitting data are highly coincident, while the error is extremely small due to superior accuracy.