A Generalized Polynomial Chaos-Based Method for Efficient Bayesian Calibration of Uncertain Computational Models

Abstract

This paper addresses the Bayesian calibration of dynamic models with parametric and structural uncertainties, in particular where the uncertain parameters are unknown/poorly known spatio-temporally varying subsystem models. Independent stationary Gaussian processes with uncertain hyper-parameters describe uncertainties of the model structure and parameters, while Karhunen-Loeve expansion is adopted to spectrally represent these Gaussian processes. The Karhunen-Loeve expansion of a prior Gaussian process is projected on a generalized Polynomial Chaos basis, whereas intrusive Galerkin projection is utilized to calculate the associated coefficients of the simulator output. Bayesian inference is used to update the prior probability distribution of the generalized Polynomial Chaos basis, which along with the chaos expansion coefficients represent the posterior probability distribution. The proposed method is demonstrated for calibration of a simulator of quasi-one-dimensional flow through a divergent nozzle with uncertain nozzle area profile. The posterior distribution of the nozzle area profile and the hyper-parameters obtained using the proposed method are found to match closely with the direct Markov Chain Monte Carlo-based implementation of the Bayesian framework. Efficacy of the proposed method is demonstrated for various choices of prior. Posterior hyper-parameters of the model structural uncertainty are shown to quantify acceptability of the simulator model.

Keywords:

• bayesian framework
• generalized polynomial chaos
• Karhunen-Loeve expansion
• hyper-parameters

AMS Subject Classifications:

• 15A29
• 33C45
• 58C40
• 60G15
• 62F15

Introduction

With the present ubiquitous use of computer simulators for scientific investigations, uncertainty quantification and calibration of the simulator models is identified as an important area of research.[1–7] Significant developments of the last decade have established the Bayesian framework as a preferred method for uncertainty quantification and calibration of computer simulators.[8–16] This paper explores a Bayesian framework for calibration and credibility assessment of a computer simulator, $\mathbf{T}(\mathbf{x},\mathbf{u})$, with possibly multivariate output. The simulator accepts deterministic control inputs $\mathbf{x}$ and uncertain parameters $\mathbf{u}$. This paper concerns simulators with uncertain spatial functions, thus $\mathbf{u}=\mathbf{u}\left({\mathbf{x}}_s\right)$, where arguments ${\mathbf{x}}_s$ may be different from $\mathbf{x}$. The uncertain functions are treated as uncertain subsystem models in this paper. Though the simulator is deterministic in a sense that each run with same arguments always provide same output, the simulator output is uncertain owing to uncertainty in these parameters. Aim of the framework is to infer $\mathbf{u}\left({\mathbf{x}}_s\right)$ using Bayesian calibration whenever new experimental observations are available. Figure 1 shows schematic of the proposed framework for a complex system consisting of physically or mechanically interconnected subsystems.

Fig. 1 Conceptual architecture of bayesian calibration for uncertain models.

The most appealing facet of the Bayesian framework is its ability to provide the complete posterior statistics. However, barring very simple cases, statistical sampling techniques are required for solution of the Bayesian calibration and uncertainty propagation problems. Markov Chain Monte Carlo (MCMC) methods [17, 18] are most widely used sampling techniques for Bayesian calibration. However, exploration of the posterior distribution using MCMC requires collection of a large number of samples for satisfactory approximation (often in the range of ${10}^3-{10}^6$), rendering the Bayesian framework computationally prohibitive for a large-scale system simulator. Thus, for a generalized application of a Bayesian framework to complex large-scale system simulators, development of a computationally efficient Bayesian calibration technique is essential.

Marzouk et al. [19, 20] have proposed a spectral projection based method for computationally efficient Bayesian calibration. The method uses a spectral expansion of a prior in generalized Polynomial Chaos (gPC) basis. Note that the gPC method is generalization of the Polynomial Chaos based spectral projection method, which is extensively investigated in the literature as a computationally efficient alternative to statistical methods for uncertainty propagation with comparable accuracy.[21] Polynomial Chaos method is based on the concept of Homogeneous Chaos introduced by Wiener [22, 23], where a random variable is spectrally expanded in terms of Hermite polynomials. Cameron and Martin [24] have shown that any non-linear functional can be expanded in terms of a series of Hermite polynomials in $L^2$ sense. Although earlier attempts at using the polynomial Chaos method (especially for turbulent fluid flow modeling) were not very successful,[25–27] the method is found to be useful for the solution of stochastic finite element [28–30] and stochastic fluid flow problems.[31, 32] Xiu and Karniadakis [33, 34] have generalized the Polynomial Chaos method for spectral projection in terms of the Askey scheme of polynomials.[35] The resultant gPC method has been applied by various researchers for uncertainty propagation through simulators of systems of engineering importance.[36–39]

The method of Marzouk et al. [19, 20] uses gPC expansion of $\mathbf{u}({\mathbf{x}}_s)$ as(1) $$\begin{array}{c}\hfill \mathbf{u}({\mathbf{x}}_s)\approx \sum _{p=1}^P{\hat{\mathbf{u}}}_p({\mathbf{x}}_s)H_p(\mathit{\xi }),\end{array}$$(1) where ${\hat{\mathbf{u}}}_p({\mathbf{x}}_s)$ are gPC expansion coefficients, $H_p(\mathit{\xi })$ are orthogonal polynomials from the Askey scheme, while $\mathit{\xi }$ is a set of independent random variables. Galerkin projection is used to calculate the corresponding gPC expansion coefficients of $\mathbf{T}(\mathbf{x},\mathbf{u})$ as(2) $$\begin{array}{c}\hfill {\mathbf{T}}_q(\mathbf{x})=\frac{〈\mathbf{T}(\mathbf{x},\sum _{p=1}^P{\hat{\mathbf{u}}}_p({\mathbf{x}}_s)H_p(\mathit{\xi })),H_q(\mathit{\xi })〉}{〈H_q^2(\mathit{\xi })〉},\end{array}$$(2) where $〈·,·〉$ denote an inner product on a Hilbert space of the Gaussian random vector $\mathit{\xi }$. Equation (2(2) $$\begin{array}{c}\hfill {\mathbf{T}}_q(\mathbf{x})=\frac{〈\mathbf{T}(\mathbf{x},\sum _{p=1}^P{\hat{\mathbf{u}}}_p({\mathbf{x}}_s)H_p(\mathit{\xi })),H_q(\mathit{\xi })〉}{〈H_q^2(\mathit{\xi })〉},\end{array}$$(2) ) can be solved using either intrusive or non-intrusive Galerkin projection. The intrusive method requires new discretization, numerical algorithms, and the code development for the reformulated governing equations using Galerkin projection, hence the nomenclature. On the contrary, non-intrusive method uses stochastic collocation approach to solve the reformulated governing equations by solving the deterministic model at appropriate collocation points. See Najm [39] and references therein for detailed discussion on intrusive and non-intrusive Galerkin projection methods. The gPC expansion given by Equation (2(2) $$\begin{array}{c}\hfill {\mathbf{T}}_q(\mathbf{x})=\frac{〈\mathbf{T}(\mathbf{x},\sum _{p=1}^P{\hat{\mathbf{u}}}_p({\mathbf{x}}_s)H_p(\mathit{\xi })),H_q(\mathit{\xi })〉}{〈H_q^2(\mathit{\xi })〉},\end{array}$$(2) ) of the simulator predictions is used to define the likelihood. On availability of the experimental observations, probability distribution of $\mathit{\xi }$ is updated using the Bayesian calibration. Note here that the Bayesian inference is reformulated as sampling from the posterior distribution of $\mathit{\xi }$.

In this paper, the method proposed by Marzouk et al. [19, 20] is extended for priors with uncertain hyper-parameters. Though a family of probability distributions to represent the prior uncertainty can be specified, associated hyper-parameters are rarely known. Realistic quantification of prior uncertainty requires specification of probability distributions for uncertain hyper-parameters. Use of the uncertain hyper-parameters is more ubiquitous in case of calibration of simulators with model structural uncertainty. Hierarchical Bayesian inference is proposed in the literature for calibration in presence of uncertain hyper-parameters.[10, 11] However, the methodology proposed by Marzouk et al. [19, 20] does not explicitly consider the effect of uncertain hyper-parameters in the formulation. To make the spectral stochastic projection based Bayesian inference more accurate, it is necessary to include uncertain hyper-parameters in the formulation.

This paper proposes an extension of the method of Marzouk et al. [19, 20] to take into consideration uncertainty in hyper-parameters of the prior distribution. A methodology is proposed to obtain a Karhunen-Loeve expansion (KL expansion) of a stochastic process in terms of functions of the hyper-parameters. The prior uncertainty in these hyper-parameters is expanded in a gPC basis. Galerkin projection is used to evaluate gPC coefficients of the resultant KL expansion terms of a stochastic process. The prior uncertainty in a subsystem model, represented in terms of gPC basis, is propagated to simulator predictions using the Galerkin projection approach.[30] The Bayesian calibration is reformulated as an MCMC sampling from the posterior distribution of the random vector $\mathit{\xi }$. Samples from the posterior distribution of uncertain parameters can then be obtained by substituting each posterior sample of $\mathit{\xi }$ into gPC expansion. Posterior distribution of the model structure defines the credibility of the simulator model. The posterior parameters are identified which quantifies acceptability of the simulator.

This research extends existing state of the art by: (a) extending the gPC expansion based method of Marzouk et al. [19,20] for priors with uncertain hyper-parameters; and (b) providing guidelines for acceptability of the simulator model using hyper-parameters of the posterior distribution. Note that a preliminary version of this work was reported in,[40] while this article is significantly expanded by including (a) model structural uncertainty; (b) substantially elaborate theoretical analysis; and (c) additional numerical results to establish computational efficiency and efficacy of the method.

The proposed method is demonstrated using a simulator of a quasi-one-dimensional flow through a nozzle. The particular choice of the application is motivated by the fact that the quasi-one-dimensional nozzle flow is well understood and can be simulated with limited computational resources. Spacial variation of the nozzle area profile, $A_{\mathit{NZ}}(x)$, is assumed uncertain, thus $\mathbf{u}({\mathbf{x}}_s)=A_{\mathit{NZ}}(x)$. Gaussian process prior is used to represent uncertainty in the nozzle area profile, while the proposed Bayesian framework is used to infer the nozzle area profile. It may be noted that though the method is demonstrated for Gaussian process prior in the present paper, the proposed method is generic in nature and can admit arbitrary prior. Further, since the proposed method updates probability distribution of $\mathit{\xi }$, posterior distribution of the nozzle area profile will be non-Gaussian.

The rest of the paper is organized as follows. Section 2 provides a statistical formulation of the problem. In section 3, the proposed method is discussed in detail. In section 4, numerical results for the calibration of a quasi-one-dimensional nozzle flow simulator are presented and finally in section 5, the paper is summarized and concluded.

Statistical formulation

Let a system be investigated by observing $M$ system responses, while, $T_j(\mathbf{x},u({\mathbf{x}}_s))$ be an available simulator of the $j\text{th}$ system response, where $\mathbf{x}\in \mathcal{X}$ is a set of deterministic control inputs and $u({\mathbf{x}}_s)$ is an uncertain subsystem model, where arguments ${\mathbf{x}}_s$ may be different from $\mathbf{x}$. For brevity, the discussion presented in this paper assumes a single uncertain subsystem model; however, the proposed method can be extended to a more generic case without any change. Let $u_t({\mathbf{x}}_s)$ denote the ‘true’ subsystem model, while ${\mathit{\zeta }}_j(\mathbf{x},u_t({\mathbf{x}}_s))$ is a model representative of reality. Note that the simulator $T_j(·,·)$ approximates the system response within limits of available knowledge. Thus conditional on the true subsystem model $u_t({\mathbf{x}}_s)$, $T_j(\mathbf{x},u_t({\mathbf{x}}_s))$ deviates from ${\mathit{\zeta }}_j(\mathbf{x},u_t({\mathbf{x}}_s))$ by(3) $$\begin{array}{c}\hfill {\mathit{\zeta }}_j(\mathbf{x},u_t({\mathbf{x}}_s))=T_j(\mathbf{x},u_t({\mathbf{x}}_s))+{\mathit{\delta }}_j(\mathbf{x}),\end{array}$$(3) where ${\mathit{\delta }}_j(\mathbf{x})$ is a discrepancy function.

Let $y_{e_j}(\mathbf{x})$ represent an experimental observation of the system at a control input setting $\mathbf{x}$, while ${\mathit{ϵ}}_j(\mathbf{x})$ denote the corresponding measurement uncertainty. Relationship between experimental observation and the true system response is given by(4) $$\begin{array}{c}\hfill y_{e_j}(\mathbf{x})={\mathit{\zeta }}_j(\mathbf{x},u_t({\mathbf{x}}_s))+{\mathit{ϵ}}_j(\mathbf{x}).\end{array}$$(4) Let experimental observations be obtained at $N$ input conditions and denote a set of  experimental observations by ${\mathbf{Y}}_e=\{{\mathbf{Y}}_{e_j};j=1,...,M\}$ where ${\mathbf{Y}}_{e_j}=\{y_{e_j}({\mathbf{x}}_i);i=1,...,N\}.$ Similarly, define ${\mathit{\delta }}_j=\{{\mathit{\delta }}_j({\mathbf{x}}_i);i=1,...,N\}$ and $\mathit{\delta }=\{{\mathit{\delta }}_j;j=1,...,M\}$. Also define ${\mathbf{u}}_t=\{u_t({\mathbf{x}}_{s,i});i=1,...,N_u\}$, where typically $N_u$ is significantly greater than $N$. Note that ${\mathbf{u}}_t$ represent a realization of a random function $u_t({\mathbf{x}}_s)$ at $N_u$ input settings. Information available from the experimental observations is used in the Bayes theorem as(5) $$\begin{array}{c}\hfill \mathcal{P}({\mathbf{u}}_t,\mathit{\delta }\mid {\mathbf{Y}}_e)\propto \mathcal{P}({\mathbf{Y}}_e\mid {\mathbf{u}}_t,\mathit{\delta })\times \mathcal{P}({\mathbf{u}}_t,\mathit{\delta }),\end{array}$$(5) where $\mathcal{P}({\mathbf{u}}_t,\mathit{\delta })$ is the prior, $\mathcal{P}({\mathbf{Y}}_e\mid {\mathbf{u}}_t,\mathit{\delta })$ is the likelihood and $\mathcal{P}({\mathbf{u}}_t,\mathit{\delta }\mid {\mathbf{Y}}_e)$ is the posterior probability distribution.

Prior probability in ${\mathbf{u}}_t$ and $\mathit{\delta }$ is specified using independent Gaussian processes. However, a complete definition of the prior requires the specification of uncertain hyper-parameters of the probability distribution. Let ${\mathit{\theta }}_u\in {\mathbf{\Theta }}_{\mathit{u}}$ and ${\mathit{\theta }}_{\mathit{\delta }}\in {\mathbf{\Theta }}_{\mathit{\delta }}$ be the hyper-parameters of subsystem model and discrepancy function respectively, where ${\mathbf{\Theta }}_{\mathit{u}}$ and ${\mathbf{\Theta }}_{\mathit{\delta }}$ are set of possible values. Thus, in the presence of the uncertain hyper-parameters, posterior probability distribution (5(5) $$\begin{array}{c}\hfill \mathcal{P}({\mathbf{u}}_t,\mathit{\delta }\mid {\mathbf{Y}}_e)\propto \mathcal{P}({\mathbf{Y}}_e\mid {\mathbf{u}}_t,\mathit{\delta })\times \mathcal{P}({\mathbf{u}}_t,\mathit{\delta }),\end{array}$$(5) ) takes the form(6) $$\begin{array}{cc}\hfill & \mathcal{P}({\mathbf{u}}_t({\mathit{\theta }}_u),\mathit{\delta }({\mathit{\theta }}_{\mathit{\delta }}),{\mathit{\theta }}_u,{\mathit{\theta }}_{\mathit{\delta }}\mid {\mathbf{Y}}_e)\propto \mathcal{P}({\mathbf{Y}}_e\mid {\mathbf{u}}_t({\mathit{\theta }}_u),\mathit{\delta }({\mathit{\theta }}_{\mathit{\delta }}),{\mathit{\theta }}_u,{\mathit{\theta }}_{\mathit{\delta }})\hfill \\ \hfill & \times \mathcal{P}({\mathbf{u}}_t({\mathit{\theta }}_u),\mathit{\delta }({\mathit{\theta }}_{\mathit{\delta }})\mid {\mathit{\theta }}_u,{\mathit{\theta }}_{\mathit{\delta }})\times \mathcal{P}({\mathit{\theta }}_u)\times \mathcal{P}({\mathit{\theta }}_{\mathit{\delta }}).\hfill \end{array}$$(6) For better readability, dependence of the uncertain subsystem model and the discrepancy function on uncertain hyper-parameters is explicitly shown henceforth.

In the present paper, Bayesian inference is developed assuming a Gaussian process prior for the discrepancy functions, as it is extensively used in the literature for specification of prior on random functions.[10, 41–43] The prior uncertainty in the subsystem model and the discrepancy function is assumed to be independent. Uncertainty in the experimental observations is specified by a zero-mean normally distributed random variable ${\mathit{ϵ}}_j$ with standard deviation ${\mathit{\sigma }}_{e_j}$. The experimental uncertainty at different input settings $\mathbf{x}$ is assumed to be uncorrelated. By marginalization over ${\mathit{\delta }}_j(\mathbf{x})$, the posterior distribution is given by(7) $$\begin{array}{cc}\hfill & \mathcal{P}({\mathbf{u}}_t({\mathit{\theta }}_u),\mathit{\delta }({\mathit{\theta }}_{\mathit{\delta }}),{\mathit{\theta }}_u,{\mathit{\theta }}_{\mathit{\delta }}\mid {\mathbf{Y}}_e)\propto \prod _{j=1}^M\frac{1}{\sqrt{\mid {\mathrm{\Sigma }}_j\mid }}\text{exp}\{-\frac{1}{2}{({\mathbf{Y}}_{e_j}-{\mathit{\mu }}_j)}^T{\mathrm{\Sigma }}_j^{-1}({\mathbf{Y}}_{e_j}-{\mathit{\mu }}_j)\}\hfill \\ \hfill & \times \mathcal{P}({\mathbf{u}}_t({\mathit{\theta }}_u)\mid {\mathit{\theta }}_u)\times \mathcal{P}({\mathit{\theta }}_u)\times \mathcal{P}({\mathit{\theta }}_{\mathit{\delta }}),\hfill \end{array}$$(7) where ${\mathit{\mu }}_j=\{T_j({\mathbf{x}}_i,u_t({\mathbf{x}}_s;{\mathit{\theta }}_u))+E({\mathit{\delta }}_j({\mathbf{x}}_i;{\mathit{\theta }}_{\mathit{\delta }}));i=1,...,N\}$ while ${\mathrm{\Sigma }}_j={\mathrm{\Sigma }}_{{\mathit{\delta }}_j}+{\mathit{\sigma }}_j^2{I}_N$, $I_N$ being the $N\times N$ identity matrix.

The Metropolis-Hastings algorithm [44, 45] is used to sample from the posterior distribution (7(7) $$\begin{array}{cc}\hfill & \mathcal{P}({\mathbf{u}}_t({\mathit{\theta }}_u),\mathit{\delta }({\mathit{\theta }}_{\mathit{\delta }}),{\mathit{\theta }}_u,{\mathit{\theta }}_{\mathit{\delta }}\mid {\mathbf{Y}}_e)\propto \prod _{j=1}^M\frac{1}{\sqrt{\mid {\mathrm{\Sigma }}_j\mid }}\text{exp}\{-\frac{1}{2}{({\mathbf{Y}}_{e_j}-{\mathit{\mu }}_j)}^T{\mathrm{\Sigma }}_j^{-1}({\mathbf{Y}}_{e_j}-{\mathit{\mu }}_j)\}\hfill \\ \hfill & \times \mathcal{P}({\mathbf{u}}_t({\mathit{\theta }}_u)\mid {\mathit{\theta }}_u)\times \mathcal{P}({\mathit{\theta }}_u)\times \mathcal{P}({\mathit{\theta }}_{\mathit{\delta }}),\hfill \end{array}$$(7) ). For each sample, evaluation of $T_j(\mathbf{x},u_t({\mathbf{x}}_s;{\mathit{\theta }}_u))$, ${\mathrm{\Sigma }}_j^{-1}$ and $\mid {\mathrm{\Sigma }}_j\mid$ impose significant computational expenses on the Bayesian framework. In the present paper, a gPC expansion based method is proposed for computationally efficient evaluations of $T_j(\mathbf{x},u_t({\mathbf{x}}_s;{\mathit{\theta }}_u))$, ${\mathrm{\Sigma }}_j^{-1}$ and $\mid {\mathrm{\Sigma }}_j\mid$.

Generalized polynomial Chaos expansion of a Gaussian process with uncertain hyper-parameters

For brevity, the proposed spectral formulation is described first for a susbsystem model, $u(\mathbf{x};{\mathit{\theta }}_u)$, which is subsequently extended for a discrepancy function. The formulation is derived for a zero-mean Gaussian process, however, it can easily be extended for non-zero mean processes. The derivation uses KL expansion of the Gaussian process with uncertain hyper-parameters, which is projected on a gPC basis using the intrusive Galerkin projection.

KL expansion

Let $u(\mathbf{x};{\mathit{\theta }}_u)$ be a zero-mean Gaussian process with a covariance function $C_u({\mathbf{x}}_1,{\mathbf{x}}_2;{\mathit{\theta }}_u)$. The covariance function can be approximated as [30](8) $$\begin{array}{c}\hfill C_u({\mathbf{x}}_1,{\mathbf{x}}_2;{\mathit{\theta }}_u)\approx \sum _{n=1}^N{\mathit{\lambda }}_n({\mathit{\theta }}_u)e_n({\mathbf{x}}_1,{\mathit{\theta }}_u)e_n({\mathbf{x}}_2,{\mathit{\theta }}_u),\end{array}$$(8) where $N$ is the number of expansion terms retained in the spectral approximation. ${\mathit{\lambda }}_n({\mathit{\theta }}_u)$ and $e_n(\mathbf{x},{\mathit{\theta }}_u)$ are eigenvalues and eigenfunctions of the covariance kernel, which are given by solution of the Fredholm’s integral equation of the second kind [46](9) $$\begin{array}{c}\hfill {\int }_{\mathcal{X}}{C}_u({\mathbf{x}}_1,{\mathbf{x}}_2;{\mathit{\theta }}_u)e_n({\mathbf{x}}_1,{\mathit{\theta }}_u)d{\mathbf{x}}_1={\mathit{\lambda }}_n({\mathit{\theta }}_u)e_n({\mathbf{x}}_2,{\mathit{\theta }}_u).\end{array}$$(9) Explicit dependence of the eigenvalues and eigenfunctions on the hyper-parameters ${\mathit{\theta }}_u$ should be noted. The resultant KL expansion using (8(8) $$\begin{array}{c}\hfill C_u({\mathbf{x}}_1,{\mathbf{x}}_2;{\mathit{\theta }}_u)\approx \sum _{n=1}^N{\mathit{\lambda }}_n({\mathit{\theta }}_u)e_n({\mathbf{x}}_1,{\mathit{\theta }}_u)e_n({\mathbf{x}}_2,{\mathit{\theta }}_u),\end{array}$$(8) ) is given by [30](10) $$\begin{array}{c}\hfill u(\mathbf{x};{\mathit{\theta }}_u)\approx \sum _{n=1}^N\sqrt{{\mathit{\lambda }}_n({\mathit{\theta }}_u)}{e}_n(\mathbf{x},{\mathit{\theta }}_u){\mathit{\chi }}_n,\end{array}$$(10) where ${\mathit{\chi }}_n$ are independent zero-mean standard normal random variables. For a given ${\mathit{\theta }}_u$, the eigenvalue problem (9(9) $$\begin{array}{c}\hfill {\int }_{\mathcal{X}}{C}_u({\mathbf{x}}_1,{\mathbf{x}}_2;{\mathit{\theta }}_u)e_n({\mathbf{x}}_1,{\mathit{\theta }}_u)d{\mathbf{x}}_1={\mathit{\lambda }}_n({\mathit{\theta }}_u)e_n({\mathbf{x}}_2,{\mathit{\theta }}_u).\end{array}$$(9) ) can be numerically solved using a Galerkin projection based approach.[47] In this paper, the approach is extended for uncertain hyper-parameters as follows.

Eigenfunctions $e_n(\mathbf{x},{\mathit{\theta }}_u)$ can be spectrally approximated as(11) $$\begin{array}{c}\hfill e_n(\mathbf{x},{\mathit{\theta }}_u)\approx \sum _{i=1}^N{d}_i^n({\mathit{\theta }}_u){\mathit{\psi }}_i(\mathbf{x}),\end{array}$$(11) where ${\mathit{\psi }}_i(\mathbf{x})$ are Legendre polynomials and $d_i^n({\mathit{\theta }}_u)$ are respective expansion coefficients.[30, 47] Use (11(11) $$\begin{array}{c}\hfill e_n(\mathbf{x},{\mathit{\theta }}_u)\approx \sum _{i=1}^N{d}_i^n({\mathit{\theta }}_u){\mathit{\psi }}_i(\mathbf{x}),\end{array}$$(11) ) in (9(9) $$\begin{array}{c}\hfill {\int }_{\mathcal{X}}{C}_u({\mathbf{x}}_1,{\mathbf{x}}_2;{\mathit{\theta }}_u)e_n({\mathbf{x}}_1,{\mathit{\theta }}_u)d{\mathbf{x}}_1={\mathit{\lambda }}_n({\mathit{\theta }}_u)e_n({\mathbf{x}}_2,{\mathit{\theta }}_u).\end{array}$$(9) ), multiply both sides by ${\mathit{\psi }}_j({\mathbf{x}}_2)$ and integrate w.r.t. $d{\mathbf{x}}_2$ to obtain(12) $$\begin{array}{cc}\hfill & \sum _{i=1}^N{d}_i^n({\mathit{\theta }}_u){\int }_{\mathcal{X}}{\int }_{\mathcal{X}}{C}_u({\mathbf{x}}_1,{\mathbf{x}}_2;{\mathit{\theta }}_u){\mathit{\psi }}_i({\mathbf{x}}_1){\mathit{\psi }}_j({\mathbf{x}}_2)d{\mathbf{x}}_{1}d{\mathbf{x}}_2\hfill \\ \hfill & ={\mathit{\lambda }}_n({\mathit{\theta }}_u)\sum _{i=1}^N{d}_i^n({\mathit{\theta }}_u){\int }_{\mathcal{X}}{\mathit{\psi }}_i({\mathbf{x}}_2){\mathit{\psi }}_j({\mathbf{x}}_2)d{\mathbf{x}}_2.\hfill \end{array}$$(12) Using$$\begin{array}{ccc}\hfill & A_{\mathit{ij}}({\mathit{\theta }}_u)={\int }_{\mathcal{X}}{\int }_{\mathcal{X}}C({\mathbf{x}}_1,{\mathbf{x}}_2;{\mathit{\theta }}_u){\mathit{\psi }}_i({\mathbf{x}}_1){\mathit{\psi }}_j({\mathbf{x}}_2)d{\mathbf{x}}_{1}d{\mathbf{x}}_2,\hfill & \hfill D_{\mathit{ij}}({\mathit{\theta }}_u)=d_i^j({\mathit{\theta }}_u),\\ \hfill & B_{\mathit{ij}}({\mathit{\theta }}_u)={\int }_{\mathcal{X}}{\mathit{\psi }}_i({\mathbf{x}}_2){\mathit{\psi }}_j({\mathbf{x}}_2)d{\mathbf{x}}_2,\hfill & \hfill {\mathrm{\Lambda }}_{\mathit{ii}}({\mathit{\theta }}_u)={\mathit{\lambda }}_i({\mathit{\theta }}_u),\end{array}$$(12(12) $$\begin{array}{cc}\hfill & \sum _{i=1}^N{d}_i^n({\mathit{\theta }}_u){\int }_{\mathcal{X}}{\int }_{\mathcal{X}}{C}_u({\mathbf{x}}_1,{\mathbf{x}}_2;{\mathit{\theta }}_u){\mathit{\psi }}_i({\mathbf{x}}_1){\mathit{\psi }}_j({\mathbf{x}}_2)d{\mathbf{x}}_{1}d{\mathbf{x}}_2\hfill \\ \hfill & ={\mathit{\lambda }}_n({\mathit{\theta }}_u)\sum _{i=1}^N{d}_i^n({\mathit{\theta }}_u){\int }_{\mathcal{X}}{\mathit{\psi }}_i({\mathbf{x}}_2){\mathit{\psi }}_j({\mathbf{x}}_2)d{\mathbf{x}}_2.\hfill \end{array}$$(12) ) can be written in a matrix form as(13) $$\begin{array}{c}\hfill A({\mathit{\theta }}_u)D({\mathit{\theta }}_u)=\mathrm{\Lambda }({\mathit{\theta }}_u)B({\mathit{\theta }}_u)D({\mathit{\theta }}_u),\end{array}$$(13) which is a generalized eigenvalue problem (GEP) that can be solved using the QZ algorithm.[30]

The matrices in (13(13) $$\begin{array}{c}\hfill A({\mathit{\theta }}_u)D({\mathit{\theta }}_u)=\mathrm{\Lambda }({\mathit{\theta }}_u)B({\mathit{\theta }}_u)D({\mathit{\theta }}_u),\end{array}$$(13) ) are functions of ${\mathit{\theta }}_u$. Thus, to solve (13(13) $$\begin{array}{c}\hfill A({\mathit{\theta }}_u)D({\mathit{\theta }}_u)=\mathrm{\Lambda }({\mathit{\theta }}_u)B({\mathit{\theta }}_u)D({\mathit{\theta }}_u),\end{array}$$(13) ), consider spectral expansion of eigenvalues and eigenfunctions as(14) $$\begin{array}{c}\hfill {\mathit{\lambda }}_n({\mathit{\theta }}_u)\approx \sum _{i=1}^L{l}_i^n{\mathit{\varphi }}_i({\mathit{\theta }}_u),d_k^n({\mathit{\theta }}_u)\approx \sum _{i=1}^L{c}_{i,k}^n{\mathit{\varphi }}_i({\mathit{\theta }}_u),\end{array}$$(14) where ${\mathit{\varphi }}_i({\mathit{\theta }}_u)$ are Legendre polynomials that form complete orthonormal basis on $L^2({\mathit{\theta }}_u)$. To ensure orthonormality of the Legendre polynomials, ${\mathit{\varphi }}_i({\mathit{\theta }}_u)$ are scaled using minimum and maximum value of ${\mathit{\theta }}_u$. Though the method assumes information about the range of hyper-parameters ${\mathit{\theta }}_u$, specification of complete prior distribution is not necessary. Thus, a conservative estimate of the range suffices at this juncture. The coefficients are given by(15) $$\begin{array}{c}\hfill l_i^n=\frac{{\int }_{{\mathit{\theta }}_u}{\mathit{\lambda }}_n({\mathit{\theta }}_u){\mathit{\varphi }}_i({\mathit{\theta }}_u)d{\mathit{\theta }}_u}{{\int }_{{\mathit{\theta }}_u}{\mathit{\varphi }}_i^2({\mathit{\theta }}_u)d{\mathit{\theta }}_u},c_{i,k}^n=\frac{{\int }_{{\mathit{\theta }}_u}{d}_k^n({\mathit{\theta }}_u){\mathit{\varphi }}_i({\mathit{\theta }}_u)d{\mathit{\theta }}_u}{{\int }_{{\mathit{\theta }}_u}{\mathit{\varphi }}_i^2({\mathit{\theta }}_u)d{\mathit{\theta }}_u}.\end{array}$$(15) Using the Gauss-Legendre quadrature, the integrals can be approximated as(16) $$\begin{array}{cc}\hfill & {\int }_{{\mathit{\theta }}_u}{\mathit{\lambda }}_i({\mathit{\theta }}_u){\mathit{\varphi }}_i({\mathit{\theta }}_u)d{\mathit{\theta }}_u\approx \sum _{q=1}^{N_q}{\mathit{\lambda }}_i\left({\mathit{\theta }}_u^q\right){\mathit{\varphi }}_i\left({\mathit{\theta }}_u^q\right)w_q,\hfill \\ \hfill & {\int }_{{\mathit{\theta }}_u}{d}_k^n({\mathit{\theta }}_u){\mathit{\varphi }}_i({\mathit{\theta }}_u)d{\mathit{\theta }}_u\approx \sum _{q=1}^{N_q}{d}_k^n\left({\mathit{\theta }}_u^q\right){\mathit{\varphi }}_i\left({\mathit{\theta }}_u^q\right)w_q\hfill \end{array}$$(16) where $N_q$ are the number of quadrature points used, ${\mathit{\theta }}_u^q$ are the quadrature nodes, while $w_q$ are the respective quadrature weights. The expansion coefficients, $l_i^n$ and $c_{i,k}^n$, are calculated by solving (13(13) $$\begin{array}{c}\hfill A({\mathit{\theta }}_u)D({\mathit{\theta }}_u)=\mathrm{\Lambda }({\mathit{\theta }}_u)B({\mathit{\theta }}_u)D({\mathit{\theta }}_u),\end{array}$$(13) ) at quadrature nodes ${\mathit{\theta }}_u^q$ and substituting the solution in (16(16) $$\begin{array}{cc}\hfill & {\int }_{{\mathit{\theta }}_u}{\mathit{\lambda }}_i({\mathit{\theta }}_u){\mathit{\varphi }}_i({\mathit{\theta }}_u)d{\mathit{\theta }}_u\approx \sum _{q=1}^{N_q}{\mathit{\lambda }}_i\left({\mathit{\theta }}_u^q\right){\mathit{\varphi }}_i\left({\mathit{\theta }}_u^q\right)w_q,\hfill \\ \hfill & {\int }_{{\mathit{\theta }}_u}{d}_k^n({\mathit{\theta }}_u){\mathit{\varphi }}_i({\mathit{\theta }}_u)d{\mathit{\theta }}_u\approx \sum _{q=1}^{N_q}{d}_k^n\left({\mathit{\theta }}_u^q\right){\mathit{\varphi }}_i\left({\mathit{\theta }}_u^q\right)w_q\hfill \end{array}$$(16) ) to evaluate integrals in (15(15) $$\begin{array}{c}\hfill l_i^n=\frac{{\int }_{{\mathit{\theta }}_u}{\mathit{\lambda }}_n({\mathit{\theta }}_u){\mathit{\varphi }}_i({\mathit{\theta }}_u)d{\mathit{\theta }}_u}{{\int }_{{\mathit{\theta }}_u}{\mathit{\varphi }}_i^2({\mathit{\theta }}_u)d{\mathit{\theta }}_u},c_{i,k}^n=\frac{{\int }_{{\mathit{\theta }}_u}{d}_k^n({\mathit{\theta }}_u){\mathit{\varphi }}_i({\mathit{\theta }}_u)d{\mathit{\theta }}_u}{{\int }_{{\mathit{\theta }}_u}{\mathit{\varphi }}_i^2({\mathit{\theta }}_u)d{\mathit{\theta }}_u}.\end{array}$$(15) ).

gPC expansion of a Gaussian process

Hyper-parameters ${\mathit{\theta }}_u$ can be expanded in gPC basis as [34](17) $$\begin{array}{c}\hfill {\mathit{\theta }}_u\approx \sum _{p=1}^P{\hat{\mathit{\theta }}}_p^u{H}_p(\mathit{\xi }),\end{array}$$(17) where $H_p(\mathit{\xi })$ are the Hermite polynomials, ${\hat{\mathit{\theta }}}_p^u$ are respective expansion coefficients, while $\mathit{\xi }$ is a vector of independent standard normal random variables. Using the Galerkin projection, gPC expansion of ${\mathit{\varphi }}_i({\mathit{\theta }}_u)$ is given by(18) $$\begin{array}{c}\hfill {\mathit{\varphi }}_i({\mathit{\theta }}_u)\approx \sum _{p=1}^P{\hat{\mathit{\varphi }}}_p^i{H}_p(\mathit{\xi }).\end{array}$$(18) Using (14(14) $$\begin{array}{c}\hfill {\mathit{\lambda }}_n({\mathit{\theta }}_u)\approx \sum _{i=1}^L{l}_i^n{\mathit{\varphi }}_i({\mathit{\theta }}_u),d_k^n({\mathit{\theta }}_u)\approx \sum _{i=1}^L{c}_{i,k}^n{\mathit{\varphi }}_i({\mathit{\theta }}_u),\end{array}$$(14) ) and (18(18) $$\begin{array}{c}\hfill {\mathit{\varphi }}_i({\mathit{\theta }}_u)\approx \sum _{p=1}^P{\hat{\mathit{\varphi }}}_p^i{H}_p(\mathit{\xi }).\end{array}$$(18) ), gPC expansion of $\sqrt{{\mathit{\lambda }}_n({\mathit{\theta }}_u)}$ is given by(19) $$\begin{array}{c}\hfill \sqrt{{\mathit{\lambda }}_n({\mathit{\theta }}_u)}\approx \sum _{i=1}^L\sum _{p=1}^P{s}_i^n{\hat{\mathit{\varphi }}}_p^i{H}_p(\mathit{\xi }),\end{array}$$(19) where $s_k^n=〈\sqrt{{\mathit{\lambda }}_n({\mathit{\theta }}_u)},{\mathit{\varphi }}_k({\mathit{\theta }}_u)〉/〈{\mathit{\varphi }}_k^2({\mathit{\theta }}_u)〉,$ and $〈·,·〉$ denote an inner product. Similarly, the gPC expansion of eigenfunctions can be obtained using (14(14) $$\begin{array}{c}\hfill {\mathit{\lambda }}_n({\mathit{\theta }}_u)\approx \sum _{i=1}^L{l}_i^n{\mathit{\varphi }}_i({\mathit{\theta }}_u),d_k^n({\mathit{\theta }}_u)\approx \sum _{i=1}^L{c}_{i,k}^n{\mathit{\varphi }}_i({\mathit{\theta }}_u),\end{array}$$(14) ) and (18(18) $$\begin{array}{c}\hfill {\mathit{\varphi }}_i({\mathit{\theta }}_u)\approx \sum _{p=1}^P{\hat{\mathit{\varphi }}}_p^i{H}_p(\mathit{\xi }).\end{array}$$(18) ) as(20) $$\begin{array}{c}\hfill e_n(\mathbf{x},{\mathit{\theta }}_u)\approx \sum _{k=1}^N\sum _{i=1}^L\sum _{p=1}^P{c}_{i,k}^n{\hat{\mathit{\varphi }}}_p^i{H}_p(\mathit{\xi }){\mathit{\psi }}_k(\mathbf{x}).\end{array}$$(20) Thus, a zero-mean Gaussian process $u(\mathbf{x};{\mathit{\theta }}_u)$ can be expanded in gPC basis as(21) $$\begin{array}{c}\hfill u(\mathbf{x};{\mathit{\theta }}_u)\approx \sum _{p=1}^P{\hat{u}}_p(\mathbf{x})H_p(\mathit{\xi }),\end{array}$$(21) where the expansion coefficients ${\hat{u}}_k(\mathbf{x})$ are given by using (10(10) $$\begin{array}{c}\hfill u(\mathbf{x};{\mathit{\theta }}_u)\approx \sum _{n=1}^N\sqrt{{\mathit{\lambda }}_n({\mathit{\theta }}_u)}{e}_n(\mathbf{x},{\mathit{\theta }}_u){\mathit{\chi }}_n,\end{array}$$(10) ), (19(19) $$\begin{array}{c}\hfill \sqrt{{\mathit{\lambda }}_n({\mathit{\theta }}_u)}\approx \sum _{i=1}^L\sum _{p=1}^P{s}_i^n{\hat{\mathit{\varphi }}}_p^i{H}_p(\mathit{\xi }),\end{array}$$(19) ) and (20(20) $$\begin{array}{c}\hfill e_n(\mathbf{x},{\mathit{\theta }}_u)\approx \sum _{k=1}^N\sum _{i=1}^L\sum _{p=1}^P{c}_{i,k}^n{\hat{\mathit{\varphi }}}_p^i{H}_p(\mathit{\xi }){\mathit{\psi }}_k(\mathbf{x}).\end{array}$$(20) ) as(22) $$\begin{array}{c}\hfill {\hat{u}}_k(\mathbf{x})=\sum _{n=1}^N\sum _{m=1}^N\sum _{i=1}^L\sum _{j=1}^L\sum _{p=1}^P\sum _{q=1}^P{s}_i^n{c}_{j,m}^n{\hat{\mathit{\varphi }}}_p^i{\hat{\mathit{\varphi }}}_q^j{\mathit{\psi }}_m(\mathbf{x}){\mathcal{Q}}_{p,q,n,k},\end{array}$$(22) where$$\begin{array}{c}\hfill {\mathcal{Q}}_{p,q,n,k}={\int }_{\mathbf{\Xi }}{H}_p(\mathit{\xi })H_q(\mathit{\xi }){\mathit{\chi }}_n{H}_k(\mathit{\xi })\exp (-\frac{1}{2}{\mathit{\xi }}^T\mathit{\xi })d\mathit{\xi }\end{array}$$is a quadruple product.

Stochastic spectral projection based Bayesian calibration

The prior uncertainty in the subsystem model, $u_t({\mathbf{x}}_s;{\mathit{\theta }}_u)$, and the discrepancy function, ${\mathit{\delta }}_j(\mathbf{x};{\mathit{\theta }}_{\mathit{\delta }})$, is given by independent Gaussian processes. Using (21(21) $$\begin{array}{c}\hfill u(\mathbf{x};{\mathit{\theta }}_u)\approx \sum _{p=1}^P{\hat{u}}_p(\mathbf{x})H_p(\mathit{\xi }),\end{array}$$(21) ), the spectral expansion of the prior is given by(23) $$\begin{array}{c}\hfill u_t({\mathbf{x}}_s;{\mathit{\theta }}_u)\approx \sum _{p=1}^P{\hat{u}}_p({\mathbf{x}}_s)H_p(\mathit{\xi });{\mathit{\delta }}_j(\mathbf{x};{\mathit{\theta }}_{\mathit{\delta }})\approx \sum _{p=1}^P{\hat{\mathit{\delta }}}_p^j(\mathbf{x})H_p(\mathit{\xi }).\end{array}$$(23) Galerkin projection approach is used to propagate the prior uncertainty in $u_t({\mathbf{x}}_s;{\mathit{\theta }}_u)$ to the simulator predictions,[37] thus,(24) $$\begin{array}{c}\hfill T_j(\mathbf{x},u_t({\mathbf{x}}_s;{\mathit{\theta }}_u))\approx \sum _{p=1}^P{\hat{T}}_p^j(\mathbf{x})H_p(\mathit{\xi }).\end{array}$$(24) Using the gPC expansion of $T_j(\mathbf{x},u_t({\mathbf{x}}_s;{\mathit{\theta }}_u))$ and ${\mathit{\delta }}_j(\mathbf{x};{\mathit{\theta }}_{\mathit{\delta }})$ in (3(3) $$\begin{array}{c}\hfill {\mathit{\zeta }}_j(\mathbf{x},u_t({\mathbf{x}}_s))=T_j(\mathbf{x},u_t({\mathbf{x}}_s))+{\mathit{\delta }}_j(\mathbf{x}),\end{array}$$(3) ) to obtain the gPC expansion of ${\mathit{\zeta }}_j(\mathbf{x},u_t({\mathbf{x}}_s;{\mathit{\theta }}_u))$ as(25) $$\begin{array}{c}\hfill {\mathit{\zeta }}_j(\mathbf{x},u_t({\mathbf{x}}_s;{\mathit{\theta }}_u))\approx \sum _{p=1}^P{\hat{\mathit{\zeta }}}_p^j(\mathbf{x})H_p(\mathit{\xi }),\end{array}$$(25) where ${\hat{\mathit{\zeta }}}_p^j(\mathbf{x})={\hat{T}}_p^j(\mathbf{x})+{\hat{\mathit{\delta }}}_p^j(\mathbf{x}).$ Use the gPC expansion (25(25) $$\begin{array}{c}\hfill {\mathit{\zeta }}_j(\mathbf{x},u_t({\mathbf{x}}_s;{\mathit{\theta }}_u))\approx \sum _{p=1}^P{\hat{\mathit{\zeta }}}_p^j(\mathbf{x})H_p(\mathit{\xi }),\end{array}$$(25) ) in (4(4) $$\begin{array}{c}\hfill y_{e_j}(\mathbf{x})={\mathit{\zeta }}_j(\mathbf{x},u_t({\mathbf{x}}_s))+{\mathit{ϵ}}_j(\mathbf{x}).\end{array}$$(4) ) to obtain(26) $$\begin{array}{c}\hfill y_{e_j}(\mathbf{x})\approx \sum _{p=1}^P{\hat{\mathit{\zeta }}}_p^j(\mathbf{x})H_p(\mathit{\xi })+{\mathit{ϵ}}_j.\end{array}$$(26) Note that $\mathit{\xi }$ are the only uncertain variables in (26(26) $$\begin{array}{c}\hfill y_{e_j}(\mathbf{x})\approx \sum _{p=1}^P{\hat{\mathit{\zeta }}}_p^j(\mathbf{x})H_p(\mathit{\xi })+{\mathit{ϵ}}_j.\end{array}$$(26) ), thus, the Bayesian inference problem is reformulated as sampling from the posterior distribution of $\mathit{\xi }$.

Let $\hat{\mathbf{T}}=\{{\hat{T}}_p^j(\mathbf{x});p=1,...,P;j=1,...,M\}$ and $\hat{\mathit{\delta }}=\{{\hat{\mathit{\delta }}}_p^j(\mathbf{x});p=1,...,P;j=1,...,M\}$ define the sets of respective gPC coefficients. Using the gPC expansion of the uncertain variables ((23(23) $$\begin{array}{c}\hfill u_t({\mathbf{x}}_s;{\mathit{\theta }}_u)\approx \sum _{p=1}^P{\hat{u}}_p({\mathbf{x}}_s)H_p(\mathit{\xi });{\mathit{\delta }}_j(\mathbf{x};{\mathit{\theta }}_{\mathit{\delta }})\approx \sum _{p=1}^P{\hat{\mathit{\delta }}}_p^j(\mathbf{x})H_p(\mathit{\xi }).\end{array}$$(23) ) – (26(26) $$\begin{array}{c}\hfill y_{e_j}(\mathbf{x})\approx \sum _{p=1}^P{\hat{\mathit{\zeta }}}_p^j(\mathbf{x})H_p(\mathit{\xi })+{\mathit{ϵ}}_j.\end{array}$$(26) )) in (7(7) $$\begin{array}{cc}\hfill & \mathcal{P}({\mathbf{u}}_t({\mathit{\theta }}_u),\mathit{\delta }({\mathit{\theta }}_{\mathit{\delta }}),{\mathit{\theta }}_u,{\mathit{\theta }}_{\mathit{\delta }}\mid {\mathbf{Y}}_e)\propto \prod _{j=1}^M\frac{1}{\sqrt{\mid {\mathrm{\Sigma }}_j\mid }}\text{exp}\{-\frac{1}{2}{({\mathbf{Y}}_{e_j}-{\mathit{\mu }}_j)}^T{\mathrm{\Sigma }}_j^{-1}({\mathbf{Y}}_{e_j}-{\mathit{\mu }}_j)\}\hfill \\ \hfill & \times \mathcal{P}({\mathbf{u}}_t({\mathit{\theta }}_u)\mid {\mathit{\theta }}_u)\times \mathcal{P}({\mathit{\theta }}_u)\times \mathcal{P}({\mathit{\theta }}_{\mathit{\delta }}),\hfill \end{array}$$(7) ), the Bayesian calibration problem can be reformulated in terms of $\mathit{\xi }$ that capture all the randomness in the system model and hyper-parameters:(27) $$\begin{array}{cc}\hfill \mathcal{P}(\mathit{\xi }\mid {\mathbf{Y}}_e,\hat{\mathbf{T}},\hat{\mathit{\delta }})& \propto \mathcal{P}({\mathbf{Y}}_e\mid \mathit{\xi },\hat{\mathbf{T}},\hat{\mathit{\delta }})\times \mathcal{P}(\mathit{\xi }).\hfill \end{array}$$(27) Since $\mathit{\xi }$ are i.i.d. standard normal random variables with independent Gaussian uncertainties in experimental data, the posterior distribution of $\mathit{\xi }$ can be obtained by:(28) $$\begin{array}{cc}\hfill & \mathcal{P}(\mathit{\xi }\mid {\mathbf{Y}}_e,\hat{\mathbf{T}},\hat{\mathit{\delta }})\propto \prod _{j=1}^M\frac{1}{\sqrt{|{\mathrm{\Sigma }}_j|}}\exp \{-\frac{1}{2}{({\mathbf{Y}}_{e_j}-{\mathit{\mu }}_j)}^T{\mathrm{\Sigma }}_j^{-1}({\mathbf{Y}}_{e_j}-{\mathit{\mu }}_j)\}\hfill \\ \hfill & \times \prod _{n=1}^{N_d}{e}^{-{\mathit{\xi }}_n^2/2},\hfill \end{array}$$(28) where ${\mathit{\mu }}_j=\{\sum _{p=1}^P{\hat{T}}_p^j({\mathbf{x}}_i)H_p(\mathit{\xi })+{\hat{\mathit{\delta }}}_0^j({\mathbf{x}}_i);i=1,\cdots ,N\}$.

Markov Chain Monte Carlo (MCMC) method is used to sample from (28(28) $$\begin{array}{cc}\hfill & \mathcal{P}(\mathit{\xi }\mid {\mathbf{Y}}_e,\hat{\mathbf{T}},\hat{\mathit{\delta }})\propto \prod _{j=1}^M\frac{1}{\sqrt{|{\mathrm{\Sigma }}_j|}}\exp \{-\frac{1}{2}{({\mathbf{Y}}_{e_j}-{\mathit{\mu }}_j)}^T{\mathrm{\Sigma }}_j^{-1}({\mathbf{Y}}_{e_j}-{\mathit{\mu }}_j)\}\hfill \\ \hfill & \times \prod _{n=1}^{N_d}{e}^{-{\mathit{\xi }}_n^2/2},\hfill \end{array}$$(28) ). Note that MCMC applied to (28(28) $$\begin{array}{cc}\hfill & \mathcal{P}(\mathit{\xi }\mid {\mathbf{Y}}_e,\hat{\mathbf{T}},\hat{\mathit{\delta }})\propto \prod _{j=1}^M\frac{1}{\sqrt{|{\mathrm{\Sigma }}_j|}}\exp \{-\frac{1}{2}{({\mathbf{Y}}_{e_j}-{\mathit{\mu }}_j)}^T{\mathrm{\Sigma }}_j^{-1}({\mathbf{Y}}_{e_j}-{\mathit{\mu }}_j)\}\hfill \\ \hfill & \times \prod _{n=1}^{N_d}{e}^{-{\mathit{\xi }}_n^2/2},\hfill \end{array}$$(28) ) does not require solution of the simulation model, $T_j(\mathbf{x},u({\mathbf{x}}_s,{\mathit{\theta }}_u))$, thus, the posterior distribution can be explored efficiently. However, solution of (28(28) $$\begin{array}{cc}\hfill & \mathcal{P}(\mathit{\xi }\mid {\mathbf{Y}}_e,\hat{\mathbf{T}},\hat{\mathit{\delta }})\propto \prod _{j=1}^M\frac{1}{\sqrt{|{\mathrm{\Sigma }}_j|}}\exp \{-\frac{1}{2}{({\mathbf{Y}}_{e_j}-{\mathit{\mu }}_j)}^T{\mathrm{\Sigma }}_j^{-1}({\mathbf{Y}}_{e_j}-{\mathit{\mu }}_j)\}\hfill \\ \hfill & \times \prod _{n=1}^{N_d}{e}^{-{\mathit{\xi }}_n^2/2},\hfill \end{array}$$(28) ) requires numerical evaluation of $|{\mathrm{\Sigma }}_j|$ and ${\mathrm{\Sigma }}_j^{-1}$, which may impose non-trivial computational cost on the MCMC sampling. In the present paper, numerical evaluation of $|{\mathrm{\Sigma }}_j|$ and ${\mathrm{\Sigma }}_j^{-1}$ is accelerated as follows.

The numerical evaluation of $|{\mathrm{\Sigma }}_j|$ and ${\mathrm{\Sigma }}_j^{-1}$ is accelerated using the gPC expansion of the individual elements of the covariance matrix ${\mathrm{\Sigma }}_j({\mathbf{x}}_1,{\mathbf{x}}_2)$ as(29) $$\begin{array}{c}\hfill {\mathrm{\Sigma }}_j({\mathbf{x}}_1,{\mathbf{x}}_2)=\sum _{p=1}^P{\hat{C}}_p^j({\mathbf{x}}_1,{\mathbf{x}}_2)H_p(\mathit{\xi }).\end{array}$$(29) The determinant $|{\mathrm{\Sigma }}_j|$ and the individual elements of the inverse ${\mathrm{\Sigma }}_j^{-1}$ can also be expanded in gPC basis as(30) $$\begin{array}{c}\hfill |{\mathrm{\Sigma }}_j|\approx \sum _{p=1}^P{\hat{D}}_p^j{H}_p(\mathit{\xi }),{\mathrm{\Sigma }}_j^{-1}({\mathbf{x}}_1,{\mathbf{x}}_2)\approx \sum _{p=1}^P{\hat{I}}_p^j({\mathbf{x}}_1,{\mathbf{x}}_2)H_p(\mathit{\xi }),\end{array}$$(30) where the gPC expansion coefficients are given by(31) $$\begin{array}{c}\hfill {\hat{I}}_k^j=\frac{〈{\mathrm{\Sigma }}_j^{-1},H_k(\mathit{\xi })〉}{〈H_k^2〉},{\hat{D}}_k^j=\frac{〈|{\mathrm{\Sigma }}_j|,H_k(\mathit{\xi })〉}{〈H_k^2〉}.\end{array}$$(31) Gaussian quadrature is used to evaluate Polynomial Chaos coefficients. $〈{\mathrm{\Sigma }}_j^{-1},H_k(\mathit{\xi })〉$ and $〈|{\mathrm{\Sigma }}_j|,H_k(\mathit{\xi })〉$ are calculated by evaluating ${\mathrm{\Sigma }}_j^{-1}$ and $|{\mathrm{\Sigma }}_j|$ at quadrature nodes, while ${\hat{I}}_k^j$ and ${\hat{D}}_k^j$ are evaluated using (31(31) $$\begin{array}{c}\hfill {\hat{I}}_k^j=\frac{〈{\mathrm{\Sigma }}_j^{-1},H_k(\mathit{\xi })〉}{〈H_k^2〉},{\hat{D}}_k^j=\frac{〈|{\mathrm{\Sigma }}_j|,H_k(\mathit{\xi })〉}{〈H_k^2〉}.\end{array}$$(31) ). The resultant gPC expansion (30(30) $$\begin{array}{c}\hfill |{\mathrm{\Sigma }}_j|\approx \sum _{p=1}^P{\hat{D}}_p^j{H}_p(\mathit{\xi }),{\mathrm{\Sigma }}_j^{-1}({\mathbf{x}}_1,{\mathbf{x}}_2)\approx \sum _{p=1}^P{\hat{I}}_p^j({\mathbf{x}}_1,{\mathbf{x}}_2)H_p(\mathit{\xi }),\end{array}$$(30) ) is used in (28(28) $$\begin{array}{cc}\hfill & \mathcal{P}(\mathit{\xi }\mid {\mathbf{Y}}_e,\hat{\mathbf{T}},\hat{\mathit{\delta }})\propto \prod _{j=1}^M\frac{1}{\sqrt{|{\mathrm{\Sigma }}_j|}}\exp \{-\frac{1}{2}{({\mathbf{Y}}_{e_j}-{\mathit{\mu }}_j)}^T{\mathrm{\Sigma }}_j^{-1}({\mathbf{Y}}_{e_j}-{\mathit{\mu }}_j)\}\hfill \\ \hfill & \times \prod _{n=1}^{N_d}{e}^{-{\mathit{\xi }}_n^2/2},\hfill \end{array}$$(28) ) for MCMC sampling.

Numerical example: calibration of quasi-one-dimensional nozzle flow simulator

 

Problem setup

A quasi-one-dimensional supersonic flow through a nozzle is considered to verify the proposed Bayesian calibration method. The quasi-one-dimensional flow is uniform across the cross-section with properties varying in $x$-direction, while, the effect of change in the cross-sectional area is considered. The flow is defined using the compressible Euler equations in conservative form(32) $$\begin{array}{c}\hfill \frac{\mathit{\partial }\mathbf{q}}{\mathit{\partial }t}+\frac{\mathit{\partial }\mathbf{f}}{\mathit{\partial }x}=\mathbf{g},\end{array}$$(32) where,$$\begin{array}{ccc}\hfill \mathbf{q}=\left(\begin{array}{c}\hfill \mathit{\rho }A\hfill \\ \hfill \mathit{\rho }vA\hfill \\ \hfill \mathit{\rho }EA\hfill \end{array}\right),& \mathbf{f}=\left(\begin{array}{c}\hfill \mathit{\rho }vA\hfill \\ \hfill \mathit{\rho }{v}^{2}A+PA\hfill \\ \hfill \mathit{\rho }vEA+PvA\hfill \end{array}\right),\hfill & \hfill \mathbf{g}=\left(\begin{array}{c}\hfill 0\hfill \\ \hfill P\frac{\mathit{\partial }A}{\mathit{\partial }x}\hfill \\ \hfill 0\hfill \end{array}\right).\end{array}$$Here, $\mathit{\rho }$ is density, $v$ is velocity, $A$ is cross-sectional area, $P$ is static pressure and $E$ is the total energy per unit mass. In the governing equations, static pressure is substituted by the total energy per unit mass using(33) $$\begin{array}{c}\hfill E=\frac{P}{(\mathit{\gamma }-1)\mathit{\rho }}+\frac{1}{2}{v}^2,\end{array}$$(33) while the ideal gas equation is used for the closure.

Variation of nozzle area profile $A$ is assumed to be uncertain, which is inferred using the Bayesian framework. Stationary Gaussian process with known mean profile is used as a prior for uncertain nozzle area profile. The procedure of obtaining the stochastic spectral formulation is summarized as follows. Define $q_1=\mathit{\rho }A$, $q_2=\mathit{\rho }vA$, $q_3=\mathit{\rho }EA$ and the corresponding gPC expansions$$q_1\approx \sum _{p=1}^P{\hat{q}}_{1,p}{H}_p(\mathit{\xi });q_2\approx \sum _{p=1}^P{\hat{q}}_{2,p}{H}_p(\mathit{\xi });q_3\approx \sum _{p=1}^P{\hat{q}}_{3,p}{H}_p(\mathit{\xi }).$$Multiplying both the sides by $H_k(\mathit{\xi })$ and taking the inner product, governing equations of a stochastic quasi-one-dimensional nozzle flow is given by(34) $$\begin{array}{c}\hfill \frac{\mathit{\partial }\mathbf{Q}}{\mathit{\partial }t}+\frac{\mathit{\partial }\mathbf{F}}{\mathit{\partial }x}=\mathbf{G},\end{array}$$(34) where$$\begin{array}{c}\hfill {\left(\mathbf{Q}\right)}_p=\left(\begin{array}{c}\hfill {\hat{q}}_{1,p}\hfill \\ \hfill {\hat{q}}_{2,p}\hfill \\ \hfill {\hat{q}}_{3,p}\hfill \end{array}\right),\end{array}$$$$\begin{array}{c}\hfill {\left(\mathbf{F}\right)}_p=\left(\begin{array}{c}\hfill {\hat{q}}_{2,p}\hfill \\ \hfill \hfill \\ \hfill {\sum }_{i=1}^P{\sum }_{j=1}^P{\sum }_{k=1}^P{\hat{q}}_{2,i}{\hat{q}}_{2,j}{\hat{r}}_k\frac{〈H_i{H}_j{H}_k{H}_p〉}{〈H_p^2〉}\hfill \\ \hfill +\frac{\mathit{\gamma }-1}{\mathit{\gamma }}[{\hat{q}}_{3,p}-{\sum }_{i=1}^P{\sum }_{j=1}^P{\sum }_{k=1}^P{\hat{q}}_{2,i}{\hat{q}}_{2,j}{\hat{r}}_k\frac{〈H_i{H}_j{H}_k{H}_p〉}{〈H_p^2〉}]\hfill \\ \hfill \hfill \\ \hfill \mathit{\gamma }{\sum }_{i=1}^{N_p}{\sum }_{j=1}^{N_p}{\sum }_{k=1}^{N_p}{\hat{q}}_{2,i}{\hat{q}}_{3,j}{\hat{r}}_k\frac{〈H_i{H}_j{H}_k{H}_p〉}{〈H_p^2〉}\hfill \\ \hfill -\frac{\mathit{\gamma }(\mathit{\gamma }-1)}{2}\left[{\sum }_{i=1}^P{\sum }_{j=1}^P{\sum }_{k=1}^P{\sum }_{l=1}^P{\sum }_{m=1}^P{\hat{q}}_{2,i}{\hat{q}}_{2,j}{\hat{q}}_{2,k}{\hat{r}}_l{\hat{r}}_m\frac{〈H_i{H}_j{H}_k{H}_l{H}_m{H}_p〉}{〈H_p^2〉}\right]\hfill \end{array}\right),\end{array}$$$$\begin{array}{c}\hfill {\left(\mathbf{G}\right)}_p=\left(\begin{array}{c}\hfill 0\hfill \\ \hfill \hfill \\ \hfill \frac{\mathit{\gamma }-1}{\mathit{\gamma }}{\sum }_{i=1}^P{\sum }_{j=1}^P{\hat{q}}_{3,i}\mathit{\partial }{A}_j\mathit{\partial }x\frac{〈H_i{H}_j{H}_k{H}_p〉}{〈H_p^2〉}\hfill \\ \hfill -\frac{\mathit{\gamma }(\mathit{\gamma }-1)}{2\mathit{\gamma }}{\sum }_{i=1}^P{\sum }_{j=1}^P{\sum }_{k=1}^P{\sum }_{l=1}^P{\hat{q}}_{2,i}{\hat{q}}_{2,j}{\hat{r}}_k\mathit{\partial }{A}_l\mathit{\partial }x\frac{〈H_i{H}_j{H}_k{H}_l{H}_p〉}{〈H_p^2〉}\hfill \\ \hfill \hfill \\ \hfill 0\hfill \end{array}\right)\end{array}$$and $r=\frac{1}{q_1}$. Note that the governing equations is a set of $3\times P$ partial differential equations. These governing equations are numerically solved using the central difference scheme in the spatial dimension and the fourth order Runge-Kutta method in the temporal dimension. A uniform grid with $\mathrm{\Delta }x=0.01$ is used in the spatial dimension, while the time step $\mathrm{\Delta }t=0.0001$ is used for time integration. The inner products involved in the governing equations are evaluated as a-priory using the Gauss-Hermite quadrature.

For demonstration purpose, the proposed method is applied using ‘hypothetical test bed nozzle’ data. The ‘hypothetical nozzle’ can be created using the nozzle simulator with sensors that can measure the nozzle responses of interest with typical accuracies. A hypothetical nozzle is specified through a set of parameters and ‘true’ area profile, all of which are treated as precisely known. Total 5 sensors are assumed to be equally placed along the x-direction of the nozzle. In all the test cases presented in this section, density $\mathit{\rho }$, velocity $v$, pressure $P$ and static temperature $T$ are used as responses of interest, which are assumed to be measured by the sensors. Steady state predictions using known nozzle area profile are used as experimental observations. The experimental uncertainty is specified by a zero-mean normally distributed random variable with standard deviation given by 1% of the mean value. For all the test cases, inflow conditions are Mach number $M_{\mathit{in}}=1.5$, static pressure $P_{\mathit{in}}=1.0$ and density $\mathit{\rho }=1.0$. Figure 2 shows a typical nozzle area profile and spatial variation of normalized responses at steady state. Note that supersonic flow through a divergent nozzle results in increased velocity, while the static pressure, temperature and density are decreased. As can be observed from Figure 2, this behaviour is captured well by the nozzle simulator.

Fig. 2 Figure shows (a) nozzle area variation and (b) steady state response predictions for deterministic case. All the steady state responses are normalized using inflow values.

 

Prior uncertainty propagation

 

Prior uncertainty in the nozzle area profile is specified using a Gaussian process with known mean and the covariance function is given by(35) $$\begin{array}{c}\hfill C(x_1,x_2)={\mathit{\sigma }}^2\exp (-\mathit{\lambda }{(x_1-x_2)}^2),\end{array}$$(35) where ${\mathit{\sigma }}^2$ is the variance and $\mathit{\lambda }$ is the correlation length. ${\mathit{\sigma }}^2$ and $\mathit{\lambda }$ are treated as uncertain. Prior uncertainty in ${\mathit{\sigma }}^2$ is specified using an inverse Gamma distribution, $\mathit{IG}(\mathit{\alpha }=9.0,\mathit{\beta }=0.5)$, while Gamma distribution, $G(\mathit{\alpha }=5.0,\mathit{\beta }=0.2)$, is used as a prior for $\mathit{\lambda }$. Here $\mathit{\alpha }$ is shape parameter and $\mathit{\beta }$ is rate parameter for Gamma distribution, while $\mathit{\alpha }$ is shape parameter and $\mathit{\beta }$ is scale parameter for inverse Gamma distribution.[48] Figure 3(a) shows prior mean nozzle area profile while Figure 3(b) shows prior probability distribution. Comparison of prior mean with ‘true’ nozzle area profile is also shown in Figure 3(a).

The Bayesian framework is implemented using the intrusive gPC expansion for propagation of the prior uncertainty to the system response. The prior uncertainty is projected on a Hermite Polynomial Chaos basis. To investigate the trade-off between computational efficiency and accuracy, results of Polynomial Chaos method for uncertainty propagation are compared with the Monte Carlo simulation. Total 100,000 samples are used for the Monte Carlo method. All the computations are performed on a desktop computer with intel i5-460 processor. CPU time is estimated using a FORTRAN intrinsic cpu_time routine. Figure 4 shows the comparison of computational time requirement as a function of number of eigenfunctions used, $N$, for different Polynomial Chaos order, $p$. As can be observed from the figure, the computational cost increases polynomially with $N$. In particular, for a system with stochasticity of order $s$, the computational cost increases as $N{s}^{\mathit{pq}}$, where $q$ is the order of non-linearity of the system and $p$ is order of the gPC basis.[34]

Figure 5 shows the $L_1$-error in the mean and the variance. The $L_1$-error is calculated with respect to the Monte Carlo method. With the increase in number of eigenfunctions used and the polynomial order, error in the mean and the variance reduces. For $N=4$ and the $2\text{nd}$ order Polynomial Chaos, method provides prediction of system response with error of the order of ${10}^{-3}$ in both mean and variance at $10$-times lower computational cost. Note that the computational cost of the proposed Bayesian calibration method is dominated by the computational time requirement for solution of forward propagation problem, while the computation cost of the MCMC sampling is negligible. Thus, the conclusions drawn from the computational cost comparison for forward propagation can be extended to solution of the inverse problem without any change.

Fig. 3 Figure shows (a) comparison of prior mean nozzle area profile with true nozzle area profile and (b) prior probability distribution of nozzle area profile.

Fig. 4 Figure shows cpu time required as a function of number of eigenfunctions used. Effect of the Polynomial Chaos order is also shown.

Fig. 5 Figure shows (a) $L_1$-error in mean and (b) $L_1$-error in variance. Effect of the Polynomial Chaos order is also shown.

 

Spectral projection-based Bayesian calibration

Baseline case

The computational cost of the proposed method is compared with the direct MCMC sampling for the Bayesian calibration. Prior for the uncertain nozzle area profile is specified using the Gaussian process with the uncertain variance and the covariance length. $\mathit{IG}(9.0,0.5)$ prior is used for the variance, while $G(5.0,0.2)$ prior is used for the correlation length. First four eigenfunctions are used in the spectral expansion, while second-order Hermite polynomials are used as the gPC basis. Model structural uncertainty is quantified by $\mathit{IG}(9.0,0.5)$ prior for the variance and $G(6.0,2.0)$ prior for the correlation length of the covariance function of the discrepancy function. The gPC expansion coefficients of the simulator output, obtained using the intrusive Galerkin projection, are used in the likelihood function to define the posterior distribution, which is explored using the MCMC sampling. Gaussian distribution centred on the present state is used as a proposal distribution for the Markov Chain. Total 100,000 samples are collected after rejecting the initial 10,000 samples. Results of the proposed method are compared with the direct implementation of MCMC for the Bayesian calibration, where the simulator output is used to define the likelihood. Total 10,000 samples are collected using the direct MCMC after the initial burnout period of 1000 samples. The total CPU time required for the implementation of the Bayesian framework using the direct MCMC sampling is 12732.65 seconds. The CPU time requirement for the proposed gPC expansion based Bayesian framework is split into the time required for the intrusive Galerkin projection of the prior uncertainty to the system response and the MCMC sampling from the posterior distribution. For the present test case, intrusive Galerkin projection is implemented in 1194.15 seconds, while the total CPU time for the MCMC sampling is 8.48 seconds, taking 1202.63 seconds for complete implementation of the proposed stochastic spectral projection based Bayesian framework.

Figure 6 shows comparison of the posterior mean nozzle area profile obtained using the direct MCMC and the proposed method. The prior mean nozzle area profile is also shown in the figure. The posterior mean nozzle area profile obtained using both the methods matches closely with the ‘true’ nozzle area profile. The comparison of the posterior distribution for the hyper-parameters of the uncertain nozzle area profile is shown in Figure 7. The probability distribution of the hyper-parameters is not updated noticeably after the Bayesian calibration, as the experimental observations of the system response are not expected to contain significant information about the covariance structure of the uncertain nozzle area profile. From the figures, it may be concluded that the proposed spectral projection based Bayesian framework provide inference of the uncertain functions with accuracy comparable to the implementation of the direct MCMC sampling at significantly lower computational cost.

Fig. 6 Comparison of the posterior mean nozzle area profile with the true nozzle area profile and prior mean. The comparison is shown for the posterior mean nozzle area profile obtained using implementation of the direct MCMC sampling and the proposed method. Figure also shows 90% confidence bound.

Fig. 7 Figure shows comparison of posterior distribution for the hyper-parameters of the uncertain nozzle area profile obtained using the direct MCMC and the proposed spectral projection based method. Figure (a) shows the comparison for the variance and the Figure (b) shows the comparison for the correlation length.

 

Choice of discrepancy function

To investigate the effect of choice of the prior for discrepancy function, Bayesian framework is implemented using different priors for variance. Prior uncertain in the nozzle area profile is specified as discussed earlier. Prior uncertainty in discrepancy function is specified using zero-mean Gaussian process with squared exponential covariance function (35(35) $$\begin{array}{c}\hfill C(x_1,x_2)={\mathit{\sigma }}^2\exp (-\mathit{\lambda }{(x_1-x_2)}^2),\end{array}$$(35) ), where the variance and the correlation length are uncertain hyper-parameters. For all the test cases presented in this section, prior in correlation length is represented using $G(6.0,2.0)$, specifying correlated discrepancy function. For variance, $\mathit{IG}(6.0,2.0)$ and $\mathit{IG}(1.5,2.0)$ priors are investigated. To simulate model structure uncertainty, hypothetical test-bed data are obtained using viscous model, whereas, simulator is defined using inviscid model. Total 5 testbed data points are used with 1% standard deviation. Figure 8 shows comparison of posterior mean nozzle area profile with true nozzle area profile and the prior mean. Posterior mean for $\mathit{IG}(6.0,2.0)$ prior matches closely with the true nozzle area profile, whereas posterior mean nozzle area profile for $\mathit{IG}(1.5,2.0)$ prior deviates from the true nozzle area profile. Figure 8 demonstrates the significant impact of the prior model structural uncertainty on the inference of the uncertain functions. Confidence on the available model is specified through the prior on the variance. In the case of high confidence on the simulator model, specified through the high value of $\mathit{\alpha }$ for $\mathit{IG}(\mathit{\alpha },\mathit{\beta })$ prior, significant amount of the information provided by the data is used to update the parameter uncertainty. If the low confidence prior is specified for the model through the lower value of $\mathit{\alpha }$ for $\mathit{IG}(\mathit{\alpha },\mathit{\beta })$ prior, the calibration process uses significant information provided by the data to improve confidence on the simulator model, whereas less information is used for updating the uncertain parameters. Note that in the case of very high confident prior on the model structure (approaching the scenario of no model structural uncertainty), calibration process attributes any remnant error in the model to the parameters. Thus, the present authors propose to avoid priors that specify either low confidence or very high confidence on the model structure. In the remaining test cases presented in this paper, $\mathit{IG}(6.0,2.0)$ prior is used for the variance of the discrepancy function.

Fig. 8 Figure shows effect of the choice of prior for discrepancy function on posterior mean of the nozzle area profile. Note that $\mathit{IG}(1.5,2.0)$ represent lower confidence on the simulator model as compared to the $\mathit{IG}(6.0,2.0)$ prior.

 

Effect of simulator model error

To investigate efficacy of the calibration process in presence of error in simulator model, artificial discrepancy is introduced in the simulator model by multiplying the gPC expansion coefficients of the static pressure by 1.5. Modified chaos coefficients are used in the Bayesian framework. The framework is implemented using following two test cases: (1) simulator model without taking into account discrepancy function (which signifies full confidence on the simulator model), and (2) $G(6.0,2.0)$ and $\mathit{IG}(6.0,2.0)$ priors for correlation length and variance of the discrepancy function (${\mathit{\lambda }}_{\mathit{\delta }}$ and ${\mathit{\sigma }}_{\mathit{\delta }}^2$, respectively). Figure 9 shows comparison of posterior nozzle area profile with the true nozzle area profile. For calibration without considering discrepancy function, Bayesian framework assumes simulator to be the ‘true’ representation of the physics. The Bayesian calibration method attributes all the difference between the test-bed data and the simulator prediction to the uncertain parameters. Thus, resultant posterior mean deviates significantly from the true nozzle area profile. However, when discrepancy function is considered in the Bayesian calibration, no significant update is observed in the nozzle area profile.

Figure 10 shows comparison of the prior and the posterior probability distribution of ${\mathit{\sigma }}_{\mathit{\delta }}^2$ and ${\mathit{\lambda }}_{\mathit{\delta }}$. The posterior probability distribution of $1/{\mathit{\sigma }}_{\mathit{\delta }}^2$ for the static pressure is shifted towards left, indicating the lower posterior expected value of $1/{\mathit{\sigma }}_{\mathit{\delta }}^2$ as compared to the prior. However, no significant change is observed for the probability distribution of ${\mathit{\lambda }}_{\mathit{\delta }}$. Note that the prior $G(6.0,2.0)$ indicated correlated discrepancy in the simulator model. Thus, the posterior probability distribution of the discrepancy function for the static pressure indicates correlated discrepancy in the simulator model with the high expected value of ${\mathit{\sigma }}_{\mathit{\delta }}^2$. As the error in the experimental observations enters the Bayesian framework as the uncorrelated uncertainty, the highly correlated discrepancy is expected to result from the error in the simulator model. Thus, the posterior distribution of the discrepancy function signals the need for verification and validation of the simulator, particularly subroutines affecting the static pressure.

Fig. 9 Figure shows the effect of error in the simulator model on the posterior nozzle area profile. Results are obtained for artificially introduced discrepancy in the simulator model. The comparison is shown for the Bayesian calibration without using the discrepancy function (no discr) and with using the discrepancy function (w/ discr).

Fig. 10 Figure shows comparison of prior and posterior distribution for (a) ${\mathit{\sigma }}_{\mathit{\delta }}^2$ and (b) ${\mathit{\lambda }}_{\mathit{\delta }}$ in presence of error in the simulator model. Results are obtained for artificially introduced discrepancy in the simulator model.

 

Effect of erroneous observations

To investigate effect of the erroneous experimental observations on the Bayesian framework, proposed method is implemented with artificially introduced discrepancy in test-bed data. The methodology is demonstrated by multiplying 1.5 to static pressure data. Figure 11 shows resultant posterior mean nozzle area profile. As observed in case of artificially introduced discrepancy in the simulator, significant deviation is observed between posterior mean and true nozzle area profile when discrepancy function is not considered in the formulation. However, when discrepancy function is considered, no significant change is observed in the nozzle area profile.

Fig. 11 Figure shows the effect of erroneous experimental observations on the posterior nozzle area profile. Results are obtained for artificially introduced errors in the experimental observations.

Figure 12 shows posterior probability distribution for ${\mathit{\sigma }}_{\mathit{\delta }}^2$ and ${\mathit{\lambda }}_{\mathit{\delta }}$. Posterior distribution of $1/{\mathit{\sigma }}_{\mathit{\delta }}^2$ for static pressure has moved towards left, indicating the high posterior discrepancy. Posterior probability distribution of ${\mathit{\lambda }}_{\mathit{\delta }}$ for static pressure is shifted towards right, indicating highly uncorrelated discrepancy. Since the uncertainty due to unknown or poorly known physics is expected to result in the correlated discrepancy, uncorrelated discrepancy may be attributed to the experimental observations. Thus, the posterior distribution indicates need for the review of the experimental observations.

Summary

The numerical test cases presented in this paper have demonstrated ability of the proposed Bayesian framework to infer uncertain spatial/temporal functions in the presence of model structural uncertainty. In addition to the inference of the uncertain functions, the proposed Bayesian framework provides useful insight into the credibility of the simulator model. The posterior probability distributions of ${\mathit{\sigma }}_{\mathit{\delta }}^2$ and ${\mathit{\lambda }}_{\mathit{\delta }}$, through hyper-parameters $\mathit{\alpha }$ and $\mathit{\beta }$, are identified as indicators of the veracity and validity of the simulator model. Based on the posterior values of the hyper-parameters $\mathit{\alpha }$ and $\mathit{\beta }$ for the posterior probability distributions of ${\mathit{\sigma }}_{\mathit{\delta }}^2$ (${\mathit{\alpha }}_{{\mathit{\sigma }}_{\mathit{\delta }}}$ and ${\mathit{\beta }}_{{\mathit{\sigma }}_{\mathit{\delta }}}$) and ${\mathit{\lambda }}_{\mathit{\delta }}$ (${\mathit{\alpha }}_{{\mathit{\lambda }}_{\mathit{\delta }}}$ and ${\mathit{\beta }}_{{\mathit{\lambda }}_{\mathit{\delta }}}$), authors provide following guidelines:

1. If the posterior ${\mathit{\alpha }}_{{\mathit{\sigma }}_{\mathit{\delta }}}$ is greater than the prior ${\mathit{\alpha }}_{{\mathit{\sigma }}_{\mathit{\delta }}}$, then the simulator model should be accepted with improved confidence.

2. ${\mathit{\alpha }}_{{\mathit{\sigma }}_{\mathit{\delta }}}<1$ : The posterior distribution does not have mode resulting in maximum probability for $\frac{1}{{\mathit{\sigma }}_{\mathit{\delta }}^2}\to 0$. For such posterior, calibrated simulator model should be rejected.

   1. If mean and mode of posterior distribution of ${\mathit{\lambda }}_{\mathit{\delta }}$ indicate strong correlation for discrepancy function (typically ${\mathit{\alpha }}_{\mathit{\lambda }}<1$ or ${\mathit{\beta }}_{\mathit{\lambda }}>{\mathit{\alpha }}_{\mathit{\lambda }}$), rigorous verification and validation process is advised with a focus on subsystem models that predict system responses for which ${\mathit{\alpha }}_{{\mathit{\sigma }}_{\mathit{\delta }}}<1$.

   2. If mean and mode of posterior distribution of ${\mathit{\lambda }}_{\mathit{\delta }}$ indicate very weak correlation for discrepancy function (typically ${\mathit{\alpha }}_{\mathit{\lambda }}>1$ or ${\mathit{\beta }}_{\mathit{\lambda }}\ll \mathit{\alpha }$), review of experimental observations is advised.

   3. For all the other cases, rigorous verification and validation of simulator model is advised with a note on review of experimental observations.

3. ${\mathit{\alpha }}_{{\mathit{\sigma }}_{\mathit{\delta }}}>1$ and ${\mathit{\beta }}_{{\mathit{\sigma }}_{\mathit{\delta }}}>{\mathit{\alpha }}_{{\mathit{\sigma }}_{\mathit{\delta }}}$: Calibrated simulator can be used, however, high uncertainty in prediction of the simulator response should be expected. As per discussion void point (2), verification and validation of simulator model and a review of experimental observations is advised.

4. ${\mathit{\alpha }}_{{\mathit{\sigma }}_{\mathit{\delta }}}\gg 1$ and ${\mathit{\beta }}_{{\mathit{\sigma }}_{\mathit{\delta }}}<{\mathit{\alpha }}_{{\mathit{\sigma }}_{\mathit{\delta }}}$: Calibrated simulator can be used with high confidence.

Fig. 12 Figure shows comparison of prior and posterior distribution for (a) ${\mathit{\sigma }}_{\mathit{\delta }}^2$ and (b) ${\mathit{\lambda }}_{\mathit{\delta }}$ in presence of erroneous experimental observations.

 

Concluding remarks

This paper has proposed a gPC based Bayesian framework for calibration of computer simulators. The proposed framework has extended the established method to the priors with uncertain hyper-parameters. Efficacy of the proposed Bayesian framework is demonstrated for calibration of a quasi-one-dimensional divergent nozzle flow simulator. Ability of the method to infer spatially/temporally varying uncertain functions is shown using the update of the nozzle area profile. The proposed method has provided accurate inference of nozzle area profile at one-tenth of a computational cost as compared to the direct implementation of the Bayesian framework. Hyper-parameters of the posterior distribution of model structure uncertainty are identified that provide information about veracity and validity of the computer simulator. Based on the hyper-parameters, guidelines have been provided for acceptability of the simulator model. Although demonstrated for a specific set of priors, the proposed method is generic in nature and can admit arbitrary priors. However, depending on the gPC basis used, higher order polynomials may be required for the satisfactory spectral approximation of the prior, incurring comparatively higher computational cost on the Bayesian framework.

Acknowledgments

This work was supported in part by Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education, Science and Technology (2010-0025484), and in part by the KI Project through KAIST Institute for Design of Complex Systems.

Notes

Note that the terminology ‘true’ parameter is not universally accepted in the open literature, with one point of view considering simulator as a black box and uncertain parameters as ‘tuning-parameters’,[10] while other authors refer to uncertain parameters as ‘true’ parameters which are being inferred.[12]

For notational convenience, ${\mathbf{x}}_s$ is replaced by $\mathbf{x}$ in this section.

Note that though the proposed method is demonstrated for $\mathit{IG}/G$ priors, the method can be applied for arbitrary choices of prior.