Multiscale Model and Informatics-Based Optimal Design of Experiments

Ind. Eng. Chem. Res. 2008, 47, 6555–6567

6555

Multiscale Model and Informatics-Based Optimal Design of Experiments: Application to the Catalytic Decomposition of Ammonia on Ruthenium Vinay Prasad and Dionisios G. Vlachos* Department of Chemical Engineering, Center for Catalytic Science and Technology, UniVersity of Delaware, Newark, Delaware 19716-3110

Fundamental multiscale models are increasingly being used to describe complex systems. Microkinetic models, which consider a detailed surface reaction mechanism containing all relevant reactions, are a prototypical multiscale model example. The computational effort in calculating all parameters of a multiscale model for real systems from first principles is prohibitive, and parameter uncertainty still limits full quantitative capabilities of these models. This motivates the development of rational model-based techniques in order to refine uncertain parameters and assess the global (in the entire experimental parameter space) model robustness. Herein we describe physics-aided methods (sensitivity, partial equilibrium, and most abundant reactive intermediate analyses) and statistics-based methods (A, D, and E optimal designs) for the design of experiments. While our methods are fairly general, we demonstrate them for the catalytic decomposition of ammonia on ruthenium to produce hydrogen. A global Monte Carlo method is used to search the operating space to generate possible optimal operating conditions for experiments. Our analysis illustrates that the D optimal and sensitivitybased designs are most promising and generate conditions that delineate important chemistry. It is shown that a standard design around the D optimal point may not be useful for highly nonlinear problems. Instead, informatics methods are proposed to identify optimal regions of the operating space. It is found that the experiments conducted within these regions have a high probability of providing useful kinetics information. It is also shown that the overall direction of the reaction (ammonia decomposition vs synthesis) and the macroenvironment (type of reactor) significantly affect the optimal design. This demonstrates for the first time the effect of macroscopic scales on microscopic ones with important implications for optimal design and product design. Introduction Empirical approaches have traditionally been used in the modeling of complex systems, with parameters being estimated from experimental data. Recently, there has been a paradigm shift toward using first-principles-based models. Many of these models are multiscale in nature, show great promise in providing physical insights, and are able to predict trends or even be quantitative.1–9 However, their use comes at the cost of increased number of parameters along with an uncertainty associated with the parameter values and the inclusion of all relevant physics. An important research direction, therefore, is the assessment of the reliability and accuracy of these complex, multiscale models.1,6,8,10–12 One way to tackle this problem is to judiciously perform experiments that allow us to refine model parameters. The design of experiments (DOE) methodology usually operates with empirical models; development of DOE methods suitable for multiscale models is in its infancy. Our thesis is that despite multiscale models being computationally intensive, they are usually less expensive than experiments and allow for a global search to identify interesting characteristics and assess the reliability of a model for system tasks, such as process optimization. Heterogeneous catalytic systems are a prototypical example of such complex, multiscale systems.7 In this work, we develop techniques for model-based optimal DOE in heterogeneous catalysis. Model-based optimal DOE usually generates conditions that provide the “best” chance (using statistical criteria) of obtaining accurate estimates of parameters, with minimal experimental effort. For fundamental multiscale models, another * To whom correspondence should be addressed. E-mail: vlachos@ udel.edu. Tel.: 302-831-2830.

approach is to use physics or chemical relevant information, inherently built into the model. A preliminary report on this approach was discussed in ref 12. We term this as physicsaided approach. Its essence lies in identification of conditions where the model’s response is qualitatively different, such as in bifurcation analysis. In reaction systems, the rate-determining step(s) (RDS), the most abundant reactive intermediate(s) (MARI) on the catalyst surface, or the distance from equilibrium (of reactions) can be some of the relevant criteria. In this work, we investigate whether purely statistical methods, such as the A, D, and E optimal designs, and physicsaided designs can provide such kinetically relevant experimental conditions. An important aspect of our approach is the introduction of informatics techniques to recognize underlying patterns in optimal regions of the operating space. We choose the ammonia decomposition reaction on γ-alumina-supported ruthenium catalyst to test our proposed scheme. There is an increasing interest in this reaction,13 since ammonia is nearly the best fuel in terms of its hydrogen storage capacity (>17% w/w). The methods developed, however, are general and may be used in the analysis of other microkinetic and general multiscale systems. Despite ammonia synthesis (the reverse reaction) being a wellstudied system due to its significance in fertilizers, there is still debate about the RDS and the MARI for this system. For example, the associative desorption of N2 or the N-H bond scission has been proposed to be the RDS in NH3 decomposition.14,15 Also, elemental hydrogen or nitrogen on the surface has been suggested or predicted, by previous models, as the MARI.15–17 A possible explanation for this apparent discrepancy could be that experiments were performed at different conditions (e.g., ultra high vacuum vs high pressure) and that the RDS and MARI

10.1021/ie800343s CCC: $40.75  2008 American Chemical Society Published on Web 07/16/2008

6556 Ind. Eng. Chem. Res., Vol. 47, No. 17, 2008

change with experimental conditions. Indeed, we have recently found out that it is possible to reconcile some of the apparent discrepancies (e.g., low vs high catalyst activity) simply by conducting multiscale simulations at drastically different experimental conditions.18 Remaining discrepancies are probably caused by differences in the catalyst itself (shape, size; these geometric features affect the fraction of active B5 sites in the case of ammonia chemistry) or support; model inadequacy cannot be ruled out either. Thus, aside from demonstrating the performance of our DOE techniques, a collateral objective of our work is to clarify the confusion with respect to the RDS and MARI of the ammonia decomposition reaction. Due to microscopic reversibility, a mechanism should be applicable to both forward and reverse directions (in our case, for ammonia decomposition and synthesis). While studies usually focus on one direction of the overall reaction, the proposed approach explores diversity within the experimental parameter space and thermodynamic consistency to assess the validity of a microkinetic model in both directions without necessarily increasing the number of experiments that need to be conducted. Microkinetic and Reactor Models In this work, we use microkinetic models,19 which have recently been developed to describe the performance of catalysts.20–22 Microkinetic models consider a detailed reaction mechanism consisting of all relevant elementary reactions, and the reactor model is solved without any assumptions made about a RDS, the MARI, or partial equilibrium conditions. The rate constant of elementary reactions is expressed in the modified Arrhenius form ki ) A0i

( ) ( ) T Tref

βi

exp -

Ei , i ) 1, · · · , m RT

(1)

where m is the number of reactions. There are two major sets of parameters in microkinetic modelssthe activation energies and the pre-exponential factors. Density functional theory (DFT) simulations are recently being conducted to calculate (from first principles) the activation energies. Examples where models using DFT-based activation energies can predict rates within 1 order of magnitude or so are promising (e.g., ref 23). In general, this parameter estimation approach provides at best semiquantitative agreement and is limited to small reaction networks due to well-known limitations (uncertainty in parameters, materials gap, concentration effects, involved calculations of the Hessian matrix for computing pre-exponentials, computational demand for complex chemistries, etc.24). An alternative estimation approach is to use semiempirical techniques, such as the bondorder conservation (BOC) method. The approach used in our group (and employed in the mechanism used in this work) is hierarchical in nature and overcomes most of the aforementioned problems (for details, see refs 12, 22, and 25). The BOC method is used to generate activation energies, using heats of chemisorption as inputs (the latter are obtained from DFT calculations or experiments). DFT is subsequently used to refine the energetics by including surface coverage effects depending on the major coverages (MARI) determined from simulations. The pre-exponential factors A0i are crudely estimated using transitionstate theory.19 To improve the quantitative capabilities of a mechanism, some of the pre-exponential factors are often estimated by fitting the model predictions to experimental data, subject to thermodynamic consistency.21 There is a need to collect meaningful data to obtain good parameter estimates and to find the limits on the information that can be extracted from

Figure 1. Schematic of catalytic system for ammonia decomposition with reaction mechanism, input, and output variables. Table 1. Ranges of Operating Variables Considered in the Design of Experiments operating variable

min

max

scaling

T (K) P (atm) residence time (s) catalyst surface area/reactor volume (1/cm) inlet mole fraction H2 inlet mole fraction NH3 inlet mole fraction N2

500 0.1 0.05 150 0 0 0

1000 10 5.0 15000 1 1 1

linear logarithmic logarithmic logarithmic linear linear linear

the experimental data. We explore these aspects using the modelbased DOE. The methods described are general enough to account for the inclusion of activation energies (or bond indices in the BOC theory) and pre-exponential factors in the development of DOE. As in our previous work, herein we concentrate on the pre-exponential factors. For the majority of this paper, we model a continuous stirred tank reactor (CSTR) at steady state. At the end we briefly compare CSTR results to those of a plug flow reactor (PFR) model. The governing equations for the CSTR and some details of our system are given in the Supporting Information S_1. The reaction mechanism, parameters, input, and output variables are summarized in Figure 1. Overall Methodology Having chosen the parameters of interest, our methodological steps include specifying the bounds of the feasible operating space, choosing the method of searching the operating space, defining a metric for evaluation of optimality, and identification of parameters that can practically be deduced from the data. Next, we describe each of these steps. Feasible Experimental Operating Space and Search Method. For ammonia decomposition on ruthenium catalyst, using realistic operating conditions, we define the range of operating conditions shown in Table 1. There are two search strategies: direct optimization (using a global-like optimization algorithm, such as simulated annealing) or a global Monte Carlo (MC) search. The former identifies a single optimal point that, as we show below, appears to be unsuitable for kinetics models. Consequently, we choose a global MC search. The method involves generating an ensemble of operating conditions that adequately span the operating space in an unbiased manner. Related stochastic, possibly more efficient methods, such as the Latin hypercube method, have not been exploited herein. The

Ind. Eng. Chem. Res., Vol. 47, No. 17, 2008 6557

Supporting Information S_2 details the generation of unbiased sample sets, and Figure S1 shows sampling in composition space. Once ranges of operating conditions have been defined and samples generated, problem-dependent constraints can be imposed. The first constraint relates to points with conversion lower than 3% that are omitted from analysis, since measurements would be unreliable (more generally, there is a detection limit of instrumentation for all species)

where Y ) [y1; . . .; yn] is the vector of outputs (observables) and P ) [p1; . . .; pm] is the vector of parameters. Possible outputs include the ammonia conversion, outlet gas mole fractions, and the surface coverage of reactive intermediates. Statistical metrics rely on the FIM, which is defined in terms of the normalized sensitivity matrix S and the measurement covariance W as

C1 : conversion > 3% (2) For the A, D, and E optimal designs,26,27 the Fisher information matrix (FIM) provides a different measure of optimality28,29 for each design (see below). Herein it is constrained to be full rank so that all parameters are identifiable

In this work, we consider all measurements to have approximately the same percentage uncertainty and be uncorrelated with other measurements. This makes W a multiple of the identity matrix (the choice of W does not change the results for the optimal criteria). Three metrics were explored in this work. The D optimal design is found by maximizing the determinant of the FIM and is equivalent to minimizing the joint confidence region and covariance of the parameter estimates. The E optimal design is found by maximizing the minimum eigenvalue of the FIM and provides the maximum distance between the parameter estimate and a point in the confidence region. It minimizes the size of the major axis of the confidence region. The A optimal design is found by minimizing the trace of the inverse of the FIM; it maximizes the diagonal elements of the FIM, minimizes the dimensions of the enclosing box around the joint confidence region, and reduces correlations between parameters. Each of these metrics attempts to provide a single measure of optimality of the operating conditions with respect to joint estimation of all the parameters in the model. This is in contrast to the physicsaided approach, where a separate optimal condition is found for each parameter (see below). The objective functions for the A, D, and E optimal criteria are summarized below

(3) C2 : rank (FIM) ) full (6 in our model) In kinetics, a key issue is to identify conditions that are kinetically relevant: not equilibrium limited and without transport limitations. To obtain kinetic information, the overall reaction should not be at equilibrium. The equilibrium factor η assesses the distance from equilibrium and should satisfy 3

C3 : η ) 1 -

1 pH2pN2 , |η| > εequil Keq p2

(4)

NH3

Here, Keq is the equilibrium constant for the overall ammonia decomposition reaction. When the system approaches equilibrium, the equilibrium factor becomes zero. This criterion is rather crude when detailed reaction mechanisms are employed. Instead, the partial equilibrium (PE) index for each reaction pair can be introduced that is normalized around 0.5 (ideal PE condition). A reaction is not in PE when C4 : PE Index )

rforward and rbackward denote the rates of the forward and backward steps of an elementary reaction, respectively. Here, we take the threshold value εPE to be 0.03. Regarding transport and possible temperature gradients, while not explicitly considered herein, one can envision several ways of approaching this. For example, the optimal conditions could be checked for transport limitations and eliminated from analysis. Alternatively, one could use more complex reactor models whenever transport phenomena or temperature gradients are important. Overall, we set C5 : transport phenomena ⁄ temperature gradients being unimportant (6) Metrics for Optimal Designs. The optimality metrics described in this section are evaluated at each of the conditions in the ensemble generated in the global MC search. Statistical Design Metrics. The literature on DOE is vast, as reviewed in refs 30-32, We have chosen local methods based on parametric sensitivity of the outputs to define optimization metrics. The normalized sensitivity matrix is defined as

[

]

A optimal: min{trace(FIM)-1} E optimal: max{λmin(FIM)}

(5)

∂ln y1 ∂ln pm l ∂ln yn ∂ln pm

(8)

D optimal: max{det[FIM]}

rforward , |PE index - 0.5| > εPE rforward + rbackward

∂ln y1 ... ∂ln p1 l S) ∂ln yn ... ∂ln p1

FIM ) STW-1S

(7)

(9)

The final step in the optimal DOE is deciding on the (sub)set of parameters that can be identified for a given set of observables. A priori identifiability considers whether or not, given a particular model and a set of experimental inputs and outputs, it is possible to uniquely determine the model parameters with noise-free data. It is tested using the correlation matrix of the parameters. Practical identifiability considers the accuracy with which model parameters can be estimated, given the covariance of the measurements. This is calculated using the FIM and is the approach adopted in this work. Practical identifiability requires nonsingularity of the FIM; given a set of measured outputs, we choose the largest subset of parameters for which the FIM is nonsingular based on computing the FIM rank. Instead of using the rank, the condition number of the FIM could be used, which must be relatively low for practical identification. Physics-Aided Design Metrics. In the physics-aided design method, we identify ensembles of operating conditions, i.e., subsets of the input space or domain in the mapping from inputs to outputs, which exhibit qualitatively different responses. As an example using the RDS as a metric, in one ensemble the RDS is the N-H bond scission and in another, is the N2 desorption. As another example with the surface coverage as a metric, for some conditions the MARI is adsorbed hydrogen H*, whereas for others, it is N*. Identification of operating conditions with qualitatively different metrics ensures that sufficiently diverse conditions have been explored, and suitable experimental techniques and measurements can assess the global robustness of the model response.


In this work, we use the highest (normalized) sensitivity of conversion for each reaction pair (known also as the RDS for the mechanism). We collect all conditions (an ensemble) for which a particular reaction (say, NH2* + * ) NH* + H*) is the RDS. The optimal condition is the one with maximum sensitivity Xoptimal ) max X

d ln(conv(X)) dln pk

(10)

When the most dominant species is in sufficiently high coverage on the catalyst, it may be detected spectroscopically (e.g., by infrared spectroscopy). Thus, another chemistry-based measure is to identify conditions at which the MARI changes. It is expected that obtaining kinetically relevant information requires that at least some of the reactions are far from equilibrium. Here, we use the distance of reactions from equilibrium (eq 5) as a metric that we relate to all other designs rather than as a direct design metric. Note that all the techniques described earlier assume that the parameter values are sufficiently accurate that the corresponding calculated sensitivity matrix and the FIM are close to the true values. If this assumption is untrue, the entire process of DOE may need to be iterated after parameter refinement. Pattern Recognition in Experiment Design Using Informatics Methods. We include principal component analysis (PCA)33,34 and data clustering35 to find underlying patterns in the most important inputs (and their ranges) for DOE. A brief explanation of PCA is given in the Supporting Information S_3. Clustering is also used to find out underlying patterns in experimental conditions. In our work, common features in inputs that have similar performance (based on an optimality metric) are identified. In other words, we seek subregions (in the input space) that have similar values of the optimality metric. This enables us to find patterns in temperature, pressure, inlet compositions, etc., that lead to optimal experimental designs. Clustering provides a means of extracting the underlying connections between the inputs and optimal designs. It serves as a relatively inexpensive method of accurately constructing the response surface (based on our metric of choice) in the optimal regions of the operating space. We use the partitioning around medoids (PAM) algorithm for clustering; this is a partition-based (as opposed to hierarchical) clustering technique that is considered to be a better-behaved alternative to the popular k-means clustering procedure.35,36 We use a combined measure based on Euclidean distance, the silhouette coefficient (see Supporting Information S_4), which characterizes cohesion within clusters and separation between clusters. The difference of the average silhouette coefficients for the original data and randomly permuted data is used as a measure to choose the optimal number of clusters. Once clustering has identified regions in the operating space that score high on the optimality metric, probability distributions of occurrence in each cluster are calculated with respect to the scaled ranges of each of the operating variables. Patterns are identified using the regions within the operating space that exhibit high probability of occurrence of optimal points. For some of the variables and some of the clusters, distributions may extend over the entire space, indicating that a pattern does not really exist for those variables and clusters. Once optimal regions have been identified, next we check whether the observed patterns in these clustered regions have been generated purely by chance. This is done by hypothesis testing: the null hypothesis is that the observed pattern (of obtaining k optimal points in the region of interest within the

cluster, out of m, the total number of optimal points) can be generated purely by sampling from a random distribution, and the alternative hypothesis is that the pattern indicates a correlation between the region and optimality. The P-value is calculated for the hypothesis using the hypergeometric distribution37,38 (Supporting Information S_4). A low P-value indicates that there is a low probability that the null hypothesis is true; in this case, the pattern observed is statistically significant. Results This section describes the results of applying the model-based DOE technique to the ammonia decomposition system described earlier. The nominal parameter values are those of ref 12. Seven inputs (temperature, pressure, residence time, catalyst surface area per unit reactor volume, and the inlet mole fractions of H2, NH3, and N2) are considered (note that the surface area per unit reactor volume and the residence time are often combined but can be varied independently as done here). Note that if a statistical DOE procedure (factorial/response surface design) was carried out without relying on a microkinetic model, a sevenfactor central composite design would involve ∼92 experiments and would not necessarily provide any kinetic insight, since it would make assumptions about the response surface being quadratic (or a similar function). Output Selection and Convergence of Optimal Designs. We conduct our studies with an ensemble consisting of 39,200 points (sets of experimental conditions). First we choose the uncorrelated observables (outputs) for subsequent analysis. The observables considered include NH3 conversion, gas-phase mole fractions, and the coverage of surface species. We eliminate the surface species NH* and NH2* as outputs, since their values are negligibly small at all conditions. In most cases, H* is the MARI, and N* and NH3* are the other significant surface species. There are cases where N* is the MARI. PCA is performed on the remaining transformed outputs (zero mean and normalized to be between 0 and 1). Based on the PCA, we conclude that the NH3 conversion, the gas mole fractions of hydrogen and ammonia, and the surface coverage of H*, N*, and NH3* are the uncorrelated outputs, and these are the observables considered for the rest of the analysis. Note that the conversion of ammonia and its outlet mole fraction are uncorrelated since the inlet mole fractions of all species are randomly chosen variables. The details of the PCA are given in Supporting Information S_3 and Table S1. In general, the stronger principal components are associated with eigenvectors that correspond to the gasphase outputs (conversion, outlet mole fractions) rather than the surface species, indicating that the former outputs are more important in calculating optimality metrics. The strongest principal component lies in the direction of the conversion, and the next strongest principal component lies in the direction of the H2 outlet mole fraction. We have performed convergence analysis to ensure that a sufficiently large ensemble of operating conditions has been generated in the global MC search. This is done by looking at the evolution of various criteria as a function of ensemble size. These include the D optimal design metric (constraints C1 and C2 are used to keep or eliminate points), the maximum sensitivity coefficient of each reaction (constraint C1 is used), the probability of each reaction being the RDS, and the number of reactions in PE. All metrics do not change after a few thousand points, indicating that the MC search has possibly converged (for an example, see Supporting Information, Figure S2).


Figure 2. (a) NSC of conversion for each reaction vs percentage conversion. Conditions of maximum sensitivity are indicated by numbers in each subplot in the order (from left to right): temperature (K), pressure (atm), residence time (s), catalyst surface area per unit reactor volume (cm-1), and inlet mole fractions of H2, NH3, and N2. (b) Representative NSC for each reaction as the RDS. The magnitude of the NSCs differs substantially among reactions.

We have found that ∼45% of the conditions have reaction R1 (H2 adsorption) as RDS, and a significant percentage of cases have reaction R4: NH2* + * ) NH* + H* as the RDS (Figure S2, Supporting Information). In addition, a large fraction of conditions has all six reactions in PE (indicative of equilibrium of the overall reaction), and there is also a significant fraction (∼70%) of conditions with at least three reactions in PE. No conditions exist where zero, one, or two reactions are in PE. Despite searching the entire parameter space, it appears that the parameters of only a few reactions can be identified due to others being very fast (in PE) under all conditions.

Optimal Designs Based on Physics-Aided and Statistical Methods. Normalized sensitivity coefficients (NSCs) were calculated for the uncorrelated outputs. The sensitivity coefficients were obtained by perturbing the pre-exponential factors of associated forward and backward reactions simultaneously by the same factor (in a thermodynamically consistent manner). Representative plots of the sensitivity of the ammonia conversion with respect to the pre-exponential factors against percentage conversion are shown in Figure 2a. Ammonia conversion was picked here as a representative element of the multidimensional output vector to visualize (in 2D) these results. The


sensitivity of the NH3 outlet mass fraction against percentage conversion exhibits similar features (Supporting Information, Figure S3). The conversion is most sensitive to reaction R4, NH2* + * ) NH* + H*, whose sensitivity is much larger than those of other reactions under most conditions; however, there are cases at which each of the other reactions becomes the RDS. Figure 2b shows representative normalized sensitivity coefficients for conditions at which each reaction is the RDS. It is clear that it is possible to find conditions where each reaction is the most significant. However, it appears that there is a hierarchy in importance (based on the magnitude of the normalized sensitivity coefficient): the N-H bond scission (reaction R4: NH2* + * ) NH* + H*) is by far the most influential, followed by the nitrogen desorption step (reaction R2: N2 + 2* ) 2N*), and then by the rest of the reactions. Given the magnitude of the sensitivity coefficient, it is expected that only the parameters of reactions R2 and R4 can practically be estimated. Our results also rationalize the apparent confusion regarding the RDS in the NH3 decomposition on Ru regarding N2 desorption or the N-H bond scission; the RDS changes drastically with operating conditions. Importantly, when performing an experiment at random, as in our MC search, the probability of a reaction being important (Supporting Information, Figure S2) is unrelated to the hierarchy of importance of the chemistry (Figure 2a). Random search of the experimental parameter space leads to H2 desorption (reaction R1) being picked most frequently as the RDS, followed by the N-H bond scission (reaction R4), which is indeed the most influential reaction. When reactions other than R4 (NH2* + * ) NH* + H*) become the RDS, the sensitivity coefficients are very small. Given the hierarchy of importance on the one hand and the low probability of “hitting” conditions revealing important reactions, it becomes clear that a random (or possibly even a statistical) DOE of a few points may lead to a poor design. For the majority of the operating conditions, H* is the most abundant reactive intermediate; however, there are a few conditions for which N* is the MARI (the distributions of these conditions are shown in Supporting Information, Figure S4). When N* is the MARI, conditions are evenly spread over the entire operating space and occur very close to complete conversion, have very low sensitivities with respect to all the reactions, and the overall reaction is close to equilibrium. This indicates that a change in the MARI (from H* to N*) for this system, although potentially an indicator of a change in chemistry, does not provide kinetically relevant information. The FIM is of full rank ()6 in the NH3 example) and therefore nonsingular for many points of the ensemble, indicating that all six parameters are practically identifiable at these conditions. However, we have noticed that Matlab reports that the matrix is invertible while the determinant of the FIM is very small (see below). What this information basically reveals is that even though many of the points pass the constraints and the parameters are recorded as practically identifiable, the confidence intervals when the parameters are estimated will be large. For example, when all six parameters are considered, the D optimal metric has a maximum value of the order of 10-12. The confidence intervals are of the order of 104-109, with the largest confidence intervals being for the parameters of the reactions in PE. If the parameters for the three reactions in PE (reactions R1, R5, and R6 in most cases) are removed from the analysis, the D optimal metric has a maximum value of O(104), and the confidence intervals of the parameters are

O(10-1)-O(102). This means that the pre-exponentials can be estimated within 2 orders of magnitude. Keeping this in mind, below we occasionally use the FIM based on the three reactions R2, R3, and R4 whenever explicitly mentioned in the figure captions. For each reaction being the RDS, we identify the conditions of maximum sensitivity. A summary of the optimal conditions is given in Table S2 in the Supporting Information. Similarly, A, D, and E optimal designs were identified, and the optimal conditions are also reported in Supporting Information, Table S3. The information is also presented in Supporting Information, Figure S5. Note that the RDS criterion provides a different optimal condition for each reaction, whereas the A, D, and E criteria each provide a single set of conditions that is optimal for identifying all pre-exponential factors together. Our results indicate that each criterion gives a different set of optimal conditions. However, some characteristics are common to more than one criterion: low temperatures (60%) with more than 3 reactions in PE, and at least 2 reactions not in PE, respectively.

indicates that the D optimal design points correspond to the most kinetically relevant conditions. Identifying Optimal Regions Using Principal Component Analysis and Clustering Methods. Experimental constraints and accuracy may not allow us to conduct experiments exactly at the optimal condition. Without knowing the nature of the response surface in the vicinity of the optimal conditions, it is impossible to predict the precision with which the predicted

Figure 5. (a) D optimal metric (logarithmic scale) against the absolute values of the normalized sensitivity coefficients of reactions R2 and R4. (b) Distribution of the absolute values of the normalized sensitivity coefficient of reaction R4 for various ranges of the D optimal metric (determinant of the FIM).

optimal conditions need to be matched experimentally for optimality to be preserved. This motivates the identification of regions rather than a single point in operating space that perform well on the optimality criterion; such an approach provides experimental flexibility. In the next section, we compare traditional design based on the single optimal point and design using regions of “nearly” optimal conditions identified herein. The latter are defined as sets of conditions having a metric close to the single optimal point. In order to determine which inputs affect mostly the optimal conditions, we performed PCA on the input and the output space for each of the RDS and A, D, and E optimal conditions. Points at which metrics higher than a threshold are observed were chosen, and PCA was performed on the input and output space for this set. For the RDS cases, the PCA indicated that conversion and outlet H2 mole fraction are the outputs associated with the strongest principal components in most cases, and the residence time and inlet H2 mole fraction are the inputs with the strongest principal components. For the A, D, and E optimal conditions, the strongest principal components lie along the same output directions, but the catalyst surface area per unit volume is one of the strong input directions in all cases. Since PCA can only detect linear correlations between input and output directions, and the strong principal components indicate variables that are more important for optimal DOE, PCA does not quantitatively predict optimal regions in operating


Figure 6. Probability of occurrence of D optimal points vs dimensionless range of operating conditions for all nearly optimal points (a) and for each cluster (b). Reactions R1, R5, and R6 have been left out from the statistical design methods.

space for conducting experiments. We employ clustering analysis to detect underlying patterns in the optimal regions in the operating space. For the D optimal design, we choose points with the determinant of the FIM greater than a threshold value (we term these as nearly optimal conditions). We then cluster these using the PAM technique and select the optimal number of clusters. The average silhouette coefficient difference between the original and randomly permuted data does not change much with the number of clusters, and we choose the optimal number of clusters to be 4. Figure 6 shows the distributions of the nearly optimal points over the dimensionless ranges of the operating variables

(pressure, temperature, residence time, catalyst surface area per unit reactor volume, inlet mole fractions of hydrogen, ammonia, and nitrogen) for the D optimal conditions in each cluster. Based on these distributions, we extract regions in the input space of high probability of occurrence of D optimal points. We then repeat this exercise (selection based on the normalized sensitivity coefficient being greater than a threshold, clustering, and pattern extraction) for the conditions based on reaction R4, NH2* + * ) NH* + H* being the RDS; the results of this clustering are shown in the Supporting Information (Figure S6). While we plot the distributions with respect to each input, the distribution of the optimal points is in fact a


Figure 8. Probability of observing a certain number of reactions in partial equilibrium for the first D optimal region listed in Table S4 (Supporting Information) and the entire ensemble. Most of the D optimal points have three reactions in PE and only a very small fraction has all (six) reactions in PE, in stark contrast to randomly searching the experimental parameter space (entire ensemble).

Figure 7. (a) Probability of occurrence of D optimal points vs dimensionless range of operating conditions for the first D optimal region listed in Table S4 (Supporting Information). (b) Interpolated response surface for section of the D optimal region.

multivariate distribution. Taking Figure 6 as an example, this means that if we were to propose an optimal region based only on cluster 2 in the D optimal analysis, the region would be specified by low pressure, low inlet NH3 mole fraction, and high residence time. Low pressure by itself is insufficient to ensure optimality. We test the statistical significance of the patterns observed within the clusters by calculating the P-value in the manner described earlier. The P-values for all the patterns observed in the clusters are O(10-4) or lower, indicating that these patterns are statistically significant. We identify patterns of high probability for each input variable in each cluster and collect them together to propose two regions in the operating space that are optimal in the D optimal sense, and two regions that are optimal with respect to the reaction NH2* + * ) NH* + H* being the RDS. These regions are presented in the Supporting Information (Table S4). Assessing the Optimal Regions and Proposing Experiments. The next question we attempt to answer is the issue of sampling within the identified optimal regions. To accomplish this, we conduct a new MC search within the optimal region (with five times as many points as in the optimal region of the original ensemble) of operating conditions. These new MC sample points can be thought of as conducting experiments. We then identify operating conditions that perform at least as well (based on the optimality criterion) as the conditions that were identified as nearly optimal in the optimal region of the original ensemble. For brevity, we focus on the D optimal design and the first region in Table S4 in the Supporting Information (the results are qualitatively similar for the second D optimal region). Figure 7a shows the probability distribution of finding optimal points within the first D optimal region (corresponding to the first column of data in Table S4). Since the probability is nearly uniform for all inputs (pressure, temperature, etc.), one can

conduct experiments anywhere within the optimal region. Importantly, not every point in the optimal region meets the optimality criterion. Specifically, only ∼35% of the points were classified as optimal. Figure 7b shows the interpolated response surface of a projection (of part of the D optimal region) onto the pressure-temperature plane. The numerous local maxima correspond to nearly optimal points. Figure 8 shows the probability of observing a certain number of reactions in PE in the new ensemble as well as in the completely random search. The majority of the conditions have only three reactions (R1, R5, and R6) in PE. By comparison, a random search leads to experiments where all six reactions are in PE. Among the 35% nearly optimal points, there is an even higher probability of having three reactions in PE (not shown), while the majority of the other 65% points have all six reactions in PE. Figure 9a shows the distribution of the D optimal metric for random search (original ensemble) and for the (first) D optimal region found via clustering. The distribution of the original ensemble is bimodal, indicating experimental conditions having three reactions (peak of high D values) and six reactions (peak of low D values) in PE. The D optimal region clearly gives higher values of the determinant of the FIM, indicating the clear advantage in identifying such a region. Since it is almost impossible to dial in the conditions of the single D optimal point, a DOE is usually conducted around this point (in a small hypercube). A question raised then is how good is the sampling in this hypercube. For simplicity, all experimental variables are varied by the same, small fraction (e.g., 1%, 10%, etc.) from their optimal values, and a MC search is subsequently performed within this hypercube. Figure 9b shows the distribution of the D metric for three distances from the single D optimal point. While deviation by 1% from the D optimal point results in a distribution that is similar to that for the D optimal region, further deviation from the D optimal point results in a substantial reduction of the D metric. This behavior deteriorates very fast as the distance from the D optimal point increases and is not much different from the random search (full ensemble) shown in Figure 9a. Our results indicate that DOE has to be carried out extremely close to the optimal point. Given instrumentation limitations, this is though not easily achievable in most experiments. For example, 1% variation in reactor temperature requires very tight changes (of the order of a few degrees).


Figure 10. (a) Distribution of number of reactions in partial equilibrium and (b) probability of each reaction being the RDS in a CSTR. Results are shown for ammonia synthesis and decomposition in randomly searching the entire ensemble. Panels c and d show the corresponding results for the PFR.

Figure 9. Distribution of the FIM determinant for the full ensemble and the first D optimal region (a) and for various distances (see text) from the D optimal point (b).

The outcome from our work underscores the difficulty in DOE of complex multiscale problems, such as kinetic experiments, when design is based on the single optimal point. Our results clearly show that conducting experiments within optimal regions (instead of a single point), identified using informatics tools, has clear merit. Given the nonlinear nature of kinetics problems, the probability of hitting a nearly optimal point is only moderate (∼35%), indicating that one has to perform a moderately large number of experiments to adequately sample the optimal region. Ammonia Synthesis. Our studies so far have focused on ammonia decomposition (because of constraint C1). We also investigate the kinetic information obtained for ammonia synthesis (negative NH3 conversion in our calculations, 90% probability; Figure 10d), and there is a much smaller probability of finding conditions where other reactions are sensitive. Overall, this


Figure 11. Simulation results in a PFR. (a) Absolute values of NSCs of reactions R4 and R5 against the D optimal metric (logarithmic scale). (b) Absolute values of NSCs of reactions R2 and R5 against the D optimal metric (logarithmic scale). (c) Distribution of the FIM determinant for the full ensemble and the first D optimal region.

means that ammonia decomposition is a better reaction to run than synthesis for obtaining different important elementary steps in the PFR. On investigating further the importance of the contribution of R5 to the determination of the optimal region, panels a and b in Figure 11 show that the correlation between the D optimal metric and the NSC for R5 is not as pronounced as that for reactions R2 and R4, implying the latter are more influential reactions. Higher values of the D optimal metric are found at higher sensitivities for reactions R2 and R4 respectively, while low values of the NSC for reaction R5 are equally probable to have high values of the D optimal metric. Finally, Figure 11c shows that, for ammonia decomposition, the distribution of the D optimal metric in the optimal region is sharper (without the long tail at lower values of the determinant of the FIM seen in Figure 9a), which means that there is a higher probability (∼70%) of conducting a nearly optimal run in the D optimal region with the PFR. Details aside, it is clear that there are significant differences between the two reactors regarding DOE. This is an important finding because for the first time there is an indication that macroscopic-scale phenomena affect the most dominant chemistry at the microscopic scale. This has obvious ramifications for DOE. The difference in the relative importance of various elementary reactions in the two macroenvironments results in differences in optimal regions and characteristics. Our results highlight the intimate coupling between microscopic and

macroscopic phenomena: reactor design affects the important chemistry even in the absence of transport limitations. Conclusions We have investigated optimal experiment designs based on physics-aided (MARI, partial equilibrium, and maximization of sensitivity of each reaction) and statistical designs (A, D, and E optimal criteria) for multiscale microkinetic models. Informatics methods (PCA and clustering based on the method of partitioning around medoids) led to the identification of patterns in the input space and provided quantitative predictions of optimal regions in which to conduct experiments. Our framework was applied to the ammonia chemistry on the ruthenium catalyst. The main findings can be summarized as follows: (1) Different criteria of optimality provided different optimal operating conditions. The D optimal design has been found to provide kinetically most relevant information and is in very good agreement with the physics-aided sensitivity analysis design. (2) The MARI has been found not to be a useful physicsaided metric since conditions that cause a change in the MARI (from H* to N*) do not correspond to highly sensitive or kinetically relevant conditions for this reaction. (3) Hydrogen abstraction from NH2* (NH2* + * ) NH* + H*) and the dissociative adsorption of nitrogen (N2 +


(4)

(5)

(6)

(7)

(8)

(9)

2* ) 2N*) have been found to be kinetically relevant and highly sensitive elementary reactions for ammonia decomposition. However, the hierarchy of importance in reactions is not always reflected in the probability of generating conditions that have it as the RDS, stressing the need for optimal designs. Design of experiments around the D optimal point, which has been the traditional design approach, is far from ideal for highly nonlinear problems, unless we stay in a very small region around the D optimal point. Instead, identification of regions of nearly optimal conditions is a superior approach to design experiments and is strongly recommended. In this context, informatics methods are valuable. The optimal regions may be sampled randomly to generate conditions for experiments, since the optimal points within these regions are fairly uniformly distributed. Due to the highly nonlinear nature of kinetic problems, the probability of conducting a good experiment in the optimal region is less than 1. As a result, a larger number of experiments need to be carried out to ensure that optimal conditions are adequately sampled. Unconstrained sampling of the inlet composition space (aside from mass conservation) can provide kinetically relevant information for both directions of the overall reaction. For the ammonia synthesis and decomposition reactions, we have found that both overall reaction directions provide similar kinetic information, as a result of microscopic reversibility. While both overall reaction directions provide similar kinetic information, the probability of being equilibrium limited is lower in one direction, and this direction depends on the type of macroenvironment (reactor). The type of reactor also modifies the probability of finding a nearly optimal point in the optimal region identified using clustering. Overall, reactor and chemistry have to go hand by hand. Unconstrained inlet composition sampling is important to account for possible product inhibition encountered at higher conversions or near the exit of a chemical reactor and for sampling both reaction directions. For example, cofeeding of products and reactants even at low conversions reveals such product inhibition. This ensures that the kinetic model works over the entire composition space and for both the forward and backward reactions. The computational cost in the search is probably fairly network independent. However, that of solving the reactor model increases with the number of species. We expect that the method will be easily applicable to fairly large reaction networks, especially since it is embarrassingly parallelizable.

Acknowledgment We acknowledge financial support from the U.S. Department of Energy (DE-FG02-06ER15795) and the National Science Foundation (CBET-0651043). We also acknowledge Dr. Rajanikanth Vadigepalli of Thomas Jefferson University for useful discussions on clustering. Supporting Information Available: Details of the methodology and further illustration of results. This information is available free of charge via the Internet at http://pubs.acs.org.

Literature Cited (1) Vlachos, D. G. A review of multiscale analysis: Examples from systems biology, materials engineering, and other fluid-surface interacting systems. AdV. Chem. Eng. 2005, 30, 1. (2) Rodgers, S. T.; Jensen, K. F. Multiscale modeling of chemical vapor deposition. J. Appl. Phys. 1998, 83 (1), 524. (3) Karakasidis, T. E.; Charitidis, C. A. Multiscale modeling in nanomaterials science. Mater. Sci. Eng. C 2007, 27, 1082. (4) Broughton, J. Q.; Abraham, F. F.; Bernstein, N.; Kaxiras, E. Concurrent coupling of length scales: Methodology and application. Phys. ReV. B 1999, 60 (4), 2391. (5) Drews, T. O.; Webb, E. G.; Ma, D. L.; Alameda, J.; Braatz, R. D.; Alkire, R. C. Coupled mesoscale-continuum simulations of copper electrodeposition in a trench. AIChE J. 2004, 50 (1), 226. (6) Braatz, R. D.; Alkire, R. C.; Seebauer, E.; Rusli, E.; Gunawan, R.; Drews, T. O.; Li, X.; He, Y. Perspectives on the design and control of multiscale systems. J. Proc. Control 2006, 16, 193–204. (7) Raimondeau, S.; Vlachos, D. G. Recent developments on multiscale, hierarchical modeling of chemical reactors. Chem. Eng. J. 2002, 90 (1-2), 3. (8) Christofides, P. D. Control of nonlinear distributed process systems: recent developments and challenges: Perspective. AIChE J. 2001, 47 (3), 514. (9) Maroudas, D. Multiscale modeling of hard materials: Challenges and opportunities for chemical engineering. AIChE J. 2000, 46 (5), 878. (10) Rusli, E.; Drews, T. O.; Braatz, R. D. Systems analysis and design of dynamically coupled multiscale reactor simulation codes. Chem. Eng. Sci. 2004, 59, 5607. (11) Braatz, R. D.; Alkire, R. C.; Seebauer, E. G.; Drews, T. O.; Rusli, E.; Karulkar, M.; Xue, F.; Qin, Y.; Jung, M. Y. L.; Gunawan, R. A multiscale systems approach to microelectronic processes. Comput. Chem. Eng. 2006, 30, 1643–1656. (12) Vlachos, D. G.; Mhadeshwar, A. B.; Kaisare, N. Hierarchical multiscale model-based design of experiments, catalysts, and reactors for fuel processing. Comput. Chem. Eng. 2006, 30, 1712. (13) Thomas, G.; Parks, G. Potential Roles of Ammonia in a Hydrogen Economy. U.S. Department of Energy, 2006. (14) Bradford, M. C. J.; Fanning, P. E.; Vannice, M. A. Kinetics of NH3 decomposition over well dispersed Ru. J. Catal. 1997, 172, 479. (15) Tsai, W.; Weinberg, W. H. Steady-state decomposition of ammonia on the Ru (001) surface. J. Phys. Chem. 1987, 91, 5302. (16) Dahl, S.; Sehested, J.; Jacobsen, C. J. H.; Tornquist, E.; Chorkendorff, I. Surface science based microkinetic analysis of ammonia synthesis over ruthenium catalysts. J. Catal. 2000, 192, 391. (17) Dahl, S.; Logadottir, A.; Jacobsen, C. J. H.; Norskov, J. K. Electronic factors in catalysis: the volcano curve and the effect of promotion in catalytic ammonia synthesis. Appl. Catal. A: Gen. 2001, 222, 19. (18) Deshmukh, S.; Vlachos, D. G. A reduced mechanism for methane and one-step rate expressions for fuel-lean catalytic combustion of small alkanes on noble metals. Combust. Flame 2007, 149 (4), 366. (19) Dumesic, I. A.; Rud, D. F.; Aparicio, L. M.; Rekoske, J. E.; Revino, A. A. The Microkinetics of Heterogeneous Catalysis; American Chemical Society: Washington, DC, 1993. (20) Mhadeshwar, A. B.; Aghalayam, P.; Papavassiliou, V.; Vlachos, D. G. Surface reaction mechanism development for platinum-catalyzed oxidation of methane. Proc. Comb. Inst. 2002, 29, 997. (21) Mhadeshwar, A. B.; Wang, H.; Vlachos, D. G. Thermodynamic consistency in microkinetic development of surface reaction mechanisms. J. Phys. Chem. B 2003, 107, 12721. (22) Mhadeshwar, A. B.; Vlachos, D. G. Hierarchical multiscale mechanism development for methane partial oxidation and reforming and for thermal decomposition of oxygenates on Rh. J. Phys. Chem. B 2005, 109, 16819. (23) Honkala, K.; Hellman, A.; Remediakis, I. N.; Logadottir, A.; Carlsson, A.; Dahl, S.; Christensen, C. H.; Norskov, J. K. Ammonia synthesis from first-principles calculations. Science 2005, 307 (5709), 555. (24) Aghalayam, P.; Park, Y. K.; Vlachos, D. G. Construction and optimization of detailed surface reaction mechanisms. AIChE J. 2000, 46 (10), 2017. (25) Mhadeshwar, A. B.; Vlachos, D. G. A catalytic reaction mechanism for methane partial oxidation at short contact times, reforming, and combustion, and for oxygenate decomposition and oxidation on Pt. Ind. Eng. Chem. Res. 2007, 46, 5310. (26) Box, M. J.; Draper, N. R. Factorial designs, the X′X criterion, and some related matters. Technometrics 1971, 13, 731. (27) Liang, Y.-z.; Fang, K.-t.; Xu, Q.-s. Uniform design and its applications in chemistry and chemical engineering. Chemom. Intell. Lab. Syst. 2001, 58, 43.

Ind. Eng. Chem. Res., Vol. 47, No. 17, 2008 6567 (28) Franceschini, G.; Macchietto, S. Validation of a model for biodiesel production through model-based experiment design. Ind. Eng. Chem. Res. 2007, 46, 220. (29) Versyck, K. J.; Claes, J. E.; Van Impe, J. F. Practical identification of unstructured growth kinetics by application of optimal experimental design. Biotechnol. Prog. 1997, 13, 524. (30) Atkinson, A. C.; Donev, A. N. Optimum Experimental Designs; Clarendon Press: New York, 1992. (31) Fedorov, V. V. Theory of Optimal Experiments. Academic Press: 1972. (32) Walter, E.; Pronzato, L. Identification of Parametric Models for Experimental Data; Springer: New York, 1997. (33) Edward, J. J. A User’s Guide to Principal Components. Wiley: New York, 1991. (34) Vajda, S.; Valko, P.; Turanyi, T. Principal component analysis of kinetic models. Int. J. Chem. Kinet. 1985, 17 (1), 55.

(35) Kaufman, L.; Rousseeuw, P. J. Finding Groups in Data: An Introduction to Cluster Analysis; Wiley: New York, 1990. (36) Pearson, R. K. Computational negative controls for cluster analysis. Technical report, Thomas Jefferson University, 2003. (37) Kurdistani, S. K.; Tavazoie, S.; Grunstein, M. Mapping global histone acetylation patterns to gene expression. Cell 2004, 117, 721. (38) Tavazoie, S.; Hughes, J. D.; Campbell, M. J.; Cho, R. J.; Church, G. M. Systematic determination of genetic network architecture. Nat. Genet. 1999, 22, 281.

ReceiVed for reView March 1, 2008 ReVised manuscript receiVed May 5, 2008 Accepted June 2, 2008 IE800343S

Multiscale Model and Informatics-Based Optimal Design of Experiments

Recommend Documents