Optimization Algorithms: Optimal Parameters Computation for Modeling the Polarization Curves of a PEFC Considering the Effect of the Relative Humidity

Abstract

The development of green energy conversion devices has been promising to face climate change and global warming challenges over the last few years. Energy applications require a confident performance prediction, especially in polymer electrolyte fuel cell (PEFC), to guarantee optimal operation. Several researchers have employed optimization algorithms (OAs) to identify operating parameters to improve the PEFC performance. In the current study, several nature-based OAs have been performed to compute the optimal parameters used to describe the polarization curves in a PEFC. Different relative humidity (RH) values, one of the most influential variables on PEFC performance, have been considered. To develop this study, experimental data have been collected from a lab-scale fuel cell test system establishing different RH percentages, from 18 to 100%. OAs like neural network algorithm (NNA), improved grey-wolf optimizer (I-GWO), ant lion optimizer (ALO), bird swarm algorithm (BSA), and multi-verse optimization (MVO) were evaluated and compared using statistical parameters as training error and time. Results gave enough information to conclude that NNA had better performance and showed better results over other highlighted OAs. Finally, it was found that sparsity and noise are more present at lower relative humidity values. At low RH, a PEFC operates under critical conditions, affecting the fitting on OAs.

1. Introduction

In the last few years, several studies related to renewable energy technologies have been performed [1,2,3,4]. The development of new energy sources with high efficiency and environmentally friendliness is one of the challenges of today’s world. The fossil fuels currently used are unsustainable and have considerable impacts on the environment. The mentioned fossil fuels are the primary source to produce the consumed energy, causing severe environmental problems like climate change and pollution. In 2005, they were estimated to produce 80% of the world’s primary energy [5]. In the same year, 16.5% of the world’s primary energy was produced by renewable energy and expected to replace 80% of fossil fuels up to 2050, i.e., in about four decades [6]. This aim is in line with the European Green Deal proposed in 2011 [7].

Polymer electrolyte fuel cells (PEFCs) are presented as one of the applications of hydrogen technology that can be used to meet energy requirements. PEFCs are electrochemical devices that transform the energy from chemical reactions into electrical energy [8]. The energy conversion is done in a single process compared to combustion engines engines, in which the energy conversion is performed in multiple steps [9]. Furthermore, the energy conversion process in a PEFC is based on three main stages: hydrogen as the fuel flow at the anode while oxygen as oxidant flows through the cathode, the gas molecules are split at the active region of the membrane electrode assembly (MEA), and, finally, free electrons are transported out to the cell producing clean electrical energy. Apart from this, it can categorized in different types [10,11]. PEFCs are one of the most promising energy conversion devices to produce electricity with zero emission. In addition, it can work at low pressure, low temperature, noiseless manner, and at different current densities with high efficiency [8].

To study the behavior of a PEFC, polarization curves are typically used. The different parameters that affect their performance like temperature, relative humidity (RH), pressure, oxidant/fuel flow, and dry conditions are usually varied [12,13,14]. The polarization curves present the voltage as a function of the current density. In the mentioned curves, three regions are clearly distinguished: activation, ohmic, and concentration region. The activation region represents the energy required to start reactions at low current densities. On the other hand, the ohmic region occurs at moderate current densities and represents the ohmic losses inside the PEFC. Finally, the region of concentration losses occurs due to the abundant generation of liquid water blocking the pores of the diffusion media (DM) at high-density currents, and this impairs the development of electrochemical reactions [15].

To accomplish the optimum operation of PEFCs, several studies have been dealing with the problem of accurately identifying its internal electrochemical and mechanical parameters by modeling and simulation [16,17]. The modeling procedures may be classified as mechanical or experimental modeling. Mechanical modeling deals with the electrochemical, heat, and mass transport parameters. The experimental modeling is performed partially from the mechanical one and formulated as empirical equations [18,19,20,21,22].

The models of PEFCs should be viable, computationally effective, and capable of predicting their performances over a wide range of operating environmental conditions [23]. Appropriate model identifications can be ensured by entering correct input parameters in the governing equations of the PEFC, which comprise chemical and physical properties. The modeling methodologies of PEFCs can either be theoretical, semi-empirical, or empirical [24,25,26].

Over the last few decades, there has been a variety of studies related to PEFCs parameters’ identifications using heuristic-based methods, such as particle swarm optimization (PSO) [27], which is used to define the parameters of steady-state PEFC models. In addition, an improved grey-wolf optimizer algorithm (I-GWO) and multiverse optimization (MVO) algorithm have been applied to accomplish the same goal. The ant–lion optimizer (ALO) algorithm has also been implemented to model PEFCs [28].

Correctly identifying the unspecified parameters leads to the actual model of a PEFC. Previous studies related to determining the optimal parameters of FCs are described as follows. For example, Priya et al. [29] utilized a simple genetic algorithm to obtain the accurate model of a PEFC by identifying the optimal parameters and used the derivative of power with respect to current as the objective function (OF). Gong et al. [30] applied a differential evolution method to obtain the optimal parameters of a PEFC. This algorithm can also include a ranking-based operator to accelerate the convergence of the solutions. On the other hand, Ye et al. [31] used PSO to identify PEFC parameters to minimize its OF, which was based on the calculated and measured stack voltage. Likewise, Fathy et al. [32] introduced an MVO methodology to determine the optimal parameters of a PEFC. Ali et al. [33] applied I-GWO to successfully identify the unknown parameters of five types of commercial PEFCs.

Applying optimization algorithms (OAs) to the PEFC performance gives excellent results, which will allow the researchers to obtain the representative equation of a PEFC model successfully. The main idea behind this is that they can escape from the local minimum, which leads them to obtaining almost the global optimum [34]. Taking into consideration the mentioned methods and their promising results in generating high-quality solutions, these meta-heuristic processes have been widely used to model PEFCs, so there is much potential for improvement in this area [35].

A few applications of heuristic-based methods to define models of PEFCs are addressed. However, due to the high nonlinearity of PEFC models, the extraction of optimal values of parameters is crucial, which means that small deviations of these values may dramatically affect the quality of such models [23]. Therefore, the nonlinear behavior of the models motivates the researchers to develop and use improved algorithms like I-GWO to define the parameters of a PEFC. According to the aforementioned, OA can deal efficiently with uni and multi-modal engineering optimization problems [36].

This paper focuses the study on the computation of the optimal parameters for polarization curves’ prediction of PEFCs by using different OAs, which are based on natural phenomena. The collected data come from a lab-scale FC test system employing different conditions of RH. It is one of the most critical variables that affect the performance of a PEFC [37,38]. Thus, the main contributions of this work are:

• Analyze the RH impact on the performance of a PEFC, the sparsity and noise of the collected data, and its influence on optimization problems;
• Compare the OA efficiency through an error metric in the computation of the optimal parameters for the mathematical model of a PEFC;
• Build prediction models for the polarization curves of a PEFC at different RH levels.

The rest of this article is organized as follows: the experimental setup, which includes the lab-scale FC test system description and the experimental parameters setup, is given in Section 2. In Section 3, the mathematical model of PEFC based on the generalized steady-state model is presented. In Section 4, the principles of employed OAs, as well as a short description about them, are given. The methodology to build the OF, the conditions for the mathematical model of the PEFC, the initial parameters for the different OAs, and the metric error are presented in Section 5. Finally, the results and conclusions are drawn in Section 6 and Section 7, respectively.

2. Experimental Setup

2.1. Lab-Scale FC Test System

The test system used to conduct this study is the 850e Multi-Range FC Test System from Scribner Associates, Inc. It allows for monitoring and evaluating the performance variables of a PEFC in real-time and in an automated way. The system is constituted by a flow controller to set the stoichiometric ratio of the inlet gasses, water humidifier tanks, and heating resistors at the anode/cathode input. A heating resistor inside the PEFC hardware and a DC load to evaluate the PEFC voltage response under a specific current are part of the equipment. The system obtains different polarization curves based on several parameter setups. The collected data were obtained with the help of the FuelCell^® software from Scribner Associates Inc, North Carolina, USA. It was installed in a peripheral device connected to the FC test system. For additional technical information, readers can refer to the manual [39].

In detail, the single-cell has the following characteristics: an effective area of 25 cm², two graphite end plates with triple serpentine channels with a depth and width of 1 mm. Furthermore, a Membrane Electrode Assembly (MEA) with five layers composed of one polymeric membrane Nafion 212^®, two GDLs woven carbon cloth type with micro-porous layers (MPL), and two catalyst layers with 60 wt% of platinum nano-particles on Vulcan carbon from the cell. A representative schematic of a single PEFC is shown in Figure 1.

2.2. Experimental Parameters Setup

Ultrapure Hydrogen (H₂) and medical Oxygen (O₂) were employed as fuel and oxidant, respectively. In addition, Nitrogen (N₂) was used as purge gas at the pipes system. The inlet pressure of the gases was established at 55 PSIG, while ASTM Type I water (with 18 M$\mathsf{\Omega }$/cm minimum resistivity) was used to humidify the inlet gases. Since obtaining different polarization curves as a function of RH was proposed, the gas temperature was initially established at 40 ${}^{°}$C, increasing in steps of 10 ${}^{°}$C until reaching 80 ${}^{°}$C. Then, the fixed temperature of the MEA was configured at 80 ${}^{°}$C under ambient pressure conditions. Note that 100% of RH is reached when the gas temperature and the fixed temperature of the MEA have the same value. The RH was computed using empirical correlations provided by Espinoza-Andaluz et al. [38]. Furthermore, a current sweep was applied in the experiment with an increasing rate of 0.25 A/min, from 0 to a maximum value of 50 A. Finally, the stoichiometric ratio of the H₂ and O₂ flows were settled at 1.2x and 2.5x, respectively, where x represents the current load applied in the PEFC according to [40].

3. Mathematical Model of PEFC

Several mathematical models to predict the performance of PEFCs can be found in the literature [41,42,43,44]. In this work, that proposed by Maan [43] is implemented since it has been widely used in previous studies [18,45]. It has been shown to perform well for modeling the polarization curves considering all polarization losses. Herein, the cell output voltage is represented by the subtraction between the Nernst voltage and the polarization losses. The polarization losses are represented by three voltage drop regions (activation, ohmic, and concentration). Additionally, the expression can be multiplied by $n_{cells}$, which represent the number of cells. In this case, as it corresponds to a single cell, $n_{cells}$ has the value of 1. The resulting expression is given as follows:

$$V_{FC}=n_{cells}·\left(E_{Nernst}-V_{act}-V_{\mathsf{\Omega }}-V_{con}\right)$$

(1)

Initially, the cell voltage drops rapidly and exponentially due to the over-potential of activation necessary to begin the reactions defined by $V_{act}$. Then, the voltage drops gradually and linearly as it enters the ohmic zone ($V_{\mathsf{\Omega }}$). Finally, at high current densities, the voltage drops quickly again due to the effects of concentration ($V_{con}$).

$E_{Nernst}$ is defined as a thermodynamic potential of the cell, also known as reversible voltage. Considering the reference temperature at 25 ${}^{°}$C, the Nernst voltage can be determined as follows:

$$\begin{array}{c}{E}_{Nernst}=1.229-0.85·{10}^{-3}·(T_{fc}-298.15)\hfill \\ +4.3085·{10}^{-5}·T_{fc}·\ln \left(P_{H_2}·\sqrt{P_{O_2}}\right)\hfill \end{array}$$

(2)

where $T_{fc}$ is the PEFC operating cell temperature (K), while $P_{H_2}$ and $P_{O_2}$ are the partial pressures (atm) of hydrogen and oxygen, respectively, defined as:

$$P_{H_2}=\frac{R{H}_a·P_{H_{2}O}}{2}\left(\frac{1}{\frac{R{H}_a·P_{H_{2}O}}{P_a}·e^{\frac{1.635·\frac{I_{fc}}{A_m}}{T_{fc}^{1.334}}}}-1\right)$$

(3)

$$P_{O_2}=R{H}_c·P_{H_{2}O}\left(\frac{1}{\frac{R{H}_c·P_{H_{2}O}}{P_c}·e^{\frac{4.192·\frac{I_{fc}}{A_m}}{T_{fc}^{1.334}}}}-1\right)$$

(4)

where $A_m$ is the membrane area (cm²), $R{H}_a$ and $R{H}_c$ are the RHs of the gases flowing at the anode and cathode, respectively. $I_{fc}$ is the operating current (A) of the PEFC, and $P_{H_{2}O}$ is the saturation pressure of water (atm) defined as:

$$\begin{array}{c}\hfill P_{H_{2}O}=2.95·{10}^{-2}·(T_{fc}-273.15)-9.18·{10}^{-5}·{(T_{fc}-273.15)}^2\\ \hfill +1.44·{10}^{-7}·{(T_{fc}-273.15)}^3-2.18\end{array}$$

(5)

Following the describing equations, the activation losses $V_{act}$ can be determined as follows:

$$V_{act}=-\left[{\xi }_1+{\xi }_2·T_{fc}+{\xi }_3·T_{fc}·\ln C_{O_2}+{\xi }_4·T_{fc}·\ln I_{fc}\right]$$

(6)

where ${\xi }_1$, ${\xi }_2$, ${\xi }_3$, and ${\xi }_4$ are semi-empirical coefficients that describe the electrochemical processes in the activation over-potential zone. The oxygen concentration is represented by $C_{O_2}$ (mol/cm³). Both $C_{O_2}$ and ${\xi }_2$ are calculable parameters defined in Equations (7) and (8), respectively, whereas ${\xi }_1$, ${\xi }_3$, ${\xi }_4$ are unknown in the first instance:

$$C_{O_2}=\frac{P_{O_2}}{5.08·{10}^6·e^{-\frac{498}{T_{fc}}}}$$

(7)

$${\xi }_2=2.86·{10}^{-3}+2·{10}^{-4}·\ln A_m+4.3·{10}^{-5}·\ln C_{H_2}$$

(8)

Furthermore, $C_{H_2}$ is the hydrogen concentration (mol/cm³), which is computed as:

$$C_{H_2}=\frac{P_{H_2}}{1.09·{10}^{-6}·e^{\frac{77}{T_{fc}}}}$$

(9)

Next, the concentration losses are represented by $V_{con}$, and can be computed as follows:

$$V_{con}=-\beta ·\ln \frac{J_{max}-J}{J_{max}}$$

(10)

where $\beta$ is a temperature-dependent parametric coefficient that proportionally modifies the concentration factor. In this mathematical model, it is presented as an unknown coefficient. J is the current density (A/cm²), and $J_{max}$ is the maximum value for J.

The last polarization loss to describe is ohmic losses, $V_{\mathsf{\Omega }}$, defined as:

$$V_{\mathsf{\Omega }}=I_{fc}·(R_m+R_c)$$

(11)

where $R_c$ ($\mathsf{\Omega }$) is the connection resistance, i.e., the resistance to the electron flow due to all conductive elements. It should be constant when the PEFC works between 40 to 90 ${}^{°}$C. In general, it is difficult to predict its value, so $R_c$ is initially unknown. On the other hand, $R_m$ represents the membrane resistance ($\mathsf{\Omega }$) and can be computed as:

$$R_m=\frac{{\rho }_m·l}{A_m}$$

(12)

where l is the thickness of the membrane (cm) and ${\rho }_m$ is the membrane resistivity ($\mathsf{\Omega }$$·$cm), defined as:

$${\rho }_m=\frac{181.6·\left[1+0.03·\left(\frac{I_{fc}}{A_m}\right)+0.062·{\left(\frac{T_{fc}}{303}\right)}^2·{\left(\frac{I_{fc}}{A_m}\right)}^{2.5}\right]}{\left[\lambda -0.634-3·\left(\frac{I_{fc}}{A_m}\right)\right]·e^{4.18·\left(\frac{T_{fc}-303}{T_{fc}}\right)}}$$

(13)

where $\lambda$ is a modifiable parameter representing the degree of the polymeric membrane performance and is influenced by the feed gases’ relative humidity and stoichiometric ratio. In this work, lambda has a maximum possible value of 23 to assure better fitting according to [43] and a minimum value of 13 under ideal conditions.

Eventually, it is clear that, in the above equations, there are six unknown parameters required to be defined: ${\xi }_1$, ${\xi }_3$, ${\xi }_4$, $\lambda$, $\beta$, and $R_c$. They were briefly described as the necessary parameters to complete the polarization curves modeling. For more detailed information and deductions of the mentioned parameters, readers are referred to the following studies [43,44]. Generally, these parameters are obtained from experimental tests with the use of thermodynamics-electrochemical constants. Hence, there is not a manual where values of these parameters are stipulated. The unknown parameters demand great accuracy in order to effectively model the PEFC performance. For this, several studies have defined some limits to these parameters intending to establish bounds over the computation works [33,45,46,47,48,49]. Upper and lower limits for the mentioned parameters are presented in

Table 1. Upper and lower limitsfor unknown parameters of a PEFC mathematical model.

Limits	${\mathit{\xi }}_1$	${\mathit{\xi }}_3$ (1 × ${10}^{-5}$)	${\mathit{\xi }}_4$ (1 × ${10}^{-5}$)	$\mathit{\lambda }$	${\mathit{R}}_{\mathit{c}}$ (m$\mathsf{\Omega }$)	$\mathit{\beta }$
Upper	−0.8532	9.80	−9.54	23.00	0.80	0.5000
Lower	−1.1997	3.60	−26.00	13.00	0.10	0.0136

.

4. Optimization Algorithms

4.1. Principles

Mathematically, optimization aims to search a decision space supported by a set of decision variables that would produce the maximum acceptability considering the defined goal. The maximum acceptability is based on the objective function (OF) initially defined [50]. Thus, computational intelligence (CI) has been able to make meaningful contributions to solving such optimization problems (OPs) thanks to the development of several efficient and competitive search algorithms.

CI-based OAs are centered around their inspiration sources. Thus, based on their origins, OAs could be classified as nature-based and non-nature-based. Nature-inspired computing (NIC) is a CI branch dedicated to studying any nature-based strategy used in optimization. Many phenomena in nature have been a source of inspiration for developing new problem-solving techniques [51]. Biology-based algorithms are some of the most popular in nature-based OAs. Figure 2 shows the general flowchart that an OA follows in order to achieve its aim, where it can be highlighted that $p_1$, $p_2$, $p_3, \dots ,$ $p_n$ are the initial conditions or parameters with any OA that start to work.

Biology-based OAs like grey-wolf, ant–lion, and bird behavior are considered in the current study. At the same time, astronomy-based OAs like multi-verse and neural-network-based OA are also considered. The mentioned OAs are described in order to give an overview of them and their operation.

4.2. Multi-Verse Optimization (MVO)

The multiverse concept comprehends the fact that universes interact and even collide with each other. The multi-verse theory inspires the MVO algorithm based on three main concepts of cosmology: white holes, black holes, and wormholes. The search process of this population-based algorithm is divided into two phases: exploration and exploitation. White and black holes in MVO are applied to explore search spaces by the algorithm. Wormholes help the algorithm to work with the exploitation of the search spaces. We assume each solution is equal to a universe and each variable in this solution represents an object of this universe [32].

4.3. Improved-Grey Wolf Optimizer (I-GWO)

The GWO algorithm imitates the leadership hierarchy and hunting mechanism of grey wolves in nature. The hunting technique and the social hierarchy of grey wolves are employed to design a mathematical model, so the optimization process will take place using the GWO [36]. Several GWO variants have been developed to avoid the local optima and accelerate convergence speed by making specific modifications to the mechanism of GWO. The improvement in the algorithm is supposed to alleviate the lack of population diversity, the imbalance between exploitation and exploration, and will also avoid the premature convergence in GWO [52].

4.4. Ant–Lion Optimizer (ALO)

ALO mimics the hunting behavior of ant–lions in the early growth stage or larvae stage, so ALO mathematically models ant–lions hunting behavior. The ants play the role of search agents and move over the decision space while ant–lions can hunt them and become fitter [53]. The new ant positions will be constantly evaluated and compared with those of ant–lions. If an ant becomes fitter than its ant–lion, its position will be taken by the ant–lion. If the best ant–lion reached in the current iteration becomes fitter than the elite, then the elite will be updated. This process will continue until it finds the best solution.

4.5. Bird Swarm Algorithm (BSA)

Many bird species are gregarious, such as finches. They may roost in communes and fly in flocks [54]. These behaviors are considered emergent behaviors and are based on simple rules such as separation, alignment, and cohesion. Swarm behaviors develop complex motions and interactions through simple social interaction. Thus, based on these birds’ behavior and specific rules that can summarize the whole iteration process, the BSA can be developed. For more information about BSA, we recommend the reader to check [55], an article that simplifies the understanding of social behaviors of birds’ swarm in a few steps. Once the previously mentioned rules are defined, the proposed algorithm can be developed.

4.6. Neural Network Algorithm (NNA)

Artificial neural networks (ANN) are constructed out of computational models inspired by the natural nervous systems. The ANN structure comprises a few artificial neurons associated with each other to form the neural network (NN). NNs can be classified based on the structure into two kinds: feedforward and recurrent. Recurrent NNs have either a local or global feedback structure dependent on its sort. As with all the metaheuristic algorithms, the NNA begins its process by creating an initial population entitled example arrangements inside the search space [18].

5. Methodology

5.1. OF Construction

An OF, or a loss function, is the used function to minimize or maximize the cost function to find a feasible solution, which is considered the optimal solution. During this process, the aim is to find the global maximum or minimum since it guarantees that the cost function has entirely been optimized, instead of a local minimum or maximum that says loss function has been optimized but not wholly.

The total of the squared errors (TSE) between the measured (V) and computed or predicted ($\widehat{V}$) voltage points integrates the defined loss function [18], which in this case is represented by Equation (14). The predicted voltage points are computed in terms of the six aforementioned unknown parameters:

$$OF=minimized\left(TSE\right)=minimized\left(\sum _{i=1}^N{\left[V_i-{\widehat{V}}_i\right]}^2\right)$$

(14)

where i represents each voltage point, and N is the total of voltage points.

Generally, an OF is subjected to a set of inequality constraints since it helps to reduce and adjust the search space where the OA tries to find the optimal solution. For this case, the used constraints are displayed in Table 1 where, for each unknown variable, an upper and lower bound is defined.

5.2. Boundary and Initial Conditions for the Mathematical Model of the PEFC

As it was mentioned in Section 2.2, five RH levels were computed and used in the experiments, according to the different gases’ temperatures. Following the description of the mathematical model displayed in Section 3, specifically in Equation (10), there is a variable called $J_{max}$, which indicates that the maximum value of the current density. This value is found when the PEFC voltage of the polarization curve reaches zero.

Thus, with the collected datasets and by using the Curve Fitting tool from MATLAB, $J_{max}$ values are computed and showed in

Table 2. RH and $J_{max}$ values for the different gases temperature tested in the experiments.

Gas Temperature (${}^{°}$C)	RH (%)	${\mathit{J}}_{\mathit{max}}$ (A/cm2)
40	18.00	1.0277
50	27.70	1.6372
60	42.62	2.1890
70	65.57	2.3779
80	100.00	2.4529

.

5.3. Input Parameters for OAs

Every OA is characterized by having input parameters that begin the optimization process and search for the optimal solution. The most common parameters are the maximum number of iterations and the total population. Similar to other heuristic-based techniques, the choice of parameters is achieved based on trial and error methodologies over a lot of independent experiments and observing the performance reached of the OA under analysis. Sometimes, equations are adequately defined by OAs, which help to find the best or the most recommended values for its input parameters.

The numerical simulations of this work are implemented by using MATLAB environment (version R2091b-9.7.0.1190202) under Windows 10 Pro for Workstations 64-bit, with a processor Intel(R) Core(TM) i7-1065G7 CPU @ 1.30GHz 1.50 GHz and 12 GB of RAM.

For the five OAs implemented in this paper,

Table 3. Input parameters for OAs (* means NOT APPLICABLE).

OA	M	p	d	${\mathit{c}}_1$	${\mathit{c}}_2$	${\mathit{a}}_1$	${\mathit{a}}_2$	FQ
NNA [56]	500	50	6	*	*	*	*	*
BSA [55]	500	50	6	1.5	1.5	1	1	3
ALO [53]	500	50	6	*	*	*	*	*
I-GWO [52]	500	50	6	*	*	*	*	*
MVO [57]	500	50	6	*	*	*	*	*

shows the used input parameters. As mentioned, the values were selected based on trial and error and according to authors’ recommendations of the original papers.

M is the maximum number of iterations until reaching the convergence, p is the population size which is an important parameter, and it directly influences the ability to search for an optimum solution in the defined search space, and d is the dimension of the problem to be solved, i.e., the number of unknown parameters/variables.

On the other hand, for BSA, other additional parameters need to be defined:

• $c_1$ and $c_2$ are two positive constants which can be respectively called cognitive and social accelerated coefficients into the birds’ foraging behavior.
• $a_1$ and $a_2$ are two positive constants’ values between 0 and 2, related to the birds’ vigilance behavior.
• $FQ$ is the frequency in which each bird flies to another place.

For more details related to the BSA parameters, please refer to [55].

5.4. Error Metric for Predictions

To validate the different results computed through the OA, it is intended to compare the voltage points experimentally obtained with the predicted voltage points. To achieve this, root mean square error or RMSE (Equation (15)) is used as a common error metric in regression problems:

$$RMSE=\sqrt{\frac{1}{N}·{\sum }_{i=1}^N{\left[V_i-{\widehat{V}}_i\right]}^2}$$

(15)

where V is the measured voltage point, $\widehat{V}$ is the predicted voltage point, and N defines the number of data points considered in the training of the OA.

5.5. Noise Metric for Experimental Data

When data are noisy, it can also be interpreted as there are outliers in itself, i.e., values that break down the smoothness or continuity in most cases. Particularly in this kind of experimentally collected data, these two effects are typical.

Thus, to analyze the noise in data, a vector with all values can be considered, and taking each successive point, the absolute difference is calculated. Assuming that data describe a smooth curve, these differences must be so small. However, if an outlier is detected, the difference is larger and more significant. This concept is based on total variation (TV) metric [58], which is applied as a basic distance norm in bounded data into a mathematics area. Equation (16) shows the expression to compute the average total variation (aTV):

$$aTV=\frac{1}{K-1}·\sum _{t=1}^{K-1}|V_{t+1}-V_t|$$

(16)

where K is the length of the vector of values, $V_{t+1}$ is the voltage point taken at time $t+1$, and $V_t$ is the voltage point taken at time t.

6. Results and Discussion

6.1. RH Impact on Performance of a PEFC

The polarization curve is a characteristic performance curve of a PEFC that shows the operating range according to the independent variables settled. Each polarization curve is composed of three defined polarization regions: activation, ohmic, and concentration. The activation polarization zone is characterized by a rapid voltage drop at the beginning of the curve. It is produced by irreversibilities that PEFC has to overcome to operate in a stable form. The irreversibilities are the back diffusion and crossover fuel, which take place in the polymeric membrane. The next region is the ohmic polarization zone. Its name is attributed to ohm’s law; this zone is characterized as the fuel cell voltage is proportional to the current density, and its slope depicts the cell’s resistance. Finally, it is the concentration polarization zone; to enter this region, the polarization curve suffers a second rapid drop. It is produced commonly due to the lack of fuel at the catalyst layer. At high current density, the water generation within the PEFC increases, obstructing the porous media. Thus, the electrochemical reactions can not hold, and the performance drops sharply.

In Figure 3, the polarization curves obtained in the present study are shown for different RH of the reactant gases. RH is one of the most important independent variables that affect the performance of a PEFC, as shown in studies [37,38,59]. It is caused mainly for the polymeric membrane used in PEFC. The membrane performance is higher when the RH of the gases is close to 100%, while the performance is lower when the RH falls in low values. This is the reason because the reactant gases have to enter humidified inside the PEFC. However, there are applications for PEFC that need to operate at low RH; for this, studies to analyze and quantify the effect of this variable are required.

As shown in Figure 3, the polarization curve with the best performance is tagged with 808080, i.e., 100% RH. It is reached when the temperature of the reactant gases is equal to the temperature of the PEFC. In this case, the maximum current density for the PEFC is around 2.0 A/cm². However, the recommended range to operate is between 0.2 up to 1.8 A/cm², which depicts the ohmic polarization zone, i.e., the more stable region for operating in real applications. Then, it can be noted that the PEFC performance decreases when the RH decreases too. The lower RH presented in this study is when the temperature of the reactant gases is 40 ${}^{°}$C, while the cell’s temperature is maintained at 80 ${}^{°}$C. In this regard, a 16% RH value is obtained. The maximum current density with previously settled conditions is around 0.9 A/cm². At the same time, its ohmic polarization zone is a little unstable and has a reduced operating range concerning the other curves. This will be discussed in the later section. In addition, there is a difference of 1.4 A/cm², comparing the maximum current densities for 100% RH and 16% RH. This indicates the significant influence of HR on the performance of a PEFC.

6.2. Sparsity and Noise of Experimentally Collected Data

All polarization curves were obtained using the same experiment configuration, where the sample rate to collect data was 0.25 A/s. The test range goes from 0 to 50 A, which is the current load maximum limit of the system. Data are obtained until the experiment drops on its own due to the performance based on conditions settled. In addition, the threshold voltage set up is 0.20 V for all cases. This setup of parameters is according to [40]. When the RH has low values, data are obtained in lower amounts due to decreased operating efficiency, which is considered sparsity in data. Hence, at low RH, the threshold voltage is probably reached for low current density, making the test stop because PEFC would be working in critical conditions.

On the other hand, concerning noise in data,

Table 4. Noise metric applied to measured voltage points in each polarization curve.

RH (%)	18.00	27.70	42.62	65.57	100.00
aTV (1 × 10−3)	7.550	4.351	3.055	2.565	2.028

shows the aTV computed for all collected subsets, where each one is associated with a RH.

At first glance, in Figure 3, not much could be inferred about the data dispersion. However, the aTV parameter shows that there is greater dispersion for curves with lower RH than when working with higher RH. This means that, with a higher value of RH, the experiment becomes more efficient and can reach 2.0 A/cm² without stopping due to a voltage drop, which is very usual at low RH values. Furthermore, the prediction of the OAs at a low value of RH could cause less efficient results than at a high value of RH. As there are less data and more noise, it could originate trends not so fitted with the experimental curve compared to the predicted curves at low and high RH.

6.3. Computed Unknown Parameters

Before computing the parameters, it is important to mention that the unknown parameters can be influenced by temperature, relative humidity, the stoichiometric ratio of the feed gases, and the composition or microstructural configuration of the MEA. For each configuration of the aforementioned variables, the unknown parameters will have their corresponding associated values. Despite this, the microstructural configuration does not change in the current study since the MEA used was the same in all tests. Furthermore, the stoichiometric ratio remained constant in all tests. Hence, in the current study, the primary variable that determines the parameters is the relative humidity.

With the described parameters in Section 5.3, OAs were performed and the unknown parameters are obtained according to the different RH employed in each experiment. To make this methodology reproducible and more statistically significant, the OAs were replicated 100 times. It was similarly done in the described tests in the articles corresponding to each different OA. Next,

Table 5. Unknown parameters computed and performance information for the modeled polarization curve at 18% RH.

Parameters	NNA	BSA	ALO	I-GWO	MVO
${\xi }_1$	−1.1149	−1.1149	−1.1994	−1.1154	−1.1997
${\xi }_3\left(1\times {10}^{-5}\right)$	9.7999	9.7999	7.9154	9.7918	7.9091
${\xi }_4\left(1\times {10}^{-5}\right)$	−25.9999	−25.9999	−25.9999	−25.9999	−25.9999
$\lambda$	13.0000	13.0000	13.0000	13.0004	13.0000
$R_c$ (m$\mathsf{\Omega }$)	0.8000	0.8000	0.8000	0.7997	0.8000
$\beta$	0.0645	0.0645	0.0645	0.06440	0.0645
$ϵ\left(1\times {10}^{-2}\right)$	1.8664	1.8664	1.8665	1.8664	1.8665
$t_{training}$ (s)	35.8161	28.1427	34.8917	87.3449	30.7480
$G{O}_{best}\left(1\times {10}^{-2}\right)$	3.1351	3.1351	3.1353	3.1352	3.1353
$G{O}_{worst}\left(1\times {10}^{-2}\right)$	3.1353	3.4024	3.6350	3.1376	3.3613
$G{O}_{std}\left(1\times {10}^{-7}\right)$	7.2980	9504.9800	11,392.5770	48.9820	5320.4450

,

Table 6. Unknown parameters computed and performance information for the modeled polarization curve at 28% RH.

Parameters	NNA	BSA	ALO	I-GWO	MVO
${\xi }_1$	−1.0192	−1.0320	−1.1828	−1.1194	−1.1997
${\xi }_3\left(1\times {10}^{-5}\right)$	9.7997	9.5129	6.1490	7.5609	5.7752
${\xi }_4\left(1\times {10}^{-5}\right)$	−24.4553	−24.4632	−24.4484	−24.4732	−24.4214
$\lambda$	14.2717	14.2110	22.2903	14.2064	14.1801
$R_c$ (m$\mathsf{\Omega }$)	0.8000	0.7994	0.8000	0.7993	0.8000
$\beta$	0.0582	0.0578	0.0587	0.05779	0.0579
$ϵ\left(1\times {10}^{-3}\right)$	5.8480	5.8490	5.8490	5.8490	5.8500
$t_{training}$ (s)	68.9683	88.9019	82.6906	201.3102	80.6157
$G{O}_{best}\left(1\times {10}^{-2}\right)$	4.9934	4.9944	4.9949	4.9955	4.9959
$G{O}_{worst}\left(1\times {10}^{-2}\right)$	5.2459	6.3568	6.6148	5.1094	6.4502
$G{O}_{std}\left(1\times {10}^{-7}\right)$	6.5355	4079.0520	4193.5260	163.7320	4009.5830

,

Table 7. Unknown parameters computed and performance information for the modeled polarization curve at 43% RH.

Parameters	NNA	BSA	ALO	I-GWO	MVO
${\xi }_1$	−0.9891	−0.9891	−1.1605	−1.0420	−1.1997
${\xi }_3\left(1\times {10}^{-5}\right)$	9.7999	9.8000	5.9719	8.6177	5.0981
${\xi }_4\left(1\times {10}^{-5}\right)$	−20.3072	−20.3074	−20.3087	−20.3098	−20.2981
$\lambda$	22.9999	22.9999	23.0000	22.9994	22.6254
$R_c$ (m$\mathsf{\Omega }$)	0.7999	0.7999	0.8000	0.7996	0.7997
$\beta$	0.0753	0.0753	0.0753	0.07520	0.07435
$ϵ\left(1\times {10}^{-2}\right)$	7.9120	7.9120	7.9130	7.9130	7.9190
$t_{training}$ (s)	191.0309	135.0847	111.0385	385.0889	108.358416
$G{O}_{best}\left(1\times {10}^{-2}\right)$	1.1456	1.1455	1.1458	1.1459	1.1476
$G{O}_{worst}\left(1\times {10}^{-2}\right)$	1.3451	1.4892	1.4107	1.1553	1.6837
$G{O}_{std}\left(1\times {10}^{-7}\right)$	6011.3840	7731.1360	7497.3040	142.4000	9409.8590

,

Table 8. Unknown parameters computed and performance information for the modeled polarization curve at 66% RH.

Parameters	NNA	BSA	ALO	I-GWO	MVO
${\xi }_1$	−0.9696	−0.9692	−1.0620	−1.0076	−1.1997
${\xi }_3\left(1\times {10}^{-5}\right)$	9.7916	9.7999	7.7240	8.9422	4.6376
${\xi }_4\left(1\times {10}^{-5}\right)$	−15.0015	−15.0028	−14.9989	−14.9991	−15.0440
$\lambda$	22.9999	22.9999	23.0000	22.9785	23.0000
$R_c$ (m$\mathsf{\Omega }$)	0.7999	0.7999	0.8000	0.7997	0.8000
$\beta$	0.0657	0.0657	0.0658	0.0657	0.0655
$ϵ\left(1\times {10}^{-2}\right)$	5.7980	5.7980	5.7990	5.8000	5.8000
$t_{training}$ (s)	298.3339	141.1235	88.8788	168.1751	130.6152
$G{O}_{best}\left(1\times {10}^{-2}\right)$	0.6724	0.6724	0.6725	0.6728	0.6729
$G{O}_{worst}\left(1\times {10}^{-2}\right)$	0.9459	1.2206	9.0701	0.6787	1.1342
$G{O}_{std}\left(1\times {10}^{-7}\right)$	6883.3250	10,564.2350	5710.6560	116.0200	8616.5400

and

Table 9. Unknown parameters computed and performance information for the modeled polarization curve at 100% RH.

Parameters	NNA	BSA	ALO	I-GWO	MVO
${\xi }_1$	−0.9843	−0.9843	−1.1993	−0.9906	−1.1997
${\xi }_3\left(1\times {10}^{-5}\right)$	9.8000	9.7999	4.9772	9.6620	4.9669
${\xi }_4\left(1\times {10}^{-5}\right)$	−9.5400	−9.5400	−9.5400	−9.5413	−9.5400
$\lambda$	23.0000	23.0000	23.0000	22.9963	23.0000
$R_c$ (m$\mathsf{\Omega }$)	0.8000	0.7999	0.8000	0.7988	0.8000
$\beta$	0.0703	0.0703	0.0703	0.0701	0.0705
$ϵ\left(1\times {10}^{-3}\right)$	3.5840	3.5840	3.5850	3.5870	3.5850
$t_{training}$ (s)	95.4813	75.8022	159.2261	230.6826	110.6548
$G{O}_{best}\left(1\times {10}^{-2}\right)$	2.3894	2.3894	2.3903	2.3931	2.3907
$G{O}_{worst}\left(1\times {10}^{-2}\right)$	3.8519	6.7273	3.9850	2.5145	5.4944
$G{O}_{std}\left(1\times {10}^{-7}\right)$	3668.1860	7889.4520	3649.9570	173.5280	5206.5160

resume the computed unknown parameters, the calculated prediction errors for each OA, which are defined by the RMSE metric, the information about the training time elapsed by the OAs, and the computed global optima (GO), during those training stages.

As observed, for all OAs, the error is computed at different RH percentages. Comparing the algorithm, the efficiency of the OA is contrasted in both training error and training time. As it has been mentioned, NNA is one of the most recent proposals in the OA field, which has shown a better accuracy over other OAs, considering variables of analysis like the training error and training time.

Following that hypothesis, in this work, NNA once again shows an advantage in general over BSA, ALO, I-GWO, and MVO. For example, for RH values such as 18%, 66%, and 100%, the training error is the lowest, while, in other RH values like 28% and 43%, the training error of I-GWO or BSA is lower than NNA. If the training time is analyzed, in the two most recently mentioned RH percentages, NNA is better since it takes less time to train the model. Thus, in a nutshell, NNA applied to this optimization problem shows excellent performance, considering that these experiments have been developed under different operation conditions.

Another point of analysis of OAs is about the convergence curves, which indicate how the OA has optimized the unknown parameters to reach the global minimum or maximum of its cost function through the training in the established epochs or iterations. In this case, Figure 4 shows the convergence curves for all trained OAs, considering the performed experiment under 100% RH, where it is clear that all curves reach the convergence along the iterations; nevertheless, the way in which each OA operates in its search space is different, and it varies along the iterations as it can be seen in Figure 4.

Furthermore, for performance information given in the global optima reached through OAs, it is important to understand the statistical significance and level of reproducibility of the models. For instance, in the aforementioned tables, the global optima obtained along the trained OA models are similar, though the standard deviation could be high or low, which indicates that, if the standard deviation is high, it is highly likely that the model is not reproducible since, each time, the results can present large variations. Instead, if the standard deviation is low, it is highly likely that the model is reproducible, since the obtained results each time are very similar and reliable. Following this analysis, NNA shows a better advantage over the rest of the OAs in the $G{O}_{std}$ computed in all the trained models.

Next, Figure 5, Figure 6, Figure 7, Figure 8 and Figure 9 show the prediction curves for the best OA, i.e., in this case NNA, at different RH values.

On the other hand, in Figure 5, Figure 6, Figure 7, Figure 8 and Figure 9, the modeled polarization curves are presented according to the OA with the best performance, where it is also clear that data points do not trace a smooth curve in low RH percentages. By contrast, there is a bit of noise in the data. As the RH values increment, this noise vanishing and data points describe curves more smoothly. To justify this, it is crucial to analyze the operating conditions under the lab-scale FC test system works from a technical point of view, as it is described in Section 6.2.

The results have shown how the OAs help obtain the best parameters missing for modeling the polarization curves with a considerably reduced time since the experimental methods for obtaining the unknown parameters involve extensive time and resources consumption. Furthermore, as previously mentioned, the NNA had the best performance compared with other natural OAs, even in tests with noise (low relative humidity). Moreover, the models run in the current study yield as final results modeling polarization curves with higher accuracy that can be used in the future for designing, controlling, and simulating complex engineering problems of PEFC.

7. Conclusions

Some different OAs have been applied to this engineering problem to define the optimal values of PEFCs’ unknown parameters so correct modeling, control, and simulation can be carried out for a specific experiment related to some of the PEFCs’ applications as polarization curves. The tests are performed at different RH values to show the effectiveness of OAs for every dataset that was experimentally obtained. According to the numerical results, several conclusions can be made as shown next:

• In all cases, for the five different OAs (NNA, MVO, BSA, ALO, and I-GWO) applied in this PEFC optimization problem, well fittings between measured and predicted voltage points are reached when using the optimal values of the unknown parameters for PEFCs.
• Statistical performance measures have been made to evaluate the efficiency and competency of the five algorithms used to carry out this experiment, concluding that NNA proves to give the best results for the optimal values of PEFCs’ unknown parameters in almost all scenarios.
• NNA and MVO show a better response than the other three algorithms when only referring to the training time. Although NNA optimal values are better, focusing on the metric used to measure the error between measured and predicted voltage points.
• The optimal values for PEFCs’ unknown parameters were obtained at different RH percentages. The NNA optimizer performed the best training in three out of five scenarios, as, at RH values of 28% and 43%, I-GWO and BSA showed more accurate results when focusing on the statistical performance measures.
• The comparisons that are properly detailed in this paper give the authors enough information to confirm and conclude that the NNA optimizer has a better performance and shows the best results, comparing it with other highlighted OAs.

Likewise, about the experimental section, a few conclusions have also been performed, as follows:

• The polarization curves obtained show the great influence that the RH has on the performance of a PEFC, obtaining a significant decrease in range operation of 1.4 A/cm² when the RH is modified from 100% RH to 16% RH.
• Another RH impact on the performance of a PEFC is the increased presence of data sparsity and noise at low RH due to the PEFC works under critical operating conditions. Therefore, it also affects predicted curves that have a lower fit.
• To sum up, the data sparsity and noise also affect OAs, since it represents a higher difficulty in finding the optimal parameters that allow for reaching the best fitting of the polarization curves to data through the PEFC mathematical model.