Proton Exchange Membrane Fuel Cell Stack Design Optimization Using an Improved Jaya Algorithm

Fuel cell stack configuration optimization is known to be a problem that, in addition to presenting engineering challenges, is computationally hard. This paper presents an improved computational heuristic for solving the problem. The problem addressed in this paper is one of constrained optimization, where the goal is to seek optimal (or near-optimal) values of (i) the number of proton exchange membrane fuel cells (PEMFCs) to be connected in series to form a group, (ii) the number of such groups to be connected in parallel, and (iii) the cell area, such that the PEMFC assembly delivers the rated voltage at the rated power while the cost of building the assembly is as low as possible. Simulation results show that the proposed method outperforms four of the best-known methods in the literature. The improvement in performance afforded by the proposed algorithm is validated with statistical tests of significance.


Introduction
Energy has been identified as "humanity's number one problem for the next 50 years" [1]. Fuel cells [2][3][4][5][6], which are electrochemical devices that convert chemical energy to electrical energy, offer a viable alternative to fossil-fuel-based sources of energy. Of the various types of fuel cells (see, for example, references [6,7] for an overview), proton exchange membrane (or polymer electrolyte membrane) fuel cells, or PEMFCs, provide a relatively low-cost, low-temperature, high-efficiency and near-zero-pollution energy source for both stationary and portable applications. This paper uses a machine learning approach to address a PEM fuel cell problem of practical interest.
Fuel cells are complex systems whose electrochemistry, thermodynamics and engineering have not yet been fully understood. The complexity and non-linearity of the underlying processes in fuel cells (and many other areas of energy research) are often too hard for physics-, chemistry-, or engineering-based methods to build and interpret models that offer a balance between theoretical rigor and practical usefulness. It is here that data-driven or machine learning approaches come in handy, providing an alternative route to predictive analysis and parameter optimization.
The literature on the use of machine learning in energy research is vast; a few representative examples are mentioned below. A recent example of the use of deep learning (e.g., [8]) in nuclear fusion energy is found in [9] where a method has been developed for predicting disruptive instabilities in controlled fusion plasmas in magnetic-confinement tokamak reactors. Related work on the same type of problem has used machine learning strategies such as neural network [10][11][12], fuzzy logic and regression trees [13], support vector machine classification [14], and genetic algorithms [15]. Deep neural network is shown to outperform linear regression and (shallow) neural network for a short-term natural gas load forecasting application [16]. A problem of allocating optical links for connecting automatic circuit breakers in a utility power grid has been solved using a multi-objective

Theoretical Background
Part of this subsection is taken from [46,47]. The reversible thermodynamic potential or equilibrium voltage or open-circuit electromotive force (EMF) of the fuel cell is given by the Nernst equation, which is generally considered to be the cornerstone of fuel cell thermodynamics [2,46,47]: where E 0 is the reference (standard) EMF at unit activity and atmospheric pressure, i and j are the numbers of reactant and product species, a represents the activity, c i is the stoichiometric coefficient of species i, R is the universal gas constant, F is Faraday's constant, n is the number of electrons transferred for each molecule of the fuel participating in the reaction, and T is the temperature [46,47]. For a hydrogen-oxygen fuel cell (e.g., SOFC or PEMFC), the reactants are hydrogen and oxygen, and water (steam) is the product. (This paper uses many of the notations of [42,46]. A nomenclature is provided at the end of the paper.) The reference EMF, E 0 , depends on temperature T: where E 0 0 is the standard EMF at temperature T 0 , and ∆s is the change in entropy. The activity a of an ideal gas is given in terms of its pressure (or partial pressure) p: where p 0 is the standard-state pressure (1 atm). When the fuel cell is operated below 100 • C, so that liquid water is produced (as in proton exchange membrane fuel cells), the activity of water can be taken to be unity (a H 2 O = 1). In that case, the Nernst equation takes the form [46,47] where use has been made of the fact that n = 2 for a hydrogen fuel cell. The terminal (output) voltage is generally obtained by subtracting from E Nernst the following types of losses (or "irreversibilities") [46,47]: Ohmic loss • Fuel crossover and internal current loss Defining the current density, i den , as where A cell represents the cell active area, we can express the activation loss, caused by the slowness of the electrochemical reactions taking place on the surface of the electrode [6], as where α is the electron transfer coefficient, i den is the current density, and i 0,den the exchange current density (with i den > i 0,den , so that the logarithm is positive). Since the Tafel constant (Tafel slope) A [2] is given by this loss can be expressed as Concentration loss or mass transport loss results from the decrease in concentration of the reactants at the triple-phase-boundaries as the fuel is used: where B is a parametric coefficient (V), and i limit,den (> i den ) is the limiting current density. Ohmic loss is caused by the electrical resistance of the electrodes, the polymer membrane, and the conducting resistance between the membrane and the electrodes [2]: where r area is the area-specific resistance.
Losses also occur because of the leakage of fuel through the electrolyte and because of internal currents [48]. This type of loss is usually modeled by an additional component of current, called the fuel crossover and internal current.
As mentioned above, the terminal voltage of a single cell is obtained by subtracting the four losses from E Nernst . This simple, lumped, zero-dimensional [46] model leaves out the issues of, for example, parasitic power, aspect ratio, heat diffusion, water blockage, and flow channel effects.

Problem Statement
When a number of (identical) fuel cells are connected in series, the equivalent voltage of the group is the sum of the voltages of the individual cells, while the same amount of current flows through each cell. For a parallel connection of the cells, on the other hand, the total current is given by the sum of the individual currents, while the voltage of the group is the same as that of an individual cell.
Consider a group of N s cells connected in series and N p such groups connected in parallel. Let the entire assembly be called a stack. If i load,den represents the load current density for the entire assembly (stack), the stack terminal voltage is given by [42] V stack = N s E Nernst − A ln i load,den /N p + i n,den i 0,den + B ln 1 − i load,den /N p + i n,den i limit,den − (i load,den /N p + i n,den )r area , where i n,den is an additional component of current density brought into the equation to account for the combined effect of fuel crossover and internal current. The fuel crossover and internal current (the "leakage" or "lost" current) does not "flow" through the same path that the "regular" current follows [48], and therefore the ohmic loss for a cell should be modeled by the expression i load,den /N p × r area . For conformity with the treatment in Ref. [42] (this is a primary reference against which comparative results are studied later in this paper), however, the expression for V stack in Equation (9) retains the inexact inclusion [48] of i n,den in the calculation of the ohmic loss term. Table 1 shows the values of the cell parameters used in this problem. Given the cell parameters (except for the cell area) and the lower and upper bounds on the number of cells and also on the cell area, the problem is to produce an optimal stack design that minimizes the number of cells in series in each group, the number of such groups in parallel, and the cell area such that the rated load voltage at the maximum power point of the stack is 12 V and the maximum power is at least 200 W (corresponding to 730 kWh per year) [41,42]. These requirements come from a "research project aimed to design a power supply system to provide dc electricity to a single dwelling in a remote area of a developing country" [41]. As in references [41][42][43], the rated load voltage and the rated power come from the problem statement.

Previous Work on This Problem
A simple selection-crossover-mutation genetic algorithm was applied [41] to find (near-)optimal values of the three variables N s , N p and A cell . The linear ranking selection used in that work favored trial solutions with a low difference between the rated voltage (12 V) and the trial solution's maximum-power-point voltage (the sampling method accompanying the selection strategy was stochastic universal sampling). A penalty was applied to trial solutions not meeting the stipulated power of 200 W. The exact definitions of the objective function and the penalty function, however, are not provided in that paper.
A stochastic heuristic was developed in [42] to solve this problem. That method searched for the three variables N s , N p and A cell by minimizing an objective function that linearly combined several terms, favoring low values for all the three variables, in addition to favoring small absolute differences between the problem-specified rated voltage and the trial solution's output voltage at the maximum power point. An additive component of the objective function was a penalty term applicable to trial solutions producing a maximum power below the stipulated power (200 W). Two different penalty schemes were used: static (stationary) and dynamic (time-variant). For varying the relative importance of the different components of the objective function, different numerical constants were used, with values chosen heuristically. Unlike the genetic algorithm or any other member of the evolutionary computation family, that heuristic was point-based, not population-based.
Techniques from qualimetric engineering and extremal analysis were employed [43] to attack this problem by minimizing the total stack area N s N p A cell . Unlike the previous two methods (and unlike the present paper's approach), that approach used the data provided in the problem statement on the rated voltage and the rated power to simplify (reduce) the problem to one where the load current at the maximum power point was calculated simply as 200 W/12 V ≈ 16.667 A and was held fixed at this value for the rest of the computation (never revised to obtain the true current at the maximum power point). A second major simplification was achieved by using the Taguchi method [49] to set N p to its minimum possible value, namely unity. Finally, values of N s and A cell were determined by using a quasi-analytical approach involving a mix of numerical estimation, approximation and curve-fitting.
The authors of [41,42] obtained a trial solution's maximum-power-point voltage numerically from Equation (9), by iterating over current. That iterative computation was avoided in [43] through the use of the fixed value of 16.667 A as an estimate of the current corresponding to the maximum power point.

Cost Function
Following reference [42], the objective or cost function to be minimized is given by where V load,rated = 12 V represents the rated output terminal voltage of the PEMFC stack; V load,maxpp stands for the output voltage at the maximum power point of the stack; P load,rated = 200 W is the rated output power of the stack; P load,max is the maximum output power of the stack; K num = 0.5, K vdiff = 10, and K area = 0.001 are (heuristically chosen) positive constants [42] that allow us to vary the relative weights of the different components of the cost function; and the penalty term is zero if P load,max is at least P load,rated but positive otherwise. P load,max and V load,maxpp are obtained numerically from Equation (9), via iterations over current (see Section 4.2).
The cost function can be considered an anti-fitness function where the fitness function is traditionally maximized in the evolutionary computation literature. For the values in Tables 1 and 2, the cost function is positive. Fewer cells and smaller cell areas lead to lower costs. A larger deviation (in either direction) of the maximum-power-point voltage from the rated voltage makes the cost higher (worse). For a penalty-free solution vector (N p , N s , A cell ) (i.e., one for which P load,max ≥ P load,rated ), the above cost function is some measure of the monetary cost of building the stack. The present paper uses the adaptive penalty method (Section 3.6 of [42]) along with the coefficients/constants (Table 3 of [42]) of [42]. Table 2. Bounds of the design variables (from references [41,42]).

Variable Lower Bound Upper Bound
Number of cells in series in a group 1 50 Number of parallel groups 1 50 Cell area (cm 2 ) 10 400

The Improved Algorithm
The proposed algorithm is based on the general concept of evolutionary algorithms (e.g., [44]) and particularly on the Jaya algorithm [45]. The pseudocode of the original Jaya [45] is shown in Algorithm 1. The improved method modifies the Jaya algorithm in two ways. The first of these two modifications causes a major change in the algorithm design, while the second is relatively minor. The outline pseudocode of the improved algorithm is presented in Algorithm 2. The first modification introduces a new policy for updating (and using) the best and the worst individuals. In a population of size N, the individuals (trial solutions) can be thought of as occupying N consecutive positions (slots or locations) in an array or some other data structure. Each individual is a vector of three variables (parameters): an integer representing the number of cells in series in a group (N s ), a second integer for the number of parallel groups (N p ), and a floating-point number representing the cell area (A cell ). Population initialization is done by choosing values for each of the three problem parameters (N s , N p , A cell ) uniformly randomly from the interval defined by their respective lower and upper bounds ( Table 2).
A new individual x new is produced from the current individual x current at generation g (see Line 9 in Algorithm 2) as follows: for i = 1 to 3: where x i , i = 1 to 3, represent the three variables (N s , N p , A cell ) to be optimized; r g,i,1 and r g,i,2 are random numbers between 0.0 and 1.0; and x bestPosition and x last represent, respectively, the best and the worst individual in the population at the time of the creation of x new from x current . If x new i happens to fall outside its bounds, it is clamped at the appropriate (lower or upper) bound.
In the standard Jaya algorithm, the population's best individual and the worst individual are determined once at every new generation and are not updated during the course of a generation; that is, they are determined again at the next generation.
In the improved algorithm, however, whenever a new individual replaces an existing individual, we check whether the new individual is better (less costly) than the current best individual, and if it is, it (the new individual) becomes the current best, outsmarting the previous (most recent) best individual. By (conditionally) updating the best individual as soon a new individual appears in the population, our method is able to utilize the new (hitherto unavailable) information at the earliest possible opportunity, without having to wait for the entire (new) population to be created. This demonstrably leads to a better search for the (near-)optimal solution.
A similar strategy can in principle be used to continually update the worst individual in the population. That, however, turns out to be computationally expensive because finding the worst individual is an O(N) operation (where N is the population size), which would entail a total additional cost of O(N 2 ) per generation. We bypass this problem by putting, at each generation, the population's worst individual at a location (position in the population) that we know will be the last one to be accessed during the course of any given generation. Since, at any generation, the individuals in a population are processed sequentially, starting with the first and ending with the last, the last position is a sensible choice for holding the worst individual. At each generation, just before the beginning of the loop over all individuals, we swap the population's worst individual with the one at the last location (position) in the population.
The proposed algorithm's policy of continual update of the best and the worst individuals in the population renders the concept of the generation unnecessary or irrelevant. The proposed method thus implements a "steady-state" [50,51] mechanism of sorts. The reasons we still keep the generation loop are twofold: (i) we need to keep the worst individual at the last position; and (ii) the Jaya algorithm resets the random parameters (r g,i,1 , etc.) exactly once every generation, and therefore for a head-to-head comparison with Jaya, the proposed algorithm needs to do the same.
The second modification consists in altering the acceptance criterion for the new individual. The standard Jaya algorithm accepts the new individual (the child) only if it is better than the original individual (the parent). The proposed method accepts the new individual if it is better than or equal in cost to the original. This keeps the search process moving even on a flat landscape ( [31,52]).
Note that inside the for-loop spanning Lines 8-22 in Algorithm 2, we do not need to check if the "incoming individual" is worse than the current individual because the "incoming individual" is guaranteed to be at least as good as the one that it replaces. Note also that the if-statement on Lines 25-27 is required because, after the completion of the for-loop (on Lines 8-22), the best individual in the entire population may reside at the "last position" and thus take part in the swap operation (on Line 24). This is essentially the same reason why the finding of the best individual and the initialization of bestPosition to the location (position) of the best individual (Line 4) must follow, not precede, the swap statement on Line 3 of the pseudocode.

Computational Complexity
The improved method (Algorithm 2) incurs nominal extra computational cost compared to Jaya. The additional computation involves: For an entire generation, therefore, the total additional cost is no more than linear in the population size, or O(N).

Simulation Results
The relative effectiveness of five methods in optimizing the stack design are studied in this section. We present simulation results of head-to-head comparisons between the proposed algorithm and Jaya, between the proposed algorithm and the stochastic heuristic of [42], between the proposed algorithm and the genetic algorithm of [41], between the proposed algorithm and the quasi-analytical method of [43], and between Jaya and the heuristic of [42].

Proposed Algorithm vs. Jaya
The proposed algorithm and Jaya were coded in Java, and for each of the thirteen configurations in Table 3, 100 independent runs of either method were executed. The best, mean, and standard deviation of the best-of-run costs from 100 independent runs of either algorithm are presented in Table 4 for each of the thirteen configurations. Since this is a problem of minimization, numerically lower values of the cost function indicate better performance.
The proposed method outperforms Jaya in the majority of cases on each of the following three metrics: (a) the mean of the 100 best-of-run costs; (b) the best of the 100 best-of-run costs; and (c) the standard deviation of the hundred best-of-run costs.
Results of large-sample unpaired tests between the two means of the best-of-run costs in Table 4 are shown in Table 5, where the difference was obtained by subtracting the proposed method's cost from the corresponding cost produced by the Jaya algorithm. The sample size is 100 for either algorithm for each of the thirteen configurations. For all configurations but one, the proposed method is apparently better, judging by the sign of the z-score. The p-values (not shown in this table) corresponding to many of the z-scores are small, but it is arguable whether they are small enough for the difference to be statistically significant (much depends on the choice, arguably arbitrary, of the level of significance). For the lone negative z-score, the relatively high value of p implies that the proposed method is certainly not inferior to Jaya. Table 4. Comparison of the solution quality between Jaya and the proposed method (each row corresponds to 100 independent runs).  Table 5. Large-sample unpaired tests between the two means of best-of-run costs in Table 4. That each configuration in Table 3 represents a certain fixed set of algorithm parameter values implies that a performance comparison confined to any one configuration presents only a fragmented picture of the capabilities of the competing algorithms. For a more meaningful analysis, therefore, a configuration ensemble must be studied. This is what is done in Table 6 where Wilcoxon signed-rank tests are performed on the data of Table 4. The thirteen configurations yield as many independent samples of each algorithm's performance in Table 6 where two performance metrics are considered separately: the average of the 100 best-of-run costs and the best of the 100 best-of-run costs. Table 6. Wilcoxon signed rank tests on the data in Table 4 (sample size = 13).

Test Description Average of Best-of-Run Costs Best of Best-of-Run Costs
Null hypothesis Jaya mean-Proposed mean = 0 Jaya best-Proposed best = 0 Alternative hypothesis Jaya mean-Proposed mean > 0 Jaya best-Proposed best > 0 A comparison of the means of the best-of-run costs of the two methods (the "Mean" column under Jaya versus the "Mean" column under "Proposed Algorithm" in Table 4) is shown in the second column of Table 6 where a two-sample paired test was performed to test the null hypothesis Jaya mean -Proposed algorithm mean = 0 against the one-sided alternative Jaya mean > Proposed algorithm mean.
In this table, n denotes the effective sample size, obtained by ignoring the zero differences, if any. The W-statistic was obtained as min(W+, W−). The critical W was obtained from standard tables [53]. The z-score was obtained (using the normal distribution approximation) as where mean = n(n + 1) 4 , and std dev = n(n + 1)(2n + 1) 24 .
The p-value was calculated from standard tables of the unit normal distribution. The data in the second column of Table 6 show that, for the given sample size and the given significance level, the W-statistic (=13) is less than the critical W (=21), a fact that allows us to reject the null hypothesis in favor of the alternate hypothesis. Again, the p-value is less than 0.05; thus, the difference, Jaya mean-Proposed algorithm mean, is significant at the 5% significance level.
The last column of Table 6 shows the comparison between the two best values. The effective number of samples, n, is 12 in this case. Unlike the comparison in the second column of the table, the W statistic (=24) in the last column is larger than the corresponding critical W (=17), which means that the null hypothesis cannot be rejected at the 5% level. The z-score and the p-value, however, do not establish any superiority of Jaya over the proposed method at the 5% level.
The results in Tables 4-6 are all about the quality (cost) of the solution without any reference to the computation time needed to obtain that quality. In the next three tables, we present a comparative study of the time (as measured by the number of cost function evaluations) taken by the algorithms to obtain solutions of a given quality. Specifically, we set a cut-off or target cost of 13.5 (following [42]) and record, for each run of the algorithm, the number of cost evaluations needed to produce, for the very first time in the run, a solution better (that is, numerically smaller) than or equal to that target cost.
Let us denote by firstHitEvals the number of evaluations needed by an algorithm to hit (meet or beat) the target for the very first time in a given run (execution) of the algorithm. Given that the algorithms discussed here are really heuristics, once an algorithm hits the target for the very first time in a run, it is generally likely (but never guaranteed) that it will produce better solutions than the target value if the run is allowed to proceed further, that is, if the run consumes more evaluations than firstHitEvals. Note also that it is possible for an algorithm to never hit the target in a given (finite) number of evaluations.
Comparative results on the firstHitEvals metric are presented in Table 7 where each row was obtained from 100 independent runs of either algorithm. For each of the two competing algorithms, the #Success column shows the number of runs (out of the 100) that did hit the target within the maximum quota of maxEvals evaluations (recall that the maxEvals values for different configurations are given in Table 3). The proposed method's success rate is better than its competitor's in all thirteen cases. Again, the mean of firstHitEvals is smaller (better) for the proposed method in the majority of configurations; the pattern with the standard deviation of firstHitEvals is similar. Table 8 shows the results of large-sample unpaired tests on the two means of firstHitEvals in Table 7. Unlike Table 5, which has a sample size of 100 for either algorithm for each of the 13 configurations, Table 8 does not in general have the same sample size for the two algorithms for a given configuration, and the sample sizes are not always 100 (this is simply because not all runs in Table 7 were successful). For example, the sample sizes for the second configuration in Table 7 (and also in Table 8) are 34 and 42.
The use of the sample standard deviation in lieu of the population standard deviation in the calculation of the z-statistic in Table 8 is not inappropriate because the sample sizes, with the exception of the first configuration, are greater than 30, so that the central limit theorem may be used to justify the normal distribution approximation [54]. With a few exceptions, the p-values corresponding to the z-scores in Table 8 cannot be used to show the superiority (in a statistically significant way) of the proposed method for individual configurations. As discussed earlier, the overall algorithm behavior across different configurations can be captured by the configuration ensemble. This is done in Table 9 where Wilcoxon signed rank tests conducted on the data in Table 7 are presented. The results in Table 9 show that not only does the proposed algorithm succeed in hitting the target significantly more often than its competitor, it also needs significantly fewer evaluations to achieve that feat. Table 9. Wilcoxon signed rank tests on the data in Table 7 (sample size = 13).

Test Description Average of firstHitEvals Number of Successful Runs
Null hypothesis Jaya mean-Proposed mean = 0 Jaya successes-Proposed successes = 0 Alternative hypothesis Jaya mean-Proposed mean > 0 Jaya successes-Proposed successes < 0 Significance level α 0.05 0.05 No. of zero differences 0 5 Ignore zero diff. Include zero diff. The next three tables show the performance of the algorithms on the metric of firstHitCost, which is the cost of the solution obtained at firstHitEvals; clearly, firstHitCost is less than or equal to the target cost. The mean and standard deviation of firstHitCost values obtained from 100 independent runs (these runs are the same as the ones used to obtain firstHitEvals in Table 7) are shown for each of the two competing methods in Table 10. Table 10. Mean and standard deviation of firstHitCost (each row is obtained from the same set of 100 independent runs used in Table 7 for each algorithm).  Table 11 shows results of large-sample unpaired tests on the two means of firstHitCost in Table 10. As before, the one-tailed p-values (not shown in Table 11) produced by the z-scores fail to establish a statistically significant outcome with respect to individual configurations, a fact that leads us to the statistical analysis of the configuration ensemble in Table 12, which presents the results of the Wilcoxon signed rank test on the set of 13 samples in Table 10. Table 12 shows that the average cost of the solution produced after consuming firstHitEvals evaluations is statistically significantly better (at the 5% level) for the proposed method. Table 11. Large-sample unpaired tests on the two means of firstHitCost in Table 10. The best stack designs produced by Jaya and the proposed algorithm, given by the solution vectors (N s , N p , A cell ) and the corresponding P load,max , V load,maxpp and cost, are shown in Table 13 where each row corresponds to a given configuration (popSize-maxGen combination). A solution in this table represents, for the corresponding algorithm and the corresponding configuration, the best of all the solutions produced in the 100-run suite. While the proposed method's best solution produces a better (lower) cost in the majority of cases in Table 13, neither algorithm is statistically significantly better than the other on this particular metric (as seen in the results of Table 6). The mean best solution of the proposed method, however, was found to be statistically significantly better than that of Jaya (Table 6).

Effect of Current
Step Size on Numerical Calculations of P load,max and V load,maxpp The maximum power point corresponding to a given stack configuration vector (N s , N p , A cell ) and a given set of cell-parameter values (i 0,den , etc.) cannot be obtained analytically. In other words, the problem of finding argmax with V stack given by Equation (9), is not analytically solvable. A numerical, iterative method has been used in this paper, where a loop over current values is executed in order to numerically find an approximation to the maximum power point. Since the maximum power point needs to be determined for every single cost (objective function) evaluation, the total computation time to complete all such iterations for all runs in a test suite may not be trivial. It is, therefore, good to be able to use a relatively large step size in incrementing the current value in the loop. A large step size, however, reduces the accuracy of the computed P load,max and V load,maxpp . The results presented so far in this paper used a current step of 50 mA. To study whether the step size affects the conclusions about the relative performance of the algorithms, we obtain a new set of results using a step size of 1 mA. These results, arguably more accurate than their 50 mA counterparts, are presented in Tables 14-16.
A target cost of 13.62 is used in Table 15 because the previously used target of 13.5 was met by no run in the new set of results. Results of Wilcoxon signed-rank tests on the new set of data are shown in Table 17 which establishes the proposed method as significantly better than Jaya on the success rate as well as on the mean of best-of-run costs. Jaya is significantly better on the best of best-of-run costs metric, while the methods are statistically tied (the one-tail p-value is close to 0.5) on the remaining two metrics. Given the vagaries of chance, for stochastic heuristics, the average performance rather than the single best performance is generally considered to be a reliable indicator of a method's performance.
The new data, therefore, do not offer a convincing reason to argue that Jaya performs better than the proposed method. Combining the new results with the old ones for all five performance metrics, we conclude that the proposed heuristic meets or beats the Jaya algorithm.  Tables 14-16.

Proposed Algorithm vs. Point-Based Stochastic Heuristic
We now present performance comparisons of the proposed method against the heuristic of [42]. Each of the four maxEvals values reported in Table 9 of [42] corresponds to two different configurations in Table 3 of the present paper (see Table 18). The best-cost solution vector from Table 4 on p. 536 of reference [42] is shown in Table 19 where the values of N s , N p and A cell are copied from [42], P load,max and V load,maxpp are computed by the iterative numerical method described in Section 4.2 (using two different values of the current step size), and the cost is computed from Equation (10). The first row of Table 19 shows that the current step of 50 mA produces a close agreement of P load,max , V load,maxpp and cost values in this table with the corresponding values in Table 4 of [42]. For the rest of this subsection, therefore, results corresponding to only the 50 mA step size are used. Table 9 of [42] and Tables 4 and 10 of the present paper show that the proposed algorithm outperforms the algorithm of [42] in all of the following metrics: (i) best of the 100 best-of-run costs; (ii) mean of the 100 best-of-run costs; (iii) standard deviation of the 100 best-of-run costs; (iv) count (out of the one hundred) of successful runs; and (v) mean firstHitCost obtained from the successful runs (the mean firstHitCost of [42] is slightly better for one of the two cases-Configuration 11-of maxEvals = 4000).
To investigate the statistical significance, if any, of the difference between the means of the best-of-run costs of the two algorithms, we performed unpaired t-tests, the results of which are presented in Table 20. We tested the null hypothesis µ 1 − µ 2 = 0 against the one-sided alternative µ 1 − µ 2 > 0, where µ 1 and µ 2 are the (population) means of the method in [42] and the proposed method, respectively. We chose a level of significance of 0.05. Since the sample variances differ by a factor of about 10,000, the standard two-sample t-test cannot be used. We, therefore, used the Smith-Satterthwaite test [32,54] corresponding to unequal variances. The test statistic is given by and its sampling distribution can be approximated by the t-distribution with s 2 1 n 1 + s 2 2 n 2 2 (s 2 1 /n 1 ) 2 n 1 − 1 + (s 2 2 /n 2 ) 2 n 2 − 1 degrees of freedom (rounded down to the nearest integer), wherex 1 andx 2 represent the two sample means, s 1 and s 2 are the two sample standard deviations, and n 1 and n 2 are the two sample sizes (100 each). Table 20 presents, for each configuration, the following: the critical t value (obtained from standard tables of t-distribution) at 95% (right-tail probability of 0.05) and for the specified degrees of freedom; and • 95% confidence interval for the difference of the two means,  The results show that, for each case, the t-statistic exceeds t 0.05 for the relevant degrees of freedom. Therefore, the null hypothesis is rejected in all of these cases. Furthermore, none of the confidence intervals contains zero. We thus conclude that the improvements produced by the proposed method are statistically significant.
Smith-Satterthwaite tests were also performed for comparing [42] against Jaya, and the results (see Table 21) establish Jaya as significantly better than [42].

Proposed Algorithm vs. Genetic Algorithm
The best solution produced by the genetic algorithm in [41] recommends the following stack design: N s = 21, N p = 1, A cell = 156.25 cm 2 . Plugged into Equations (9) and (10), this design yields the maximum power point value and the corresponding cost shown in Table 22 where two sets of calculations (for the two step sizes for current) were used. The proposed method's best solutions outperform both cases in Table 22 not only on the objective function (cost) but also on the output voltage at the maximum power point: the V load,maxpp values in Table 22 are less than the rated voltage of 12 V.

Proposed Algorithm vs. Quasi-Analytical Approach
The best solution produced by the quasi-analytical method in [43] yields the following stack design: N s = 22, N p = 1, A cell = 151.4 cm 2 . Plugging these values into Equations (9) and (10) gives the maximum power point values and the corresponding costs in Table 23. All of the best solutions of the proposed method in Tables 13 and 16 have a better (lower) cost than those in Table 23. The objective function minimized in [43] is the total stack area given by N s N p A cell , and yet the proposed method's best solutions produce better (smaller) total areas, while meeting the requirements of the rated voltage and rated power.
The characterization in [43] of the design vector (22, 1, 154.16) as the best design in [42] does not seem to be correct. The vector (22, 1, 154.16) is mentioned in [42] as a "typical" design, not the best or optimal design. In fact, several of the designs in Table 4 of [42] have cell areas that are smaller than 154.16 cm 2 , with two of them even smaller than the "optimal" area of 151.4 cm 2 in [43]. The design (22,1,149.597) from Ref. [42] outperforms the "optimal" design of (22, 1, 151.4) of [43] on the cell area metric as well as on the objective (cost) function for both the step sizes (see Tables 19 and 23).  Figure 1 shows the polarization curve for the stack design produced by the best of best-of-run solutions of the proposed algorithm for Configuration 13 (the last row of Table 16). The corresponding power characteristics are plotted in Figure 2. The nature of these plots agrees well with that of polarization and power characteristics of typical PEM fuel cells in the literature.

Conclusions
A PEMFC stack design problem of practical importance has been addressed in this paper. This constrained optimization problem requires finding optimal values of three parameters of the stack configuration (namely, the number of cells in series, the number of groups in parallel, and the cell area) such that the configuration delivers the rated voltage at the rated power (the load voltage at the maximum power point is to be 12 V and the maximum power is to be at least 200 W), while keeping the total cost at a minimum. A new approach has been developed, based on the Jaya algorithm and ideas from evolutionary computation, and the new approach has been empirically (that is, via numerical simulation) compared with the Jaya algorithm and with the methods in [41][42][43], using the following performance measures: (a) the best of best-of-run costs from 100 independent runs; (b) the average of the 100 best-of-run costs; (c) the standard deviation of the 100 best-of-run costs; (d) the number of cost function evaluations needed to reach a pre-determined target cost for the very first time in a run of the algorithm; and (e) the cost of the very first (earliest) solution in a run that meets (or beats) the target cost. The proposed method has been shown to outperform the three methods in [41][42][43] and is competitive with Jaya. The improvement in performance provided by the proposed algorithm has been substantiated with statistical tests of significance.
Following the authors of [41][42][43], this paper considered only a stand-alone PEM fuel cell stack. Auxiliary systems, such as heat exchangers, compressors, etc. that often significantly affect the mass, size and cost of the entire system, were not considered; the membrane type, the type of the bipolar plate material, and many operating and constructive parameters were not considered, either. Inclusion of some of these issues would be in the agenda for future research.