Comparison of Sampling Strategies for Gaussian Process Models

Ind. Eng. Chem. Res. 2011, 50, 1379–1388

1379

Comparison of Sampling Strategies for Gaussian Process Models, with Application to Nanoparticle Dynamics Andres F. Hernandez and Martha A. Grover* School of Chemical & Biomolecular Engineering, Georgia Institute of Technology, Atlanta, Georgia 30332-0100, United States

Gaussian process modeling (GPM) is widely used in engineering design as a surrogate model for more complex and computationally demanding simulations. However, few fundamental studies focus on GPM as a surrogate for stochastic or high-dimensional dynamic systems. This paper provides a critical evaluation of selection strategies for sample collection in GPM from a stochastic dynamic multivariable simulation. Since the simulation may have high computational demands, it is important that the finite number of sample locations be carefully designed to be highly informative for building the GPM. Because the simulations are stochastic, one must determine whether repetitions at each point are necessary to accurately capture the noise behavior, or rather if the points should be spread out over the sampling space. A case study based on nanoparticle synthesis is used here to motivate and illustrate the GPM approach and several options and challenges for sample selection in the context of dynamic systems modeling. Our results show that the use of more repetitions does not provide a clear benefit when the total number of simulation samples is fixed. Our work also shows that including trend functions in GPM improves the prediction accuracy over a single discrete time step. However, when GPM is used recursively to predict dynamic trajectories, the trend functions may instead lead to extrapolation of the dynamic trajectory creating larger errors. 1. Introduction Design of experiments (DOE) encompasses an array of statistical methods for selecting experiments to be performed. Some common design objectives in DOE include accurate identification of parameters in a model, discrimination between multiple models, and location of optimal operating points.1 The newer field of design and analysis of computer experiments (DACE)2-4 bears many resemblances to DOE. In both cases, the goal is to select a set of sample points, but in DACE these samples are collected from a computer simulation, not a physical experiment. Typically a DACE sample set is used to construct an empirical model of the simulation, which approximates the input-output mapping of the system, but with a simpler underlying mathematical structure requiring reduced computation.5 DACE sample sets have also been used for process optimization.6 However, there are also key differences between DOE and DACE. DACE has been developed primarily for deterministic computer simulations, so there is no need to perform replicates of a particular sample point, as in DOE. (Each replicate would yield an identical result, providing no new information.) Even though each sample from a computer simulation may have significant cost of time and resources, it is generally assumed that each computer sample is less costly than a physical experiment. As a result, a common approach in deterministic DACE is to perform a space-filling design, where the sample points are selected to cover the input sampling space for accurate mean prediction. The sample points in space-filling designs can be chosen according to a mathematical criterion such as maximizing the minimum distance between points, or alternatively by following a uniform grid spacing or a randomized sampling approach like Latin hypercube. Space-filling designs typically require significantly more samples than would be practical for DOE. DOE methods are often used to identify * To whom correspondence should be addressed. E-mail: martha. [email protected].

parameters in a simple polynomial model of the response, while a space-filling design enables black-box modeling of more complicated nonlinear behavior, via, for example, higher order regressors or nonparameteric models.2 Gaussian process modeling (GPM) is a nonparametric approach that is commonly used for metamodeling in DACE. It was originally developed in the geostatistics community for empirical modeling of actual sampled data from the field, under the name “kriging,” but it has a number of features that make it attractive for DACE. The main characteristic of GPM for DACE is its nature as an exact interpolator of the data, a desirable property when the sampled data is deterministic. GPM is based on linear regression and thus has a clear statistical foundation. It also combines aspects of global and local modeling. Thus, it is flexible enough to describe any smooth function via interpolation, given a sufficient amount of data. However, when fewer samples are available, GPM is also able to use regressors to predict overall global trends in a more parsimonius manner. Although GPM was originally created to model stochastic geostatistics data, it has been developed for metamodeling and DACE in the context of deterministic data. However, computer simulations may also be stochastic, and Kleijnen and co-workers have created implementations for GPM metamodeling with stochastic simulations.7 Multiple realizations of the same sample point (replicates) are collected and averaged, and then this sample mean is treated as a deterministic outcome to build a GPM metamodel. Using this approach, sequential sampling was performed on a queueing simulation, by selecting each new sample at the point with maximum prediction variance.8 This method was shown to be more effective than a standard spacefilling design using the Latin hypercube method. Kleijnen and co-workers have also used GPM for constrained optimization of a stochastic simulation.9 The variance due to random uncorrelated noise has been included within the local variance function of the GPM in kriging, via the “nugget” term, but a

10.1021/ie1007954  2011 American Chemical Society Published on Web 12/28/2010

1380

Ind. Eng. Chem. Res., Vol. 50, No. 3, 2011

Figure 1. (a) Platinum nanoparticle synthesis under thermal decomposition in an inert atmosphere. (b) Transmission electron microscopy image from a sc-CO2 process for Pt nanoparticles on carbon nanotubes (Image by Dr. Galit Levitin, Georgia Institute of Technology). The scale bar is 10 nm.

rigorous approach to model uncorrelated noise via the correlation function had not been applied in DACE prior to our recent work.10 Only recently has the GPM approach been applied to the field of dynamic systems modeling. Black box modeling of nonlinear dynamics systems is an extremely challenging but important problem. An alternative approach to black-box modeling is provided by methods such as equation-free computing and in situ adaptive tabulation. Instead of building an explicit model, simulations can be run as needed to evaluate derivatives11 or stored in a table for future reuse.12,13 The equation-free approach has been applied for a number of tasks, including the prediction of distributions via moment approximations14 and the dynamic propagation of uncertainty due to parameter uncertainty.15 Despite the flexibility of equation-free computing, there are a number of advantages to building an explicit model. This model can be built ahead of time, such that no additional simulations are performed at the prediction stage. Additionally, the metamodel once built can be validated and verified, so that it can be used with confidence in control and design. Volterra kernel modeling can be effective for local nonlinear dynamics with fading memory,16 but is not suited to batch processes in which the system state changes significantly over time.17 Neural nets have been employed for black-box dynamic modeling, but as global function approximators they can also be unreliable or difficult to train.18 ARMAX time-series models provide a simple, effective, and compact method for black-box dynamic modeling, often following a simple linear structure with several delay terms in the input and output. The ARMAX approach has been extended to GPM by Rasmussen19,20simplementing a recursive formulation of the GPM to approximate the discretetime dynamics of the state variable. Applications include gas-liquid separation,21 building temperature control using weather forecasting data,22 and filtering for estimation.23 All these dynamic studies use the GPM to predict the evolution of a single dynamic variable, such as temperature in the forecasting application.22 Although the actual system may be of significantly higher dimension (e.g., weather dynamics), the delay terms can be used as a system “memory,” to capture higher-order dynamics via a multistep model. An alternative approach is considered here, in which a higher state dimension is used, with no delay terms. Thus, the resulting GPM has the Markov property, such that the current state is sufficient to predict the future dynamics. While multioutput GPM’s have been studied in the static case, it is still not clear how to select the functional formation for the covariance matrix between the multiple outputs.24 A GPM was used recently to model the dynamics of a two-state system, the inverted pendulum. This

model was then used to predict a single time step, for use in control via neurodynamic programming (NDP).25 Here Deisenroth and co-workers developed a GPM for multivariable state prediction. However, it is only implemented as a single-step prediction in the NDP, and therefore its prediction accuracy is not evaluated. The work presented in this paper is a critical evaluation of sampling strategies for the creation of the GPM for dynamic systems modeling. Specifically we address three key challenges and open questions in multidimensional dynamic GPM for metamodeling of stochastic simulations. Our case study on the synthesis of platinum nanoparticle catalysts on carbon nanotubes requires multiple states representing several species, as well as the nanoparticle size distribution. Owing to the expense of each sample, it is not practical to sample over a large multidimensional space, much of which may not be reachable from the initial condition set via the dynamics. The region of the state space relevant for prediction is not known a priori in the dynamics context, and this is a new challenge for DACE. We also evaluate the trade-off between performing multiple realizations (replicates), versus spreading out points over the design space. Finally, we consider the effect of using global regression functions in the dynamics context. 2. Application in Stochastic Dynamic Predictions for Platinum Nanoparticle Synthesis This section describes a case study on nanoparticle synthesis dynamics. The dynamic process model here is a mass action kinetics model of the supercritical fluid phase, coupled to a stochastic simulation for the surface kinetics. While a number of simplifying assumptions were made in constructing the model, it enables us to demonstrate and highlight many of the challenges and opportunities for GPM as an approximate model for nanoscale phenomena. The process under consideration here is the deposition of platinum nanoparticles on carbon nanotubes in supercritical carbon dioxide (sc-CO2). Large-scale production of nanoparticles of controlled size are needed for a variety of applications, from fuel cells to drug delivery. The sc-CO2 process has been investigated recently for nanoparticle synthesis, with advantages including the high solubility of precursor in sc-CO2, the lack of organic solvents which are environmental hazards, the scalability to three-dimensional processing via porous supports or powders, and the ease of separation of the final product from the solvent (by lowering the pressure).26 A brief description of the process is provided in Figure 1a, as motivated by the procedure of Bayrakceken and co-workers.27


This thermal reduction process involves two stages: the adsorption of the platinum precursor onto the surface, and the subsequent growth of platinum nanoparticles on the surface. Specifically, the organometallic platinum precursor dimethyl(1,5-cyclooctadiene)platinum(II) (PtOL) is first solubilized in a high pressure sc-CO2 environment, so that it can be adsorbed by the functionalized multiwalled carbon nanotubes (CNT) which are also present in the chamber. Once the adsorption step is complete, the pressure is released, and N2 gas flows into the chamber at an elevated temperature. The platinum precursor PtOL undergoes a thermal reduction in an inert N2 atmosphere, leaving isolated platinum atoms (Pt1) on the surface of the carbon nanotubes, and releasing the organic ligands (OL) back into the fluid phase. The platinum atoms then form nuclei on the surface which grow into larger nanoparticles. Many other variations on this process also exist, such as the use of hydrogen in a chemical reduction of the platinum precursor.28 Figure 1b is a transmission electron microscopy (TEM) image of platinum nanoparticles deposited on multiwalled carbon nanotubes by a sc-CO2 process. The typical nanoparticle size is 5-10 nm, with a visible distribution of sizes. From a practical standpoint, it is often not desirable to have a wide nanoparticle size distribution, since many applications demand a specific size for optimal performance. However, in practice there will always be some distribution of sizes, and a key manufacturing question is to determine what distribution will provide adequate performance while minimizing processing cost. 2.1. Description of Mechanistic Model. A coupled deterministic-stochastic model is presented here, motivated by the previous experiments27 and previous modeling studies.29-31 A system of ordinary differential equations (ODE) is constructed for the adsorption step and a kinetic Monte Carlo (kMC) simulation is used for the platinum reduction, nucleation, and nanoparticle growth. Owing to the sequential nature of the process, as shown in Figure 1a, the final conditions from the ODE adsorption simulation are used as initial conditions for the kMC growth simulation. The main assumption in the ODEkMC coupled model is the independence between the process steps, that is, nanoparticles are not formed during the adsorption step and platinum precursor is not adsorbed in the growth step. The kinetic and thermodynamic parameters in this model are listed in Table 1. The exact nanoscale mechanisms occurring during the process are not fully understood or quantified, and this is generally a challenge for the robust optimal processing of nanomaterials. For example, in this case study the release of the organic ligands might also occur during the adsorption phase, and there could also be nanoparticle nucleation occurring during adsorption. Despite the limitations of the existing process models, their construction and validation in conjunction with experiments encodes current understanding of the process and thus are needed for process engineering.26 The simulations begin with the specification of the precursor mass (mPtOL) and carbon nanotube mass (mCNT), in the ranges listed in Table 1, while the rest of the parameters in the table are constant. Before the ODE adsorption process is simulated, the solubility of the platinum precursor in sc-CO2 is evaluated, to determine the initial concentration of platinum in the supercritical fluid phase. Using reported experimental information for the solubility of the precursor in sc-CO2,33-35 a Chrastil model36 was fit and then used to compute the initial precursor concentration. The fitted Chrastil model is ln S ) 5.13 ln(Fsc-CO2(T, P)) -

9240 + 31.6 T

(1)

1381

Table 1. Model Parameters for Platinum Nanoparticles on Carbon Nanotubes Using sc-CO2 description

variable

value

Operating Conditions pressure27 temperature27 reactor volume32 precursor mass carbon nanotube mass adsorption time growth time

P T V mPtOL mCNT tads tgr

24.2 MPa 343 K 54 cm3 130-170 mg 130-170 mg 2h 2h

Adsorption Isotherm Langmuir adsorption constant27 Keq adsorption capacity27 Q0

0.299 g sc-CO2/mg PtOL 334 mg PtOL/g CNT

Kinetic Parameters adsorption desorption reduction nucleation growth

kads kdes kred knuc kgr

4.11 × 101 cm3/(mol PtOL · s) 3 × 10-4 1/s 1 × 10-3 1/s 1 × 101 cm3/(mol Pt1 · s) 1 × 106 cm3/(mol Pt1 · s)

Material Constants molecular weight of precursor molecular weight of active site

MWPtOL 333.34 g/mol MWR 45.0 g/mol

where S is the solubility limit at the process conditions in units of (mg PtOL)/(cm3 sc-CO2), and Fsc-CO2 is the density of the sc-CO2. The sc-CO2 density depends on the temperature and pressure and is computed using an equation of state.37 If the precursor density is lower than S, all platinum precursor can be solubilized in the fluid phase and the initial concentration of platinum is mPtOL/V, otherwise the maximum platinum precursor is limited by the solubility and the initial concentration is defined according to this limit. The adsorption-desorption dynamics of the precursor on the carbon nanotubes in the sc-CO2 fluid can be modeled as PtOL(f) + R(s) h PtOL-R(s) where PtOL(f) is the platinum precursor in the fluid phase, R(s) is an active site on the carbon nanotube support, and PtOLR(s) is the adsorbed platinum precursor on a carbon nanotube active site. The forward reaction is adsorption and the reverse reaction is desorption, with rate constants kads and kdes, respectively. Their ratio determines the Langmuir adsorption constant Keq that is reported in Table 1. The initial conditions for the platinum precursor concentration and active sites concentration must also be specified. After the solubility check described previously, the initial concentration of the platinum precursor in the sc-CO2 fluid is CPtOL,0 )

mPtOL MWPtOLV

(2)

where MWPtOL is the molecular weight of the platinum precursor and V is the reactor volume. The initial concentration of the available active sites is related to the initial amount of carbon nanotube mCNT in the system. Different functional groups (such as the carboxyl group-COOH) have been incorporated onto the CNTs to enhance their reactivity.38 We assume here that the adsorption of platinum precursor only occurs at these functionalized sites and not directly on other sites on the carbon nanotube surface; however, this is in fact another open mechanistic question in the literature.38 The initial concentration of active sites is thus dependent on the extent of functionalization of the CNT, and here we use a value of 4.5 w/w % of the carbon nanotube mass,

1382


based on typical commercially available values. Using this information, we can calculate the initial concentration of active sites on the CNT as follows: CR,0 )

0.045mCNT MWRV

(3)

where MWR is the molecular weight of the functional group, which is a carboxyl group here. With the initial conditions CPtOL,0, CR, 0 and the adsorption time tads, the ODE adsorption model simulates the process dynamics and in the end predicts the final concentration of adsorbed platinum precursor on the carbon nanotube surface, CPtOLR(t ) tads). This value for the adsorbed precursor is then used as the initial condition for the kMC growth model. The kMC growth model is evaluated using a stochastic simulationsin particular, the Gillespie algorithm39 is used to simulate individual chemical reactions as probabilistic events that occur at known rates. The kMC simulation includes three surface reactions on the CNT surface: thermal reduction of the adsorbed platinum precursor PtOL to release the organic ligands from the platinum atom, nucleation of a nanoparticle from two individual Pt atoms, and subsequent growth of platinum nanoparticles by the incorporation of additional Pt atoms. Specifically, these reactions take the form kred

PtOL-R(s)

98

2Pt1(s)

98

Ptn(s) + Pt1(s)

98

knuc

kgr

Pt1(s) + OL(f) + R(s) Pt2(s) Ptn+1(s),

(4) n ) 2, 3, ...

where Pt1(s) is an isolated platinum atom on the CNT surface, Ptn(s) is a nanoparticle with n platinum atoms, and OL(f) is the organic ligand released to the nitrogen fluid phase from the adsorbed platinum precursor. The kMC growth model assumes that platinum is released in its elemental state; two platinum atoms form a stable nucleus; a size-independent growth rate exists for the nanoparticles, and no aggregation between particles contributes to nanoparticle growth.31 The number of nanoparticles of each size (Ptn) thus defines the nanoparticle size distribution. The rate constants for these reactions are listed in Table 1 and they define the typical rates of reaction for the probabilistic events. Particular features of the distribution like mean and variance are related to the ratio between nucleation and growth rate,29 and their numerical values in Table 1 were estimated to fit the measured nanoparticle size distribution.27 Each kMC simulation proceeds until no single platinum atoms are on the CNT surface or until the growth time tgr is reached. The kMC system size, defined by the number of platinum atoms used to simulate the overall process, was 500 000 atoms. 3. Gaussian Process Models and Dynamic Systems Modeling 3.1. Mathematical Form of Gaussian Process Models. Gaussian process modeling (GPM) is an empirical modeling approach for generalized linear regression models. The key component of Gaussian process modeling is the explicit modeling of the distance-based correlation between the residuals of the output and its subsequent use in the prediction.40,41 This concept is not usually employed in system identification and

model reduction, and enables interpolation of the original training samples when they are noise-free. Here we use the Gaussian correlation function υ:Rd × Rd f R such that

[

V(si, sj) ) σc2 exp -

]

d (si,q - sj,q)2 1 + σu2δij 2 2 q)1 lq

∑

(5)

where si, sj ∈ Rd. This function includes a local correlation term with variance σc2, as well as an uncorrelated term with variance σu2 to model simulation noise. δij is the Kronecker delta, and θ ) [σu2, σc2, l12, ..., ld2,]T are the GPM parameters that control the features of the correlation. We chose this commonly used correlation function for its flexibility in identifying the local correlations in each input dimension, and furthermore because it is an infinitely differentiable function, allowing us to quantify the uncertainty in the GPM parameters once they are estimated (see section 3.2).20 A correlation matrix V(A, B) is defined element-wise according to Vij ) V(Ai, Bj), where Ai, Bj ∈ Rd×1 are columns of the matrices A and B. The GPM is an empirical model, and thus it must be constructed from sampled data. Here ns is the number of samples, and the matrix X ∈ Rd×ns contains the values of the input vector x ∈ Rd at each of the ns sample points. The corresponding output matrix is Y ∈ Rdy×ns. In most GPM applications, dy ) 1, corresponding to a scalar output, while in the multidimensional dynamic setting considered in section 3.2, dy ) d. We define the entire collection of input and output points as D ) {X, Y}. The typical case of dy ) 1 is presented first. Using the covariance function V in eq 5, a GPM predictive distribution for a new input x is expressed under the best linear unbiased estimator42 y ∼ N ( yˆ(x, D), σy2(x, D))

(6)

yˆ(x, D) ≡ E[y|x, D] ) hT(x)βˆ + VT(x, X)V-1(X, X)(YT - H(X)βˆ )

(7) σy2(x, D) ≡ Var[y|x, D] ) V(x, x) - [hT(x) VT(x, X)] × -1 0 h(x) HT(X) V(x, X) H(X) V(X, X)

[

][

]

(8)

where V(x, X) ∈ Rns×1 contains the correlation values between x and each of the columns in the matrix X in D, according to eq 5. The GPM combines local modeling via a correlation function with global modeling, based on p regressors with unknown coefficients βˆ ∈ Rp. The regressors are specified via the function h : Rd f Rp, and the matrix H ∈ Rns×p represents the settings in the sample set D, such that Hij ) hj(Xi). The p regression coefficients are estimated as βˆ (θ) ) (HTV-1(X, X)H)-1HTV-1(X, X)YT

(9)

which is the generalized least-squares estimator of the regression coefficients. Based on this mathematical formulation, three elements are needed to build a GPM: (1) The location of the ns sample points in the set D. Space-filling methods are typically used for DACE,2 although in a dynamic context, the region to be “filled” is not clearly defined. (2) A basis set of regression functions h(x). The most common regression function is a constant with a single βˆ coefficient, although other researchers have used polynomial regression functions in GPM models,7,43 showing improvements in the prediction compared to the constant GPM model. (3)


1383

matrix of the likelihood function, along with uncertainty estimates from the previous state covariance matrix via eq 15. Thanks to the continuous infinite-differentiable property of the Gaussian correlation function, we can analytically derive the expressions for the Fisher information matrix that otherwise would have not been possible. Similar to Girard and MurraySmith,47 the corrected GPM mean and variance prediction for a particular state variable i, zî, c[k + 1] and Σzi, c[k+1], are Figure 2. Flowchart of the recursive one-step-ahead prediction scheme using a GPM for dynamic systems modeling.

A parameter estimation method for θˆ in the covariance function. For the results presented here, GPM parameters are obtained with a restricted maximum likelihood estimator (RMLE),44 which is a likelihood function that accounts for the degrees of freedom associated with the regression functions. 3.2. Application of GPM in Dynamic Systems Modeling. An approximation for a dynamic system can be written based on the GPM structure presented in the previous section. A GPM can be used as an autoregressive model where the state prediction information is fed back and used as the input to predict the new state.21,25,45 This discrete-time model is based on a time step of ∆t, with t ) k∆t, k ∈ Z: z[k + 1] ) f(z[k])

k ) 0, 1, ...

(10)

where f: Rd f Rd. The corresponding approximate version for the mean state prediction is zˆ[k + 1] ) fˆ(zˆ[k], D)

k ) 0, 1, ...

(11)

where fˆ: Rd f Rd comprises d GPMs. The sample set becomes D ) {X, Y}, X, Y ∈ Rd×ns, with dy ) d such that Yi ) f(Xi)

(12)

where f is the original simulation in eq 10. Similar to in situ adaptive tabulation,12 the information in D is used to approximate the dynamics at the new value of z in a one-step-ahead prediction (Figure 2). Usually the recursive formulation uses the GPM mean prediction in eq 7 as the feedback information, but alternatively a stochastic realization from the predictive distribution in eq 6 could be used to generate an approximate stochastic realization. To implement the recursive formulation in eq 11 for multivariate systems, a GPM is constructed for each of the d state variables. Using eqs 6-8, we assume the distribution of the next state z[k + 1] to be a multinormal distribution with a state mean vector zˆ[k + 1] ∈ Rd and a state covariance matrix Σz(k+1) ∈ Rd×d as z[k + 1] ∼ MN(zˆ[k + 1], Σz[k+1]) zî[k + 1] ) ˆfi(zˆ[k], D) Σz[k+1],i,i ) σzi[k+1]2(z[k], D)

i ) 1, ..., d i ) 1, ..., d

(13) (14) (15)

Both the mean and the covariance are defined element-wise according the d scalar GPMs, with the covariance having all offdiagonal elements equal to zero. Fully multivariable methods for modeling the covariance have been explored under the name “cokriging,”24 but it is not clear what if any covariance function is suitable, in analogy to eq 5. It is possible to formulate a correction in the state prediction using a truncated Taylor series expansion46 that combines error estimation in the GPM parameters via the Fisher information

zî,c[k + 1] ) zî[k + 1] + 1 2

1 2

d+2 d+2

∑∑

Cov(θp, θq)

p)1 q)1 d d

∑ ∑ Cov(z , z ) p

q

p)1 q)1

∂2(zî[k + 1]) ∂zp∂zq

d+2 d+2

Σzc[k+1],i,i ) Σz[k+1],i,i +

∑∑

1 2 p)1 d

Cov(θp, θq)

q)1

d

∑ ∑ Cov(z , z )

1 2 p)1

p

q

q)1

∂2(zî[k + 1]) ∂θp∂θq

|

|

+

θ)θˆ

(16) z)zˆ[k]

∂2(Σz[k+1],i,i) ∂θp∂θq

∂2(Σz[k+1],i,i) ∂zp∂zq

|

|

+ θ)θˆ

(17) z)zˆ[k]

The first terms on the right-hand sides of eqs 16 and 17 are the GPM mean and variance predictions from eqs 14 and 15. The second term corrects the GPM prediction using a parameter covariance matrix Cov(θi, θj), which is calculated from the observed Fisher Information Matrix. The final term in eqs 16 and 17 is the correction for the state uncertainty. Here the state covariance matrix Cov(zp, zq) is simply Σz[k] from eq 15. 3.3. Implementation in Model Reduction of Nanoparticle Dynamics. In the nanoparticle case study considered here, the coupled ODE-kMC simulation is approximated by an ODEGPM simulation. The kMC simulation defines the “original” simulation according to eq 10, with ∆t ) 300 s. Five state variables are used in z to describe the growth simulation, that is, d ) 5. The first two states (z1 and z2) are the volumetric concentrations of platinum atoms (Pt1) and adsorbed platinum precursor (PtOL-R) on the CNT surface. The remaining three state variables describe the nanoparticle size distribution Ptn, using a reduced state created from a moments approximation.48,15 The basic equation for the moments of a distribution is ∞

mq )

∑n N q

n

(18)

n)2

where mq is the qth moment of the nanoparticle size distribution, and Nn is the number of nanoparticles of size n in the kMC simulation. The dynamics of the population balance defined by eq 4 are infinite dimensional, while a GPM must be built for each finite state that is modeled; therefore some approximation of the size distribution is required. If the nanoparticle distribution followed a normal distribution across the state space at any discrete time k, we would only need the first three moments to completely describe the distribution. Specifically, the moments can be used to compute the mean µ and variance σ2 of a normal distribution: m1 m0

(19)

m2 - µ2 m0

(20)

µ)

σ2 )

Here, we use these three moments to approximate the nanoparticle size distribution zˆ3-5[k] ) [m0, m1, m2] in eqs 18-20 as a normal distribution, based on our observation of

1384


Figure 3. Approach for approximate dynamic modeling using precollected information. (a) Creation of the dynamic region based on precollected dynamic trajectories. (b) Approximated dynamic trajectory using the information set D and the GPM in a recursive one-step-ahead prediction.

Figure 4. Approximated nanoparticle dynamic trajectory using a GPM reduced order model. mpre ) 155 mg, mCNT ) 140 mg: (a) concentration of elemental platinum on the CNT surface; (b) zero moment of the nanoparticle size distribution m0.

typical distributions such as the one shown in Figure 5b. Clearly the true distribution cannot be exactly normal, since negative nanoparticle sizes are not possible, so some error will be incurred with this reduction. At early times such as that shown in Figure 5a, the Gaussian distribution must be significantly truncated during sampling (n g 2), since here the spread of the distribution is large relative to the mean. The initial sample set D is constructed from trajectories of the original simulation, based on a prespecified design region for the initial conditions. Specifically, the initial precursor mass mPtOL and carbon nanotube mass mCNT may vary in the ranges described in Table 1, creating a rectangular design region for the initial conditions. The ODE-kMC simulation is then run for 7200 s, corresponding to 25 discrete time steps. Five precollected dynamic trajectories are simulated at the four corners and the center of the operating condition region. As illustrated by Figure 3, the set of operating conditions may be well-defined, but in a dynamic context that does not guarantee that dynamic region is well-defined, since the trajectories of this nonlinear stochastic system are not known ahead of time. An example of a predicted trajectory is shown in Figure 4. Two of the five dynamic states are plotted versus time, beginning from an initial state defined by a precursor mass of mpre ) 155 mg and a carbon nanotube mass of mCNT ) 140 mg. Since the reactor volume V is fixed, these two species along with the parameters in Table 1 define the initial state of the system and the dynamics. The trajectories in Figure 4 compare the results of five individual stochastic realizations of the ODE-kMC model (blue) to the prediction by the ODE-GPM model (red). Since the GPM predicts both the expected trajectory and the prediction variance, it is possible to plot both the mean value and the uncertainty. The dashed red line denotes yˆ ( z95(σy2)1/2 where z95 is the normal score for a 95% confidence interval, and yˆ is

the GPM output as in eq 6. In the particular trajectory of Figure 4a, the elemental surface platinum Pt1 exhibits little variability in the stochastic simulations and it is well predicted by the GPM, including the small variation among realizations. However, the prediction of the zero moment of the distribution in Figure 4b exhibits offset in the mean prediction at steady state. The nanoparticle size distribution is plotted in Figure 5, at two representative times. The blue bars are histograms coming directly from a single stochastic realization of the original simulation, using the same precursor and nanotube masses as in Figure 4. In contrast, the red bars are based on the reconstructed predictions from the GPM. Because the GPM only predicts the first three moments of the nanoparticle size distribution, these moments are then used to construct a Gaussian distribution. The red bars in Figure 5 are histograms resulting from stochastic samples from this Gaussian distribution, until it reaches the kMC system size of 81000 atoms, which is the mean value of adsorbed platinum atoms over the different operating conditions. Figure 5a corresponds to a time early in the process, and neither histogram appears smooth due to the small number of particles in the simulation domain. At the later time in Figure 5b, both histograms appear smoother and consistent with a Gaussian distribution. At both times, the mean and variance of the GPM histogram appears consistent with that of the original simulation. The plots in Figures 4 and 5 are illustrative examples, but to understand and quantify the overall accuracy of the GPM approach and its many options for implementation, we compute and compare the prediction error associated with all possible combinations shown in Table 2. All computations were performed in MATLAB. The average time to estimate the parameters in the GPM is 5-10 s for all d GPMs that are trained. This short parameter estimation time is


1385

Figure 5. A realization of the nanoparticle size distributions at two different process times: (a) t ) 900 s (b) t ) 7200 s. The GPM realization is drawn from a continuous distribution with mean and variance according to eqs 19 and 20. Table 2. Summary of Potential Combinations to Create Approximated Models, Evaluated for Local and Global Prediction Analysis settings

options

number of sample points repetition level sampling strategy

30, 45, 60, 75, 90, 105, 120 1, 3, 5 1. random selection of precollected points in dynamic trajectories 2. random walk in the state 1. constant function 2. linear functions of the state variables

regression functions

due in part to the small number of sample points that are used (Table 2 shows that in our work we are using ns ) 120 sample points in each GPM compare with 400 sample points in other GPM implementations.25) The restricted maximum likelihood estimation was performed using MATLAB’s fmincon function, with each range parameter constrained according to 10-5 < lq < 6 in the normalized coordinates. The initial guess for lq is 1, and there is no guarantee that the global minimum has been found, as is typical for GPM.49 The inversion of V(X, X) in the restricted maximum likelihood is typically a bottleneck in GPM, and here was implemented using a Cholesky decomposition.20 Since this work does not use a sequential DACE approach to improve the GPM models,10 we train the GPMs once for each sample set D from the ODE-kMC simulation. Thus, V-1(X, X) is only calculated once, prior to the prediction of all dynamic trajectories. The overall training time of the methodology was therefore limited by the data collection from the ODE-kMC (dynamic trajectories and additional repetitions) and not by the GPM calculations. 4. Results 4.1. Local Prediction Analysis. Before computing the prediction error associated with dynamic trajectories, we first consider the prediction error associated with a single evaluation of the GPM, which corresponds to the prediction for a single time step. For this purpose we define the local error metric LEM using normalized variables to eliminate scaling effects around the five state variables nm

nt

d

∑ ∑ [ ∑ (fj (z [0]) - ˆf (z [0], D )) ]

1 LEM ) nmnt m)1

j

i

j)1

j

i

2 0.5

m

i)1

(21) At nt test points for z, the mean prediction of the GPM, fˆ(z), is compared to a sample mean f(z) computed from nr realizations

Figure 6. Local error analysis LEMy for different repetition levels using points from precollected dynamic trajectories. Regression functions: (s) constant RF, (---) linear RF. RL ) repetition level.

of the original stochastic simulations. For any finite sample size nr, f(z) will fluctuate from the true mean, and this will of course limit the minimum value of the LEM that can be achieved. However, by using a fixed value of nr, we will still be able to use the LEM to compare the performance of various implementations of the GPM, as described in Table 2. The summation over m in eq 21 accounts for stochastic elements in the construction of a GPM. Because the samples in the set Dm are chosen randomly, the resulting GPM will also vary. Thus, nm otherwise identical GPMs are constructed for evaluation of the local error via eq 21. Each of these nm GPMs is then evaluated at each of the nt different test points for z. These are the two outer summations in eq 21. Figures 6 and 7 contain the LEM for the various implementation options given in Table 2. In all cases, nr ) 20 realizations are performed at each of the nt ) 960 values of z. The nt test samples used to calculate the LEM were generated from 40 dynamic trajectories, with initial conditions selected using a Latin hypercube design. These trajectories are different from those in the training set D that are used to construct the GPM. A total of nm ) 10 different GPMs are constructed, and parameter estimation for θ is performed using restricted maximum likelihood estimation (RMLE). For the comparison in Figure 6, the sample points in D used to construct the GPM are selected randomly, directly from the set of five precollected trajectories, while in Figure 7 the points in D are sampled randomly from a convex hull constructed around the set of precollected trajectories. Several trends can be seen in Figure 6. Overall the LEM decreases as the number of samples ns is increased, as expectedssince the GPM performs local function approxima-

1386


Figure 7. Local error analysis LEMy for different repetition levels using random points selected in the convex hull. Regression functions: (s) constant RF, (---) linear RF. RL ) repetition level.

tion, the accuracy should improve as more samples are included in D. The lowest value, LEM ) 0.1, was achieved using the linear regression functions and a repetition level RL ) 1. This value of the LEM represents a 2% error in the local mean prediction of the nanoparticles synthesis model in normalized variables. In Figure 6 notice that all LEM lines at high values of ns seems to have similar values, which suggests that neither the repetition level nor the regression trend functions should be critical once enough samples are included in D. Both sample sets are created from dynamic trajectories beginning the same operating region, so both sample sets should be located in the same portion of the state space. Therefore, the distancebased correlation component of the GPM will strongly dominate the predictions in this region. The strength of this correlation term in the GPM is quantified by the correlated residual parameter σc2. If the magnitude of σc2 is significantly greater than the uncorrelated noise parameter σu2, then the GPM performance will be strongly influence by the location of the sample points in the model. On the other hand, if σu2 is similar to or greater than σc2, it will be difficult to accurately estimate a distance-based correlation between the pairs of sample points in D. For the nanoparticle synthesis model, the estimated value of σu2 is 2 orders of magnitude smaller than the estimated σc2 parameter, indicating a low-tomoderate noise level in the kMC simulation. The same comparisons are made in Figure 7 as in Figure 6, except that the GPMs in Figure 6 were constructed from a set D of points randomly selected from the initial five trajectories, while in Figure 7, the points in D are sampled randomly from the convex hull bounding the initial five trajectories. Since the initial trajectories began at the corners of the operating condition space (for precursor mass and CNT mass), they are expected to also bound (or nearly bound) the trajectories for intermediate values of precursor and CNT mass. In both Figures 6 and 7, the same nt test points as in Figure 6 are used to compute the LEM. Our rationale for sampling from the convex hull is that it could provide better spacing and interpolation for new query points, compared to sampling only from existing trajectories. This might also be accomplished using fewer sample points ns in D. Comparing Figure 6 to Figure 7 we conclude that sampling from the convex hull yields less accurate predictions, compared to sampling from the trajectories. The lowest value for the LEM of 0.25 was similarly achieved using linear regression functions and a repetition level RL ) 1, and this represents a 4.5% error in the local mean prediction of the nanoparticles synthesis model in normalized variables. The reason for a higher LEM, compared to the previous case, is that the nanoparticle dynamics occurs in a smaller region of the constructed convex hull and so many

of the ns sample points may not be useful for predicting the dynamics. Therefore, by building the GPM using points far away from where the nt test points are, the GPM prediction of the nt test points is worse. This indicates that the convex hull is not a minimal representation of the nonconvex region bounding the trajectories. Since not all points in the convex hull sampling scenario may be nearby to the trajectories, the number of effective samples in D may be reduced, such that the details of the particular GPM implementation are now more significant. Although increasing the number of repetitions used in D might improve the overall accuracy, by enabling an independent estimation of the simulation noise σu2 and the correlated residual σc2, that trend is not seen in Figure 7. Instead, a greater benefit is realized by spreading out the sample points over the space. For example, the LEM for the constant regression function is lower when ns ) 120 different values of z are used in D, as opposed to 24 different values of z with RL ) 5 repetitions of each. This conclusion is consistent with our previous observations, since by spreading out the points, we can cover more effectively the region where the nt sample points are located. Additionally, the use of regression functions in GPM provides greater accuracy, at all values of ns and for all repetition levels. In GPM, regression functions act as global functions that describe trends in the data across the state space, with the remaining residual modeled using the distance-based correlation. In the case of the convex hull, the linear regression functions have added a global trend component that aids in the prediction. Overall, the GPM constructed from the precollected dynamic trajectories has the lower LEM, compared to the GPM constructed from samples in the convex hull. However, this may not be a completely fair comparison, since the test set is selected from trajectories, which are more similar to the sample set in Figure 6. In either case, a lower number of repetitions provides the lowest error, with no repetitions (RL ) 1) being best. The use of the linear regression functions provides increased accuracy at each of the different repetition levels, compared with the constant regression function. 4.2. Global Prediction Analysis. With an understanding of the one-step-ahead error, we now compute the error associated with the recursive use of the GPM for a dynamic trajectory. In this case the error is propagated from one time step to the next, since only the initial condition is known for certain. nm

DEM )

nt

nk

[

d

∑ ∑ ∑ ∑ (fj (z [0]) -

1 nmntnk m)1

j)1 k)1

k i

j

i)1

]

0.5

{fˆ k(zj[0], Dm)}i)2

(22)

To quantify the dynamic error, we define the dynamic error metric (DEM), which is identical to the LEM except that the equation is now summed over all points in the trajectory described by the discrete time index k. We expect that the mean error per step will be larger than that in the LEM, due to the additional propagation of error along the trajectory. This uncertainty is modeled via the prediction variance σy2 in eq 8, with the propagation computed via eqs 16 and 17. In the results plotted in Figures 8 and 9, we again compare the various sampling strategies in Table 2. The number of kMC realizations for each particular initial condition is nr ) 10, and the test set now consists of nt ) 40 different initial conditions for the precursor mass mpre and carbon nanotube mass mCNT. These initial conditions are selected from the ranges defined in Table 1, using a stochastic spacing filling design via the Latin


1387

2

Figure 8. Global prediction analysis DEM for different repetition levels using points from precollected dynamic trajectories. Regression functions: (s) constant RF, (---) linear RF. RL ) repetition level.

σc from a GPM with constant functions is greater than the estimated σc2 from a GPM with trend regression functions), it avoids the potential extrapolation problems associated with the trend functions. Even in the case of a poor coverage of the dynamic region, the constant regression function will predict a constant value of the known dynamic information across the state space when no sample points are nearby. The best overall performance for dynamic predictions is seen using RL ) 1, which spreads out the sample points over the dynamic region instead of performing replicates. For the lowest DEM value of 0.2, the error per state is approximately 4%. The implementation with linear trend functions exhibits lower LEM values compared to an implementation with constant functions but it lacks robustness and reliability when dynamic trajectories are performed. The use of the convex hull does not yield any obvious benefit, compared to the simpler sampling-fromtrajectories method. The accuracy of the GPM predictions are ultimately limited by the noise level in the kMC simulations. 5. Conclusions

Figure 9. Global prediction analysis DEM for different repetition levels using random points selected in the convex hull. Regression functions: (s) constant RF, (---) linear RF. RL ) repetition level.

hypercube method.50 The trajectories have the same length as the initial trajectory set, with ∆t ) 300 s and nk ) 25. The same nm ) 10 different GPMs that were constructed for the LEM calculations are used again here. The coupled ODE-GPM model executes a prediction of a single dynamic trajectory in 7 s compared to 58 s using the original coupled ODE-kMC model, which represents a 90% decrease in the computational load (or 99% compared to nr ) 10 kMC repetitions for averaging). Figures 8 and 9 exhibit a similar trend in the repetition level RL, as for the LEM. The advantage gained by spreading the points over the region is greater, compared to the advantage of using repetitions at a smaller number of distinct points to better estimate σu2. A new feature seen in all plots is that the GPMs constructed using linear regressions functions do not provide reliable predictions, due to the extrapolation of the trend functions into regions that do not have a large number of samples. When a particular random set of sample points Dm is generated, there is a possibility that those selected samples might not effectively cover the nanoparticle dynamic region. In such cases, the recursive formulation of the GPM could drive the dynamic prediction to regions far away from the sample points, such that the prediction relies only on the regression functions. In the case of the sampling-from-trajectories method in Figure 8, the linear regression function approach does provide a low error once enough samples are obtained, but in general that implementation may be less desirable, due to the danger of extrapolation and the lack of any improved accuracy. In contrast, using a constant regression function is a more conservative and robust implementation of the GPM. Because this GPM prediction relies purely on the distance-based correlation between the pairs of sample points (i.e., the estimated

Gaussian process modeling is a flexible framework for empirical modeling, which has been applied here to a multivariable stochastic simulation of nanoparticle synthesis dynamics. The GPM approach combines aspects of regression modeling and interpolation, in a statistically rigorous framework. Such features are not present in other black-box modeling approaches, such as response surface modeling and neural networks. Additionally, the GPM enables a calculation of the uncertainty in the mean prediction, which is very useful in blackbox empirical modeling. The nanoparticle case study has shown that repetitions in the sample set may be an inefficient use of samples in GPM, and instead the samples should be distributed throughout the dynamic region. The convex hull sampling method did not provide a clear benefit compared to direct sampling from the initial trajectory set, and thus may not justify the added computation and complication associated with constructing it and sampling from it. The regression functions improved the accuracy of the single-step error, but were detrimental when used recursively for trajectories due to extrapolation. Here we analyzed the accuracy of the GPM prediction for various options in the initial sampling of the stochastic simulations. However, the prediction variance in the GPM opens the option for even more efficient sampling. Instead of selecting the samples in a purely random fashion, their locations could be further optimized, for example to maximize their separation. The sample set used to build the GPM could also be progressively updated, via sequential sampling at locations of the state space that are most uncertain, or most likely to be needed for optimized trajectories. Acknowledgment The authors gratefully acknowledge funding from the Air Force Office of Scientific Research (Grant FA9550-07-0161) and the National Science Foundation (CBET-0933403), and helpful input from Prof. J. C. Lu, Prof. Dennis Hess, and Dr. Galit Levitin. Literature Cited (1) Montgomery, D. C. Design and Analysis of Experiments, 6th ed.; John Wiley & Sons, Inc.: Hoboken, NJ, 2005. (2) Sacks, J.; Welch, W. J.; Mitchell, T. J.; Wynn, H. P. Design and analysis of computer experiments. Stat. Sci. 1989, 4, 409–435.

1388


(3) Kleijnen, J. P. C.; Sanchez, S. M.; Lucas, T. W.; Cioppa, T. M. State-of-the-art review: A user’s guide to the brave new world of designing simulation experiments. INFORMS J. Comput. 2005, 17, 263–289. (4) Chen, V. C. P.; Tsui, K.-L.; Barton, R. R.; Meckesheimer, M. A review on design, modeling and applications of computer experiments. IIE Trans. 2006, 38, 273–291. (5) Hu, W.; Li, E. Y.; Li, G. Y.; Zhong, Z. H. Development of metamodeling based optimization system for high nonlinear engineering problems. AdV. Eng. Software 2008, 39, 629–645. (6) Davis, E.; Ierapetritou, M. A kriging method for the solution of nonlinear programs with black-box functions. AIChE J. 2007, 53, 2001– 2012. (7) van Beers, W. C. M.; Kleijnen, J. P. C. Kriging for interpolation in random simulation. J. Operat. Res. Soc. 2003, 54, 255–262. (8) van Beers, W. C. M.; Kleijnen, J. P. C. Customized sequential designs for random simulation experiments: Kriging metamodeling and bootstrapping. Eur. J. Operat. Res. 2008, 186, 1099–1113. (9) Kleijnen, J. P. C.; van Beers, W.; van Nieuwenhuyse, I. Constrained optimization in expensive simulation: Novel approach. Eur. J. Operat. Res. 2010, 202, 164–174. (10) Hernandez, A. F.; Grover, M. A. Stochastic dynamic predictions using Gaussian process models for nanoparticle synthesis. Comput. Chem. Eng. 2010, 34, 1953–1961. (11) Kevrekidis, I. G.; Gear, C. W.; Hummer, G. Equation-free: The computer-aided analysis of complex multiscale systems. AIChE J. 2004, 50, 1346–1355. (12) Pope, S. B. Computationally efficient implementation of combustion chemistry using in-situ adaptive tabulation. Combust. Theory Modell. 1997, 1, 41–63. (13) Varshney, A.; Armaou, A. Multiscale optimization using hybrid PDE/kMC process systems with application to thin film growth. Chem. Eng. Sci. 2005, 60, 6780–6794. (14) Zou, Y.; Kevrekids, I. G. Uncertainty quantification for atomistic reaction models: An equation-free stochastic simulation algorithm example. Multiscale Model. Simul. 2008, 6, 1217–1233. (15) Alexander, F. J.; Johnson, G.; Eyink, G. L.; Kevrekidis, I. G. Equation-free implementation of statistical moment closures. Phys. ReV. E 2008, 77, 026701. (16) Doyle, F. J., Pearson, R. K., Ogunnaike, B. A. Identification and Control using Volterra Models; Springer: Berlin, 2002. (17) Boyd, S.; Chua, L. O. Fading memory and the problem of approximating nonlinear operators with Volterra series. IEEE Trans. Circuits Syst. 1985, 32, 1150–1161. (18) Himmelblau, D. Accounts of experiences in the application of artificial neural networks in chemical engineering. Ind. Eng. Chem. Res. 2008, 47, 5782–5796. (19) Williams, C. K. I., Rasmussen, C. E. In AdVances in Neural Information Processing Systems 8; Touretzky, D. S., Mozer, M. C., Hasselmo, M. E., Eds.; The MIT Press: Cambridge, MA, 1996; pp 514520. (20) Rasmussen, C. E., Williams, C. K. I. Gaussian Processes for Machine Learning; The MIT Press: Cambridge, MA, 2006. (21) Kocijan, J.; Likar, B. Gas-liquid separator modeling and simulation with Gaussian-process models. Simul. Model. Pract. Theory 2008, 16, 910– 922. (22) Zavala, V. M.; Constantinescu, E. M.; Krause, T.; Anitescu, M. On-line economic optimization of energy systems using weather forecast information. J. Process Control 2009, 19, 1725–1736. (23) Ko, J.; Fox, D. GP-BayesFilters: Bayesian filtering using Gaussian process prediction and observation models. Auton. Robots 2009, 27, 75– 90. (24) ver Hoef, J. M.; Cressie, N. Multivariable spatial prediction. Math. Geol. 1993, 25, 219–240. (25) Deisenroth, M. P.; Rasmussen, C. E.; Peters, J. Gaussian process dynamic programming. Neurocomputing 2009, 72, 1508–1524. (26) Erkey, C. Preparation of metallic supported nanoparticles and films using super-critical fluid deposition. J. Supercrit. Fluids 2007, 47, 517– 522.

(27) Bayrakceken, A.; Kitkamthorn, U.; Aindow, M.; Erkey, C. Decoration of multiwall carbon nanotubes with platinum nanoparticles using supercritical deposition with thermodynamic control of metal loading. Scr. Mater. 2007, 56, 101–103. (28) Lin, Y.; Cui, X.; Yen, C.; Wai, C. M. Platinum/carbon nanotube nanocom-posite synthesized in supercritical fluid as electrocatalysts for lowtemperature fuel cells. J. Phys. Chem. B 2005, 109, 14410–14415. (29) Finney, E. E.; Finke, R. G. Nanocluster nucleation and growth kinetic and mechanistic studies: A review emphasizing transition-metal nanoclusters. J. Colloid Interface Sci. 2008, 317, 351–374. (30) Gummalla, M.; Tsapatsis, M.; Watkins, J. J.; Vlachos, D. G. Multiscale hybrid modeling of film deposition within porous substrates. AIChE J. 2004, 50, 684–695. (31) Venables, J. A. Rate equation approaches to thin-film nucleation kinetics. Philos. Mag. 1973, 27, 697–738. (32) Zhang, Y.; Kang, D.; Saquing, C.; Aindow, M.; Erkey, C. Supported platinum nanoparticles by supercritical deposition. Ind. Eng. Chem. Res. 2005, 44, 4161–4164. (33) Yoda, S.; Mizuno, Y.; Furuya, T.; Takebayashi, Y.; Otake, K.; Tsuji, T.; Hiaki, T. Solubility measurements of noble metal acetylacetonates in supercritical carbon dioxide by high performance liquid chromatography (HPLC). J. Supercrit. Fluids 2008, 44, 139–147. (34) Aschenbrenner, O.; Kemper, S.; Dahmen, N.; Schaber, K.; Dinjus, E. Solubility of β-diketonates, cyclopentadienyls, and cyclooctadiene complexes with various metals in supercritical carbon dioxide. J. Supercrit. Fluids 2007, 41, 179–186. (35) Aschenbrenner, O.; Dahmen, N.; Schaber, K.; Dinjus, E. Adsorption of dimethyl(1,5-cyclooctadiene)platinum on porous supports in cupercritical carbon dioxide. Ind. Eng. Chem. Res. 2008, 47, 3150–3155. (36) Chrastil, J. Solubility of solids and liquids in supercritical gases. J. Phys. Chem. 1982, 86, 3016–3021. (37) Span, R.; Wagner, W. A new equation of state for carbon dioxide covering the fluid region from the triple-point temperature to 1100 K at pressures up to 800 MPa. J. Phys. Chem. Ref. Data 1996, 25, 1509–1596. (38) Ye, X.-R.; Lin, Y.; Wang, C.; Engelhard, M. H.; Wang, Y.; Wai, C. M. Supercritical fluid synthesis and characterization of catalytic metal nanoparticles on carbon nanotubes. J. Mater. Chem. 2004, 14, 908–913. (39) Gillespie, D. T. A general method for numerically simulating the stochastic time evolution of coupled chemical reactions. J. Comput. Phys. 1976, 22, 403–434. (40) Cressie, N. Statistics for Spatial Data, 3rd ed.; Wiley Interscience: New York, 1993. (41) Koehler, J. R., Owen, A. B. In Handbook of Statistics; Ghosh, S., Rao, C. R., Eds.; Elsevier Science: New York, 1996; pp 261-308. (42) Goldberger, A. S. Best linear unbiased prediction in the generalized linear regression model. J. Am. Stat. Assoc. 1962, 57, 369–375. (43) Joseph, V. R.; Hung, Y.; Sudjianto, A. Blind kriging: A new method for developing metamodels. J. Mech. Des. 2008, 130, 031102. (44) Kitanidis, P. K.; Shen, K.-F. Geostatistical interpolation of chemical concentration. AdV. Water Resour. 1996, 19, 369–378. (45) Azman, K.; Kocijan, J. Application of Gaussian processes for blackbox modelling of biosystems. ISA Trans. 2007, 46, 443–457. (46) Todini, E.; Ferraresi, M. Influence of parameter estimation uncertainty in kriging. J. Hydrol. 1996, 175, 555–566. (47) Girard, A., Murray-Smith, R. In Switching and Learning in Feedback Systems (Lecture Notes in Computer Science 3355); MurraySmith, R., Girard, A., Eds.; Springer: Berlin, 2005; pp 158-184. (48) Ramkrishna, D. Population Balances: Theory and Applications to Particulate Systems in Engineering; Academic Press: Sang Diego, CA, 2000. (49) Martin, J. D.; Simpson, T. W. Use of kriging models to approximate deterministic computer models. AIAA J. 2005, 12, 115–125. (50) Loh, W.-L. On Latin hypercube sampling. Ann. Stat. 1956, 24, 2058–2080.

ReceiVed for reView April 1, 2010 ReVised manuscript receiVed October 15, 2010 Accepted December 3, 2010 IE1007954

Comparison of Sampling Strategies for Gaussian Process Models

Recommend Documents