Tree Growth Algorithm for Parameter Identification of Proton Exchange Membrane Fuel Cell Models

T alarming increase in environmental pollution and the growing shortage in conventional energy sources, lead to searching for alternative environmentally friendly sources. Photovoltaic, wind and fuel cells (FC) are considered the most promising sources that are applied in small and large scales [1]-[5]. Thanks to its high efficiency, simple operating principle, low exhausts, durability and reliability, the fuel cell has attracted the attention of researches and decision-makers in the last decades [6]. Particularly, the proton exchange membrane fuel cell (PEMFC) has been demonstrated as a suitable solution for various applications, because of its high efficiency, short startup time and its good performance under low temperatures [7]-[9]. The operating temperature of PEMFC occurs in the range of 70-85°C. The operation of the PEMFC can be well understood if an efficient reliable model is established. Accordingly, the design of the PEMFC stack can be improved and optimum operation can be reached as well as the integration of the fuel cell stack into other devices and systems [5]. In the last years, many research papers tried to model the operation of PEMFC and many models have been provided in literature [10]-[15]. The semi-empirical model proposed by Amphlett et al. is considered as the most acceptable model for many researchers [10]. In this model the problem of modelling the PEMFC has been turned into an optimization problem of some unknown parameters in the parametric equations that describe the model. Meanwhile, the model includes multi-variable nonlinear equation that makes it very difficult to estimate the values of these parameters using traditional optimization techniques. Many researches have been demonstrated for extracting the optimal design parameters using optimization algorithms in both parameter and non-parameter modelling [16], [17]. The polarization curves of the PEMFC is highly nonlinear and the proposed model is based on a set of nonlinear equations that contain a set of semi-empirical adjustable parameters. The conventional optimization methods are not suitable for such complicated optimization problems, so metaheuristic optimization techniques are used for such issues [18], [19]. Based on parameter modelling, various algorithms have been proposed to solve the optimization problem of PEMFC parameters. In order to improve the parameters’ accuracy of the PEM fuel cell a hybrid genetic algorithm was proposed in [20]. The particle swarm optimization (PSO) has been proposed for PEMFC parameters estimation problem in [21]. The differential evaluation (DE), as well as the hybrid adaptive differential evaluation (HADE) algorithms, have been introduced for solving the optimization problem presented in [22], [23]. More recent optimization methods have been applied to solve the problem of PEMFC’s parameter such as: the harmony search algorithm (HAS) [24], the seeker optimization algorithm (SOA) [25], the multi-verse optimizer (MVO) [26], the adaptive RNA genetic algorithm [27], Eagle strategy based on JAYA algorithm and Nelder-Mead simplex method (JAYATree Growth Algorithm for Parameter Identification of Proton Exchange Membrane Fuel Cell Models


I. Introduction
T HE alarming increase in environmental pollution and the growing shortage in conventional energy sources, lead to searching for alternative environmentally friendly sources. Photovoltaic, wind and fuel cells (FC) are considered the most promising sources that are applied in small and large scales [1]- [5]. Thanks to its high efficiency, simple operating principle, low exhausts, durability and reliability, the fuel cell has attracted the attention of researches and decision-makers in the last decades [6]. Particularly, the proton exchange membrane fuel cell (PEMFC) has been demonstrated as a suitable solution for various applications, because of its high efficiency, short startup time and its good performance under low temperatures [7]- [9]. The operating temperature of PEMFC occurs in the range of 70-85°C.
The operation of the PEMFC can be well understood if an efficient reliable model is established. Accordingly, the design of the PEMFC stack can be improved and optimum operation can be reached as well as the integration of the fuel cell stack into other devices and systems [5]. In the last years, many research papers tried to model the operation of PEMFC and many models have been provided in literature [10]- [15]. The semi-empirical model proposed by Amphlett et al. is considered as the most acceptable model for many researchers [10]. In this model the problem of modelling the PEMFC has been turned into an optimization problem of some unknown parameters in the parametric equations that describe the model. Meanwhile, the model includes multi-variable nonlinear equation that makes it very difficult to estimate the values of these parameters using traditional optimization techniques.
Many researches have been demonstrated for extracting the optimal design parameters using optimization algorithms in both parameter and non-parameter modelling [16], [17]. The polarization curves of the PEMFC is highly nonlinear and the proposed model is based on a set of nonlinear equations that contain a set of semi-empirical adjustable parameters. The conventional optimization methods are not suitable for such complicated optimization problems, so metaheuristic optimization techniques are used for such issues [18], [19]. Based on parameter modelling, various algorithms have been proposed to solve the optimization problem of PEMFC parameters. In order to improve the parameters' accuracy of the PEM fuel cell a hybrid genetic algorithm was proposed in [20]. The particle swarm optimization (PSO) has been proposed for PEMFC parameters estimation problem in [21]. The differential evaluation (DE), as well as the hybrid adaptive differential evaluation (HADE) algorithms, have been introduced for solving the optimization problem presented in [22], [23]. More recent optimization methods have been applied to solve the problem of PEMFC's parameter such as: the harmony search algorithm (HAS) [24], the seeker optimization algorithm (SOA) [25], the multi-verse optimizer (MVO) [26], the adaptive RNA genetic algorithm [27], Eagle strategy based on JAYA algorithm and Nelder-Mead simplex method (JAYA-Tree Growth Algorithm for Parameter Identification of Proton Exchange Membrane Fuel Cell Models NM) [28], grey wolf optimizer (GWO) [29], hybrid Teaching Learning Based Optimization -Differential Evolution algorithm (TLBO-DE) [30], shark smell optimizer (SSO) [25], Cuckoo search algorithm with explosion operator (CS-EO) [31], selective hybrid stochastic strategy [32], bird mating optimizer [33], grasshopper optimizer (GHO) [34], Chaotic Harris Hawks optimization (CHHO) [35], and Modified Artificial Ecosystem Optimization (MAEO) [36].
To extract the optimal values of the unknown parameters on the proposed model, a metaheuristic technique inspired by the competition among trees in reaching the sources of food and light, Tree Growth Algorithm (TGA), is proposed to solve the optimization problems [37]. TGA is a very simple code that can be applied to various types of problems. For example, in [38], tree growth algorithm has been used to solve the localization problem in wireless sensor networks. The present work also applies TGA to the targeted problem. The main contributions of this study can be summarized as follows: • Developing the TGA for extracting the optimal best parameters of different PEMFC stacks.
• Demonstrating an effective mathematical model of PEMFC stacks which imitates the principle operation of different commercial PEMFCs through estimating the optimal values of the unknown parameters of these cells; • Studying the effect of cell temperature and reactants' pressure variations on the electrical characteristics of various PEMFC stacks; • A comprehensive comparison between the results obtained by the TGA and those obtained by other metaheuristic optimization algorithms has been provided; • Four different models of PEM fuel cells are introduced to verify the effectiveness of the TGA optimization method; • Parametric and non-parametric statistical analysis have been demonstrated to validate the goodness of the proposed metaheuristic technique; • The results obtained by the application of TGA prove its reliability and superiority in estimation of the effective parameters of PEMFC stacks.
The rest of the paper is arranged as follows: Section II introduces the model of the PEMFC as well as the objective function. The tree growth algorithm is briefly explained in section III. Section IV introduces the application of TGA for parameter estimation of different PEM fuel cell stacks under various operating conditions. Statistical analysis of the obtained results from the TGA is demonstrated in Section V. Section VI is dedicated to the conclusion.

A. PEMFC Model
A PEMFC consists of two isolated electrodes (anode and cathode) separated by a thin solid membrane able to conduct protons as described in Fig. 1 [39].
The overall chemical reaction occurring in the fuel cell can be represented as follows [39]: The voltage and current generated from a single fuel cell are too small, therefore a number of fuel cells are connected in parallel and/or series to form a fuel cell stack having a reasonable voltage and current rating. When N cells of identical fuel cells are connected in series, the resultant voltage of the stack will be calculated as follows, ( where, V cell denotes the voltage of a single fuel cell.
Because of the three types of voltage loss occurred in the fuel cell, namely, the activation loss V act , the concentration voltage loss V con , and the ohmic voltage loss V ohm , the terminal voltage of the fuel cell will be calculated as follows [40], where, E Nernst denotes the theoretic voltage of the fuel cell, which can be expressed as given in the following formula [33], where, T is the cell temperature in Kelvin; and denotes the partial pressures of the reactants (i.e. hydrogen and oxygen) at the inlet channels of the fuel cell stack (atm).
During the operation of the fuel cell, if the inputs to the PEMFC stack are hydrogen and natural air, then can be expressed as follows [40],

I fc
where, P c represents the pressure of the input channel at the cathode (atm), RH c denotes the relative humidity around the cathode, I fc is the current generated by the cell (A), A is the area of the membrane surface (cm2), sat P H O 2 represents the water vapor pressure at saturation, which is defined as follows [40],  If the inputs of the fuel cell stack are hydrogen and pure oxygen, the partial pressure of oxygen will be calculated as follows [40], In both previously mentioned conditions, the partial pressure of the hydrogen is calculated as follows, where, P a denotes the pressure at the anode side channel (atm), RH a represents the relative humidity of water vapor in the side of the anode.
The voltage loss due to the activation process V act can be calculated as, where, ξ 1 , ξ 2 , ξ 3 , ξ 4 denote semi-empirical coefficients; C O2 represents the concentration of the oxygen at the cathode (mol.cm-3) and is calculated as [40]: The second type of voltage loss occurred in the fuel cell, the ohmic voltage loss V ohm can be defined as, where, R M is the resistance of the membrane surface, RC denotes the resistance that the protons face when transferring through the membrane. The membrane resistance can be calculated as follows [40], where, ρ M denotes the specific resistance of membrane material (Ω.cm), l is the thickness of the membrane (cm). ρ M can be expressed as follows, where λ denotes an adjustable empirical parameter, which needs to be extracted.
The last type of losses that takes place in the fuel cell is the voltage loss due to concentration, which can be expressed as follows [40], where b denotes the parametric coefficient that needs to be estimated; J and Jmax are the current density and the maximum current density (A/cm2), respectively.

B. Formulation of the Objective Function
Through a closer look at the previous equations, it will be noticed that the operation of the PEMFC depends on a set of unknown adjustable parameters. For the optimization problem provided in this study, the target from the parameter estimation is to extract the optimal values of these parameters, which let the proposed model to match well with the experimentally measured data of the PEMFC. Therefore, the proposed objective function is a measure of the quality of the extracted parameters. The degree of matching between the calculated data and the data obtained from experiments can be formulated as the difference between the estimated output voltage based on the extracted parameters and the measured voltage. Accordingly, in this paper the Total Square Deviation (TSD) between the measured voltage of PEMFC and computed stack voltage is defined as the objective function (OF) [29]- [31], [41].
according to the following constraints, (16) where, X denotes a vector of the parameters that have to be estimated, V meas is the experimentally measured voltage, V stack is the computed voltage from the proposed model of the PEMFC, and N denotes the length of the data points.

III. Tree Growth Algorithm (TGA)
The TGA optimization algorithm was firstly proposed by Mostafa Hajiaghaei-Keshteli and Armin Cheraghalipour, which is inspired by the competition among trees in a certain area for absorbing the resources of light and food [37]. TGA in its process is divided into four phases. The first phase N1, which includes the best trees that have a good opportunity for acquiring food and light. Thus these trees will focus their effort on obtaining food, as the light source is already guaranteed by their height. In the second group, called the competition for the light group N 2 , some trees will move a distance by changing their angle to catch the light source through the tall elder ones. In the third group N 3 , which is called remove and replace group, the trees which do not have a good chance for growth are cut by the foresters and replaced by small new ones. In the final group called the reproduction group N 4 , the tall good trees begin to generate new small trees around the elder mother.
TGA algorithm is mathematically executed in the following manner. Firstly, an initial generation of trees N is randomly generated with respect to the predefined upper and lower boundaries and the fitness function is calculated for each individual search agent. After that, the initial generation is arranged with respect to the value of fitness function and the present best solution j T GB in the j-th iteration of the search space is determined. Then equation (17) is used to execute the local search to the individual agents in the first population N1 [37], where, θ denotes the reduction rate of trees due to their age and reduction of food sources in the surrounded area. r is a randomly distributed number between [0, 1]. T i j and T i j+1 are the current solution and the next solution in the population N 1 .
When the new solution is produced, a comparison with the previous solutions is conducted, in with the worst solutions will be ignored. The solutions of the second group N 2 will be moved towards the best solutions of the first group N 1 under different angles α. Equation (18) is used to define the distance between the selected solutions in the second group and the other solutions [38], [42], Once the distance d i has been calculated, two solutions x 1 and x 2 , having the minimum distance according to any solution, form a linear combination between the trees as described in the following expression (Fig. 2), where, λ denotes a control parameter in the range [0,1]. Fig. 2. Linear combination [37].
All possible solutions from the second group N 2 are moved between the two adjacent solutions with an angle α i as shown in Fig. 3 and formulated as follows, 2 2 j j N N i When the third subpopulation N 3 is reached, the worst solutions are removed and new solutions are generated randomly. Then a new population is created depending on the previous three populations, N = N 1 +N 2 + N 3 . The new generated subpopulation N 4 is created and modified according to the best solution in the first group N 1 using a masked operator. Then the new solutions from the randomly created subpopulation are combined with the N population.
The fitness of the created population N +N 4 is determined and the best solutions are stored and considered as the initial population in the next iteration of the tree growth algorithm search space. The roulette wheel or tournament is employed for this process. The flowchart that describes the procedure of the TGA algorithm is shown in Fig. 4.

IV. Simulation Results and Discussions
The simulation process has been performed using MATLAB software simulation package. The simulation was implemented using Intel Core i3-M370 CPU@2.40GHz and 4.00MB RAM Laptop. In order to validate the ability of TGA algorithm in identifying the unknown parameters of the PEM fuel cell, four different PEM stacks have been adopted from [30] - [33], [41]; namely BCS-500W FC, 250W PEMFC stack, SR-12-500W FC and Temasek 1kW FC stack. The estimation process has been used to identify the seven unknown parameters (ξ 1 , ξ 2 , ξ 3 , ξ 4 , β, Rc and λ). The upper and lower limits of the parameters, which have been optimized in this work are shown in Table  I [ 30], [31]. The characteristics of the four PEMFC stacks included in this case study are summarized in Table II. The TGA has been used in this study, while the following control variables have been adopted: the maximum number of iterations is adjusted at 500 iterations and the population size is 20. Due to the randomness nature of the technique used in this study, the best solution is taken as the minimum value obtained within 30 independent executions for the optimization algorithm. In addition, to validate the effectiveness of the developed algorithm, the dynamic response of the studied PEMFC stacks have been introduced. The results obtained from the proposed model is compared with the experimental data given in the datasheet of each fuel cell type.

A. Parameters' Estimation of the Proposed PEMFC Stacks
The TGA was applied for extracting the optimal values of the unknown seven parameters under the upper and lower boundaries given in Table I. The results of the minimum values of the objective function with the lowest SSE over the 30 runs as well as the optimized parameters of the four fuel cell stacks under study are given in Table  III. Once again, by a deep closer look to Table III, it can be seen that small absolute deviations of the estimated values proved the agreement between the experimental datasheet values and the computed values of voltage. Moreover, the results obtained by the application of TGA are compared with the results introduced in literature [23], [28]- [31], [41]. The minor values of the objective function point out the significance of the suggested TGA in solving the optimization problem. Fig. 5 shows the convergences curves of the SSE of the proposed optimization technique for all PEMFC stacks under consideration. It can be seen that the value of TSD is minimized after 98 iterations in the case of 250W FC stack, 220 iterations for BCS 500W stack, 310 iterations for SR-12 500W PEMFC and finally after 435 iterations in the case of Temasek 1kW fuel cell stack. It may be noticed that the convergence curves prove smooth, rapid and steadily progress to the optimal final values for all fuel cell stacks under the demonstration.
The I-V and I-P polarization curves obtained from applying the best solution obtained from TGA-based method compared with the measured data for the four types of PEMFC stacks are shown in Fig. 6(a) -Fig. 6(h). Through a closer look to Fig. 6, it can be admitted that the computed curves based on the estimated values of the unknown parameters give a good matching with experimentally measured ones, which validate the insignificant value of the TSD given in Table III [20], HADE [23], ARNA-GA [27], JAYA-NM [28], TLBO_DE [30]. It is clearly seen from Table III that the developed algorithm successes to find the minimum value of the TSD of 0.7496, which is less than the corresponding values obtained from other algorithms. The values of the optimized parameters as well as the value of TSD for BCS 500W PEMFC stack compared with that identified using GWO [29], SSO [41] and CS-EO [31] are provided in Table IV. Also, in this case, the developed algorithm provides the optimal values of the unknown parameters with an insignificant value of TSD of 0.083525, which ensures the coincidence between the polarization curves based on the proposed model and the measured ones. The estimated parameters compared with that obtained based on GWO [29], SSO [41] and CS-EO [31] for SR-12 500W FC stack are shown in Table V. Table VI shows the values of the seven unknown parameters using the TGA, GWO [29] and SSO [41] for Temasek-1kW PEMFC stack. Similarly, for SR-12 500W FC and Temasek 1kW PEMFC stack, TGA provides a superior performance over all optimization techniques as it is obviously seen from the small values of TSD. Although there are some deviations between the computed polarization characteristics and the measured ones provided in the datasheet of the manufacturers, they are acceptable in engineering standards.

B. Simulation Under Different Operating Conditions
In this section, various combinations of cell temperature and the inlet pressure of oxygen and hydrogen are proposed to demonstrate the performance of the fuel cell stacks under study. Accordingly, the polarization characteristics of the PEMFC stacks are elucidated and the behavior of the stack efficiency is demonstrated as well. Starting from the importance of the PEMFC efficiency and according to [43], the efficiency of fuel cell stack is calculated as follows, This equation is an approximation form for the exact efficiency called the voltage efficiency. v max denotes the maximum value of the output voltage generated from the PEMFC under hydrogen higher heating values, which equals to 1.48 V/cell. μ F is the utilization factor. It is assumed that the flow rate of hydrogen is controlled depending on the load condition, which leads to a constant utilization factor that equals to 95%.
Under the optimal values of the unknown parameters of the PEMFC stacks, the polarization curves (I-V and I-P characteristics) under different cell temperatures, while keeping the partial pressures of the reactants (P H2 /P O2 ) constant at the values given in datasheets, are produced to validate the effectiveness of the developed TGA-based model. To avoid repeated figures, only two of the four proposed stacks are demonstrated in this section. The I-V, I-P, and efficiency of the SR-12 PEM 500W stack at 303, 323, and 353K are described in Fig.   7(a) - Fig. 7(c), respectively. Fig. 8(a) -Fig. 8(c) show the I-V, I-P polarization characteristics as well as the efficiency of 250W FC stack at 323, 353, and 383K, respectively.
The impact of changing the pressures of the reactants in the inlet channels (P H2 /P O2 ) under constant cell temperature, described in the datasheet of the manufactures, is introduced. Fig. 9(a) -Fig. 9(c) and Fig. 10(a) -Fig. 10(c) depicted the I-V, I-P, and efficiency of the BCS 500W PEM stack and Temasek 1kW FC stack at pressures of (1/0.2075bar), (1.5/1bar) and (2.5/1.5bar), respectively.  In this section, a deep statistical analysis has been demonstrated to give a clear assessment of the developed algorithm. In addition, a sensitivity analysis is provided as a measure of the stability of the optimization algorithm proposed in this study. The comparisons between the different PEMFC stacks are based on many metrics, mainly minimum and maximum values of TSD, mean value of TSD, Median.
In addition to the previously mentioned metrics, standard deviation (SD), Relative error of the objective function (RE), Root mean square error (RMSE), Mean absolute error (MAE) and efficiency have been provided to examine the accurateness of the developed TGA-based simulation model, which are arithmetically calculated based on (22) to (26), respectively.
where, TSDi refers to the fitness function at each run. TSD min denotes the minimum fitness function obtained over the 30 executions.
TSD denotes the average value of the observed TSD over the simulation period. The summary of the studied metrics for the different PEMFC stacks is depicted in Table VII. It can be observed that the insignificant values of MAE and RMSE proved a well matching between the calculated values based on the estimated parameters and the measured ones. The values given in Table VII introduces a clear explanation of Fig. 11, in which the values of TSD in the case of Temasek 1kW stack over the 30 runs are changing in a narrow range.

VI. Conclusions
Extracting the values of seven unknown parameters of the PEMFC model is one of the most challenging points that attract the attention of many researchers. An effective PEMFC model based on TGA has been proposed in this paper, which is considered a suitable tool for simulation and performance evaluation of the PEMFC stacks under a wide range of operating scenarios. Many case studies have been performed, from which the estimated data, based on the optimal values of the unknown parameters, provide a good matching with the experimental data of different commercial fuel cell stacks. The results obtained from TGA-based model have been compared with different optimization methods. Different steady-state tests scenarios have been demonstrated to validate the effectiveness of the developed TGA-based technique. Moreover, statistical analysis has been conducted to measure the significance and precision of the optimal values obtained based on TGA method. Simulation results as well as statistical measurements emphasizes the superiority of the TGA over many optimization algorithms in extracting the parameters of proton exchange membrane fuel cells. In the future studies, the developed algorithm can be applied for simulating the dynamic behavior of PEMFC and solid oxide fuel cell (SOFC).