Stochastic Nonlinear Optimization for Robust Design of Catalysts

ARTICLE pubs.acs.org/IECR

Stochastic Nonlinear Optimization for Robust Design of Catalysts Chang Jun Lee,† Vinay Prasad,† and Jong Min Lee*,‡ † ‡

Chemical and Materials Engineering, University of Alberta, Edmonton, Alberta, Canada, T6G 2V4 School of Chemical and Biological Engineering, Seoul National University, Seoul, 151-744, Korea ABSTRACT: Computational methods for designing an optimal catalyst have recently received much attention, especially for energy-related applications. What is lacking in the previous methods is an explicit method to handle uncertainties in the complex models used, so that a robust design is achieved. This work proposes a stochastic optimization method for designing a robust catalyst. In particular, reactions involved in catalytic decomposition of ammonia are presented, and uncertainties associated with experimental determination of kinetic parameters are represented as exogenous variables with assumed probability distributions. The problem is formulated in terms of finding the optimal binding energies that maximize conversion in a microreactor. The resulting stochastic optimization problem is nonlinear, and involves the expectation operator as well as integration in the objective function. This difficult optimization problem is tackled by a population sample based approach, referred to as particle swarm optimization. The results show that the value of solving the stochastic problem is significant, and that it can provide a more robust solution compared to the certainty equivalence approach.

1. INTRODUCTION Energy concerns affect all of us in the world today, and much research is being devoted to finding more efficient and greener methods of energy generation. Catalysis has a very important role to play in terms of both improved efficiency and environmental impact, since a large majority of the processes used to convert fossil-based and renewable fuels into energy make use of catalysts. An emerging area in catalysis is the rational model-based design of catalysts, which aims to complement high-throughput experimentation for the discovery of new, multifunctional catalysts. Computational methods for optimal design of catalysts in energy processes have recently received much attention to replace intuitive or trial-and-error approaches. The current approach to computation-based rational catalyst design involves the development of deterministic models and the search for an optimum set of catalyst properties using deterministic optimization.1-3 The models are usually highly simplified, and the decision variables in the optimization are so-called model descriptors such as binding energies of elementary species involved in the chemistry. An unappreciated fact is the existence of uncertainty in the model parameters, even though the optimal parameter values are determined using experiments and parameter estimation. There are many sources of this uncertainty, for example, the inherent uncertainty in the first-principles calculation of binding energies using density functional theory (or through experiments), the uncertainty in the scaling relations (e.g., Brønsted-EvansPolanyi relations4) that relate activation energies of individual elementary reactions in the mechanism to the binding energies, and (perhaps most importantly) the uncertainty in pre-exponential factors across catalysts made from different materials. This uncertainty may lead to the deterministically identified optimal catalyst properties not being the most probable optimal catalyst. Stochastic optimization approaches optimize the expectation of the objective function of the decision variables and the r 2011 American Chemical Society

uncertainties and thus come up with a decision that will perform best on average. These techniques take advantage of the fact that probability distributions governing the uncertain parameters are known or can be assumed reasonably. The solution from stochastic optimization provides more realistic predictions of performance than that of the certainty equivalence (CE) approach, where uncertain parameters are replaced by their means. This motivates the objective of our research, which is to use stochastic optimization to identify the optimal catalyst. Stochastic optimization techniques have been applied to similar problems in computer aided molecular design and process design under uncertainty.5-8 Maranas uses chance-constrained programming for optimal polymer design under uncertainty, where a mixed integer nonlinear programming problem is converted to a deterministic equivalent linear problem.6 Shastri propose a modified L-shaped algorithm (BONUS) to convert a stochastic nonlinear problem into linear subproblems.8 However, the models used in these previous approaches are primarily linear, and highly restrictive assumptions such as convexity and differentiability of models are used in the development of solutions. Models describing the activity of catalysts do not generally have analytical solutions, and solving them involves numerical integration. Moreover, reaction systems generally exhibit nonlinearity due to the exponential term in kinetic expressions. These make it impossible to convert the stochastic problem into a deterministic counterpart, and precludes the use of these methods in a complex system such as ours. In this work, computationally efficient stochastic approaches to optimal selection of a catalyst are proposed and demonstrated. In the proposed scheme, the expectation of the objective value is Received: October 17, 2010 Accepted: January 21, 2011 Revised: January 9, 2011 Published: March 01, 2011 3938

dx.doi.org/10.1021/ie102103w | Ind. Eng. Chem. Res. 2011, 50, 3938–3946

Industrial & Engineering Chemistry Research

ARTICLE

evaluated using a sample average approximation (SAA) method, due to the simplicity of the approach. However, exponential increase in the number of samples with the number of uncertain parameters defies using all possible combinations of the uncertain parameters. In order to reduce the number of random variables (uncertain parameters), linearity analysis, together with partial least squares, is proposed, given the fact that a CE solution is optimal when the objective function has only a linear dependence on the variables. In other words, the CE approach is used for uncertain parameters whose effect on the output is found to be approximately linear, and the parameters showing nonlinear effects are treated with the SAA approach. In the optimization step, a gradient-free sample-based method, referred to as particle swarm optimization (PSO), is employed to circumvent the issues involving local optima and computational time associated with gradient-based approaches. The proposed stochastic approach is applied to the optimal design of catalysts, which involves the specification of optimal binding energies of hydrogen and nitrogen on the catalyst surface. The next step is to identify materials (perhaps bimetallic alloys) that possess these optimal binding energies. That is a separate optimization problem, and is not considered in our work. The system we use to demonstrate our approach is the catalytic decomposition of ammonia to produce hydrogen. In addition to being a well-studied model system often used for demonstrations of proof of concept in catalytic studies, it also has practical importance related to hydrogen production for fuel cells. Hydrocarbon feeds produce carbon monoxide on reforming in addition to hydrogen, and it poisons fuel cell catalysts; this problem is circumvented by the use of ammonia as the feed. The infrastructure is also well developed because of the fertilizer industry, and thus ammonia has been considered as a viable alternative in the short term as a fuel for hydrogen production for fuel cells, even though it is made from hydrogen. This paper is organized as follows. In section 2, the basic model and optimization problem that we use are reviewed. Section 3 explains the PSO briefly, and section 4 proposes an efficient sampling approach to reduce the number of random variables. Section 5 discusses the results of our stochastic model compared to the results of deterministic problem. Finally, section 6 gives concluding remarks.

Table 1. Rate Constants of Elementary Reactions at 300 K reaction type

rate constant rffiffiffiffiffiffiffiffiffiffiffiffiffi β RT T e - Ea ðθ, TÞ k ¼ s 2πMW T0

adsorption

desorption

surface

k ¼ A

β T e - Ea ðθ, TÞ T0

k ¼ A

β T e - Ea ðθ, TÞ T0

Table 2. UBI-QEP Relations of Atomic Bonding Energies to Species Binding Energiesa term

expression

QH

input

QN

input QH

Q*H

2QN

Q*N

2-

QN

k1 k2

QNH2*

Q 2 Q þ DNH2

QNH *

Q 2 Q þ DNH

k5

NH/ þ / sF Rs N/ þ H/

ðR3Þ

k7

ðR4Þ

k6

NH2 / þ / sF Rs NH/ þ H/ k9

NH3 / þ / sF Rs NH2 / þ H/ k10

k3

N2 þ 2/ sF Rs 2N/ k4

ðR1Þ

ðR2Þ

QN þ DNH3

a QH and QN are the binding energies, and the quantities with an asterisk (/) are those used within UBI-QEP theory. NH and NH2 species are modeled with the strong binding UBI-QEP equation, and NH3 is modeled with the weak one.10.

2. MICROKINETIC AND REACTOR MODEL UNDER UNCERTAINTY

H2 þ 2/ sF Rs 2H/

1 3

2

QNH3*

k8

2.1. Basic Kinetic Models. The model used in this work describes the catalytic decomposition of ammonia in a plug flow reactor (PFR), and the kinetics are described by a global rate expression that was developed to explain ammonia decomposition on ruthenium catalyst (see refs 9 and 10). The model was developed considering the following reversible elementary reactions (/ indicates a vacant site or, in conjunction with a chemical species, an adsorbate):

1 3

k11

NH3 þ / sF Rs NH3 / k12

ðR5Þ ðR6Þ

These elementary reactions are considered as six linearly dependent reactions. The rate constants of elementary reactions are expressed in the modified Arrhenius form as shown in Table 1. For details about the parameters, see ref 10. The activation energies of these reactions are calculated from the binding energies of hydrogen and nitrogen based on the unity bond index-quadratic exponential (UBI-QEP) method.11 The UBI-QEP method uses semiempirical relationships to relate the species binding energies to the atomic binding energies10 (as shown in Table 2) and bond indices (BI) to describe how the activation energy of a particular reaction relates to the binding energies (QA, QB) and product dissociation 3939

dx.doi.org/10.1021/ie102103w |Ind. Eng. Chem. Res. 2011, 50, 3938–3946


ARTICLE

energies (DAB): QA QB 2QA QB Ea ¼ þ DAB BI2 þ - 2DAB BI QA þ QB Q A þ QB Q A QB - QA - QB ð1Þ þ QA þ QB

Table 3. Nominal Values for the Elementary Reactions (Taken from ref 10) reaction

The global rate expression is developed without the assumption of a single rate-determining step among the elementary reactions, but from detailed analysis of a full microkinetic model9,10 based on many simulations under a range of operating conditions. Reactions R1, R5, and R6 are assumed to be in partial equilibrium, and the most abundant reactive intermediates are adsorbed elemental hydrogen and nitrogen. Equation 2 is the global rate expression we adopt from ref 10 for computing the conversion. r ¼

-2ðk4 ω - k3 pN2 Þ !2 rffiffiffiffiffiffiffiffiffiffiffi k11 k1 1 þ pNH3 þ pH þ ω k12 k2 2 2

ffi sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi rffiffiffiffi k3 k2 k7 k9 k11 ω¼ pN þ pNH3 pH2 -0:5 k4 2 k1 2k4 k10 k12

ð2Þ

ð3Þ

where k is the reaction constant and p is the partial pressure. Nine reaction constants are used in this equation. 2.2. Stochastic Modeling. The nominal binding energies of nitrogen and hydrogen are modeled as functions of surface coverages (θN, θH): QH ¼ QH0 þ C1 θH

ð4Þ

QN ¼ QN0 þ C2 θN

ð5Þ

where the binding energies on a clean surface, QN0 and QH0, are in the ranges of 100-175 and 45-75 kcal/mol, respectively. These ranges encompass most practical values of the binding energies. C1 and C2 describe how the binding energies change with surface coverage. Q, C1, and C2 are determined using density functional theory calculations or experimental studies. However, it is not possible to specify the exact values of C1 and C2 at each iteration of the optimization procedure, since each set of nominal binding energies represents a different catalyst surface with potentially different coverage dependence. Thus, these inherent uncertainties associated with the coverage dependence of binding energies are modeled as exogenous random variables added to the nominal values. QH ¼ QH0 þ εQH

ð6Þ

QN ¼ QN0 þ εQN

ð7Þ

where both εQH and εQN are assumed to follow uniform distribution U[-5,5], since the error of the density functional theory calculations or temperature programmed desorption experimental studies for determining binding energies of each species is approximately (5 kcal/mol and this error may lie anywhere within the range.10 There are also uncertainties in the activation energies of individual elementary reactions, and they are also modeled as exogenous random variables, added to the nominal value as

term

A [s-1], S

β

BI

H2 þ 2/ f 2H/

ðR1fÞ

S1f

0.87

0

0.5

2H/ f H2 þ 2/

ðR1bÞ

A1b

2.4 1013

0

0.5

N2 þ 2/ f 2N/

ðR2fÞ

S2f

2.2 10-7

0

0.8

2N/ f N2 þ 2/

ðR2bÞ

A2b

1.6 1010

0

0.8

NH2 / þ / f NH/ þ H/ ðR4fÞ

A4f

7.2 1012

0

0.5

NH3 / þ / f NH2 / þ H/ ðR5fÞ

A5f

1.4 1013

0

0.5

NH2 / þ H/ f NH3 / þ / ðR5bÞ

A5b

8.8 1012

-0.39

0.5

NH3 þ / f NH3 /

ðR6fÞ

S6f

0.63

0

0.5

NH3 / f NH3 þ /

ðR6bÞ

A6b

4.4 1013

0

0.5

follows: j ¼ 1, 2, 4, 5, 6

Ea, Rj ¼ Ea0, Rj þ εEa, j ,

ð8Þ

where j is the type of reaction. εEj is assumed to follow Gaussian distribution, N(0,(3/2.58)2).10 Since elementary reaction R3 is not used in the global rate expression, the total number of random variables for activation energies is 5. ERj,0 is calculated with the UBI-QEP method, using QN and QH. The N(0,(3/ 2.58)2) distribution is our estimate of the errors in the fitting of UBI-QEP. The nominal values of pre-exponential factors are estimated from extensive experimental data10 on ruthenium catalysts. The same nominal parameters are assumed to hold for all catalysts in this study. The uncertainty in pre-exponential factors across catalysts is made from the experiments or different materials and these factors can be considered as the exogenous random variables. Afj ¼ Afj, 0 þ 10

εA f

ð9Þ

εA b

ð10Þ

Abj ¼ Abj, 0 þ 10

j

j

where j = 1, 2, 4, 5, 6, and Aj,0 is the nominal value. Since the kinetic constants for the elementary reaction R3 and the backward reaction of elementary reaction R4are not used in the global rate expression, the total number of random variables for pre-exponential factors is nine. In this study, εAj is assumed to follow the uniform distribution U[-2,2]. The estimates of pre-exponential factors calculated from transition state theory are typically accurate only to 2 orders of magnitude.10 If pre-exponential factors are fitted from experiments, too, there are multiple minima. Finally, as the binding energies change as we go between different catalyst materials, the pre-exponential factors for each material are different. For ammonia decomposition, the literature shows a spread of approximately (2 orders of magnitude over the single metal catalysts studied. This is why a uniform distribution, U[-2,2], is selected. In case of adsorption reactions (R1f, R2f, and R6f in Table 3), the range of the uniform distribution is (100% of S. In total, 16 random variables are 3940



ARTICLE

involved in our stochastic model. The nominal values of preexponential factors, bond indices, and β used here are shown in Table 3.10

3. STOCHASTIC OPTIMIZATION USING SAMPLE AVERAGE APPROXIMATION AND PARTICLE SWARM OPTIMIZATION In order to design the catalyst for improving the conversion of ammonia decomposition in a PFR, the optimal atomic binding energies should be determined to maximize the conversion and the following stochastic optimization problem is formulated. max

QHO , QNO

Eε ½x

ð11Þ

subject to 1 x¼ F0

Z

V

ð - rÞ dV

ð12Þ

0

where F0 is the flow rate, V is the reactor volume, and r is given by eq 2. Equation 12 is the conversion for gas-phase reaction in a PFR.12 As described under section 2, the set of random variables, ε, is modeled as a random vector with known probability distributions. Note that solving this optimization problem is challenging due to the following reasons: • The overall reaction rate is highly nonlinear. • The conversion can only be evaluated by numerically integrating the rate expression, given the initial conditions, parameters, and random variables. • The expectation operator is also difficult to evaluate, because the components of the random variables ε ∈ R16 are independent of each other and the computational complexity increases in a faster-than-linear fashion in the number of sampled points. One computational bottleneck to solving stochastic optimization problems is the calculation of the expected objective function. Many stochastic optimization problems involve a large number of scenarios, making it prohibitively expensive to calculate the expectation of the objective function. To overcome this limitation, several researchers have suggested approximation schemes. There are two major approaches to evaluating the expectation approximately. One is based on the decomposition of sample spaces, and the other is based on sampling data points. The former class of algorithms is a successive approximation method in which decomposition of the sample space is pursued for assuring efficiency in solution.13 However, this approach is limited because the convexity, monotonicity, and differentiability of constraints should be established before the decomposition. Sample average approximation (SAA) evaluates the expected objective value (f(Q) = E[x(Q, ε)]) by generating samples based on Monte Carlo sampling techniques and averaging the corresponding objective values as follows: N

∑

^f ðQ Þ ¼ 1 xðQ , εk Þ N N k¼1

ð13Þ

where N is the number of sampled points from known probability distributions. The sample-based approaches are easy to implement even for a system such as ours, with complex constraints and the function evaluation for the optimization involving numerical integration.

It should be noted that SAA is used only for handling the expectation operator and one still has to solve the maximization problem by an appropriate numerical procedure. Since the objective function is not in an explicit form with respect to the decision and uncertain parameter vectors, it is difficult to apply gradient-based deterministic optimization methods. In this work, a heuristic optimization method, referred to as particle swarm optimization (PSO), is used due to its advantages in handling complex objective functions without evaluating derivatives. PSO is a population-based sampling optimization technique motivated by the social behavior of collections of animals.14 It starts with randomly generated swarms of particles (initial guesses for the optimal states) and remembers the best solution found. The particles move around the solution space with velocities determined by the algorithm, and they move toward the global optimal solution over the optimal procedure. The attractive features of PSO are that it does not need to evaluate derivatives of the objective function and constraints, and there are relatively a small number of parameters to adjust.15 Table 5 shows a set of parameters used in this study. The steps involved in the PSO algorithm15 are as follows: 1. Initialize the search parameters. • Niter = number of iterations i to be performed • Npt = number of particles • Nd = number of search dimensions • Np = number of parameters to be estimated • C1, C2, Wint = weight parameters • xmin and xmax = lower and upper limits of the variables to be optimized 2. Calculate the maximum particle velocities for the Np parameters. vmax ¼ j

xmax - xmin j j

2 3. Calculate the initial particle positions and velocities:

ð14Þ

þ rðxmax - xmin xip, d, j ¼ xmin j j j Þ

ð15Þ

vip, d, j ¼ vmax j ð2r - 1Þ

ð16Þ

where p, d, j, and r denote the particle, the search direction, the parameter index, and a random number in the range [0,1]. 4. Evaluate the objective function for each particle which moves in Nd search dimensions with all Np parameters. ind and Npt (vectors with dimension of Np) that 5. Update xp,j contain the best position found by each particle in the swarm. 6. Update xglo (a vector of dimension Np) that contains the best position found by the whole particle swarm. 7. When the maximum number of iterations is reached, the PSO search is terminated. Otherwise update the particle velocities and position using i i vip,þd,1j ¼ Wint vip, d, j þ C1 r1 ðxind p, j - xp, d, j Þ þ C2 r2 ðxj - xp, d, j Þ glo

xip, d, j ¼ xip, d, j þ vip,þd,1j

ð17Þ ð18Þ

8. If the absolute particle velocity is higher than the maximum velocity, then use vip,þd,1j ¼ vmax signðvip,þd,1j Þ j

ð19Þ

9. If the particle position is out of the search limits, then the particle is placed at xind p,j . 10. Increase the iteration count by 1 and return to step 4. 3941



ARTICLE

4. AN EFFICIENT SAMPLING APPROACH The main difficulty with sample-based approaches to stochastic optimization is that the number of scenarios can grow considerably with the number of uncertain parameters. Furthermore, exhaustive sampling is necessary when there is no correlation structure in the uncertainties. However, it is still an open question as to what a reasonable number of scenarios will be to obtain the optimal solution of the problem. Monte Carlo sampling (MCS) technique is one of the most popular sample-based approaches, which chooses the sample points randomly. Though MCS approaches are unlikely to scale exponentially with the number of uncertain parameters, they do not have good uniformity.16 Hammersley points can be used to uniformly sample a high-dimensional hypercube, and it is reported that a Hammersley sequence sampling technique can reduce the number of required sample points for a given accuracy of estimates to about 3-100 times on average compared to MCS or other stratified sampling techniques.17 However, it still suffers from the “curse of dimensionality” when there are a large number of uncertain parameters without any correlations. This work proposes an efficient sampling method that involves thermodynamic constraints, linearity analysis, and partial least squares (PLS) to reduce the number of uncertain parameters to be sampled. 4.1. Thermodynamic Constraints. In constructing kinetic models thermodynamic consistency at the enthalpic and entropic levels should be considered.18 For the elementary reactions in our system, the reaction rate constant is described by the Arrhenius law. We have Afj Abj

¼ eΔSj =R

j ¼ 1, :::, 6

ð20Þ

where j is the reaction index and S is entropy. Subscripts “f” and “b” indicate the forward and backward reactions. Since the reactions in our model are performed at constant temperature and pressure, ΔSj of each elementary reaction should be constant. In generating the uncertain parameters of pre-exponential factors, this thermodynamic constraint should also be satisfied. In generating samples of uncertain parameters, the following additional constraints are used: Afj

¼

Afj, 0

j ¼ 1, 2, 5, 6 ð21Þ Abj Abj, 0 Equation 21 allows us to determine the value of the preexponential factor for a reverse reaction given an estimate of its corresponding forward reaction, and vice versa. This reduces the number of uncertain parameters. 4.2. Linearity Analysis. The certainty equivalence solution, which replaces uncertain parameters with their means, is optimal when the objective is linear in the uncertainties.19 This presents significant numerical advantages in that the uncertain parameters linearly related to the objective function can be omitted in the sampling procedure. Hence, linearity analysis is employed to identify uncertain parameters whose effect on the conversion can be approximated to be linear. The output, i.e., conversion, is computed given the range of each uncertain parameter with the other uncertain parameters fixed at their mean values. A regression line is fitted to the obtained data using the least-squares criterion, and the value of R2 (percentage variation explained) is examined to

determine if there is a linear relationship between the output and the chosen parameter. 4.3. Partial Least Squares. In order to reduce the number of uncertain parameters further, a partial least squares (PLS) model is subsequently constructed to identify uncertain parameters that have negligible effects on the objective function. Even if the relationship between these uncertain parameters and the objective function is nonlinear, uncertain parameters that have negligible effects on the output can be eliminated or replaced with their means.20 The PLS model consists of outer relations (X and Y blocks individually) and an inner relation (linking the two blocks). The scaled and mean-centered X and Y matrices are decomposed into their own score and loading vectors denoted by T and U, and P and Q, respectively. The inner relationship between the two blocks of X and Y is represented as a linear algebraic relation, U = TB, between their scores. X ¼ TPT þ E

ð22Þ

Y ¼ UQ T þ F

ð23Þ

where E and F are the residual matrices for X and Y, respectively. In PLS, the loading and score vectors are determined in such a way to maximize the prediction accuracy of Y while describing a large amount of the variation in X. The PLS regression coefficient, B, represents how much input variables would affect the output. The larger the absolute value of a regression coefficient, the more the corresponding input has an effect on the output. Hence, we ignore the uncertain parameters in the sampling procedure whose coefficients in B are less than a specified threshold value. In our work, X is the set of uncertain parameters and Y is the corresponding conversion, while other parameters are fixed at their nominal values. 4.4. Uniformity in Sampling. In order to ensure uniformity in sampling, the sample points are placed along equally spaced intervals on the grid of each parameter. If four samples are generated for εEa,1, the Gaussian probability distribution function is divided into four areas having the same probability, i.e., area under the curve, and only one point is sampled in each area. The proposed sampling procedure can be summarized as follows: 1. Consider thermodynamic constraints to reduce the number of uncertain parameters. 2. Calculate the conversion according to the changes of each remaining uncertain parameter while other uncertain parameters and inputs are fixed. 3. Fit the data to a line using the least-squares method and evaluate R2. 4. Generate sampling points randomly for the remaining uncertain parameters. 5. Calculate the conversion under all possible cases. 6. Using PLS, evaluate the regression coefficient (B). 7. Identify which variables have negligible effects on the conversion. For these variables, mean values are used. 8. Determine how many samples are generated on each uncertain parameter, and span the probability distribution equally with the number of samples. 9. Generate one sample in each divided area. 10. Construct all possible scenarios of uncertain parameters using the samples. The stochastic optimal solution, which maximizes the expectation of objective function under the sampled scenarios, is searched by PSO. 3942



ARTICLE

Table 4. Operating Conditions of the Catalytic Fixed Bed Microreactor length

0.7 cm

internal diameter

0.41 cm

pressure

1 atm

flow rate

200 sccm

temperature

600 K

site density

1.66 10-9 mol/cm2

surface area/volume

8545 cm2/cm3

Table 5. PSO Parameter Values Used To Find the Optimal QN and QHa parameter

a

value

Niter

42

Npt

10

Nd

1

Np C1

2 2(0.95)i(2r - 1)

C2

2(0.95)i(2r - 1)

Wint

(0.95)i

r denotes a random variable in the range [0, 1].

Figure 1. R2 of linear regression lines between random variables and conversion. The order of random variables follows Table 6.

Table 6. Order of Uncertain Parameters 1

2

3

4

5

6

7

8

9

10 11 12 13 14 15 16

εQN εQH εEa,1 εEa,2 εEa,3 εEa,5 εEa,6 εA1f εA1b εA2f εA2b εA4f εA5f εA5b εA6f εA6b

5. COMPUTATIONAL RESULTS 5.1. Dimensionality Reduction in Sampling. In this section, the results of stochastic optimization applied to the model for ammonia decomposition is presented. The operating conditions of the catalytic fixed bed microreactor are described in Table 4. First, eq 21 was considered in generating samples of uncertain parameters for the pre-exponential factors. Since elementary

Figure 2. Absolute regression coefficients of the PLS model between random variables and conversion. The order of random variables follows Table 6.

reaction R3 and the backward reaction of elementary reaction R4 are not used in the global rate expression, the values of four pre-exponential factors of reverse reactions were determined by the values of their corresponding forward reactions. Hence, the total number of random variables is reduced to 12 with the thermodynamic constraints. Next, to assess linearity of each uncertain parameter, the conversions at various values of each uncertain parameter were computed, while the other uncertain parameters were kept at their mean values. Two thousand data points were generated for each uncertain parameter, and a regression line between the conversion and each random variable was fitted to the data. From Figure 1, it can be seen that nine uncertain parameters have linear relationships with the conversion. In this work, the uncertain parameters whose R2 values are greater than 0.999 are replaced by their mean values. The x-axis of Figure 1 follows the order of random variables in Table 6. PLS was then used to identify uncertain parameters among the remaining seven uncertain parameters with negligible effects on the conversion. Two thousand data points were generated for the uncertain parameters, and their corresponding conversions were calculated. The absolute regression coefficients (B) of the constructed PLS model are shown in Figure 2. It was found that all seven uncertain parameters have a reasonable degree of effect on the conversion. Hence, all seven uncertain parameters were sampled to calculate the stochastic solution in our model. Ten sample points were then generated for each of QH0 and QN0 and four samples were generated for the other variables, since the decision variables affect the conversion most. The probability distribution was divided into the same number of regions (10 or 4) that have the same area, and one point was sampled from each region; the total number of sampled points was 102 400. 5.2. Stochastic Optimization. The PSO was implemented to find the optimal binding energies. First, 20 random particle pairs (or swarms) of QH0 and QN0 were generated over the ranges 100-175 kcal/mol and 45-75 kcal/mol for QN0 and QH0, respectively, and were updated in each iteration 3943



ARTICLE

Table 7. Stochastic and Deterministic Optimal Solutions result

stochastic

deterministic

QH0 [kcal/mol]

64.38

64.32

QN0 [kcal/mol]

126.04

122.83

Table 8. Comparison of Stochastic and Deterministic Optimization Schemes result

stochastic solution

deterministic solution

average conversion

0.2338

0.1969

maximum conversion minimum conversion

1.0000 0.0074

1.0000 0.0105

standard deviation

0.1159

0.1236

Figure 3. Maximum conversion during iterations. Figure 5. The 3D plot of the conversion distributions with a fresh set of scenarios.

optimization are given as follows: deterministic: Q d ¼ argQ max x½Q , EðεÞ

ð24Þ

Q s ¼ argQ max E½xðQ , εÞ

ð25Þ

stochastic:

Figure 4. Values of decision variables during PSO iterations. The dotted and solid lines indicate the points of QN0 and QH0, respectively.

step. The stochastic and deterministic optimal solutions are presented in Table 7. Figure 3 and Figure 4 show the maximum conversion and the decision variables, QH0 and QN0, at each iteration. To verify the efficacy of stochastic optimization, the solutions of deterministic and stochastic problems were tested again with a fresh set of four16 sample points uniformly sampled in the uncertain parameter space for comparison. This means that four data points were generated for each uncertain parameter. The objective functions for the stochastic and deterministic

As presented in Table 8, the average conversion of the stochastic optimal solution is greater than that of the deterministic optimal solution. Figures 5 and 6 show the results of stochastic and deterministic optimization formulations under a fresh set of scenarios. The deterministic problem ignores the uncertainties, and the maximum conversion is biased downward. The benefit of the stochastic solution can be quantified using the concept of “value of stochastic solution” (VSS): VSS ¼ max Q Eε ½xðQ , εÞ - Eε ½xðQ d , εÞ ¼ 0:2338 - 0:1969 ¼ 0:0369 A VSS value of 0.0369 implies that the average performance based on the stochastic solution is improved by about 20% compared to that of the deterministic solution. To confirm whether VSS represents a meaningful difference between two solutions, a t-test was performed. The t-test assesses whether the means of two groups are statistically different from each other. 3944



ARTICLE

Figure 6. Contour plot of the conversion distributions with a fresh set of scenarios.

6. CONCLUSIONS In this work, an efficient stochastic optimization method is proposed to design robust catalyst properties based on microkinetic models for improving the conversion of ammonia decomposition. Uncertainties associated with experimental determination or kinetic parameters are represented as exogenous uncertain parameters with assumed probability distributions. Since it is difficult to implement conventional derivative-based optimizers due to the complex constraints and numerical integration needed for function evaluation, a sample average approximation (SAA) scheme is proposed to evaluate the expected conversion. PSO was also adopted to solve the maximization problem for the optimal binding energies. These optimal binding energy values can be further mapped to the actual design of the catalyst material at a later stage. The large number of uncertain parameters can render the sample-based approach computationally unmanageable. Moreover, uniformity in sampling should be ensured. To overcome these difficulties, an efficient sampling approach based on thermodynamic constraints, linearity analysis, PLS, and constrained sampling is developed. In this work, 7 out of 16 uncertain parameters were sampled in obtaining the stochastic solution, while the values of four pre-exponential factors of reverse reactions were determined by thermodynamic constraints. Five variables, which have linear or negligible effect on the conversion, were replaced with their mean values. In order to verify the proposed approach, a fresh set of scenarios including all uncertain parameters in the model was tested. The results show that the solution of stochastic optimization provides more robust results than that of the deterministic one. Rigorous statistical testing showed that a significant improvement was achieved with stochastic optimization. ’ AUTHOR INFORMATION Corresponding Author

Figure 7. Difference between the results of stochastic and deterministic solutions under a fresh set of scenarios. The positive values of the y-axis indicate the improvement using stochastic optimization.

In case of paired and large samples like this case, the paired t-test should be used. The t-value used is given by t ¼

VSS - H0 pffiffiffiffi < t0:01=2 sd = N df ¼ n - 1

ð26Þ ð27Þ

where H0 is the hypothesis value and sd is the standard deviation of difference between the results on the stochastic and deterministic solutions. N and df are the number of scenarios and degrees of freedom. The T-score t0.01/2 is 2.58 at 99% confidence level and df = 416 - 1. To have an acceptable hypothesis value, the computed t-value should be smaller than t0.01/2. The minimum value of H0 under the given conditions is 0.0369-2.742 10-6. Because 2.742 10-6 is almost 0, it can be concluded that the VSS represents a significant improvement using the stochastic optimization. It should be noted that a large number of scenarios were used for evaluating the expectation in the testing stage, and thus, the comparison is comprehensive. Figure 7 shows the difference between the conversions of stochastic and deterministic cases over all scenarios.

*E-mail: [email protected].

’ ACKNOWLEDGMENT This work was supported by a National Research Foundation of Korea grant funded by the Korean Government [NRF-2009352-D00049]. ’ REFERENCES (1) Andersson, M. P.; Bligaard, T.; Kustov, A.; Larsen, K. E.; Greeley, J.; Johannessen, T.; Christensen, C. H.; Nørskov, J. K. Toward computational screening in heterogeneous catalysis: Pareto-optimal methanation catalysts. J. Catal. 2006, 239, 501–506. (2) Jacobsen, C. J. H.; Dahl, S.; Boisen, A.; Clausen, B. S.; Topsøe, H.; Logadottir, A.; Nørskov, J. K. Optimal catalyst curves: connecting density functional theory calculations with industrial reactor design and catalyst selection. J. Catal. 2002, 25, 382–387. (3) Prasad, V.; Karim, A. M.; Ulissi, Z.; Zagrobelny, M.; Vlachos, D. G. High throughput multiscale modeling for design of experiments, catalysts, and reactors: Application to hydrogen production from ammonia. Chem. Eng. Sci. 2010, 65, 240–246. (4) Bligaard, T.; Nørskov, J. K.; Dahl, S.; Matthiesen, J.; Christensen, C. H.; Sehested, J. The Brønsted-Evans-Polanyi relation and the volcano curve in heterogeneous catalysis. J. Catal. 2004, 224, 206–217. (5) Halemane, K. P.; Grossmann, I. E. Optimal process design under uncertainty. AIChE J. 1983, 29, 425–433. 3945



ARTICLE

(6) Maranas, C. D. Optimal molecular design under property prediction uncertainty. AIChE J. 1997, 43, 1250–1264. (7) Tayal, M.; Diwekar, U. Novel sampling approach to optimal molecular design under uncertainty. AIChE J. 2001, 47, 609–628. (8) Shastri, Y.; Diwekar, U. An efficient algorithm for large scale stochastic nonlinear programming problems. Comput. Chem. Eng. 2006, 30, 864–877. (9) Deshmukh, S. R.; Mhadeshwar, A. B.; Vlachos, D. G. Microreactor modeling for hydrogen production from ammonia decomposition on ruthenium. Ind. Eng. Chem. Res. 2004, 43, 2986–2999. (10) Prasad, V.; Karim, A. M.; Arya, A.; Vlachos, D. G. Assessment of overall rate expressions and multiscale, microkinetic model uniqueness via experimental data injection: Ammonia decomposition on Ru/γAl2O3 for hydrogen production. Ind. Eng. Chem. Res. 2009, 48, 5255–5265. (11) Sellers, H.; Shustorovich, E. Intrinsic activation barriers and coadsorption effects for reactions on metal surfaces: unified formalism within the UBI-QEP approach. Surf. Sci. 2002, 504, 167–182. (12) Fogler, H. S. Elements of Chemical Reaction Engineering, 2nd ed.; Prentice-Hall: Upper Saddle River, NJ, 1999. (13) Casey, M. S.; Sen, S. The scenario generation algorithm for multistage stochastic linear programming. Math. Oper. Res. 2005, 30, 615–631. (14) Kennedy, J.; Eberhart, R. Particle swarm optimization. Proceedings of IEEE International Conference on Neural Networks; IEEE: Piscataway, NJ, 1995. (15) Schwaab, M.; Biscaia, E. C.; Monteiro, J. L.; Pinto, J. C. Nonlinear parameter estimation through particle swarm optimization. Chem. Eng. Sci. 2008, 63, 1542–1552. (16) Fu, Y.; Diwekar, U. M. An efficient sampling approach to multiobjective optimization. Ann. Oper. Res. 2004, 132, 109–134. (17) Kalagnanam, J. R.; Diwekar, U. M. An efficient sampling technique for off-line quality control. Technometrics 1997, 39, 308–319. (18) Mhadeshwar, A. B.; Wang, H.; Vlachos, D. G. Thermodynamic consistency during microkinetic development of surface reaction mechanisms. J. Phys. Chem. B 2003, 107, 12721–12733. (19) Kleywegt, A. J.; Shapiro, A. Stochastic optimization. In Handbook of Industrial Engineering, 3rd ed.; Salvendy, G., Ed.; John Wiley & Sons: New York, 2001; pp 2625-2649. (20) Zabalza-Mezghani, I.; Manceau, E.; Feraille, M.; Jourdan, A. Uncertainty management: From geological scenarios to production scheme optimization. J.Pet. Sci. Eng. 2004, 44, 11–25.

3946


Stochastic Nonlinear Optimization for Robust Design of Catalysts

Recommend Documents