Design optimization and parameter estimation of a PEMFC using nature-inspired algorithms

Numerical simulation of proton-exchange membrane fuel cells (PEMFCs) requires an adequate model and precise parameters for reproducing their operational performance quantified by the polarization curve. Bioinspired algorithms are well-suited for optimization. The simulator is stressed by inputting thousands of randomly generated parameters, and hence, a robust numerical model is required. Once the proper model and parameters reproduce the experimental data, they can be used for design improvement. This article proposes a reformulation of a macrohomogeneous mathematical model to provide higher numerical stability to the solutions. We introduce optimization problems for parameter estimation and design optimization by applying three bioinspired algorithms to maximize its performance and minimize the platinum mass loading mPt\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$m_\mathrm{{Pt}}$$\end{document}. The results are validated by comparing the experimental polarization curves with those simulated from the estimated parameters. We compare a base design’s performance with the optimized design for maximum performance. We also compare a base design with the optimized design for minimum mPt\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$m_\mathrm{{Pt}}$$\end{document}. The results show that the particle swarm optimization requires the lowest computational cost and performs the best in most cases, fitting the experimental data with errors lesser than 10-17\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^{-17}$$\end{document}. The minimization of mPt\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$m_\mathrm{{Pt}}$$\end{document} reduces the amount by 42% compared to the base case.


Introduction
Fuel cells are devices of great interest since they have high efficiency in converting hydrogen's and oxygen's chemical energy directly into electricity by electrocatalytic reactions. In recent decades, substantial research efforts have been devoted to its development at a commercial level. There are various applications of interest, such as stationary power generation, trains, buses, cars, or other transports, including airplanes, toward achieve a hydrogen-based economy (Wang et al. 2018).
A proton exchange membrane fuel cell (PEMFC) is an electrochemical device that consists of an ionomer membrane in contact on either side with a catalytic layer, a diffusion layer, and a terminal plate with flow channels for the reactants. Figure 1 shows the basic structure of a PEMFC. At the anode, hydrogen flows into the channels through the diffusion layer (DL) to the catalytic layer (CL), usually formed by the Pt/C catalyst. At Pt catalytic active sites, H 2 is oxidized to produce protons and electrons. Protons move to the cathodic catalytic layer through the membrane, while elec-  (Outeiro et al. 2008;Sun et al. 2020). Anode reaction: Cathode reaction: Global reaction: The electrode-reactions rate depends on the CL structural parameters, such as the composition, distribution of catalytic materials, and manufacturing process. Mathematical models provide useful information for studying PEMFC performance and improving its design. Different studies related to the modeling and simulation of a PEMFC have been carried out (Outeiro et al. 2008;Secanell et al. 2007;Song et al. 2004;Khajeh-Hosseini-Dalasm et al. 2010;Marr and Li 1999;Heidary et al. 2016). However, mathematical models require a set of parameters that must be estimated because they are very difficult to measure experimentally.
In the literature, three principal classes of mathematical modeling for the cathode catalyst layer have been studied: interface models (Lu et al. 2002;Berning and Djilali 2003a, b), macrohomogeneous models (Tiedemann and Newman 1975;Marr and Li 1999;You and Liu 2001;Khajeh-Hosseini-Dalasm et al. 2010;Song et al. 2004;Heidary et al. 2016), and agglomerate models (Secanell et al. 2007;Shah et al. 2007;Wang et al. 2007). The macrohomogeneous models are of intermediate complexity and have the advantage of providing the expected results at low and medium current densities (Khajeh-Hosseini-Dalasm et al. 2010). This model stands out for being an intermediate alternative between the interface model and the agglomerated model. It provides adequate information to study the relationships between the catalyst layer's physical parameters and the performance of the PEMFC (Secanell et al. 2007;Khajeh-Hosseini-Dalasm et al. 2010). In addition, it permits the analysis of the performance of the fuel cell through the polarization curve by using the transport and diffusion of species through the catalytic layer. A macrohomogeneous model consists of a set of coupled nonlinear differential equations that describe oxygen diffusion, overpotential, and current density (Khajeh-Hosseini-Dalasm et al. 2010;Song et al. 2004).
In mathematical modeling of PEMFC, often, the proposals privileged simple and intuitive governing equations to implement numerical modeling software (Secanell et al. 2017). As a consequence of the simplicity, there is a possible pay-off to the application scope of the models. In the macrohomogeneous model, the applicability is limited by numerical instabilities when solving the differential equations associated with the modeling. As the current densities grow, the solution profiles abruptly change their gradient in a small region and reach a behavior similar to singular perturbed differential equations (Kopteva and O'Riordan 2010). Thus, this model is used in PEMFC for medium current densities, where ohmic losses occur. This topic is addressed in the numerical stability discussion part of this article where we provide directions to improve such stability.
Numerical simulation of a mathematical model implies a set of parameters, but frequently they are not directly obtained. Therefore, estimating parameters is necessary to analyze or describe the operation and design of a PEMFC. There are a variety of methods and techniques for parameter estimation. In particular, algorithms inspired by nature and with stochastic character are widely used because they are search methods that can deal with multiple minima. In addition, they can optimize parameters to improve the energetic performance in a PEMFC. The parameter estimation problem has been approached by nature-inspired algorithms mainly for the semi-empirical model (Priya et al. 2018;Chen and Wang 2019;Askarzadeh and dos Santos Coelho 2014;Blanco-Cocom et al. 2021a). For instance, parameters such as maximum current density, active area, equivalent contact resistance to electron conduction and other parametric coefficients have been estimated using a genetic algorithm in Outeiro et al. (2008), Ohenoja and Leiviskä (2010), Moo et al. (2006) and Zhang and Wang (2013), simulated annealing in Outeiro et al. (2008), Outeiro et al. (2009), particle swarm optimization in Ye et al. (2009), Mirjalili (2020, artificial bee swarm in Askarzadeh and Rezazadeh (2011a), artificial immune system in Askarzadeh and Rezazadeh (2011b), seeker optimization algorithm in Dai et al. (2011), harmonic search in Askarzadeh and Rezazadeh (2011c), optimization algorithms based on p-systems in Zhang and Sanderson (2009), and differential evolution algorithm in Turgut and Coban (2016).
There exists a diversity of research works related to parameter estimation. In particular, modifications have been explored to improve the standard algorithms for PEMFC optimization using different objective functions. For instance, the minimum of the mean of the sum of squared error was used in an adaptation of the cuckoo search algorithm with onlooker bee search (ObCS) in Zhu and Wang (2019). The squared sum of the errors was used with the improved fluid search optimization algorithm (IFSO) to estimate the parameters of a PEMFC in Qin et al. (2020). This proposal simulates Bernoulli's principle in fluid mechanics. Another bioinspired algorithm is the chaotic binary shark smell optimization studied in Han et al. (2019), which imitates the ability of the sharks to find their prey incorporating the chaotic movements of the shark in the water with the objective function given by the minimum of the sum of squared error. Algorithms based on the behavior of coyotes in the environment are analyzed in Yuan et al. (2020). This proposal is called the Coyote optimization algorithm, which uses the objective function consisting of the sum of the absolute error between the experimental data and the simulation of the polarization curve of the PEMFC. Other algorithms based on social behavior between species, such as the imperialist competitive algorithm, shuffled frog-leaping algorithm, and firefly optimization algorithm, are explored in Kandidayeni et al. (2019). The fitness function is based on minimizing the sum of squared errors (SSE) from the output voltage of a PEMFC stack and the estimated voltage given by the semi-empirical model. Recent studies on semi-empirical mathematical models for PEMFCs have applied different meta-heuristics. The Grasshopper optimization algorithm (GOA) ) is implemented with the sum of square errors (SSE) between the experimental data and the corresponding estimated results (El-Fergany 2017). The total of square deviations (TSDs) between the actual and calculated results is used in Mirjalili et al. (2017) as objective function in the Salp swarm optimizer (SSO) in the semi-empirical model in two test case studies of commercial stacked PEMFCs, namely NedStack PS6 and BCS 500-W PEM generator. In Xu et al. (2019), the authors applied the JAYA-NM algorithm to esti-mate 7 parameters of the semi-empirical model, and they used the sum of the squared error (SSE) between the output voltage of the actual PEMFC stack and the model output voltage. The Chaos Owl Search Algorithm (COSA) to estimate parameters of the semi-empirical model of the PEMFC stacks by minimizing the sum of squared error (SSE) between the estimated and the measured output voltage for two different case studies was analyzed in Zhang et al. (2020). Moth-flame Optimization (Ben Messaoud et al. 2021) was applied to estimate 6, 7 and 11 parameters of the semi-empirical model and used three commercial PEMFCs to obtain experimental data (250 W, NedSstack PS6 and BCS 500-W). Farmland Fertility optimization algorithm (FFA) was applied in Menesy et al. (2021) to study four different PEMFC stacks with a semiempirical model. The author used the sum of the squared errors between the experimental data and the simulated by the algorithm. The Univariate Marginal Distribution Algorithm (UMDA G ) was applied in Blanco-Cocom et al. (2021b). In this study, the authors applied three different formulations for the objective function to estimate parameters in the model and compared the results with those other algorithms in the literature. In Yakout et al. (2022), the authors studied four different evolutionary optimizers, namely Gorilla Troops Optimizer (GTO), Honey Badger Algorithm (HBA), African Vultures Optimization Algorithm (AVOA), and Chimp optimization algorithm (ChOA), which were applied for parameter estimation in a mathematical model for a solid oxide fuel cell (SOFC). The authors used the mean-square error (MSE) between the calculated and measured terminal voltages. A review that analyzes various mentioned algorithms with the corresponding objective functions is presented in Yang et al. (2020).
Fuel cells are part of the renewable energy economy system nested in the energy storage and conversion subsystem (Wang et al. 2018). Therefore, mathematical modeling and simulation can reduce the design costs and assembly times by reducing experimentation in the laboratory. In this sense, to optimize the fuel cell performance, it is necessary to have the optimal configuration of the reaction parameters, such as platinum mass, diffusion coefficients, and pressures, and design parameters such as cell dimensions and membrane thickness (Secanell et al. 2007;Marr and Li 1999;Mansouri et al. 2020). Hence, these parameters can be optimized to obtain lower manufacturing costs (Wang et al. 2018) or to improve performance (El-Fergany 2018; Amadane et al. 2018). For instance, in Marr and Li (1999), the authors use the macrohomogeneous model to infer that the effective use of platinum catalyst decreases when the current density increases. In Song et al. (2004), the problem of maximizing the current density in a PEMFC at a given electrode potential subject to porosity is approached. In Khajeh-Hosseini- Dalasm et al. (2010), the minimization of activation overpotential for a given current density is pre-sented, and other authors have analyzed the geometry in the flux channels of the gas diffusion layer (Lu et al. 2002;Berning and Djilali 2003a, b;Wang et al. 2018).
In this article, we introduce a unified approach for the numerical simulation and design optimization of a PEMFC. We select the macrohomogeneous model because of its reliability, precision, and scale of detail. The macrohomogeneous model had not been studied from the perspective of parameter estimation or optimization through meta-heuristics. The model is applied for two purposes: (a) estimating a set of unknown physical parameters and (b) optimizing a design for improving its performance or minimizing its cost. The objective function is nonlinear and nondifferentiable, as it is modeled as the sum of squared errors of a discrete polarization curve. Each point in the polarization curve corresponds to the solution of a system of differential equations. To the best of our knowledge, there is no previous research on optimizing the cell design either in terms of performance or cost using numerical simulations. As the objective function presents multiple local minima, evolutionary algorithms (EAs) are the appropriate optimizers, considering that they can escape from local minima, do not require derivatives and are suited for dealing with global optimization problems. Evaluating a candidate solution, that is, an individual in the context of EAs, requires a PEMFC numerical simulation. We note that the simulation software receives random inputs that could lead to numerical failures. In this sense, a contribution of this work is to present a robust simulation algorithm that is extensively stressed by inputting thousands of random configurations that are successfully treated. Another contribution is the comparison of evolutionary algorithms from different families for estimating parameters and optimizing the PEMFC performance. In addition, we introduce objective function definitions for three goals: estimating a set of unknown parameters for reproducing experimental data, improving a base design by maximizing its performance, and reducing the cost via the minimization of the amount of platinum, without a significant detriment of a base performance. The method is validated by comparing the estimated parameters with the true ones and by comparing the estimations with reported estimates in the specialized literature. Note that if the parameters are the same as those reported in the literature and the numerical model reproduces the polarization curve, then the numerical model is valid. Finally, performance optimization and cost minimization are evidenced by comparing the optimized design with a base design.

The macrohomogeneous mathematical model
The macrohomogeneous model considers a catalyst layer with a finite thickness in which electrochemical reactions take place (see Fig. 2). It is integrated by catalyst clusters of carbon and platinum, a solid GDL phase, a Nafion ionomer phase, and void spaces. Through CL, the diffusion of oxygen, overpotentials, and current occur (Khajeh-Hosseini-Dalasm et al. 2010). This model describes the catalyst layer as a onedimensional system in a steady state. The membrane and the GDL phases penetrate into the CL, and the void spaces are fully flooded by water (Khajeh-Hosseini-Dalasm et al. 2010;Heidary et al. 2016). The standard mathematical model is described by the following nonlinear boundary value problem (Marr and Li 1999;Khajeh-Hosseini-Dalasm et al. 2010;Song et al. 2004;Heidary et al. 2016), using the notation shown in Table 1: on Ω = [0, l c ], with l c being the thickness of the CL, subject to the boundary conditions: where the transport due to diffusion of oxygen in the cathode CL following Fick's law is described in Eq. (4), the activation overpotential caused by ohmic losses in Eq. (5), and the electrochemical reaction rate following the Butler-Volmer equation in Eq. (6). In Table 2, we summarize the key algebraic expressions of the system (4)-(6).     Song et al. (2004) and Marr and Li (1999)

Polarization curve
The polarization curve is a graph associated with the performance of a fuel cell and presents three characteristic regions (see Fig. 3): • Region 1: Activation losses. The region localized at low current densities is due to the activation losses or kinetic activation of the electrochemical reaction that occurs at the catalytic layers. Its shape is due to the exponential voltage drop for low current densities. • Region 2: Ohmic losses. It is observed in the zone of moderate current densities and has a linear trend due to electronic resistances or ionic mass transport losses. • Region 3: Concentration losses. The region observed at high current densities is due to the resistance to the mass transfer of oxygen from the bulk to the active sites localized in the catalytic layer. Its shape is characterized by the rapid drop in voltage until it reaches the value of 0 V.
The polarization curve for the macrohomogeneous mathematical model is described by (7) according to Khajeh-Hosseini-Dalasm et al. (2010). In this paper, it is used for characterizing the cell. The parameters are estimated to approximate an experimental polarization curve, and they are optimized to increase the polarization curve area or to maintain a similar polarization curve while minimizing the platinum mass loading.
where R Ohmic is the specific ohmic resistance of the cell and E Nernts is the reversible cell voltage (Song et al. 2004;Khajeh-Hosseini-Dalasm et al. 2010), where E 0 = 1.229 V, T represents the cell temperature (K), and P H 2 and P O 2 are the partial pressures (atm) of hydrogen and oxygen, respectively.

Reformulation of the macrohomogeneous model as singularly perturbed ODE system
A drawback of the macrohomogeneous model is the numerical instabilities, mainly in the concentration loss region (Fig. 3). On the other hand, a requirement for applying bioinspired algorithms is to have a reliable numerical simulator. This is so since the simulator is stressed by inputting thousands of random parameters, and it must deliver precise numerical solutions for each input. Hence, in this section, we revisit the nonlinear coupled system of the PEMFC to improve its numerical stability. The most common coupled system of nonlinear ordinary differential equations consists of the three equations of the system (4), (5), (6). This system of ODEs formulated as an initial value problem is usually solved using standard numerical techniques such as the Runge-Kutta method (Song et al. 2004;Marr and Li 1999;Khajeh-Hosseini-Dalasm et al. 2010;Heidary et al. 2016). With these methods, the system presents numerical stability if I δ belongs to Region 1 and Region 2 of the polarization curve. However, the solvers become unstable if I δ belongs to Region 3 (concentration losses, see Fig. 3).
To obtain more realistic solutions, the macrohomogeneous system should be formulated as a boundary value problem (Khajeh-Hosseini-Dalasm et al. 2010;Heidary et al. 2016), where we need to estimate the initial value condition for Eq. (5) (overpotential).
As described in Song et al. (2004); Khajeh-Hosseini-Dalasm et al. (2010); Marr and Li (1999) Eq. (4) is obtained from: with the boundary values, From physical experiments in PEMFC, it has been estab- Shen et al. (2011) and Chan et al. (2012). Then, Eq. (9) is a second-order singularly perturbed differential equation, and its numerical implementation is important to obtain regularized solutions. To address this issue, we introduce two approaches for improving the numerical stability of the system. The first is rewriting Eq. (4) using the change of variable given by subsystem of Eq. (11), The second approach is given by using D eff O in the second equation, rewriting Eq. (11 as Eq. (12), On each subsystem, we have a first-order singularly perturbed ODE system. First, the subsystem of Eq. (11) is more numerically stable than the standard Eq. (4). To clarify this, we identify the slow variable as N O 2 and the fast variable as O 2 . Now, we consider the change of variable Then, we have a new formulation: With this change of variable, called homogenization (Pavliotis and Stuart 2008;Kadalbajoo and Patidar 2002), the diffusive coefficient disappears. However, the equations become completely determined by the slopes N O 2 and di dτ . It is important to note that these profiles have fast growth in a small subdomain Ω f of Ω, which implies that there exists a dependency of the solutions with the correct approximations of N O 2 and di dτ . With the second proposed subsystem of Eq. (12), we use a similar analysis, and we use the same change of variable to transform it into In this case, we observe that the coefficient D eff O 2 appears in the two equations, which compensates for the effect of the fast growth of the slope N O 2 and provides more numerical stability for the state variable O 2 . With this last formulation, we obtain the highest numerical stability among the three models. The reformulation of the system of equations is justified since a robust and stable numerical method is necessary due to the complexity of the combinations of values in the parameters that are evaluated when applying the optimization methods. The system of ODEs in this work is solved using the Runge-Kutta with an implicit Lobatto technique (Kierzenka and Shampine 2008).
To demonstrate the numerical stability of our proposals, we carried out numerical experiments with the three systems according to the following: • Mod 1 is the formulation of the system solved in the following order: Eqs. (4), (5), and (6). This is the most widely used system in the literature (Khajeh-Hosseini-Dalasm et al. 2010;Song et al. 2004;Marr and Li 1999). • Mod 2 is the first formulation introduced in this article, solved in the following order: subsystem of Eqs. (11), (5), and (6). • Mod 3 is the second formulation introduced in this article, solved in the following order: subsystem of Eqs. (12), (5), and (6).
The boundary conditions for Mod 2 and Mod 3 are: In Fig. 4, we present the results of the current density i(z) using the boundary condition i(l c ) = 9672 mA cm −2 for the three systems. The expected solution curve for the current density must increase within a small interval close to zero, tending to 9672, and must stabilize at that value. Thus, the solution obtained with the system Mod 1 presents an erroneous profile. Figure 4 shows that the red dotted line (Mod 2 ) presents the expected curve profile. Hence, the result from system Mod 2 is more numerically stable than that from system Mod 1 . Nevertheless, in Fig. 5, in a close-up of the solutions, we notice that Mod 1 , the solid green line presents the largest abrupt changes, while Mod 2 presents small variations around 9672 mA cm −2 . Finally, the solution using system Mod 3 (blue solid-dotted line) avoids any numerical issue in contrast to those obtained from systems Mod 1 and Mod 2 . Then, the solution of i(z) using system (12), labeled Mod 3 , is the most stable of the three approaches and closely reproduces the expected polarization curve.

Objective functions and algorithms for parameter estimation
Nature-inspired algorithms are a type of meta-heuristic (Wong and Ming 2019) that emulate the behavior of species in nature to improve a fitness function. Regarding the family of algorithms, they emulate different natural processes, and as a consequence, the performance varies from problem to problem and from one family to another. This paper is devoted to three of the most successful ones: The genetic algorithm (GA) (Priya et al. 2015;Ohenoja and Leiviskä 2010) emulates the evolution of a population of individuals who compete for inheriting their genes. The fittest are selected for reproduction via recombination of genes. In the present setting, genes are the parameters of the PEMFC, and the fitness is the error with respect to the experimental data. Particle swarm optimization (PSO) (Ye et al. 2009) emulates swarm behavior, such as flocks of birds. The particles are updated using a velocity, which is a combination of information of the best particle in a neighborhood that is defined as a graphical model, the best global particle and inertial vectors. The swarm of particles moves through the landscape fitness. Then, they concentrate the search in the most promising region. Finally, the estimation of distribution algorithms (EDAs) (Valdez et al. 2010) was derived from statistical and probabilistic modeling of GAs. They estimate a search probability distribution and sample new candidate solutions. Even though EDAs address the same kind of problems as other nature-inspired algorithms, they present additional benefits, such as convergence to the best solution known and the possibility of incorporating problem-specific knowledge or search bias, which promotes the search around the local and possible global minimum and circumvents setting arbitrary stopping criteria for the algorithm.
In (Blanco-Cocom et al. 2021b), the authors applied different meta-heuristics to a similar mathematical model (the semi-empirical model for a PEMFC). The reported simulations of PSO, GA and UMDA G are similar to those reported with other algorithms, such as JAYA, CMAES, MPSO, ARNA-GA, HABC, GOA and CHHO. In the above work, three of the algorithms with the best improvements in solutions were GA, PSO, and UMDA G . For this reason, in this work, we select these three algorithms to compare our results.
All the algorithms require as the main input an objective function to minimize/maximize. Accordingly, we introduce three objective functions and a brief recall of the optimization algorithms used for approximating the optimum. The first objective function is used for estimating unknown parameters, the second for the performance maximization, and the third for the minimization of the platinum mass loading m Pt . The last two objective functions are part of the contributions of this article. In Fig. 6, we present the general scheme of the methodology. The upper flowchart depicts the general methodology, Fig. 6a, while the lower flowcharts depict the optimization algorithms in Fig. 6b-d. Our proposed methodology identifies the design parameters or the parameters of interest and selects an optimization process (parameter estimation, performance maximization, or minimization of m Pt ); with this selection, the objective function is defined for the optimization algorithm (GA, PSO, and UMDA G ); see Sect. 4.1. After optimization, we obtain a set of optimal design parameters under the objective function selected. Finally, we generate the optimized or fitted polarization curve.

Parameter estimation problem
Given a set of unknown parameters θ = {θ 1 , θ 2 , . . . , θ l } ∈ = [L 1 , U 1 ] × [L 2 , U 2 ] × · · · × [L l , U l ] ⊂ R l , l is the number of unknowns. We define the objective function, SSE : → R + , as the sum of the squared error between experimental data obtained from a polarization curve and the subject to where m is the number of the experimental data, V FC i and V z i ; θ are the i-th experimental datum and the voltage output of the macrohomogeneous model, respectively, at current density z i , and L h and U h are the lower and upper bounds for θ h . The parameter estimation problem consists of determining the set θ ∈ that returns a minimum value for SSE( θ).

Performance maximization problem
We consider maximizing the voltage generated by the cell according to the polarization curve. Then, this problem is equivalent to minimizing the difference of an ideal voltage to the actual voltage delivered by the numerical simulation (see Fig. 3). Thus, the second optimization problem is Eq. (17) given by: subject to where z i is a current density value and θ ∈ are the l parameters to estimate in the model for the fuel cell. Note that the ideal voltage, 1.2 V in this case means that the cell maintains a maximum constant voltage in all the regions of the polarization curve. This is impossible, but the optimization problem looks for the cell that best approximation of this ideal performance.

Minimization of the platinum mass loading
An objective function to minimize the platinum mass loading, and consequently, the cost of the cell, follows the idea of generating a polarization curve with the same area as the experimental polarization curve, this is defined with the first term of Eq. (19), where d A = (A sim − A exp )/A exp , A sim is the simulated polarization curve, and A exp is the experimental polarization curve. The second term is the variable m Pt ∈ [L m pt , U m pt ], which is the platinum mass loading: With this formulation, we preserve the performance while minimizing the platinum mass loading.

Genetic algorithm
Genetic algorithms are based on the processes of biological evolution, and they evolve a population of individuals that are a codification of candidate solutions. The population evolves via crossover, mutation, and selection operators. The crossover preserves genetic information; that is, it recombines existing partial solutions, the mutation applies small perturbations to instances, generating values that do not depend on the current population, and the selection defines the fittest individuals for reproduction, biasing the search toward the most promising regions. In each iteration, the algorithm selects a set of the best-fitted individuals according to their fitness values. These are the parents, who generate individuals with expected better qualities called children. Through the generations, the population evolves until a stopping criterion is reached. We use the GA according to Labs (2020). This version is as follows: • Generate a random population and compute their fitness based on the objective function. • Score the elements of the population following the fitness values. • Scale the fitness to obtain the expectation value.
• Select the parents based on a selection criterion and the expectation value. • Select an elite subset of the current population based on the best expectation values and pass them to the next generation. • Create a child population from the selected parent population using the mutation and crossover operators. • Replace the old population with elite and children populations. • The algorithm iterates until the stopping criteria are reached.

Particle swarm optimization
Particle swarm optimization is an algorithm inspired by the social behavior of animals, such as ants, fish, and birds. The algorithm is based on the criterion of the existence of an "invisible" method to share information between unstructured and generally dispersed biological groups. It emulates the interaction between the members of a biological group (Meng and Pian 2016). The PSO version used here (Mirjalili 2020) is based on the algorithm described in Kennedy and Eberhart (1995), using modifications suggested in Mezura-Montes and Coello (2011) and in Pedersen (2010). The algorithm works as follows: • Create an initial swarm of particles and assign initial velocities. • Evaluate the objective function at each particle location.
• Determine the best location and the best (lowest) function value. • Compute new velocities using the current velocities, the particle best location, and the best location of their neighbors. • Until the algorithm reaches a stopping criterion, update the particle positions (the new position is the old position plus the randomly scaled velocity), velocity, and neighbors.

Univariate marginal distribution algorithm with Gaussian models
The estimation of distribution algorithms (EDAs) were derived from the probabilistic modeling of genetic algorithms (Goldberg et al. 1993(Goldberg et al. , 1989. They are based on the dynamics of populations where the offspring of a generation are generated by sampling from a probabilistic distribution model, described by the parent set. Hence, the learning of these algorithms comes from estimating a parametric probability distribution from the selected set, with the goal of associating the highest probability to the most promising regions. This is different from other nature-inspired algorithms in which the evolution of the populations is carried out by crossover and mutation operators (Goldberg et al. 1993(Goldberg et al. , 1989Bandyopadhyay et al. 1998;Kargupta 1996;Kargupta and Goldberg 1997;Lobo et al. 1998;van-Kemenade 1998). The univariate marginal distribution algorithm with Gaussian models (UMDA G ) is a kind of EDA (Mühlenbein and Paaß 1996;Larrañaga and Lozano 2002;Larrañaga et al. 2000a, b;Bosman and Thierens 1999;Blanco-Cocom et al. 2021a) that assumes that the probabilistic distribution of the populations is a l-dimensional joint probability function that is factored as a product of l univariate and independent probability functions. Then, at each iteration, the algorithm computes the join probability given by Each function P k,i (θ i k,h ) is a normal distribution function with mean and variance, where p k,i (θ i k,s |S θ k ) are weights associated with the N ind individuals, for 0 < i ≤ l.
In the algorithm, at each iteration, the weights p k,i (θ i k,s |S θ k ) are computed using a threshold ν and a selected subset S θ k ⊂ θ k , which is called the truncation method (Valdez et al. 2010). The assigned probability to each θ i k,s in the iteration k for 0 < i ≤ l and θ k,s ∈ θ k is computed following the function of the threshold ν, The stop criteria are a function of the norm of the normalized variance of the optimized parameters (Blanco-Cocom et al. 2021a), This criterion provides a natural way to stop the algorithm compared with the standard criteria as a function of the number of generations.

Results
In this section, we apply the evolutionary algorithms presented in Sect. 4 to the macrohomogeneous model in different scenarios. In Sect. 5.1, we use algorithms to estimate a set of parameters to reproduce the experimental polarization curve. In Sect. 5.2, we perform optimization to increase the voltage generated by the cell and the minimization of the platinum mass loading used in the cathode.

Validation for parameter estimation
This section presents comparative optimizations with three experimental datasets: Song et al. (2004), Khajeh-Hosseini-Dalasm et al. (2010 and Salva et al. (2015). Using these datasets, we estimate parameters for reproducing the polarization curve with the three optimization algorithms: GA, PSO, and UMDA G . We carried out 30 executions of each algorithm to obtain statistical results for comparing their performance. For all simulations, we set the parameters for the GA and PSO as suggested by MATLAB (Labs 2020;Mirjalili 2020). For the three algorithms (GA, PSO, UMDA G ), we fixed the population size as 100 individuals (particles) and the number of generations as 200 nvars, where nvars is the number of variables and tolerance is 10 −6 .

Dataset 1: Song et al. Song et al. (2004)
The first dataset is obtained from the research of Song. Song et al. (2004). The experimental polarization curve for this dataset presents a whole desirable polarization curve as shown in Fig. 3. We estimate the set of parameters given by θ = { c , m Pt , ρ Pt , D eff O 2 ,GDL }. The first three parameters were selected because we consider that they are important in the modeling. D eff O 2 ,GDL appears in the initial condition for state variable O 2 and describes the effective diffusion coefficient of oxygen in GDL; more details are provided in Song et al. (2004). The values of the base parameters of the macrohomogeneous model are presented in Appendix: Table  7.
The upper and lower bounds for searching each parameter are shown in Table 3, and these intervals are set according to the results reported in Song et al. (2004). In Table 4, the statistical results obtained from the estimation of parameters are presented by applying the three evolutionary algorithms. The results are similar to the order of approximation for an SSE of less than 10 −7 . In Fig. 7, we plotted the polarization curve of the best solution for each estimate. In Fig. 7b, we observe that the best approximation to the experimental data is that delivered by the UMDA G , but the standard deviation

Dataset 2: Khajeh-Hosseini-Dalasm et al. Khajeh-Hosseini-Dalasm et al. (2010)
The second dataset is obtained from Khajeh-Hosseini-Dalasm et al. (2010) and describes regions 1 and 2 of the ideal polarization curve (see Fig. 3). Originally, the experimental data are from Ticianelli (1988). We select to estimate the set of parameters given by θ = {m Pt , ρ Pt , c }, and the upper and lower bounds for each parameter are shown in Table 3. They are selected according to the reported results in Khajeh-Hosseini-Dalasm et al. (2010) and Song et al. (2004). The values of the complete base parameters used in Khajeh-Hosseini-Dalasm et al. (2010) are presented in Appendix: Table 8.
Comparative statistical results are presented in Table 4, and the results show that the estimations have an SEE on the order of 10 −17 for each algorithm. In Fig. 8, we plotted the polarization curve of the best solution for each estimate. In this case, the UMDA G has means in the three parameters that is closer to the reported values compared to the means of each parameter obtained by GA and PSO.

Data set 3: Salva et al. Salva et al. (2015)
The third dataset is taken from Salva et al. (2015). In this case, the authors only report the experimental data without performing mathematical modeling. Therefore, in this work, we try to adjust the macrohomogeneous model to the experimental data, thus providing a model with the adjusted parameters that can reproduce it. We estimate the loads of platinum and carbon and the porosity of the CL. Then, θ = {m Pt , m c , c }.
The proposed values of the base parameters used in this work for the macrohomogeneous model are shown in 9. The upper and lower bounds for the search space for optimization are shown in Table 3.
A comparative plot of the best parameter estimations is shown in Fig. 9, and we can observe that the polarization curves are similar for the three algorithms. Table 4 shows the statistical results, observing that the parameters reach stable SSE values of approximately 0.21417 for the three algorithms with a standard deviation of less than 10 −6 .

Performance maximization and platinum mass loading minimization
In the previous subsection, we have shown that the three comparison algorithms perform adequately for the three datasets; nevertheless, the experimental dataset that is reproduced with the highest accuracy is dataset 2 of Khajeh-Hosseini- Dalasm et al. (2010). Hence, in this section, we optimize the cell performance according to Eq. (17) and the platinum mass loading according to Eq. (19). In both cases, we use the dataset of Khajeh-Hosseini- Dalasm et al. (2010) because the model closely reproduces the experimental data. We apply the three comparison algorithms considering that all of them are statistically consistent and deliver a low error.

Performance maximization
We aim to optimize the performance required for the PEMFC to deliver a voltage as close as possible to the ideal voltage of (a) (b) Fig. 9 a Comparative plots of the polarization curves with the best set of estimated parameters for the UMDA G , GA, and PSO for data from (Salva et al. 2015); b zoom of the approximations at i = 3600A m −2  Table  3 (Dataset 2).
The parameter values after optimization for the three algorithms are shown in Table 5. Figure 10  For the sake of completeness, the voltage for each current density value is increased at least 0.42% and up to 0.44%. The total voltage gain is 4.26% with respect to the base polarization curve.
The results validate that the same methodology and numerical modeling with a different objective function are successfully applied for both purposes, namely estimating unknown parameters and optimizing the design performance.

Minimization of platinum loading m Pt
To minimize the platinum loading, we use the objective function given by Eq. (19). The search limits are set based on the ranges of values reported in the specialized literature, with the purpose of searching on the largest feasible range. The lower bound of the search limits for m pt is set to 0.2, as suggested in Qi and Kaufman (2003). The upper bound is set to 0.5, which is the average of the values suggested in Khajeh-Hosseini-Dalasm et al. (2010). The lower bound for c is 0.2, which is the suggested value in Song et al. (2004), and the upper bound is 0.5, which is the suggested value in Khajeh-Hosseini-Dalasm et al. (2010). Finally, the search range for p pt is Wang et al. (2007) and Zhang and Wang (2013) in order to have sufficient flexibility to compensate for the possible reduction of m pt . The optimum approximations are presented in Table 6. The GA and PSO provide very similar sets of values. The best results reduce the quantity of platinum by 42%, while the voltage difference is on average less than 1.5%, in contrast to the base design. That is, the optimized cell delivers almost the same voltage as the base design, while the most expensive component has been significantly reduced. Figure 11 compares the polarization curves with respect to the base design, and the polarization curves are very similar for the three algorithms. The reported values for a low-cost cell design are approximately m pt = 0.2, p pt = 24, and c = 0.5. An important result is that the optimized value of m Pt belongs to the range of 0.20 ± 0.05 mg cm 2 , which is congruent with the experimentally reported results in Qi and Kaufman (2003).   Section 5.1 demonstrates that our methodology is capable of fitting experimental data. In the first two cases, in Figs. 7 and 8, the datasets and mathematical models are reported in the literature. Hence, we a priori know that those datasets can be well-fitted by the model. In contrast, the third data set in Fig. 9 is reported as experimental data without an accompanying model, and it is the estimated model with the proposed methodology that fits the data but not as well as the others. For instance, the curvature is more pronounced in the data than in the fitted polarization curve. There are several possible explanations for this deviation. The first is the model's capacity to recreate the curve. However, in Fig. 7, the model is capable of changing the curvature's sign to fit the data. Hence, another possible explanation is related to the nonestimated parameters due to the constraint of the family of solutions. This is the most plausible explanation because for the third dataset, we do not have the model information that we have for the other two cases. In other words, for the third case, some material properties or conditions given by nonestimated parameters could differ from those used. A solution to alleviate this problem is to estimate a larger set of parameters. However, this would increase the complexity of the problem by increasing the space of possible solutions. This would increase the computational cost, and the numerical instability of the model without the guarantee of improving the current estimation. Furthermore, even if the fitting in the third case is not as precise as in the other two cases, it is close enough to reproduce the general cell performance, and all the estimated parameters are physically plausible.
Regarding the performance and computational Time of the optimization algorithms, all of them report a low error, with a precision lower than what is required for the physical specifications. According to Tables 4, 5 and 6, the best performing are GA and PSO; nevertheless, the PSO required less time than the GA. Hence, in general, we recommend the last for the optimization process. Furthermore, PSO reports the lowest variance in the parameters and objective function values. Table 4 shows the average execution times in hours of a single parameter estimation process achieved by each algorithm. The times increase due to the number of experimental points that are used, since for each experimental point, the system of differential equations proposed in this work (Mod 3 ) must be solved. If the estimation or optimization uses a large amount of experimental data, the time required is increased. In addition, the numerical complexity increases if the experimental point is in the concentration loss zone, where a finer mesh is required to find the appropriate solution compared to the experimental data corresponding to the activation zone (Blanco-Cocom et al. 2022). The most expensive operation is the objective function evaluation, which requires a numerical simulation of a cell, but it is the same for all algorithms.
Thus, the computational complexity of the meta-heuristics depends on the population size and the number of generations, that is, O(n pop n gen ). These values are different for each algorithm and, in general, depend on the empirical convergence of the meta-heuristics and are a priori unknown. Hence, a fair computational complexity approximation is, actually, the computational time reported in Table 4.
Once a model reproduces the cell performance and the estimated parameters are physically plausible, we argue that the model reproduces the performance even if the cell conditions change. Therefore, the model and estimated parameters can be used to simulate various scenarios and select the most favorable one. In this sense, Sect. 5.2 presents two experiments of conflicting problems. The first, in Sect. 5.2.1, maximizes the performance, but it requires increasing the cell cost, and the improvement in voltage is quite small. In contrast, Sect. 5.2.2 presents the minimization of the platinum load, which has a significant impact on cost reduction, while the performance decreases only slightly. These cases are practical applications of parameter estimation and demonstrate the benefits of model fitting to improve the cell design.

Conclusion
In this article, we present various proposals to solve the parameter estimation problem for the macrohomogeneous model using three bioinspired algorithms: GA, PSO, and UMDA G . The first contribution of this work is a modified macrohomogeneous model based on its relationship with the singularly perturbed differential equations. This manner of solving these kinds of equations provides more numerical stability, even when meta-heuristic algorithms used for parameter estimation require hundreds of evaluations of random parameters. The second contribution is the use of three different objective functions depending for the three problems: parameter estimation, performance maximization, or platinum mass loading minimization. The third contribution is the use of three metaheuristic algorithms: GA, PSO, and UMDA G . Our results demonstrate that all of them provide satisfactory results. Nevertheless, the best are GA and PSO. In particular, PSO consistently delivers the best objective values and requires less computational effort. Nevertheless, the nature of the probabilistic construction of UMDA G produces a region of confidence; that is, for this algorithm, we know that the results are the best from a normal probability distribution with a small variance.
The proposed methodology can support fuel cell experts in simulating different conditions with more numerical stability. Additionally, PEMFC models very often include nonphysical parameters or parameters that integrate a set of physical variables that are not directly measurable. Hence, estimating these quantities is necessary even if one has access to the physical cell. Even though the performance maximization is not the most promising application because the improvement is not significant, the computational time to achieve it is in the order of hours. Hence, it is worth applying it. In addition, a similar procedure of minimizing the cost of an expensive component, such as platinum mass loading, produces a greater reward with a significant decrement of the cost, namely greater than 42%.
We conclude that a first step for applying stochastic optimization methods to the parameter estimation and design optimization of fuel cells is to develop a numerically stable simulation method. Another conclusion is that a similar assembly of an optimizer and simulator can be applied to a variety of problems, for instance, for parameter estimation or performance optimization, if the optimization problem is modeled adequately.
There are promising directions for future work. For instance, one could apply similar numerical techniques, based on perturbed differential equations, to more complex models such as those required to model fuel cells in 2D and 3D. Regarding the meta-heuristics, this work demonstrates that all of them deliver solutions with the required numerical precision. Hence, the niche of opportunity is to find another competitive metaheuristic to achieve computational cost reduction. Regarding the performance optimization, in this work, we do not find a significant improvement, which could be different if the model includes the cell geometry and the optimization algorithm proposes geometry modifications as well as other parameters required for 2D and 3D simulations. Finally, our results suggest a conflict between platinum mass loading minimization and performance maximization, and future work must contemplate exploring this conflict as a multiobjective problem. Table 7 Parameter values used in the baseline calculation used in (Song et al. 2004)  Expressions in Table 2 f , A s , i ref 0 Reference values to estimate m Pt = 0.35 (mg cm −2 ) ρ Pt = 21.4 (g cm −3 ) c = 0.4 Table 9 Parameter values used in this work to adjust dataset 3 obtained from Salva et al. (2015) Base Expressions in Table 2 f , A s , i ref 0 Parameters to estimate m Pt (mg cm −2 ), m C (mg cm −2 ), c