Adjoint Method for the Optimization of the Catalyst Distribution in Proton Exchange Membrane Fuel Cells

In this article we present an adjoint method for the optimization of the catalyst distribution in proton exchange membrane fuel cells (PEMFCs). By using the theory of functional analysis we derive analytical equations for the sensitivity functions of the cell voltage with respect to the catalyst distribution in a very general framework, independent on the transport model used to simulate the PEMFC. Then we present an efﬁcient numerical algorithm to calculate the sensitivity functions using the adjoint method. The adjoint method has the advantage that it can be applied to the optimization of systems with a large ( > 10 4 ) number of optimization variables that are computed simultaneously and independently to maximize the objective function. Finally, we apply the method to the optimization of 2-D platinum distribution in PEMFCs. We show that the optimum platinum distribution varies with the operating conditions, position of landings and openings, cell geometry

Proton exchange membrane fuel cells (PEMFCs) have attracted the attention of the research community because of their high power density, energy efficiency, and environmentally friendly characteristics. Operating at low temperatures, PEMFCs are suitable for a variety of applications ranging from automotive applications to stationary and portable applications. In the area of automotive applications, PEM-FCs have already been introduced commercially or for demonstration purposes by all the major car manufacturing companies; in addition, they are increasingly being used in busses, motorcycles, boats, submarines, and airplanes. 1,2 In the area of stationary and portable applications, PEMFCs are employed as emergency power systems, such as uninterrupted power supplies, and have the potential to be used for energy storage in electric grid systems. 3 Since the technical feasibility of PEMFCs has already been established, a significant amount of current research focuses on increasing the durability and decreasing the total manufacturing cost, which, at large manufacturing volumes, is mainly given by the cost of the catalyst layers (CLs), bipolar plates, and membrane. [4][5][6][7][8] (For instance, the US Department of Energy estimates that 45% of the cost of fuel cell stack fabricated at a volume of 500,000 systems/year is due to the cost of platinum, 27% to bipolar plates, 10% to the cost of membrane 9 ). In this article, we address the problem of decreasing the manufacturing cost of the CLs by presenting a large-scale optimization technique to optimize the platinum deposition in the CLs of PEMFCs. The technique can also be applied to the optimization of other two-dimensional (2-D) and three-dimensional (3-D) field variables, such as porosity distribution and ionomer content, in which the number of optimization parameters is very large (e.g. more than 10 4 ).
There is much work on the numerical optimization of catalyst layers in fuel cells. The existing optimization techniques can be divided in heuristic and non-heuristic techniques. To search the optimization space, heuristic techniques usually use stochastic approaches. Examples of such techniques include the genetic algorithm, particle swarm optimization, differential evolution, simulating annealing, harmony search, memetic algorithms, and bird mating optimization. [10][11][12][13][14] Heuristic techniques are regarded as computationally too expensive for the optimization of two-dimensional (2-D) or three-dimensional * Electrochemical Society Member. z E-mail: pandrei@eng.fsu.edu (3-D) catalyst distributions in finite element simulations. For this reason, heuristic techniques can be applied to optimize only a small number of variables (called design variables) and often need to be used with simplified models, such as compact, reduced-order models, or equivalent circuit models. For example, Ebrahimi et al. 15 assume that the catalyst density varies only along the x-direction (along the length of the cell) and according to a predefined polynomial function, with unknown coefficients. In this way the number of design variables is reduced to the number of coefficients of the polynomial function and the optimization problem is solved using the Nelder-Mead Simplex technique. Song et al. 16 optimize the catalyst distribution and the Nafion content in the CL by partitioning the CL into several layers along the length of the cell and treating the platinum loading and Nafion content constant in each layer. In this way, the number of design variables is significantly reduced again and is equal to the number of layers in which the CL is partitioned.
Non-heuristic techniques are deterministic and often require evaluating the gradient of the objective function. [17][18][19][20][21] In general, these techniques have local convergence and their rate of convergence can be increased, using second-order derivatives such as the Hessian matrix, or taking advantage of the particular functional form of the objective function. Depending on whether the gradient function can be evaluated of not, as well as on the accuracy with which the first and the higher order derivatives can be computed, non-heuristic techniques usually require between a few and 100 iterations to optimize the objective function. Due to this low number of iterations, non-heuristic techniques are almost always preferred in optimizing PEMFCs described by finite element models. For example, Mathieu-Potvin and Gusselin 22 optimize the catalyst distribution by dividing the CL into 5 by 4 rectangular regions and interpolating the catalyst distribution using Lagrange interpolation. The 20 independent variables are then optimized using a gradient-approach. Since the sensitivity of the objective function to the catalyst loading in each region and the Hessian matrix are computed independently the algorithm requires around 1,400 device simulations (the relatively large number of devices simulations is due to the fact that the sensitivity of the objective function to the catalyst loading in each region is evaluated independently, and not by using the adjoint method). Secanell et al. 23 use a gradient-based approach to optimize the catalyst loading and three other cell parameters. Since the sensitivities of the objective function with respect to the optimization variables are computed analytically and evaluated E3233 numerically independently, the optimization study was applied to only four design variables.
In addition to the above numerical optimization studies in which the catalyst loading was assumed nonuniform and its optimal functional form was computed using finite element or finite difference models,  derived analytical expressions for the optimum catalyst deposition along and across the CL by solving the oxygen diffusion equation in the cell analytically. Although his approach requires many approximations related to the water vapor diffusion, electron and proton conductivity, and neglects crossover effects, the results obtained are significantly important because they predict the optimum functional form of the catalyst distribution analytically and in qualitative agreement with the existing finite element simulations.
Finally, it is also worth mentioning that PEMFCs with nonuniform catalyst loadings have also been studied experimentally by many groups. In general, the experimental work has focus on introducing a gradient of the platinum lading in the cathode CL in the direction along the length of the cell usually by introducing adjacent layers with different platinum loadings. 27 In this article we use the theory of functional analysis to derive a formula for the sensitivity of cell voltage and current density to the catalyst distribution and present a gradient-based technique for the optimization of the catalyst profile in PEMFCs using finite elements. The optimization is performed by allowing the catalyst amount to vary independently from location to location inside the CL, i.e. the amounts of catalyst at each node of the finite element discretization are independent variables. In this way, the number of optimization variables is equal to the total number of nodes in the catalyst layers, which, in this work, is larger than 10 4 . In order to compute the optimum catalyst distribution in the CL, we define the catalyst sensitivity functions of the cell voltage and current density. To make the method computationally tractable, the catalyst sensitivity functions are evaluated using the adjoint method, which is appropriate for large-scale systems.
It should be noted that adjoint methods have been used before independently by the applied mathematics community to solve 1-D problems in fluid dynamics, climate, and heat transfer, 28,29 by the aerodynamic community to perform sensitivity analysis studies, [30][31][32][33] and by our group to estimate parameter variability and optimize the doping profiles in 2-D and 3-D semiconductor devices. [34][35][36][37][38] In the area of electrochemical energy storage, Carraro et al. 39 and Kapadia et al. 40 were the first to use adjoint method to study the sensitivity of cell voltage with respect to cell parameters for solid oxide fuel cells, however, they limited their analysis to only a few parameters such as average porosity, reaction rate, tortuosity, and mean pore radius, and have not performed a large-scale sensitivity analysis of the cells. It was only recently that a widely spread commercial simulation software, COMSOL, 41 has introduced the adjoint method for the calculation of sensitivities of platinum distribution in PEMFCs and the adjoint method become more popular for optimization purposes. The first to use the adjoint method to optimize the 2-D catalyst distribution was Cetinbas et al. 42,43 who used COMSOL to simulate the cathode and calculate the sensitivities of the cell current with respect to the catalyst and ionomer distributions (their finite element simulation did not include the membrane and anode). Then, they used the "fmincon" solver in Matlab to optimize the catalyst iteratively and obtain the bidirectionally-graded catalyst and ionomer profiles. In contrast to the above work, this article derives analytical expressions for the catalyst sensitivity functions of cell voltage and introduces an efficient numerical implementation for the computation of these sensitivities in a more general framework that allows the optimization of multiple field variables under constraints. In addition, the adjoint method presented in our work also uses an improved optimization technique and can be implemented numerically together with the solver and optimizer in stand-alone finite element simulators for fuel cells. 44 Finally, we simulate the entire cell including the membrane and the anode and cathode CLs and gas diffusion layers (GDLs) in order to better capture the drift, diffusion, and kinetics of hydrogen gas and protons.
The remainder of this article is organized as follows. In the next section we define the catalyst sensitivity functions using the formalism of variational analysis and present the technique that we propose for the evaluation of these functions. Since this article focuses on the mathematical derivation and numerical computation of the sensitivity functions, the presentation is made without considering any particular form of the transport and kinetic models in the PEMFCs. Then, we present the numerical algorithm for the optimization of the catalyst distribution and describe its numerical implementation in the case of finite element simulations. Finally, we present sample numerical examples and draw conclusions in the last section.

Catalyst Sensitivity Functions
Definition of the catalyst sensitivity functions.-The optimization algorithm presented in this work is based on the evaluation of the catalyst sensitivity functions, which describe the sensitivity of the objective function with respect to the amount of catalyst at each location inside the cell. Since the objective function is either the discharge current, I, or the cell voltage, V, we distinguish two catalyst sensitivity functions: (a) the catalyst sensitivity function of the cell voltage, γ V (r), and (b) the catalyst sensitivity function of the cell current, γ I (r). The catalyst sensitivity function of cell voltage, γ V (r) is defined as how much the cell voltage is changing if we add one unit amount of platinum (usually defined as 1 μg of Pt) at location r inside the catalyst layer. The catalyst sensitivity function of cell current, γ I (r) is defined as how much the cell current is changing if we add one unit amount of platinum at location r. The units of the catalyst sensitivity functions are V/g Pt in the case of γ V (r) and A/g Pt in the case of γ I (r).
To define the catalyst sensitivity functions mathematically, let us denote the catalyst density, defined as the mass of the catalyst per unit volume of the CL, by μ Pt (r). The total mass of the catalyst (measured in grams) can be computed by integration [1] where denotes the total region occupied by two CLs. Mathematically, the catalyst sensitivity functions are related to the Gâteaux derivatives of the terminal voltages and currents in the direction of the catalyst variation. It can be shown that, by modifying the catalyst distribution by an infinitesimally small amount δμ Pt (r), the cell voltage changes under constant current by [2] If the measurements are made under constant cell voltage, the cell current changes by Equations 2 and 3 can be derived using the Riesz representation theorem but, in this work, they will be considered as definitions of the catalyst sensitivity functions. These equations can be interpreted geometrically easily by looking at the polarization curves for two "infinitesimally" close platinum distributions, as shown in Fig. 1.
In the previous equations we have explicitly indicated that the mass of the catalyst per unit volume of the cathode are functions of r, however, V [μ Pt (r)] and I [μ Pt (r)] should be regarded as functionals of μ Pt and not as functions of r.
Since the catalyst sensitivity functions are related to the gradients of the objective function with respect to the amount of catalyst, these functions are instrumental for the design and optimization of fuel cells. Next, we introduce the finite element model that we use to simulate the PEMFCs in a very general form (see Eq. 4) and describe the technique for computing the sensitivity functions. One should notice that the numerical algorithm for the computation of γ I (r) and γ V (r) requires the solution of only one sparse system of linear equations and can be applied to any general transport model that can be discretized numerically as a nonlinear system of equations.
Transport and kinetic equations.-In this work we assume that, for any given cell current I, the transport of oxygen, hydrogen, protons, liquid water, and water vapors is governed by a system of nonlinear equations that can be represented in vector form as where x is the vector of state variables (state vector) and μ Pt is a vector whose components are the values of the catalyst density at each mesh point. In general, system (4) is obtained by writing the transport equations on a finite element grid and expressing the equation at each mesh point (i.e. node or vertex) of the grid in discretized form. This procedure is relatively straightforward in the case of models that can be written as a system of partial differential equations (such as the fluid dynamics equations) but it can become more complicated or sometimes even impossible in the case of models based on particle tracking or other iterative methods. Therefore, the optimization method presented in this work can be applied to continuum models or to quantum models in which the state variables are computed using the Schrödinger equation, the effective-mass Schrödinger equation, or functional theory (like in the case of density functional theory), however, it might be more difficult to apply it to kinetic Monte-Carlo or molecular dynamics models. The components of the state variable vector depend on the particular transport model used in the description of the fuel cell. For instance, if the PEMFC is described by the mass transport equations, the components of the state vector consists of the oxygen, hydrogen, and water vapor concentrations, electron and proton potentials at each mesh point of the finite element discretization. Both x and μ Pt are column vectors. Vector function F denotes the discretized transport equations (including discretized boundary conductions) for the chemical species involved in the transport. The number of equations in system (4) is equal to the dimension of the state variable vector and the discretized system can be solved iteratively, for instance using the Newton method or other techniques.
Once the state variable vector is computed, the cell voltage can be evaluated using the following equation: where v is a function of the state variable that can be derived by discretizing the analytical equations for the cell voltage. Depending on the transport model used to simulate the PEMFC, function v usually returns the difference between the electrostatic potentials of electrons at the cathode and anode contacts.
Computation of the catalyst sensitivity functions.-The catalyst sensitivity functions depend on the current profile of the catalyst distribution and could in principle be computed as follows: (a) solve discretized transport Equations 4 using the catalyst distribution for which we want to compute γ I (r) or γ V (r) (b) compute the value of the discharge current or cell voltage (whichever is of interest) (c) add a small amount of catalyst at location r and solve discretized transport Equations 4 again (d) compute the value of the discharge current or cell voltage using the new catalyst distribution (e) compute the catalyst sensitivity function at location r and repeat steps (c) and (d) for each location inside the CL The above algorithm suffers from two deficiencies. First, the computed values of the catalyst sensitivity function are subject to large numerical errors because the computation involves subtracting two large but close quantities. Second, the algorithm is computationally very expensive because it requires solving the transport equations for many times (once for each node of the finite element discretization). Because of its large computation cost, the above algorithm is practically impossible to use for the computation of sensitivity functions unless it is implemented in very large parallel computer clusters (with over 10 4 processors in the case of 2-D simulations or 10 6 processors in the case of 3-D simulations). Next, we present a new approach based on small-signal perturbations that avoids both deficiencies.
Let us consider that there is an infinitesimally small variation δμ Pt . This variation results in a new distribution of the catalyst profile, μ Pt = μ Pt 0 + δμ Pt and of the state variable, x = x 0 + δx. The variation of the state variable, δx, can be computed using small-signal perturbations where F x and F μ Pt are the Jacobian matrices of the discretized transport equations with respect to the state variable and parameters μ Pt , respectively. The last equation can be solved to compute where F −1 x is the inverse matrix of F x evaluated at (x 0 ,μ Pt 0 ). Since we assume that the cell is operated under constant d.c. current, the variation of the cell voltage due to δμ Pt can be computed using small-signal perturbations as where v x denotes the gradient vector (a row vector) of the cell voltage with respect to the state variable. Equations 7 and 8 imply The last equation shows that the small-signal variation of the cell voltage is proportional to the variation of the catalyst profile. By comparing Eq. 9 with 2 we conclude that coefficients [10] give the value of the catalyst sensitivity function of the cell voltage (i.e. γ V ) at mesh point n. In the above equation (...) n denotes n-th component of the enclosed vector and V n is the volume of the Voronoi cell at mesh point n inside the fuel cell. Eq. 9 also shows that in order to compute the variation in the response function we need to compute the inverse matrix, F −1 x , which usually requires a very large computational overhead. In order to reduce this overhead let us introduce row vector g V = v x F −1 x , which can be computed by solving the following linear system of equations [11] where g V and v x denote the transpose vectors of g V and v x . With this notation, the mesh point values of the catalyst sensitivity function of the cell voltage can be computed as Notice that the computational cost in computing the discretized values of the sensitivity function is given by the time required to solve linear system (11) and the time to perform the vector-matrix multiplications in by Eq. 12.
The sensitivity function of the discharge current can be computed using the same line of reasoning as above. We obtain γ I,n = − g I F μ Pt n V n [13] where g I is the solution of the following linear system of equations [14] where i x is the transpose of the gradient of the discharge current.

Optimization Algorithm
In this section we present the algorithm for the optimization of the catalyst distribution. This algorithm uses the catalyst sensitivity functions computed in the previous section and is based on steepest descents. We consider the case when the cell current is constant and maximize the cell voltage by solving numerically the following two optimization problems: Problem 2: optimization with constraints max μ Pt (r) V (μ pt ) such that m Pt is constant (e.g. is specified by the manufacturer) [16] In the first problem we optimize the catalyst distribution such that the discharge voltage is maximum. In general, we expect the cell voltage to increase when we increase the amount of the catalyst, however, if the amount of the catalyst becomes too large, the volume fractions of the Nafion and pores decrease and result in a decrease of the proton conductivity and H 2 and O 2 diffusion coefficients. Hence, there should be an optimum catalyst profile that maximizes the power of the cell.
For low amounts of Pt, the power of the PEMFC increases abruptly with the amount of catalyst in the catalyst layer, however, this increase becomes relatively small for moderate and large amounts of Pt. For this reason, it becomes financially more advantageous to impose an upper limit on the amount of catalyst used in the cell. In this case one should solve Problem 2.
Next, we focus on solving problem (16); the unconstrained maximization problem (15) is discussed at the end of the section. Using the formalism of variational (functional) analysis let us consider that we start with an initial platinum distribution μ Pt0 , which is modified by an infinitesimally small amount δμ Pt . Problem (16) becomes max δμ Pt (r) To control the amount of platinum that changes we impose the additional constraint |δμ Pt (r)| 2 d r ≤ δω 2 [18] where δω is a small quantity measured in g Pt m −3/2 . With this constraint the optimization problem becomes max δμ Pt (r) V (μ Pt0 + δμ Pt ) , such that |δμ Pt (r)| 2 d r ≤ δω 2 and δμ Pt (r) d r = 0 [19] where ω is a small quantity. The last integral in Eq. 19 has been obtained by applying Eq. 1. Problem (19) can be solved by defining [20] where λ and λ m are the corresponding Lagrange multipliers and taking the first-order (Karuch-Kuhn-Tucker) necessary conditions where δ δ(δμ Pt ) L(δμ Pt , λ, λ m ) is the functional (variational) derivative of the Lagrangean with respected to the catalyst variation and The first two equations can be solved to obtain [27] Integrating Eq. 27 over and using (26) we obtain that λ m is the average of the catalyst sensitivity function [28] where V is the volume of the catalyst layers (if the optimization is performed only over the cathode CL, V is the volume of the cathode CL). Once Lagrange multiplier λ m is computed, it can be inserted in (27) to compute the variation in the catalyst profile. Equation 27 shows that the catalyst density should be decreased in regions where γ V (r) > λ m and increased in regions where γ V (r) < λ m . In the case of unconstrained maximization (i.e. Problem 1), the catalyst profile should be modified using the following equation [29] where δω is again a small quantity proportional to the amount of change requested by the optimizer. As expected, in the case of unconstrained optimization the catalyst density should be decreased in regions where γ V (r) > 0 and increased in regions where γ V (r) < 0. Equations 27-29 can be extended to the optimization of other two-dimensional or three-dimensional field parameters in PEMFCs such as the porosity distribution in the CLs or GDLs or the ionomer content. For instance, let us consider the optimization of the porosity distribution in the CLs and GDLs layers. The problem can be formulated mathematically as compute the optimum porosity distribution ε(r) that maximizes the cell voltage under constant current density, while keeping the total amount of the platinum loading constant. In this case it is convenient to introduce the porosity sensitivity function  [30] In this case the porosity needs to be changed iteratively according to Equation 31 δε [31] where λ ε m is a new Lagrange multiplier, which is computed as

Numerical Implementation
The numerical algorithm for solving optimization problems (15) and (16) is summarized in Fig. 2. More details about the numerical implementation and computational complexity of the algorithm are presented below.

Start with an initial catalyst distribution, μ (0)
Pt (r), that can be obtained by distributing the given amount of platinum uniformly inside the CLs (in the case of constrained optimization) or by starting with a small amount of platinum and distributing it uniformly inside the CLs (in the case of unconstrained optimization).

Solve transport Equations 4 to compute the state variables and cell
voltage for the initial catalyst distribution under the cell current for which the optimization needs to be performed. The solution of the transport equations can be computed by applying the Newton technique to the discretized transport Equations 4, that is by The transport equations can also be solved using commercially available software packages such as COMSOL or Matlab. 3. Assemble vector v x . 4. Compute vector g V by solving sparse linear system (11). Notice that the Jacobian matrix F x is already known at step 1. 5. Assemble matrix F μ Pt , which gives the derivatives of the transport equations with respect to the catalyst density. 6. Compute the discretized values of the catalyst sensitivity function γ V using Eqs. 12. 7. Update the catalyst profile using equation μ Pt (r) is computed using Eq. 29 in the case of unconstrained optimization (Problem 1) or Eq. 27 in the case constrained optimization (Problem 2). If the optimization has converged (for instance, δμ The most computationally expensive part of the algorithm is to solve the discretized transport equations at step 2. This step requires between 3-10 Newton iterations and each iteration requires solving one sparse linear system of equations. Computing the catalyst sensitivity functions requires the solution of only one sparse linear system of equations at step 4, which has the same numerical complexity as performing one Newton iteration. All the other steps involve sparse matrix-vector multiplications and can be performed efficiently on most computers. The total number of iterations, j, in the above algorithm varies between 5-30, depending on the initial guess and on how the "step size" δω is computed in Eqs. 27 and 29. To decrease the number of iterations parameter δω should be re-evaluated each time the catalyst is updated. Most step-size selection algorithms that do not involve higher (than one) order derivatives can be applied for this purpose (for instance, Nocedal et al. 45 ). In our implementation we prefer to start with an initial step that results in a 10% modification in the magnitude of the catalyst density and, then, adjust the step depending on how close the projected voltage increase is to simulations.
The above numerical algorithm can be easily extended for the optimization of the cell current under constant cell voltage. In this case, steps 3-6 should be modified as follows: at step 3, v x should be replaced b i x ; at step 4 g V should be computed by solving linear system (14); and at step 6 the discretized values of the catalyst sensitivity function γ I should be computed using Eqs. 13.

Optimization Study and Sensitivity Analysis
The above algorithm has been implemented in our device and network simulation platform, RandFlux. 46 RandFlux is a mixed-mode finite element simulator that allows simulating networks of one or multiple finite element devices. The software can import finite element meshes generated by external mesh generators such as Gmsh (http://geuz.org/gmsh) and export the final data for visualization in VisIt (https://wci.llnl.gov), Paraview (http://www.paraview.org), or MATLAB. So far, RandFlux includes libraries for simulating and optimizing semiconductor devices, Li-ion and Li-air batteries, and PEMFCs. The Jacobian matrices and the gradients with respect to the optimization quantities (such as catalyst density) are computed using automatic differentiation. This approach simplifies the computer code significantly, particularly when the functional form of the transport equations is complex.
In this section we present an optimization study for a fuel cell in which the thickness of the GLDs is 250 μm, the thickness of CLs is 15 μm, the thickness of the membrane is 25 μm, and the channel and land widths are 1 mm. The cell is assumed to operate at 80 • C and the inlet pressures at the anode and cathode are 1.8 atm. The fractions of the void and electron conductive materials are 0.5 in both gas diffusion layers. The fraction of the void in the catalyst layers is 0.4, the fraction of the electron conductive material is 0.3, and the fraction of the Nafion is 0.3 (binder-free cell). The PEMFC is modeled by solving the mass continuity equations for oxygen hydrogen, and water vapors in all five regions of the fuel cell including the anode GDL, anode CL, membrane, cathode CL, and anode GDL, and the current continuity equations for protons and electrons in Nafion and in the electron conductive material, respectively (see Fig. 3). More details about the model used in our work are given in the Appendix.
The values of all the material dependent parameters have been taken from the literature. [47][48][49][50][51] The effective diffusion coefficients and electric conductivities are computed using standard Bruggeman relationships with β= 1.5. Reaction rates are assumed to follow Butler-Volmer kinetics in uniformly distributed spherical agglomerates of platinum particles, electron conductive material, and Nafion of radius r agg = 130 nm. The "bulk" diffusion coefficient of hydrogen in the pores D H 2 = 0.69905 cm 2 s −1 and in Nafion D N H 2 = 1.2 × 10 −6 cm 2 s −1 , the "bulk" diffusion coefficient oxygen in the pores D O 2 = 0.69905 cm 2 s −1 and in Nafion D N O 2 = 8.46 × 10 −6 cm 2 s −1 , and the "bulk" diffusion coefficient of water vapors in pores D H 2 O = 0.258 cm 2 s −1 . In addition, the Henry constant of hydrogen gas at the pore-Nafion interface H H 2 ,N = 6.69 × 10 10 Pa cm 3 mol −1 , the Henry's constant of oxygen gas at the pore-Nafion interface H O 2 ,N = 3.1664 × 10 10 Pa cm 3 mol −1 , and the thickness of the thin film on the Nafion surface t N = 10 nm. The conductivity of the membrane was 52 σ m = (0.00326 + 0.005139λ) exp[1268( 1 300 − 1 T )] (S cm −1 ), where λ = 0.043 + 17.81a− 38.85a 2 + 36 if a ≤ 1 and λ = 12.6 + 1.4a if a > 1, where λ is the water content in Nafion and a is the ratio of the number of water molecules to the number of (SO 3 )H + sites. The saturation pressure of water vapor was 52 P sat H 2 O = −2.1794+0.02953T −9.1837×10 −5 T 2 +1.4454×10 −7 T 3 . The exchange current densities for the anode reaction was i 0A = 1.7 A cm −2 while for cathode was given by i 0C = 0.729 A cm −2 if the overpotential was below 0.8 V and 0.007906 A cm −2 if the overpotential was above 0.8 V. The anodic charge transfer coefficient at the anode was 0.5, while the cathodic charge transfer coefficient at the cathode was 0.39 if the overpotential was below 0.8 V and 0.007906 if the overpotential as above 0.8 V. 47 Finally the reaction voltage at the cathode was computed using U 0 = 1.229 − 0.000846(T − 300) (V).
Since the amount of catalyst in the anode CL is usually much smaller than in the cathode CL (because of the much faster hydrogen reaction), we have optimized only the catalyst distribution in the cathode CL, however the transport equations are solved throughout the entire cell. Therefore, the optimization domain denoted by in Equations 1-3 is the volume of the CL at the cathode. Fig. 4 shows the geometry and the 2-D finite element mesh that was used to simulate the cell. The finite element grid contains 6,342 nodes and 11,374 elements and was obtained by refining the grid until the computed values of the cell voltage did not change by more than 1 mV when the density of the mesh points was doubled. Fig. 4 also shows the region where the catalyst distribution and sensitivity functions are represented in the contor plots that are presented later in this section.
Unconstrained optimization.-To verify the numerical implementation of our algorithm we start by considering the case of unconstrained optimization, in which we compute the optimum distribution of the catalyst without any constraints on the total amount of Pt used. The convergence of the algorithm was checked by staring with different initial guesses for the catalyst loading (m Pt /A, where A is the cross-sectional are of the cell) in the cathode. The following two r Initial Guess 2: m Pt /A = 5 mg/cm 2 , which corresponds to ρ Pt (r) = 3.57 mg cm −3 .
The results of the unconstrained optimization are presented in Fig.  5, which shows the increase in the cell voltage at different current densities with respect to the iteration number. The results of the optimization in the case when the initial platinum loading is "Initial Guess 1" are represented by continuous lines, while dash lines show the results when the initial platinum loading is "Initial Guess 2". The optimization in Fig. 5a is performed in order to increase the cell voltage at a discharge current of 0.5 A/cm 2 , while the optimization in Fig.  5b is performed in to increase the cell voltage at a discharge current of 3 A/cm 2 . Notice that the final (optimized) cell voltage is the same for both guesses. Fig. 6 presents the catalyst density and sensitivity function of the cell voltage obtained at different iterations when the initial guess of the platinum loading is "Initial Guess 1". Fig. 7 presents the catalyst density and sensitivity function of the cell voltage obtained at different iterations when the initial guess of the platinum loading is "Initial Guess 2". In both cases the optimization is performed at a current density of 0.5 A/cm 2 . As expected, in the case of "Initial Guess 1" the catalyst density is monotonically increasing during the optimization problem, while in the case "Initial Guess 2" the catalyst density is monotonically decreasing. However, the final catalyst distribution is the same for both initial guesses and corresponds to an average catalyst loading of approximately 2 mg Pt /cm 2 . It is interesting to observe that in the case of "Initial Guess 1", the catalyst density increases close to the membrane interface during the first iterations and toward the GDL during the final iterations. As we will confirm later, this fact suggests that, in the case of low platinum loading it is more advantageous to use more catalyst at the membrane interface than at the GDL interface. In the case of high platinum loadings it is more advantageous to use less catalyst at the membrane interface than at the GDL interface. Fig. 8 shows the catalyst density and sensitivity functions of the cell voltage obtained at different iterations when the optimization is performed for a current density of 3 A/cm 2 and the initial guess of the platinum loading is "Initial Guess 1" ("Initial Guess 2", not shown here, results in the same final profile of the catalyst distribution). Due to the large value of the final catalyst loading but also to the fact that the optimization is performed at a high value of the current density, when the water concentration is very high, the final catalyst distribution is larger at the GDL side. Fig. 9 shows the polarization curves computed by using the optimized catalyst distributions obtained in the optimizations presented in Figs 4 and 6, respectively. There is relatively insignificant difference between the two polarization curves because the catalyst profiles are relatively similar an correspond to a relatively large platinum loading of 2 mg Pt /cm 2 .
Constrained optimization.-Next, let us focus on the case of constrained optimization, in which we compute the optimum distribution of the catalyst by fixing the total amount of Pt in the catalyst layer to smaller amounts, such as the ones used practical applications. We present optimization results for two platinum distributions with average loadings of 0.1 mg Pt /cm 2 and 0.2 mg Pt /cm 2 , respectively. The platinum loading was initially distributed uniformly throughout the catalyst layer. Then, the catalyst sensitivity function of the cell voltage was computed and the catalyst was redistributed according to Eq. 27. The cell current was kept constant and equal to 0.5 A/cm 2 in both cases. Fig. 10 presents the increase of the cell voltage as a function of the iteration number. Notice that the cell voltage increases by approximately 9 mV when the average catalyst loading is 0.1 mg Pt /cm 2 and by 35 mV when the 0.2 mg Pt /cm 2 by optimizing the cell. Similar values were obtained by Ebrahimi et al. 15 who computed the optimum catalyst distribution using the Nelder-Mead Simplex technique and obtained an increase in the power density of 7% and by Cetinbas et al. 43 who obtained an increase of 10-14%. Figures 11a and 11b present the optimized catalyst distributions for the two average Pt loadings, respectively. In both cases the power of the cell increases when the density of platinum particle is higher at the membrane side than at the GDL side. In fact, our numerical simulations predict that more than half of the platinum particles should be placed in the CL in a region of width of 5 μm adjacent to the membrane, while the remaining particles can be distributed non-uniformly throughout the remaining volume of the CL. According to these simulations, the optimum concentration of platinum particles is 4-8 times larger at the membrane side with respect to GDL side. These results are in good agreement with the existing optimization studies presented in the literature. For instance, by assuming that the catalyst distribution varies only in the x-direction, Ebrahimi et al. 15 obtain that the catalyst density, denoted by h(x) in their study should be 6-8 times larger at the membrane than at the GDL layer. Cetinbas et al. 43 obtain that the   catalyst density should by only 3.5 times larger at the membrane than at the GLD layer, while Song et al. 53 obtain an increase of over 4 times. Due to the low cell current and relatively small diameters of the agglomerate particles the platinum should be deposited relatively uniformly in the y-direction (along the membrane interface). Since the two catalyst distributions depend approximately only on the x-variable, Fig. 11c presents cross-sections made at y = 0.5 mm. These plots are good qualitative agreement with the analytical results obtained by Kulikovsky. 24 Figure 9. Polarization curve obtained using the two optimized catalyst distributions obtained in Fig. 6 and Fig. 8, respectively.

Conclusions
We have introduced a numerical algorithm for the computation of the optimum catalyst distribution in PEMFCs that maximizes the cell voltage of these cells under a given current dentiy. The computed catalyst distribution is optimal in the sense that any variation from this distribution is guaranteed to reduce or keep constant the power density of the cell. Although we do not prove that the computed catalyst distribution is globally optimum, our simulation results suggest that the solution is unique. In agreement with previous experimental results  and theoretical predictions we obtain that, depending on the structure and dimensions of the PEMFC, the cell voltage can be increased between 2-10% by keeping the mass of the catalyst constant but placing the platinum particles non-uniformly throughout the CLs.
The numerical approach presented in the article takes advantage of an adjoint space technique that can be used to evaluate the first-order derivatives of the objective function with respect to the optimization parameters efficiently, by solving only one sparse linear system of equations. This technique has been used before by the applied mathematical and computation electronics communities and is adjusted in this article to optimize the platinum distribution inside the CLs. We show that the optimized platinum distribution profiles depend on the discharge conditions under which the optimization is performed (for instance, on the discharge current) and on the basic geometry of the cell. In particular, the optimized profiles depend on the positions of the landings and gas channel openings, and on the thickness and microstructure of the catalyst layers. In general, one should place more catalyst under the gas channel than under landings and toward the membrane interface. In this way, the electrostatic potential drop in the electrolyte is smaller and the oxygen concentration at the reaction sites is larger.
The numerical algorithm presented in this article requires a relatively small computational overhead compared to the computation of polarization curves. Indeed, computing the variation of the catalyst distribution at each iteration requires solving only one system of linear equations of the same dimension and sparsity as the transport system of equations. Since solving the transport equations requires approximately 5-10 Newton iterations and assuming that each polarization curve is given by 20 points, it means that the computation time spent on computing the optimum catalyst distribution is between 100-200 times smaller than the time required to compute the polarization curves.
Future work will focus on adding other field parameters in the optimization including the porosity, ionomer content, and particle size distributions. 54,55 The computational cost for computing the sensitivity function for these distribution is comparable to the computational cost for the calculation of the catalyst sensitivity function and require only a few additional matrix-vector multiplications similar to Eqs. 12 and 13. Additional parameters, such as operating temperature and inlet pressure can also be optimized using the adjoint method and will be discussed in future work. equations give similar values for A S at low platinum concentrations, we found that Eq. A6 approximates the catalyst surface area more accurately at high platinum concentrations. At the anode reaction we use 50 A S (μ Pt ) = 7.401(Pt|C) 4 − 18.11(Pt|C) 3 + 15.45(Pt|C) 2 −6.453 (Pt|C) + 2.045] × 10 6 [A7]