Non-parametric correlative uncertainty quantification and sensitivity analysis: Application to a Langmuir bimolecular adsorption model

Uncertainty quantification (UQ) is a critical tool to assess the predictive ability of a model. This information can be used to map the envelop of model predictions, to improve a model via error reduction methods, and to inform control and optimization strategies in system tasks.^1–9 Sensitivity analysis (SA) is one of the most effective tools for identifying influential model parameters.^8,10–12 The SA approaches are typically classified as local or global methods. Local sensitivity analysis (LSA) computes variability in model predictions due to infinitesimal perturbations in the model parameters.⁸ The resulting local sensitivity indices (LSIs) include gradient⁸ and information-based methods.¹³ LSIs have also been applied to system optimization and model calibration problems.^14,15 Global sensitivity analysis (GSA) determines variability in model predictions over a range of parameters due to uncertainty in those parameters. GSA techniques include analytical, regression, screening and variance-based methods.^8,16–27

While significant progress has been achieved in UQ and SA methods over the years, physical systems often exhibit correlated parameters. In recent work, we introduced such a mathematical framework and demonstrated the impact of correlations on a model predictive ability for a complex reaction network.²² Our ability to understand and improve methods relies on developing simple but physically sound models that we can analyze mathematically.

In this paper we propose non-parametric methods for the GSA of chemical reaction models with correlated parameter dependencies. A Langmuir bimolecular hydrogen/oxygen competitive adsorption model is employed as a benchmark to motivate and concretely illustrate the derivation and algorithmic aspects of the proposed method. This system describes the competitive adsorption of hydrogen and oxygen on a Pt(111) surface. Such systems are encountered in catalytic oxidation, such as emissions abatement, small scale power generation, fuel cells and batteries. Here, parameter correlations stem from correlated quantum-scale computational data calculated using Density Functional Theory (DFT). These correlations are assimilated into the model as an informed prior distribution for the model’s parameters. The use of non-parametric methods in modeling parameter uncertainty and understanding global sensitivity is necessitated by the limited availability of quantum-scale data. The proposed mathematical methodology employs gradient-based methods to compute correlative local/global sensitivity indices (LSI/GSI) to illustrate the relative effects of parameter perturbations, errors and uncertainties in model parameters. We show that the ranking of influential parameters depends critically on whether or not correlations between parameters are taken into account. The impact of uncertainty in the correlation on the LSI/GSI is also demonstrated. Finally, we show the necessity of the proposed non-parametric perspective by comparing with a parametric approach.

The paper is organized as follows: the general form of the studied predictive model and LSI/GSI definitions are briefly presented in the “Background on sensitivity analysis” section. A detailed description of the impact of parameter correlations on LSA/GSA is addressed in the “Parameter correlation” section. Then, we apply the proposed CLSA/CGSA in the case of correlated hydrogen/oxygen binding energies in the “Langmuir bimolecular-adsorption model” section. Finally, we discuss methods to implement the proposed correlative LSA/GSA methods, other than direct numerical integration or Monte Carlo, namely generalized Polynomial Chaos methods that account for non-parametric noise models and parameter correlations.

In this section an appropriate mathematical framework is discussed for sensitivity analysis. First, we consider an ensemble of models of the general form

where Π(x|λ) denotes the predictive forward mathematical models, i.e. the probability distribution function (PDF) of state X = x for fixed K dimensional model parameters $λ=[λ1λ2⋯λK]T∈Λ$⁠, and Λ presents the parameter space. The term p(λ) denotes the PDF of λ which contains knowledge of uncertainty in the model, i.e. once we have p(λ) we can generate ensembles of X’s for each λ. Note that X may represent a static random variable, a snapshot of the system at some fixed time, or an entire time-series for dynamics and λ may denote the model parameters or indexing of different models.

In our specific model, λ corresponds to the binding energy of atomic oxygen and hydrogen on a given metal surface. We look at the uncertainty in coverage (Π) given the binding energy and its associated uncertainties and correlations. Coverage can also depend on other quantities that could be represented by λ, such as binding site, inert species, surface defects, surface impurities, and surface temperature.^28–30 The formalism is though general and beyond the binding energy and the isotherm. Other physical systems that follow the Π(x|λ)p(λ) relationship include the dependency of molecular frequency (Π) on coverage (λ)³¹ and forecasted temperature changes (Π), with CO₂, methane, and other greenhouse gases (λ).^32–34 

The system observable can be defined over all possible realizations of the state

where h(x) denotes a desired quantity. The correlations in the parameter vector λ are also included in p(λ) and are propagated into the state X through the predictive forward model of Π(x|λ). Finally, the averaged observable for the model can be defined by

Consider a general class of nonlinear models of the form

where f is an arbitrary scalar function. The (relative) LSI of f with respect to λ of (4) at the nominal value of $λ*$ is

where

The LSI of (5) supplies useful sensitivity information in the case of almost certain parameters, i.e. for a relatively tight range of parameter values. To incorporate the knowledge of uncertain parameter distributions and provide sensitivity information over the entire range, we determine the relevant GSI by employing the partial derivative of the LSI as a basic building block to integrate the local sensitivities over the total range of parameter changes

where q denotes the type of required index (q = 1: improved Morris index, q = 2: asymptotic limit of the standard Morris index).

Note that the PDF, which incorporates knowledge of the λ distribution in the GSI of (7), must be identified subject to available experimental and/or simulation-based data. The possible correlations between the system parameters, which may be discovered during regression of the data in statistical models, must be encoded in p(λ). Such correlations play a deciding role in sensitivity analysis and their effects are quantified in the following sections.

Previously, LSA and GSA were treated for the case of independent parameters. To extend sensitivity analysis to models with correlated parameters, we partition the vector of parameters into two,

where $λ=[λ1λ2]T$⁠, Λ = Λ₁ ⊕ Λ₂, λ₁ contains all independent parameters, and λ₂ contains all dependent parameters. Parameters can be classified through their correlations, which are identified using experiments and/or computational tools for specific case studies, and by applying causality statistical methods. When λ₁ and λ₂ are correlated, perturbations in one parameter affect the other. Proper mathematical tools are needed to quantify parameter correlations and their impact on model reliability.

The correlation between λ₁ and λ₂ can be described by their joint probability distribution,

For the marginal distribution of λ₁,

Identifying p(λ₁) and p(λ₂|λ₁) in a systematic way is an essential step in our CGSA.

The joint probability distribution of p(λ₁, λ₂) can either be built directly, with a sufficiently large ensemble of experimental and/or simulation-based data,²² or computed sequentially by marginalization according to (9). The latter approach requires building PDFs with data for both p(λ₂|λ₁) and p(λ₁), followed by Monte Carlo sampling to calculate correlative indices. There are various linear regression (LR) methods that can identify the conditional PDF of p(λ₂|λ₁); including deterministic (DLR), stochastic (SLR), and Bayesian (BLR).^35,36 DLR yields a deterministic linear model whose LCSIs can be computed exactly, while its GSIs depend on the choice of p(λ₁). SLR uses a least squares model for p(λ₂|λ₁), typically in a Gaussian form.³⁶ Usually, there is not enough data for adequate fitting of a least squares model. BLR can bypass this shortcoming by putting a prior on the parameters in the linear fit.³⁵ 

Hierarchical or empirical Bayesian methods can identify the marginal PDF of p(λ₁) via deterministic linear or stochastic nonlinear fitting to the data. Bootstrapping does not require fitting but instead relies on simulation,³⁶ which is particularly appropriate for problems with little data, by creating synthetic new samples using random re-sampling according to the actual distribution of the data.³⁶ In this way we can enrich the histogram with more data points and obtain a Gaussian distribution.³⁶ If histograms are too sparse, “smoothed bootstrapping” can be used. This method applies a kernel to data from a standard histogram. A Bayesian approach can be used to fit data to a well known distribution that “looks like” the histogram. “Looks like” means that we pick a family of well known distributions and match the first few moments with the corresponding moments of the data’s histogram, i.e. mean, variance, skewness, etc.^36,37

We perform the correlative sensitivity analysis by focusing on λ₁, while still accounting for the correlations with λ₂

The correlative local sensitivity index (CLSI) at the nominal point $λ1*$ is obtained similarly to (5) by direct differentiation,

The CGSI can then be formulated

by employing the CLSI of (12) as building blocks where q = 1 or q = 2.

For deterministic correlation where λ₂ = g(λ₁), we can simplify the λ₁-marginal PDF of (10) by considering p(λ₂|λ₁) = δ(g(λ₁) − λ₂),

where δ(⋅) denotes the standard Dirac function. Therefore, from (11), we have

and the CLSI can be simplified to the following form

The additional second term in the CLSI differs from uncorrelative LSI in that the derivative with respect to the parameter λ₂ comes directly into play. The CGSI formulation for such a simplified case is

The implementation of the sampling algorithm used to compute the correlative local/global sensitivity index (CLSI/CGSI) is described in Section S5 of the supplementary material.

We consider a Langmuir bimolecular adsorption model which describes competitive dissociative adsorption of hydrogen (H₂) and oxygen (O₂) on a catalyst surface,

where H₂ and O₂ denote the hydrogen and oxygen molecules in the gas phase, $2*$ are two active sites on the metal surface, and $H*$ and $O*$ represent the adsorbed hydrogen and oxygen atoms on the surface, respectively. A schematic of this adsorption process is illustrated in Figure 1. The physical system is related to hydrogen oxidation in fuel cells and batteries.^38–42 

The coverages dynamics can be formulated by the following set of ordinary differential equations

where $θH*0$ and $θO*0$ represent the initial hydrogen and oxygen coverages, respectively. $PH2$ and $PO2$ are the partial pressures of the gas phase species,⁴³ and we set $PH2=1.01325×10−10N/m2$⁠, $PO2=1.01325×10−50N/m2$ in this paper.

The hydrogen and oxygen coverages at equilibrium are

where P is partial pressure and $K=kadskdes$ is the equilibrium constant.^44,45K is determined from DFT calculations.

By focusing on variations of binding energies and fixing the other parameters, the coverages are

where ΔE_H and ΔE_O are the binding energies of atomic hydrogen and oxygen to the surface. It is the effect that uncertainty and correlations in the binding energies have on the coverages that we explore below. The detailed formulas of $θ^H*$ and $θ^O*$ is derived in Section S1 of the supplementary material.

The Langmuir adsorption isotherm is strictly valid at low coverage with adsorption at a single site, which is our system of interest. For dissociative adsorption, which we study here, the Langmuir model requires adjacent empty sites on the catalyst surface. An ab-initio molecular dynamics (AIMD) study of hydrogen on Pd(100) showed that regardless of coverage, only two adjacent catalyst sites are necessary to dissociate hydrogen.⁴⁶ Although applications of AIMD to heterogeneous catalysis are rapidly advancing, the computational cost is still prohibitive for it to be used in generating adsorption isotherms.⁴⁷ Less computational intensive methods, such as Monte Carlo⁴⁸ and molecular dynamics with force fields,⁴⁹ are used instead to generate an isotherm.

Seller et al. have shown that, when combined with the Bragg-Williams coverage model, the Langmuir adsorption isotherm accurately recreates experimental isotherms for several systems.⁴⁸ Furthermore, the same study found that the Langmuir adsorption isotherm with mean field treatment compares favorably with coverages predicted from lattice based grand canonical Monte Carlo (GCMC) simulations under certain conditions. A force field based molecular dynamics simulation of dimethyl methylphosphonate (DMMP) also supports the validity of the Langmuir model.⁴⁹ 

Electronic contributions to adsorption enthalpies are calculated with DFT, using the Vienna ab initio Simulation Package (VASP), version 5.3.2,^50–53 with the plane wave basis set, PAW pseudopotentials,^54,55 and periodic boundary conditions. Simulation parameters are similar to those used in our previous work.⁵⁶ All VASP input files are created using the Atomistic Simulation Environment, an open-source python-based software program.⁵⁷ We employ the PBE exchange-correlation functional with D3 dispersion corrections.^58,59 Spin-polarized calculations are performed for molecules in a vacuum and systems containing Ni and Co. The first Brillouin zone is sampled using the Monkhurst-Pack (3x3x1) mesh.⁶⁰ For the purposes of this work, the level of accuracy achieved using this mesh size was sufficient. Ionic force cut-off for all calculations is set to 0.05 eV/Å. In slab calculations, we use the (4x4) supercell containing four layers of atoms, with the positions of the bottom two layers fixed. We use an adsorbate coverage of 1/16 monolayers in all calculations.

The DFT dissociative adsorption energy for molecular oxygen on a metal surface is defined in Equation 22. A similar relation holds for dissociative adsorption of molecular hydrogen.

such that

In Equation 22, $EO2$ is the DFT energy of an O₂ molecule in a vacuum, $E*$ is the DFT energy of the pristine metal slab, $EO*$ is the energy of the adsorbate-metal system, D₀ is the gas-phase bond dissociation energy at 0 K, and E_O is DFT adsorption energy of an oxygen atom. The calculated energies for a variety of metal surfaces and the resulting scaling relationship between electronic contributions of H and O dissociative adsorption energies are shown in Figure 2. We obtain a linear correlation between hydrogen and oxygen adsorption energies with a R² value of 0.87, i.e.

where a = 2.51, b = −2.46 (eV).

Vibrational contributions and other temperature effects to adsorption enthalpy and entropy are accounted for in calculations of adsorption equilibrium constants (see Section S1 of the supplementary material for details). Zero point energy (ZPE) corrections for gas phase H₂ and O₂ are calculated using their experimental vibrational frequencies.^61,62 Kinetic energy loss upon adsorption is accounted for by using the ideal gas value of $32RT$⁠. Harmonic and rigid rotor approximations were utilized to account for vibrational and rotational degrees of freedom, respectively.⁶³ Hessian matrices are constructed using 0.015 Å displacements in x, y, and z directions from adsorbate equilibrium positions.

In order to validate our computational setup and provide an error estimate, we compare the calculated and experimental adsorption enthalpies of oxygen and hydrogen on platinum in Table I.^64–67 The DFT calculations reproduce the experimental data well.

Figure 3 highlights the differences in our model, defined in (1) for $(θ^H*,θ^O*)$⁠, resulting from correlations between ΔE_H and ΔE_O. Consider

where $Π(θ^H*,θ^O*|ΔEH,ΔEO)=δ(θ^H*(ΔEH,ΔEO),θ^O*(ΔEH,ΔEO))$⁠, δ is the standard Dirac function and both $θ^H*(ΔEH,ΔEO)$ and $θ^O*(ΔEH,ΔEO)$ are given by (21).

Then, for the uncorrelated case (subscript uc), we can assume that

where p(ΔE_H) and p(ΔE_O) are defined as density functions of Gamma distributions, given by (30) and (31) in next section.

In the correlated case (subscript c), we assume that

where p(ΔE_H) is still given by a Gamma distribution but p(ΔE_O|ΔE_H) comes from a normal distribution with mean aΔE_H + b and variance determined by the data in Section VII B, which gives p_c(ΔE_H, ΔE_O) a lower variance than p_uc(ΔE_H, ΔE_O).

We can use changing of variables such that

where J is the Jacobian of the inverse of coverage function $θ^(ΔEH,ΔEO)$ from (8), evaluated at $(θ^H*,θ^O*)$⁠. Figure 3 shows the density function contours $p(θ^H*,θ^O*)$⁠, in the uncorrelated case using (26), and correlated case using (27). Note that the correlation of ΔE_H and ΔE_O reduces the variance of our model.

In this section, we compute the uncorrelated local and global sensitivity index defined in Section II B for the coverages $θ^H*$ and $θ^O*$ with respect to ΔE_H, and will turn to the correlated cases in the next three sections. In section VI A, we compute the LSIs according to (21); and in section VI B, we construct an uncorrelated prior distribution for ΔE_H and ΔE_O, and then use this distribution and LSIs to compute GSIs.

Using the binding energies of adsorbed hydrogen and oxygen from Figure 2, we can analyze the relative LSIs for $θ^H*$ and $θ^O*$ with respect to ΔE_H. $SHH$ and $SHO$ are identified using (5) and the model given in (21) (detailed calculations are presented in Section S2 of the supplementary material). For Pt, ΔE_H = 2.6581(eV), ΔE_O = 3.6604(eV), and the sensitivity of H and O coverages with respect to the H binding energy are $SHH=38.9080$ and $SHO=−0.0138$⁠. As expected, the H binding energy has a large effect on its coverage and a slight effect on the O coverage (some coupling is expected due to the competitive nature of adsorption).

To compute the corresponding GSIs on Pt by (7), we need to construct the distribution of our parameters, p(ΔE_H, ΔE_O). In the uncorrelated case, we have

where p(ΔE_H) is the prior information for ΔE_H on Pt, and p(ΔE_O|ΔE_H) = p(ΔE_O) since ΔE_H and ΔE_O are independent, assuming no correlation.

Using the experimental and DFT data shown in Table I, we construct an informative prior for $ΔEH∈EH$⁠. Let x_H = 2.63 (eV), x_O = 3.71 (eV), y_H = 2.69 (eV) and y_O = 3.68 (eV), where x_i are the values given by experiment and y_i are given by DFT, i = H, O. To quantify uncertainty from DFT error, we assume that ΔE_H on Pt follows a gamma distribution with mean x_H and the standard deviation given by the difference between experiment and DFT, (x_i − y_i). We can construct the distribution for $ΔEO∈EO$ in the same way under the uncorrelated assumption. The explicit density functions are shown below,

where $ai=xi2/(xi−yi)2$ and $bi=(xi−yi)2/xi$⁠, i = H, O.

The GSIs, $ξHH$ and $ξHO$⁠, are formulated by (7) with q = 2 and p(λ) = p(ΔE_H)p(ΔE_O). Computing the integral in (7) numerically gives $ξHH=1509.8865$ and $ξHO=0.2194$⁠. Again, the H binding energy has a major effect only on the H coverage and a slight effect on the O coverage.

The following sections cover the CLSI. In section VII A, we consider the simplest correlation model, both deterministic and linear, to compute the CLSIs defined in Section III. In section VII B, we construct parametric models for p(ΔE_O|ΔE_H) using the data shown in Figure 2, and compute the corresponding CLSIs. Section VIII covers non-parametric models.

To calculate the CLSI for $θ^H*$ and $θ^O*$ with respect to ΔE_H, with the formula defined in (12), we use the conditional probability p(ΔE_O|ΔE_H). In the deterministic case, we can assume the conditional distribution of $EO$ has a mean of g(ΔE_H) and zero variance whose PDF can be described by the standard Dirac function

as presented in Section III. Then, using the data from Figure 2, one can determine the function g(ΔE_H) with different fitting models, like polynomials and smoothing splines. In this paper, we use the linear function and set g(ΔE_H) = a ΔE_H + b, as shown in Figure 2. Then, the conditional distribution can be written as

For brevity, we only consider the CLSI formulation, as defined in (12) for the hydrogen coverage $θ^H*$⁠, with respect to adsorbed hydrogen binding energy on the surface. The rest of CLSIs can be formulated by following the same procedure.

The CLSI for $θ^H*$ with respect to ΔE_H at the nominal hydrogen binding energy of $ΔEH^$ according to (12) takes the following form

where

The CLSI derivations in the presence of deterministic linear correlation are briefly described in Section S2 of the supplementary material and the results of $SH,corrH$ and $SH,corrO$ are shown in FIG. 1. The numerical results for Pt are $SH,corrH=38.9021$ and $SH,corrO=97.7442$⁠. The corresponding uncorrelated LSI indices from (5) are $SHH=38.9080$⁠, and $SHO=−0.0138$⁠; hence the correlation between ΔE_H and ΔE_O does not affect the sensitivity of $θ^H*$ with respect to ΔE_H, but does impact the sensitivity of $θ^O*$ with respect to ΔE_H. The LSI changes from slightly negative to highly positive. This is rationalized from the slope of the correlation depicted in Figure 2. Specifically, an increase in the binding energy of H leads a much higher increase in the binding energy of O and thus to an increase of the O coverage.

The above is a perfect linear model (deterministic) and ignores the variation around the linear fit of the binding energies. To capture correlations from the linear model, we set up a linear probabilistic model for ΔE_H and ΔE_O by introducing a random variable, ω, in the correlation,⁶⁸ 

To determine the distribution of ω, we can fit the data or adjusted data (to match the required domain of some distribution) using parametric models, like normal or gamma. Here we choose the normal distribution and fit the parameters using MATLAB by the Maximum Likelihood Estimation (MLE) method.⁶⁸ The result is shown in Figure 4.

Using (12) to compute the CLSI of $θ^H*$ with respect to ΔE_H, we consider

where p(ω) is the PDF for ω in (36). Instead of using the Monte Carlo method (discussed in Section S5 of the supplementary material), we can also use numerical integration to approximate the integral in (37).

The CLSI at the nominal point ΔE_H$*$ can then be obtained by direct differentiation

The gradient of $ln⁡FHH(ΔEH)$ is commonly estimated, such that

The CLSI for $θ^O*$⁠, $SH,corrO(ΔEH*)$⁠, is computed similarly. The numerical results for Pt are $SH,corrH=35.9874$ and $SH,corrO=9.4965$⁠. Compared to the deterministic model results in the previous subsection, we find that uncertainty significantly impacts $SH,corrO$⁠. This is a rather interesting result because the correlation in the data (linear) results in the H binding energy having a significant effect on the O coverage but uncertainty significantly diminishes this effect. We give results from other parametric models in Section S3 of the supplementary material.

For small data sets, such as ours, parametric models are not usually adequate. Instead, we consider non-parametric methods.⁶⁹ A possible non-parametric method is the empirical distribution function,

where I is the identity function. With this method, $EP^[f]$ for some function f can be approximated via bootstrapping.⁶⁹ 

For categorical distributions, the bootstrap distribution is close to the posterior distribution with a non-informative symmetric Dirichlet prior according to Bayes method. It also has the same support, mean, and nearly the same covariance matrix as the data in the histogram. The bootstrap distribution is obtained without specifying either the prior or sampling from the posterior distribution.⁷⁰ 

We can also use curve estimation for our model.⁶⁹ A simple density estimator is a histogram, which is a piece-wise constant function where the height of the function is proportional to number of observations in each bin

where B₁, …, B_n are the histogram bins, h = 1/n is the bin-width, and ν_i is the number of observations in B_i, as shown in Figure 4.

Smoother estimators, called kernel density estimators,⁶⁹ converge faster to the true density than fitting from histograms because histograms are discontinuous

where h > 0 is the bandwidth and K is the kernel, defined to be any smooth function satisfying K(x) ≥ 0, ∫ K(x)dx = 1, ∫ xK(x)dx = 0 and $σK2=∫x2K(x)dx>0$⁠.

In the main text of this work we use the histogram to approximate the distribution of ω, and use

to compute $SH,corrH(ΔEH*)$ and $SH,corrO(ΔEH*)$ using Equation 38. The numerical results for adsorption on Pt are $SH,corrH=35.2196$ and $SH,corrO=12.4873$⁠. Results from the kernel density estimators with uniform and standard normal kernel, $N(0,1)$⁠, are given in Section S4 of the supplementary material.

Figure 5 summarizes all local sensitivity analysis results (the magnitude is plotted so a semi-log scale can be used). Correlations play a significant role as demonstrated in our earlier work.²² Clearly, the uncertainty in the correlations must properly by accounted for and, given the limited number of data we have for physical models, non-parametric models of the uncertainty are essential. For a large sample size, both (parametric and non-parametric) models should converge to the real distribution.^36,69 Because we only have a few samples here, the non-parametric models approximate the noise term better, as shown in Figure 4.

In this section we compute the CGSIs for the correlation model previously used to determine the CLSIs. According to (13), the CGSIs are

As discussed in Section VII, we assume the prior distribution of ΔE_H on Pt satisfies

with a PDF of

using the data in Table I according to (30). Then, from $SH,corrH$ and $SH,corrO$⁠, we numerically calculate the CGSIs according to (44) and (45). The results are shown in Figure 6. Correlations have only a slight effect on the H coverage as we consider the H binding energy as an independent variable and the O binding energy as the dependent parameter.