Proton Exchange Membrane Fuel Cells Modeling Using Chaos Game Optimization Technique

For the precise simulation performance, the accuracy of fuel cell modeling is important. Therefore, this paper presents a developed optimization method called Chaos Game Optimization Algorithm (CGO). The developed method provides the ability to accurately model the proton exchange membrane fuel cell (PEMFC). The accuracy of the model is tested by comparing the simulation results with the practical measurements of several standard PEMFCs such as Ballard Mark V, AVISTA SR-12.5 kW, and 6 kW of the Nedstack PS6 stacks. The complexity of the studied problem stems from the nonlinearity of the PEMFC polarization curve that leads to a nonlinear optimization problem, which must be solved to determine the seven PEMFC design variables. The objective function is formulated mathematically as the total error squared between the laboratory measured terminal voltage of PEMFC and the estimated terminal voltage yields from the simulation results using the developed model. The CGO is used to find the best way to fulfill the preset requirements of the objective function. The results of the simulation are tested under different temperature and pressure conditions. Moreover, the results of the proposed CGO simulations are compared with alternative optimization methods showing higher accuracy.


Introduction
In addition to the environmental impact of fossil fuels, they are declining and unsustainable sources of energy. For this reason, more focus has been given to develop clean renewable energy resources [1,2]. FCs are among promising energy resources. In FCs, electricity is generated by chemical reactions from chemical energy [3][4][5]. PEMFC is one of the FCs that have been widely used. The major characteristics and benefits of the PEMFCs are: (1) low operating pressure and temperature, (2) high performance ratio, and (3) no waste materials [6]. The fuel cell temperature requirements vary from 50 to 100 • C and the output power is between 30 and 60 percent under these temperature operating conditions [6,7]. In [8], a switched reluctance motor powered by the PEMFC's stack is used to optimize the operating performance. The performance is examined in terms of: (1) torque-ampere ratio, (2) torque smoothness, and (3) starting torque. The work in [9] applies an algorithm to a step-up converter in order to improve the performance of PEMFCs. A comprehensive review of the benefits of different design and optimization algorithms for PEMFCs is provided in [10]. In [11], a model is introduced for the PEMFC to save the undesirable thermal energy generated during the operation. The model uses the unoccupied capacity of the PEMFC to produce/store hydrogen.
Many researchers have recently drawn attention to the construction of a PEMFC model. A hybrid optimization approach is introduced in [12] for estimating the parameters of the PEMFC. The analysis in this work is performed on a 250 W stack, NedStack PS6, BCS 500-W, and SR-12 PEM 500 W fuel cells. To obtain results which are very similar to experimental results, the researchers strive to find a trustworthy PEMFC model with software programs and modeling that reduce substantial time and effort [13]. The PEMFC modeling is an important and promising research field. Both conventional and metaheuristic methods can be used to extract the parameters of the fuel cell. For example, the work in [14] introduces an evaluation of parameters' effects on the FC performance and model validation besides identification of optimal parameters. Taleb [15] presents monitoring of the operation of the PEMFC using a model-based approach. Experimental validation is done to the model using a frequency identification method. Then, time series identification is conducted using a least square method. In [16], a semi-empirical model is proposed based on the data of a three-cell SOFC stack, which predicts the performance of real stacks. The effects of key operation parameters were studied. The generalized reduced gradient technique is introduced in [17]. Moreover, in that context, a large number of complex analysis techniques were implemented to obtain an accurate fuel cell model [18][19][20].
Traditional optimization methods are unable to achieve sufficient results over a fair period. There are also drawbacks of traditional optimization methods which includes: (i) sensitivity for the initial problem estimate, (ii) reliance on solution exactness on the differential equation solver, and (iii) the possibility of getting trapped in a local solution. The probability of becoming trapped in a local optimal solution increases when dealing with highly nonlinear problems.
To efficiently deal with complex optimization problems, various algorithms have been developed. One of these algorithms is the CGO algorithm, which was introduced by Talatahari, Siamak, and Mahdi Azizi. [21]. The new optimization approach can be chosen to solve various problems in the field of engineering and other areas such as manufacturing applications and scientific models. The key idea of the CGO algorithm consists of some theories of chaos, which consider the construction of fractals according to the chaos game and the fractal problems of self-similarity. The findings have shown the superiority of the CGO in most cases. The theory of chaos shows that a few slight modifications in the initial guess of a complex system led to drastic variations in the conditions of these systems as they depend on the initial conditions.
In this paper, a new application of the CGO algorithm is presented to achieve the optimal PEMFC model design parameters. The objective of this study is to develop a PEMFC theoretical model. The results of the developed model are tested using measurements data provided by practical experiments. The experimentally measured data should yield close results to those found by the theoretical model. Considering that the data provided by the manufacturers are insufficient, and because of the nonlinear properties of the PEMFC as well as the seven design variables, the current-voltage characteristics of the PEMFC must be modeled. Thus, a NLOP is mathematically formulated in this study. In this formulated problem, the objective function is to achieve the minimum value of the square errors between the experimentally measured PEMFC voltage and the determined voltages by the developed model. The CGO algorithm is implemented explicitly to minimize the FF. The theoretical model's simulation effects are also tested at different temperatures and pressures.
The contribution of this research study can be summarized as: (1) the static theoretical models of the different types of PEMFCs are introduced and the estimated output voltages of the PEMFCs are compared with the measured voltages to evaluate the adequateness of the model, and (2) the robustness of the proposed CGO algorithm in minimizing the objective function is verified by comparing the results with other optimization methods.

Modeling of The PEMFCS Methods
To explain PEMFC equations and to assess their simulation effectiveness, the model given by Mann et al. [22] is chosen. This model is often used as a reference for the research analysis in more than one literature paper. The number of series cells (N cells ) is expressed by (1) [22][23][24] in the mathematical equation of terminal voltage from the PEM fuel cell (V stack ).
where E nernst is the Nernst voltage of a cell, v act is the overpotential of cell activation, v Ω is the ohmic voltage drop in the cell, and v conc is the over-potential of the concentration. If the temperature is below 100 • C, E nernst is calculated using Equation (2). The voltage drops v act , v Ω , and v conc can be calculated using (3)-(5), respectively [22,24]. where v where ρ m = 181.6 1 + 0.03 where T f c represents the PEMFC temperature in Kelvin (K), P H and P O are the hydrogen (H 2 ) pressure and oxygen (O 2 ) pressure, respectively, I fc is the fuel cell current (A), M A is the membrane surface area (cm 2 ), C O 2 is the concentration of O 2 (mol/cm 3 ), ξ 1 − ξ 4 are semi-empirical factors, R m denotes the membrane resistance (Ω), R c is the leads' resistances (Ω), l is the membrane thickness (cm), ρ m denotes the membrane resistivity (Ω.cm), λ is a tunable design variable, β represents a constant of an empirical range, J refers to the current density (A/cm 2 ), and J max is the maximum current density (A/cm 2 ). ξ 1 − ξ 4 are design variables which are established according to the geometrical dimensioning and the material resources of the manufactured fuel cells [23]. The fact that the design variable λ is dependent on different factors such as relative humidity, the stoichiometric ratio, makes its design a complex task. Hence, in this paper, the value of λ is assumed to be varied from 13 to 23. As λ becomes less, this indicates higher cell relative humidity. On the other hand, as λ becomes greater, it is an indication of the supersaturated conditions [7,22,23]. β is a design variable for Equation (5) which can be computed using Equation (6).
where is the ideal gas constant, F refers to the Faraday's constant, and α is the charge transfer coefficient.
It can be noted from Equations (5) and (6) that the drop in the concentration voltage follows a linear relationship with the temperature and the current density.
There are seven design variables (i.e., ξ 1 − ξ 4 , λ, Rc, and β) to be determined when constructing the model for the PEM fuel cells to utilize it in any further electric power system investigations.

Problem Formulation
The problem introduced in this article is solved using the newly developed CGO algorithm. The PEMFC model introduces nonlinear characteristics including multiple unknown design variables due to a shortage in manufacturing data about PEMFCs. Therefore, a detailed model is a challenge that must be addressed. To achieve the unique, dependable model, there are seven design variables to be tuned. These variables must be accurately calculated in order to minimize the FF. The FF is the sum of the square of the errors which is the difference between the terminal voltage calculated from the PEMFC and the terminal voltage observed by measurements as provided in Equation (7) [21,23]. The proposed fitness function is minimized by using the CGO algorithm. The solution of the optimization problem yields the design value for the aforementioned design variables: ξ 1 − ξ 4 , λ, Rc, and β [25].
where N samples refers to the number of the voltages measured, m defines the iteration number, V FC,exp is the PEMFC measured voltage (V), and V FC,est is the PEMFC computed voltage by the introduced model (V).

The Chaos Game Optimization Algorithm
Studying the use of chaotic systems attracted many researchers as the theme of chaos achieves a high interdisciplinary level [26]. Employing the chaotic systems in the optimization algorithms is an example, a new metaheuristic Chaos Game Optimization is introduced. The CGO algorithm's key concept is based on the chaos theory in which fractals are configured through the technique of the chaos game. The theory of chaos game is used as the key principle of the algorithm in the CGO method.

Inspiration
Chaos theory focuses on the basic characteristics of dynamically sensitive systems to the initial conditions. In view of the randomness of these complex processes, the principle of chaos refers to the presence of certain fundamental patterns.
The method of constructing fractals using the initial form and a random initial point is the Chaos Game. The main goal is to construct an iterative sequence of points to enter a sketch of similar form in various scales. Firstly, it is important to correctly place a triangle vertex. A random, initial point for the formation of fractal is chosen. The next item is then calculated as a part of the space between the first point and the random vertex. A fractal is generated while repeating this step with regard to the initial point and the vertex selected. Initially, the principal fractal form, which is triangular, is formed by three vertices. Dice with two, red, blue and green faces are used. The fractal, which is a seed, is chosen as a first random point. The first seed is pushed towards the corresponding vertex in the middle of the gap between the seed and the vertex while the dice are rolled, according to which color is selected. In the next iteration, the new location of the seed is used. The Sierpinski triangle is obtained by rolling the dice several times. The schematic explanation of such an approach is shown in Figure 1.

Mathematical Model
A number of solution candidates (X) representing the seeds is assumed by the CGO algorithm. Each candidate solution (X i ) consists of some tuneable variables (X i,j ) representing the seed's location. The triangle of Sierpinski is considered the exploration field looking for solutions. The following is the mathematical representation: where n represents the population size of the exploration field while d refers to the problem dimension. The initialization of the seeds is random as follows: where x j i (0) is the initial population, x j i,min and x j i,max refer to the min. and max. values for the jth design variable, and rand is a random number between 0 and 1.
To finalize the structure of a triangle, the mathematical model generates different seeds within their upper and lower limits. A temporary triangle with three seeds is drawn as follows:

•
The Global Best (GB) found so far; • The Mean Group (MGi) which is the mean values of random seeds; • The ith solution candidate (X i ).
The vertices of the triangle are GB, MGi and X i . A temporary triangle is generated for each initial seed to generate new seeds to complete a new triangle. The basic schemes of the formation of temporary triangles are shown in Figure 2.
The ith iteration contains the three corners of a Sierpinski triangle and the seeds obtained in the previous iteration. For new seeds, the temporary triangle is used. The dice are rolled. The seed in X i is transferred to the GB or to the MGi depending on the resulting color. This is modeled by a random integer function. It generates two integers, either 0 or 1, to allow green or red faces to be selected. According to the chaos game, seed movement should be minimal. Some generated factorials are used in this regard. The defined process for the first seed is presented mathematically as follows: where X i refers to the ith candidate, α i defines the random factorial which limits the seeds' movement, and β i and γ i are random integers of 0 or 1 representing the probability of rolling some dice. Dice of three blue and red faces are used for the second seed (GB). The dice are rolled, and the GB is shifted to X i or MGi, according to the resulting color. The seed moves towards the X i when a blue face is raised; if the red face is raised, the seed moves towards the MGi. The second seed goes towards a point on the connecting lines between X i and MGi. The following equation describes this mathematically: Regarding the third seed, the dice are rolled. According to the resulting color, the seed moves towards the X i or GB, which produces only two integrals as 0 and 1. According to some random factors, the seed will travel towards the connected lines between the X i and GB as expressed in Equation (12).
Another method is also used to create the fourth seed to introduce a mutation step in the location updates of the seeds. The location modification is based on certain changes in the selected decisions. The mathematical model of the process mentioned is: where k represents a random integer between 1 and d, and R represents random number ranges from 0 to 1.
Four different formulations for α i which regulate the motion limitations of the grain are presented to monitor and change the exploration and exploitation rate of the CGO algorithm as follows: where Rand denotes random number ranges from 0 to 1, and δ and ε refer to random factors between 0 and 1.
The fitness of new candidates is compared with the present ones and those of high fitness are kept, meanwhile, the candidates with lower fitness are rejected.
The pseudo-code of the CGO is expressed as follows:

Simulation Results
Three test cases of PEMFCs are investigated in this work. The considered PEMFCs are Ballard Mark V, AVISTA SR-12 500 W, and the 6 kW Nedstack PS6 PEMFCs. Experimental measurements are observed for the currents, voltages, and their characteristics. These data are used in the assessment of the efficacy of newly developed CGO by calculating the unknown seven PEM fuel cell design variables. The details of the upper and lower limits of the design variables for the PEMFCs are found in [23,[27][28][29]. Since the design variables can be calculated offline, the time needed for simulation is not the highest priority. However, the simulation time of the proposed CGO is also contrasted with the times of other approaches. These simulations are tested on a Microsoft Windows 10 Pro 64-bit version: 'R2014b' MATLAB platform (10.0, Build 18363). The MATLAB m-files are run on a Toshiba 2.00 GHz AMD A8-6410 CPU, fitted with 4 GB of RAM. In the following sections, the results of more detailed simulation are shown with respect to the SSE convergence and the I/V polarization curves for the three studied cases of PEMFCs.

Case One: The Ballard Mark V
The characteristics of the PEMFC used in this case are illustrated in Table 1. The first case fuel cell used in this experiment contains 35 cells connected in series. The current flowing through the fuel cell under test should not exceed 70 A due to thermal considerations.   Figure 4. It shows a comparison between the estimated voltage by the model and 13 samples of voltages observed experimentally. The blue curve, which is the calculated voltage by the proposed PEMFC model, best fits the 13 voltage readings, which are shown as red crosses, measured experimentally in the lab. The current vs voltage curves show that as the voltage of the fuel cell decreases upon the increase in the drawn current, the appropriate fitting between the estimated and observed voltages is clear. Additionally, the I/P relationship is shown in Figure 5. The blue curve, which is the calculated output power by the proposed PEMFC model, best fits the 13 power readings, which are shown as red crosses, calculated from the measured voltages and currents. As the current increases, more power can be provided by the fuel cell. It validates the exactness of the model and its capability to estimate very close results to those observed in the laboratory.    Table 2. It presents the best candidate solutions of the seven design variables (ξ 1 − ξ 4 , λ, R c , and β) obtained by each optimization method. This table also shows the settings of the optimization which are: the population size, the number of iterations, and the time taken by each optimization method in seconds to finish the simulation. Finally, the best result of the objective function achieved by each optimization method is shown in the last row. There are internal losses (voltage drop) in the voltage of the PEMFC such as V ohmic , V act , and V conc . The ohmic losses (V ohmic ) increases linearly upon the current increase. V act has the highest effect on the total voltage drop (losses). Meanwhile, V conc has the lowest effect on the total V losses . The plot of the relation between the internal voltage losses and the current is shown in Figure 6. The simulations are done more than once under different temperatures (30 • C, 50 • C, and 70 • C) at constant partial pressures (P H 2 /P O 2 = 1/1 atm) to investigate the effect of tuning the temperature on the voltage. This investigation and results are shown in Figure 7. It is noted that the PEMFC output voltage increases when the temperature rises at the same pressure. The difference between the output voltage values at different temperatures increases slightly at higher values of current. In Figure 8, the output power of the fuel cell is plotted against the current. It is observed that the effect of changing the temperature on the output power increases more at higher values of current. Meanwhile, at lower values of current the effect of temperature can hardly be observed. Besides studying the impact of temperature, the simulations are also done under the different pressures at constant temperature to measure the effect of tuning the pressure on the PEMFC voltage. The resulted I/V curves are shown in Figure 9 which show that the voltage increases when the pressure increases. In Figure 10, the output power of the fuel cell is plotted against the current. It is observed that the effect of changing the pressure on the output power increases at higher values of current. However, when comparing the effect of changing the temperature at constant pressure and vice versa, it is observed that the effect of the temperature on the output power is much higher. Therefore, changing the temperature leads to a wide range of output power from the fuel cell.

Case Two: The AVISTA SR-12 500 W
The test in the second case of this research study is conducted on 48 cells connected in series, AVISTA SR-12 500 W PEMFC. The specifications of the AVISTA SR-12 500 W PEM fuel cell are found in [25,28,[30][31][32]. Figure 11 shows the convergence curve of the objective function for 2000 iterations. Moreover, Figure 12 presents the I-V characteristics in which the voltage measured for 20 samples is compared with the voltage of the developed model. It can be clearly observed that both voltages are close in value. Regarding the fuel cell of case two, the blue curve of the calculated voltages follows the measured voltage-red marks-much better than that in the case one fuel cell. The relationship between the fuel cell current and power are provided in Figure 13. The results shown in these figures assure the robustness of the model. The values of the measured powers, marked in red, almost satisfy the blue curve of the estimated output power of the fuel cell.  Further numerical comparisons between the newly developed CGO method and other methods in extracting the design variables of the same PEMFC type of case two is given in Table 3. It provides the values of the seven design variables, the best solution, and the simulation time taken by each optimization method. The internal voltage losses of the PEMFC of case two are plotted against the fuel cell current in Figure 14. In contrast to the case one fuel cell, the effect of V conc and V act on the total voltage losses is more observable. V conc increases nonlinearly at higher values of the fuel cell current. The ohmic losses has the minimum share in the total voltage losses. The simulations are repeated under different values of temperatures (30 • C, 50 • C, and 70 • C) at constant partial pressures (P H 2 /P O 2 = 1.47/0.2 atm). The effects of changing the temperature on the output voltage and the output power are shown in Figures 15 and 16, respectively. It can be noted that the output voltage of the fuel cell increases upon the temperature rise.  The simulations are further repeated at different pressures under constant temperature. This is done to test the effect of changing the pressure on the voltage. The effects of changing the pressure on the output voltage and the output power are shown in Figures 17 and 18. The voltage increases upon the increase in pressure.

Case Three: The 6 kW Nedstack PS6
The third considered PEMFC in this study is the 6 kW Nedstack PS6 PEMFC. The specifications of this PEMFC are provided in Table 4. The test is performed on 65 cells connected in series. The current passing through an area of 1 cm 2 should not exceed 5 A. The dimensions of the fuel cells under test are shown in the second and third rows in the table.    Similar to the previous two cases, in Figure 19, the change in the values of the FF through the iterations can be seen. For further analysis, the I/V characteristics are provided in Figure 20. It shows a comparison between the estimated voltage and 29 samples of measured voltages. The fitting between the estimated and measured voltages is clear. The relationship between the current and the power is also provided in Figure 21.   Similar to results shown in Tables 2 and 3 of the first and second fuel cell cases, further comparisons between the newly developed CGO method versus WOA [29] and GHO [23] are presented in Table 5. They present the best candidate solutions of the design variables obtained by each optimization method, the best solution which corresponds to the minimum sum of the square of the errors between the calculated output voltages and the measured output voltages of the fuel cell, and the simulation time taken by each optimization method to finish the simulation and reach such results measured in seconds. The plot of the internal voltage losses is presented in Figure 22. As seen in this figure, V act has the maximum effect on the total voltage losses. V ohmic has the second highest share in the total voltage losses. It increases almost linearly with the fuel cell current. V conc has a very small value, and it is an almost constant value which is independent of the fuel cell current. The simulations are done more than once under different temperatures (30 • C, 50 • C, and 70 • C) at constant partial pressures (P H 2 /P O 2 = 1/1 atm). The results of the investigation on the effect of changing the temperature are shown in Figures 23 and 24. The effect is that the voltage increases when the temperature rises. The difference in output voltage increases at higher values of current. The simulations are also repeated under different pressures at constant temperature to measure the effect of tuning the pressure on the PEMFC voltage. The I/V curves are shown in Figures 25 and 26. From these two figures, it can be seen that the voltage increases when the pressure increases. After observing the results at different temperatures at constant pressure and vice versa, it is concluded that the effect of changing the temperature on the output voltage and the output power is higher than the effect of changing the pressure.

Robustness Statistical Analysis
Many runs were conducted to further check the robustness of the CGO. The robustness of an optimization algorithm means that many runs should be performed of the same problem to assure that this method can reach such a solution at every run, and that the best solution is not obtained by chance at a single run. For each case, 30 independent runs were performed. The results of case one are shown in detail in Table 6. The values of the standard deviation and the variance are low. This means that the results achieved by each run are very close to each other, indicating the robustness of the proposed optimization process.

Conclusions
This article proposed a novel attempt by the CGO to extract the design variables of the PEM fuel cell model. The goal is to achieve an effective PEMFC model that provides comprehensive simulation results for FCs. The FF is formulated to minimize the squaring error between the FC measured voltage and the calculated value. The CGO was introduced to minimize the FF under the problem constraints. The accuracy of the PEMFC model is proved under various temperature and pressure conditions. The efficiency of the constructed model is assessed by comparing the numerical model performance with the experimental results of the three commercial PEMFCs stacks tested. For all case studies, the simulation results coincided with the experimental outcomes. In comparison, the CGO results provided better results than the results based on alternative optimization approaches. This results in a high dominance in the literature of the CGO-based model over other models based on other optimization methods. The proposed CGO algorithm shall be used to solve several other engineering problems such as operation and control of microgrids and smart grids. In future work, the dynamic model of the PEMFCs can be considered instead of the static model. The aging factor of the fuel cell can be also considered.