OPTIMAL DESIGN OF OPERATION PARAMETERS FOR PEMFC

(1)

OPTIMAL DESIGN OF OPERATION PARAMETERS FOR PEMFC

Ying-Pin Chang^1,* Mu-Sheng Chiang²

1Department of Electrical and Information Engineering Nan Kai University of Technology

Nantou, Taiwan 542, R.O.C.

2Department of Mechanical Engineering Nan Kai University of Technology

Nantou, Taiwan 542, R.O.C.

Key Words: neural-network, orthogonal arrays, response surface methodology, PEMFC, operation parameters.

ABSTRACT

This paper presents a method for combining sequential neural-network approximation and orthogonal arrays (SNAOA) in determining the major operation and design parameters which affect the performance of proton exchange membrane fuel cells (PEMFC). An orthogonal array was first conducted to obtain the initial solution set. The results obtained from the orthogonal array were then used as the experimental data for response surface methodology (RSM) that could predict the operation parameters at discrete levels. The set was then treated as the initial training sample and a back- propagation sequential neural network was trained to simulate the feasible domain for seeking optimal operation parameters of PEMFC. With this method, the size of the training sample was greatly reduced due to the use of the orthogonal array. In addition, a restart strategy was also incorporated into the SNAOA so that the searching process could have a better opportunity to reach a near global optimum with the objective of reaching maximum output power of the PEMFC, which has a separate flow field in the cathode. The major parameters harnessed in this study include operating temperature, humidification temperature, reactant flow rate, split point, and split flow rate.

According to this novel methodology, the optimal parameters with a maximum power output were: operating temperature 78C, anode humidification temperature 72C, anode flow rate 296 sccm, cathode flow rate 295 sccm, split flow rate 145 sccm and split point 44%.

I. INTRODUCTION

A fuel cell (FC) is a device that converts the chemical energy from fuel into electricity through an electrochemical reaction with oxygen as the oxidizing agent. Hydrogen is the most common fuel, but hydrocarbons such as natural gas and alcohols such as methanol are sometimes used. An FC is different from a battery in that it requires a constant source of fuel and oxygen to run, but they can produce electricity

continuously as long as these inputs are supplied. In the archetypal design, a proton-conducting polymer membrane, the electrolyte, separates the anode and cathode sides. In a PEMFC stack, the cells are electrically connected in series and the polarization curves of the individual cells can be obtained by measuring the current of the entire stack and the voltages of individual cells [1].

For a given set of operating parameters such as system pressure, temperature, humidification condition, and gas stoi-

*Corresponding author: Ying-Pin Chang, e-mail: cyp@nkut.edu.tw

(2)

18 Journal of Technology, Vol. 33, No. 1

chiometry, the fuel cell voltage is essentially determined according to the operating current [2-7]. In [2], the Nexa power module is evaluated during membrane-electrode- assembly (MEA) and at the stack levels. The I-V Curves of the Nexa PEMFC system are measured through periodic current interruption to maintain the isothermal stack temperature and the uniformity analysis is mainly performed on the load of 800 W for all MEAs. Hydrogen for current polymer electronic membranes (PEM) and alkaline FCs must be supplied with no more than a few tens of ppm of CO or CO2, respectively [3]. If hydrogen is generated, as it is used, it must be produced efficiently over a broad power demand range, and the load changes on the order of seconds. This paper generated hydrogen for a broad variety of demands from 1.09/l molar mix of methanol/water using a commercial water-gas shift catalyst and a membrane reactor. In [4], Catalyst layers for the PEMFC were prepared by spraying and sputtering to deposit Pt amount of 0.1 and 0.01 mg cm², respectively. These a Pt layers were then assembled to fabricate an MEA having either single or double-layered catalysts. Experimental results indicate that Pt loading in state-of-the-art PEMFCs could be reduced by approximately 50% with no performance loss using both spraying and sputtering methods in the MEA fabrication process. In [5], the water transport in PEMFCs has been investigated by measurements of the effective or net drag coefficient. Results are presented for a wide range of operating conditions as well as for different types of membrane-electrode-assemblies. The understanding of current and temperature distributions along with the variation of gas composition in the cell of the PEMFC are crucial for designing cell components such as the flow field plate and the membrane-electrode assembly. The in- fluences of flooding and stoichiometry variation of the feed gas were discussed in terms of the rate of electrochemical reaction, from the measured distributions of local currents in a segmented single cell [6]. In [7], the spatial current den- sity distributions in a single PEMFC with three serpentine flow channels are measured using a segmented bipolar plate and printed circuit board technique. The effects of key operating conditions such as stoichiometry ratios, inlet humidity levels, cell pressure, and temperature on the local current den- sity distributions for co, counter, and cross-flow arrangements are examined.

Recently, numerical modeling and computer simulation have achieved a great progress for estimating the FC operation parameters [8-11]. In [8], a robust design analysis based on RSM was performed on a simplified fuel cell stack in order to identify the effect of assembly parameters on MEA pres-

sure distribution. The assembly pressure and bolt position were considered as randomly varying parameters with a given probabilistic property and acted as the design variables. In [9], a technique combining the use of a genetic algorithm neural networks model and the Taguchi method is employed to estimate the output voltage of PEMFC. In [10], a model of an electrochemical-based PEMFC mechanism suitable for engineering optimization is developed, and an effectively informed adaptive particle swarm optimization which is based on parameter-identification-techniques for this model is presented. In [11], artificial intelligence is utilized in optimiz- ing a hybrid PV/PEMFC energy system concerning the control and management. Also, simulation studies are carried out by using a real climate data and practical load profile.

With the fact that orthogonal arrays are widely employed in the experimental process for solving practical problems encountered, such arrays are geometrically balanced in their coverage of the experimental region with only a few representative experiments being required to be implemented.

A sliding-level orthogonal differential evolution algorithm with a two-level orthogonal array is proposed for solving worst- case tolerance design problems [12]. The report in [13]

has also presented a method of combined sequential neural network approximation and orthogonal arrays (SNAOA) for the planning of large-scale harmonic filters.

In [8], the planning of operation parameters on the performance of PEMFC with discrete variables in the power system belonging to the constrained combinatorial optimization problems is conducted. Such planning is difficult to solve by conventional methods, such as exhaustive searching.

Owing to the ANN’s excellent ability in dealing with the combinatorial optimization problems and the merit of orthogonal arrays that can systemically reduce the number of trails in the experimental process, in this study, a SNAOA is adopted to find the optimal design of discrete-value operation parameters in a PEMFC.

II. MODELING AND OPTIMIZATION APPROACH

Sequential approximation method:

1. Model Representation

A combinatorial optimization problem can be formulated as Eq. (1):

MaximumM X (1) ( )

subject to: g X_j( )0 , j !1, ,n_g and ^g^b ^X ¹ ^,

(3)

1, , _b b ! n ,

where M(X) is the objective function of discrete variable vector X, and X = [X1, X2, , Xj, , XD]^t. Each element of X belongs to an individual solution set, i.e.,

^,1^, ^,2^, ^,3^, ^, ^, ^, ^{1, 2, 3,} ^,

i i i i i mD

X X X X " X i ! D (2)

gj(X): inequality constraint, and gb(X): “pass-fail” binary constraints.

The feasible domain of the optimization model in Eq.

(1) can be simulated using a few representative data points to form the approximated model as Eq. (3):

MaximumM X (3) ( )

subject to: NN(X) = 1

where the binary constraint NN(X) = 1 approximates the feasible domain. If NN(X) = 1, the discrete point X is feasible; and if NN(X) = 0, the discrete point X is infeasible.

The discrete variable constraints are embedded in the format of the input nodes of the neural network. Throughout this report, the optimization model will be denoted Mreal, and the approximate model will be denoted MNN.

2. Response Surface Methodology

The methodology employed in this study aims to find the optimal settings for an automated RSM procedure when there is very little information about the objective function.

A framework of the RSM procedures for finding optimal solutions is presented which emphasize the use of both stop- ping rules and restart procedures. It is expected that the proposed settings would achieve considerable improvement.

A second-order model is commonly used for the multi- disciplinary design in RSM. In general, the response model can be expressed as follows:

2

1 1

k k

Y o i i iX i iiXi i j ijX Xi j

⁽⁴⁾

where Y is the response variable, X is the independent va- riable, is the reaction of the observed noise or error, and ^o^,ⁱ^,ⁱⁱ^,^ij are coefficients of the second-order regression equation. Fig. 1 shows the experimental flowchart using RSM

III. SYSTEM UNDER STUDY

A fuel cell (FC) is an electrochemical energy converter.

Constructan experimental system and the objective function

Set input A-F six factors within reasonable experimental range Using two-level factorial design

Apply to the path of steepest ascent

Evaluate the objective function of the optimal area.

Analyze the impact and effect of six factors

Conduct the second regression analysis

Find the optimal combination of factors, and then obtain the

best objective function Perform 27 sets of

central composite design

Use a model by RSM and find (Np-1) random initial vectors by SNAOA

Optimal area ? N Y

End

Fig. 1 Experimental flowchart using RSM

Load

H⁺ H⁺ e⁻ e⁻

o2 H₂O H2

Fuel in Oxidant in

Anode

Polymer electrolyte

Cathode

Depleted oxidant and product gases out Depleted fuel and

product gases out

H⁺

Fig. 2 PEMFC principle

It transforms the chemical energy of the fuel into electrical energy by two separated electrochemical reactions. In a hydrogen-fueled PEMFC, hydrogen is oxidized to protons and electrons at the anode. Protons migrate through the membrane electrolyte toward the cathode. As the membrane is an electric insulator, electrons are forced to flow through an external circuit. At the cathode, oxygen reacts with protons to produce water, which is the only waste product from a hydrogen-operated PEMFC. Fig. 2 shows the operating principle of PEMFC. In this single cell, a proton exchange membrane is used and the structure is divided into seven

(4)

Fig. 3 Performance test equipment of PEMFC

layers. The observation sequence is as follows: first observe the anode terminal, then the anode gas flow channel, the anode gas diffusion layer, and an anode catalyst layer, proton exchange membrane, the cathode catalyst layer, the cathode gas diffusion layer, and the cathode gas flow channel. Usu- ally, a researcher employs a five layer membrane electrode assembly MEA which contains an anode gas diffusion layer, anode catalyst layer, proton exchange membrane, cathode catalyst layer and cathode gas diffusion layer.

The control module of the performance test equipment of the PEMFC used in this study includes the following five components: temperature and humidifier control modules, gas supply equipment, and electronic load system and flow control system, as shown in Fig. 3. The experiment was conducted on a machine called the Beam^® test platform.

The electrical specifications are as follows: the maximum power is 600 W; the maximum current upper limit is 80 A, and the maximum voltage upper limit is 8 V. The test gas is hydrogen in the anode and the cathode has two different gas supplies: oxygen and air. In addition to these gases, nitrogen is used as the cleaning flow medium. This set of experimental test platform is well suited for the study of single cells and stacks. The BEAM^® test platform also connects to a computer and the test data will be stored into the com- puters through the digital transmission system.

The bipolar plate in the FC stack that acts as an anode for one cell and a cathode for the adjacent cell is shown in Fig. 4. The plate may be made of metal or conductive polymer (which may be a carbon-filled composite). It usually

Fig. 4 Bipolar plate in an FC

Fig. 5 FC flow channel plate with gas hole

incorporates flow channels for the fluid feeds and may also provide passage for heat transfer, as well as electrical con- duction. The plates used in this study have an additional gas hole in the flow channel plate of the fuel as shown in Fig. 5. That is, the flow channel of the bipolar panel opens another hole to allow flow not only from the inlet port for the reaction gas, but also from the middle of the channel input to reduce the concentration variation in different locations.

IV. OPTIMAL APPROACH OF OPERATION PARAMETERS ON PEMFC

Experiments were conducted at different temperature and flow rate setting with a split flow field in the PEMFC.

The prediction model of major operation parameters influ- encing the performance of cell reaction was also developed using SNAOA.

1. Objective Function

The objective of this study was to implement an operation parameter design for a PEMFC so that its output power can be significantly increased. The constraints of the operation parameters were also considered. The calculation of operation parameters can be formulated as a combinatorial optimization problem as follows:

1 2 3 4 5 6

Maximum ( )

Maximum ( )( , , , , , _o) M X

M X X X X X X X P

(5)

where M is the objective function of six factors for operation

(5)

Table 1 L27 (3⁶) orthogonal array Exp. X1 X2 X3 X4 X5 X6

1 1 1 1 1 1 1 2 1 1 1 1 2 2 3 1 1 1 1 3 3 4 1 2 2 2 1 1 5 1 2 2 2 2 2 6 1 2 2 2 3 3 7 1 3 3 3 1 1 8 1 3 3 3 2 2 9 1 3 3 3 3 3 10 2 1 2 3 1 2 11 2 1 2 3 2 3 12 2 1 2 3 3 1 13 2 2 3 1 1 2 14 2 2 3 1 2 3 15 2 2 3 1 3 1 16 2 3 1 2 1 2 17 2 3 1 2 2 3 18 2 3 1 2 3 1 19 3 1 3 2 1 3 20 3 1 3 2 2 1 21 3 1 3 2 3 2 22 3 2 1 3 1 3 23 3 2 1 3 2 1 24 3 2 1 3 3 2 25 3 3 2 1 1 3 26 3 3 2 1 2 1 27 3 3 2 1 3 2

parameters of the PEMFC, X1, X2, X3, X4, X5, and X6 de- notes the six factors of operation parameters, respectively and Po represents the output power of the PEMFC. Table 1 show that operation parameters of the PEMFC and their level values.

2. Constraints

i. Limits of FC Temperature (^C)

max min

1 1

X X X (6)

where X₁^maxand X₁^min are the upper and lower FC operation temperature respectively.

ii. Limits of FC Humidification Temperature (^oC)

max min

2 2 2

X X X (7)

where X₂^maxand X₂^min are the upper and lower FC humidification temperature respectively.

iii. Limits of Anode Flow Rate (sccm)

max min

3 3 3

X X X (8)

where X₃^max and X₃^min are the upper and lower anode flow rate limits respectively.

iv. Limits of Cathode Flow Rate (sccm)

max min

4 4 4

X X X (9)

where X₄^max and X₄^min are the upper and lower cathode flow rate limits respectively.

v. Limits of Cathode Split Flow Rate (sccm)

max min

5 5 5

X X X (10)

where X₅^max and X₅^min are the upper and lower cathode split flow rate limits respectively.

vi. Limits of Location of the Split Point (%)

max min

6 6 6

X X X (11)

where X₆^max and X₆^min are the upper and lower location of the split point limits respectively.

V. RESULTS AND DISCUSSION

The response surface plot of the simulation result is shown in Fig. 6. Fig. 6 (a) shows the residual plot of this method. It is used to detect the normal probability. The va- lue of the Y axis represents the percentage and the X axis is the residual. The Y axis data ranges between 99% and 1%, while the distribution of the X axis data is between 0.4.

A straight line of Y = 196X 50 can be used to mimic the variation with R² = 0.98, indicating the data has the charac- teristic of normal distribution. The histogram of the residual values is plotted in Fig. 6 (b). It is used to detect several peaks, outliers and non-normalities. The random variation frequency is at the range of 0-34 with the maximum value of 34 for residual 0. It is also shown that the frequency of

0.01 residual is about 10 while 0.02 and 0.03 is 1.

Therefore, it is close to the bell-shaped curve and represents normal distribution. Fig. 6 (c) shows the residual fit figure.

The Y axis is represented the value of residual, and the X

(6)

0.50 0.25 0.00 -0.25 -0.50 99 90

50

10 1

20 18 16 14 12 0.50 0.25 0.00 -0.25 -0.50

0.4 0.2 0.0 -0.2 -0.4 30

20

10

0

50 45 40 35 30 25 20 15 10 5 1 0.50 0.25 0.00 -0.25 -0.50

Residual Plots for P

Percent Residual

Frequency Residual

Normal Probability Plot Versus Fits

Histogram Versus Order

(a) Residual

(b) Residual

(c) Fitted Value

(d) Observation Order

Fig. 6 Statistical analysis of split flow setting using RSM on PEMFC: (a) linear plot; (b) frequency histogram;

(c) residual of fits scatter plot; (d) the residual line chart

Fig. 7 Toolbox main screen of the graphical user interface (GUI): Sequential Neural-Network Approximation/

Update neural-network search for the optimal point

axis is the FC output power P (W). This figure is to de- tect non-uniform variance peaks and missing higher-order terms with outliers P. The range of P is between 11 W-22 W. When the residual value is zero, P has its maximum and smallest value simultaneously. The resulting residuals should be randomly distributed in the range between 0.5.

Fig. 6 (d) also shows the residual line chart which is to detect the time- dependence of residuals. The Y axis represents the value of residuals and the X axis represents the obser- vation order. When the order of the observation order is 10, the residual has the maximum value of 0.4. On the con- trary, when the order of observation is the first, the residual has the minimum value of -0.42.

This paper employs the graphical user interface (GUI).

It not only can build powerful computing capabilities, but also makes it easier for users to make graphical represen- tations of the diversity and integration of all relevant processing programs into a single processing system. Fig. 7 shows the toolbox main screen of the GUI for sequential neural-network approximation and update neural-network searching of the optimal point. The following shows the parameter values used in this study: number of neu-

Table 2 Operation parameters of the PEMFC and their level values

Factor Parameter Level 1 Level 2 Level 3 X1 FC temperature (C) 50 65 80 X2 FC humidification temperature (C) 55 70 85 X3 anode flow rate (sccm) 200 250 300 X4 cathode flow rate (sccm) 200 250 300 X5 cathode split flow rate (sccm) 50 100 150 X6 location of the split point (%) 35 45 55

rons in the hidden layer = 12; net.trainparam.show = 100;

net.trainparam.epochs = 1000; net.trainparam.goal = 1e-6 and threshold value = 0.25.

1. Defining Set of Initial Training Data Using Tri-Level Orthogonal Arrays

In this study, there are 6 design variables, each with 31 discrete values. Thus, a total of 6 31 neurons are used in the input layer. Each neuron in the input layer has a value of 1 or 0 to represent the discrete value corresponding to each variable. There is only a single neuron in the output layer to represent the feasibility or infeasibility of this design point. The number of neurons in the hidden layer depends on the number of neurons in the input layer. There are 12 neurons in the hidden layer in this practical case illustration.

The transfer functions used in the hidden and output layer of the network are both log-sigmoid functions. The neuron in the output layer has a value range of [0, 1]. After the training is completed, a threshold value is applied to the output layer when simulating the boundary of the feasible domain. In other words, given a discrete design point in the search domain, the network always outputs 0 or 1 to indicate the feasibility or infeasibility of the point.

This training point is shown in Fig. 8, where an empty circle represents a “0” in the node, and a black circle represents a “1.” For illustrative purposes, the 186 input nodes are arranged in six rows to represent six design variables.

The first twenty-eight nodes of the first row have a value of 1 to represent that the first variable X1 has the discrete value 77. The first 17 nodes of the second row also have the value 1 to denote that the second variable X2 has the discrete value 71. The first 31 nodes of the third row also have the value 1 to denote that the third variable X3 has the discrete value 300. The first nodes of the fourth row also have the value 1 to denote that the fourth variable X4 has the discrete value 200. The first 31 nodes of the fifth row also

(7)

X1

55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 51 52

Input layer

Output layer

(infeasible) 53 54

50

X2

60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 56 57 58 59

55 85

X3

220 222 226 230 234 238 240 242 244 246 250 254258 262266 270 274 278 280 282 284 286 288 292 296 300 204 208 212 216

200

X4

X5

70 74 78 82 86 90 92 94 96 98 100 104108 112116 120 124 126 130 134 138 140 142 144 146 150 54 58 62 66

50

X6

37.5 38 38.5 39 39.5 40 41 42.5 43 44.5 45 45.5 46 46.5 47 47.5 48 50 51.5 52 52.5 53 53.5 54 54.5 35.5 36 36.5 37

35 55

220 222 226 230 234 238 240 242 244 246 250 254258 262266 270 274 278 280 282 284 286 288 292 296 300 204 208 212 216

200

Fig. 8 Representation of an infeasible training point [77, 71, 300, 200,150, 36]

X1

55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 51 52

Input layer

Output layer

(infeasible) 53 54

50

X2

60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 56 57 58 59

55 85

X3

220 222 226 230 234 238 240 242 244 246 250 254258 262266 270 274 278 280 282 284 286 288 292 296 300 204 208 212 216

200

X4

X5

70 74 78 82 86 90 92 94 96 98 100 104108 112116 120 124 126 130 134 138 140 142 144 146 150 54 58 62 66

50

X6

37.5 38 38.5 39 39.5 40 41 42.5 43 44.5 45 45.5 46 46.5 47 47.5 48 50 51.5 52 52.5 53 53.5 54 54.5 35.5 36 36.5 37

35 55

220 222 226 230 234 238 240 242 244 246 250 254258 262266 270 274 278 280 282 284 286 288 292 296 300 204 208 212 216

200

Fig. 9 Representation of a feasible training point [65, 70, 300, 200,100, 55]

have the value 1 to denote that the fifth variable X5 has the discrete value 150. The first 3 nodes of the sixth row also have the value 1 to denote that the sixth variable X6 has the discrete value 36. By checking with Table 2, one can de- termine that [X1, X2, X3, X4, X5, X6] = [77, 71, 300, 200, 150, 36] is an infeasible design point. Hence, the single node in

the output layer has the value of 1 to denote that this point is infeasible. Similarly, Fig. 9 shows another initial train- ing point [X1, X2, X3, X4, X5, X6] = [65, 70, 300, 200, 100, 55], which is a feasible design point. There must be at least one feasible design point in the set of initial training points, since the search algorithm described later must start from a

(8)

Table 3 Arrangement of orthogonal arrays

Exp. X₁ X₂ X₃ X₄ X₅ X₆

1 50 55 200 200 50 35 2 50 55 200 200 100 45 3 50 55 200 200 150 55 4 50 70 250 250 50 35 5 50 70 250 250 100 45 6 50 70 250 250 150 55 7 50 85 300 300 50 35 8 50 85 300 300 100 45 9 50 85 300 300 150 55 10 65 55 250 300 50 45 11 65 55 250 300 100 55 12 65 55 250 300 150 35 13 65 70 300 200 50 45 14 65 70 300 200 100 55 15 65 70 300 200 150 35 16 65 85 200 250 50 45 17 65 85 200 250 100 55 18 65 85 200 250 150 35 19 80 55 300 250 50 55 20 80 55 300 250 100 35 21 80 55 300 250 150 45 22 80 70 200 300 50 55 23 80 70 200 300 100 35 24 80 70 200 300 150 45 25 80 85 250 200 50 55 26 80 85 250 200 100 35 27 80 85 250 200 150 45

0 10 20 30 40 50 60

10^-7 10^-6 10^-5 10^-4 10^-3 10^-2 10^-1 10⁰

63 Epochs

Training-Blue Goal-Black

Performance is 6.24856e-007, Goal is 1e-006

Stop Training

Fig. 10 Training performance in terms of RMSE values

feasible design point.

In this study, the L27 (3⁶) array is used. Table 3 shows the arrangement of orthogonal arrays. The three levels are determined by both end points and one middle point of each variable within the design domain. Therefore the twenty- eight initial training points are shown. Among these points, [65, 70, 300, 200, 100, 55] for starting point and [80, 70, 200, 300, 100, 35] for the restart point are feasible to conduct the search of optimum values.

In order to perform a supervised training we need a way of evaluating the SNAOA output error between the actual and the expected output. A popular measure is the root mean squared error (RMSE) as shown in Fig. 10. It depicts training performance versus training epoch number in terms of

(9)

Table 4 Iteration history of operation parameters

Point X₁ X₂ X₃ X₄ X₅ X₆ P_o (W)

a 65 70 300 200 100 55 18.7 b 77 71 300 200 150 36 19.5 c 75 71 300 200 140 35 20.6 d 80 71 300 200 150 35 20.7 e 76 70 295 204 146 45 23.1 f 50 70 250 250 150 55 23.3 g 78 72 296 295 145 44 24.0 g1 80 70 200 300 100 35 19.6 h 78 72 296 295 145 44 24.0

M(W) ^OPTIMUM

POINT

RESTART POINT

0 1 2 3 4 5 6 7 G

24

23 22

21

20

19 a

b

c ^d

e f

g h

g1 ( :FEASIBLE POINT) ( :INFEASIBLEPOINT)

Fig. 11 Iteration history of study

RMSE values.

2. Iteration History of the Study

The algorithm described above is to search for the solution point by the L9 initial training points. There are two feasible design points [65, 70, 300, 200, 100, 55] and [80, 70, 200, 300, 100, 35] in the set of twenty-seven initial training points. The feasible design point [65, 70, 300, 200, 100, 55] has a lower objective value and is first used as the starting point in the search algorithm.

Fig. 11 shows the iteration history of the search, in which the points represented by black circles are the feasible design points. Table 4 also shows that the iteration history of operation parameters forms an “a” start point to an “h”

optimum point. The search algorithm is started from the

“a” point [65, 70, 300, 200, 100, 55] at an objective value of 18.7 W to seek a solution point at an initial feasible point and terminates at another new design point from the

“g” point [78, 72, 296, 295, 145, 44] at an objective value of 24 W. This search path proceeds from the “a” point [65, 70, 300, 200, 100, 55], through the “b” point [77, 71,

300, 200, 150, 36], “c” point [75, 71, 300, 200, 140, 35],

“d” point [80, 71, 300, 200, 150, 35], “e” point [76, 70, 295, 204, 146, 45], “f” point [50, 70, 250, 250, 150, 55] to the “g” point [78, 72, 296, 295, 145, 44].

To ensure a better chance for reaching a near global optimum, this search process is restarted from another feasible “g1” point [80, 70, 200, 300, 100, 35] in the set of initial training data. In this study, the search process terminates at the same feasible point [78, 72, 296, 295, 145, 44]

in one iteration.

Note that a total of 34 points from the 31 31 = 961 possible combinatorial combinations at an objective value of 24 W are evaluated to obtain this optimum design.

VI. CONCLUSIONS

This report presents a method of combined SNAOA and RSM for determining discrete-value of operation parameters which affect the output power of the PEMFC. Orthogo- nal arrays have been employed to obtain the local optimum solution and consume less time. The results obtained from the orthogonal arrays were then used as the experimental data for RSM that could predict the operation parameters at discrete levels. A back-propagation neural network has also been used to approximate the feasible domain of constraints. A clear 0-1 binary pattern is used in the input- output layers of the neural network. Therefore, the com- putational cost of the training of the neural network is small when compared to the evaluation of the constraints. The results show that the optimal operation temperature for the present PEMFC is 78C, anode humidification temperature 72C, anode flow rate 296 sccm, cathode flow rate 295 sccm, split flow rate 145 sccm, and split point 44%.

(10)

NOTATION

M(X) objective function of discrete variable vector X.

gj(X) inequality constraint

gb(X) “pass-fail” binary constraints

Po the output power of the PEMFC

LA(B^D) the orthogonal arrays

T temperature (C)

V fluid velocity (m/s)

YPT annual savings ($/year)

{X^o} the initial training data

0

Xh the lowest objective value

G

Xb the new solution point

o, i, ii, ij coefficients of the second-order regression equation

X1, X2, X3, X4, X5, X6 the six factors of operation parameters on the PEMFC

G

Ui a usable solution

X the independent variable

Y the response variable

the reaction of the observed noise or error

REFERENCES

1. Lee, S. K., G. Mogi, Z. Li, K. S. Hui, S. K. Lee, K. N.

Hui, S. Y. Park, Y. J. Ha, and J. W. Kim. 2011. “Meas- uring the Relative Efficiency of Hydrogen Energy Technologies for Implementing the Hydrogen Economy:

An Integrated Fuzzy AHP/DEA Approach.” International Journal of Hydrogen Energy 36 (20): 12655-12663.

doi: 10.1016/j.ijhydene.2011.06.135.

2. Zhu, W. H., R. U. Payne, D. R. Cahela, and B. J.

Tatarchuk. 2004. “Uniformity Analysis at MEA and Stack Levels for a Nexa PEM Fuel Cell System.”

Journal of Power Sources 128 (2): 231-238. doi: 10.1016/

j.jpowsour.2003.09.074.

3. Buxbaum, R., and H. Lei. 2003. “Power Output and Load Following in a Fuel Cell Fueled by Membrane Reactor Hydrogen.” Journal of Power Sources 123 (1):

43-47. doi: 10.1016/S0378-7753(03)00359-8.

4. Cho, Y. H., S. J. Yoo, Y. H. Cho, H. S. Park, I. S. Park, J. K. Lee, and Y. E. Sung. 2008. “Enhanced Perform- ance and Improved Interfacial Properties of Polymer

Electrolyte Membrane Fuel Cells Fabricated Using Sputter-deposited Pt Thin Layers.” Electrochimica Acta 53 (21): 6111-6116. doi: 10.1016/j.electacta.2007.12.008.

5. Janssen, G. J. M., and M. L. J. Overvelde. 2001. “Water Transport in the Proton-exchange-membrane Fuel Cell:

Measurements of The Effective Drag Coefficient.” Journal of Power Sources 101 (1): 117-125. doi: 10.1016/

S0378-7753(01)00708-X.

6. Yoon, Y. G., W. Y. Lee, T. H. Yang, G. G. Park, and C. S.

Kim. 2003. “Current Distribution in a Single Cell of PEMFC.” Journal of Power Sources 118 (1-2): 193-199.

doi: 10.1016/S0378-7753(03)00093-4.

7. Alaefour, I., G. Karimi, K. Jiao, and X. Li. 2012.

“Measurement of Current Distribution in a Proton Ex- change Membrane Fuel Cell with Various Flow Ar- rangements–A Parametric Study.” Applied Energy 93:

80-89. doi: 10.1016/j.apenergy.2011.05.033.

8. Liu, D., X. Lai, J. Ni, L. Peng, S. Lan, and Z. Lin. 2007.

“Robust Design of Assembly Parameters on Membrane Electrode Assembly Pressure Distribution.” Journal of Power Sources 172 (2): 760-767. doi: 10.1016/

j.jpowsour.2007.05.066.

9. Chang, K. Y. 2011. “The Optimal Design for PEMFC Modeling Based on Taguchi Method and Genetic Al- gorithm Neural Networks.” International Journal of Hydrogen Energy 36 (21): 13683-13694. doi: 10.1016/

j.ijhydene.2011.07.094.

10. Li, Q., W. Chen, Y. Wang, S. Liu, and J. Jia. 2011.

“Parameter Identification for PEM Fuel-Cell Mecha- nism Model Based on Effective Informed Adaptive Particle Swarm Optimization.” IEEE Transactions on Industrial Electronics, 58 (6): 2410-2419. doi: 10.1109/

TIE.2010.2060456.

11. Hattia, M., A. Meharrar, and M. Tioursi. 2011. “Power Management Strategy in the Alternative Energy Photo- voltaic/PEM Fuel Cell Hybrid System.” Renewable and Sustainable Energy Reviews 15 (9): 5104-5110. doi:

10.1016/j.rser.2011.07.046.

12. Tsai, J. T. 2012. “An Evolutionary Approach for Worst- case Tolerance Design.” Engineering Applications of Artificial Intelligence 25 (5): 917-925. doi: 10.1016/

j.engappai.2012.03.015.

13. Chang, Y. P., C. Low, and C. J. Wu. 2007. “Optimal

(11)

Design of Discrete-value Passive Harmonic Filters Using Sequential Neural-network Approximation and Orthogonal Array.” IEEE Transactions on Power Delivery 22 (3):

1813-1821. doi: 10.1109/TPWRD.2007.899780.

Manuscript Received: Jan. 20, 2015 First Revision Received: Jan. 22, 2015 Second Revision Received: Apr. 24, 2017 and Accepted: May. 02, 2017