Inverse Uncertainty Quantification by Hierarchical Bayesian Modeling and Application in Nuclear System Thermal-Hydraulics Codes

Figure 1 shows the major components used in the IUQ framework, where x represents the control parameters, such as boundary conditions, initial conditions, etc. $\bm{\theta}$ represents calibration parameters, in this work it specifically represents the selected physical model parameters in the closure models of TH codes.

The input deck serves as the input for thermal-hydraulics (TH) codes, defining a specific system. It includes the geometrical configuration (i.e., nodalization), the materials and fluids involved, the initial and boundary conditions, and potentially the settings for the numerical solver. Control parameters x are selected from those specifications. The TH simulation codes TRACE consists of 6 conservation equations, which are then closed with additional set of closure models $M_{i}(\bm{x,\theta,y^{M}})$ (Wicaksono (2018)). The output of the simulation code $\bm{y}^{M}$ can be combined with corresponding experimental data $\bm{y}^{E}$ in the Bayes’ rule, resulting in a joint posterior probability density distribution (PDF) for the selected inputs.

The detailed procedures of the IUQ framework can be organized into the following stages:

1. 1.

   Define the problem under investigation and select appropriate experimental data and simulation codes. In this study, the BFBT (Boiling Water Reactor Full-size Fine-mesh Bundle Test) benchmark is used as the basis, with corresponding TRACE models developed.

2. 2.

   Determine influential input parameters within the chosen TH codes.

3. 3.

   Conduct a two-step sensitivity analysis to identify significant input parameters. Initially, implement a relatively straightforward perturbation method to all parameters, selecting those that are active in the model. Subsequently, employ a more precise sensitivity analysis technique, Sobol indices, to identify the influential parameters. Surrogate models may be used in this stage.

4. 4.

   Develop a hierarchical Bayesian inference model to quantify the posteriors of uncertain inputs. This stage integrates experimental data with simulation outcomes to establish input distributions that result in improvement agreement with experimental data. Surrogate model and MCMC sampling algorithms are used in this stage.

5. 5.

   Validate the derived posteriors.

A key assumption in the Bayesian IUQ framework is the model updating equation. Following the work of Kennedy and O’Hagan (2001), we represent the relationship between the computer model outputs $\bm{y}^{M}(\bm{x},\bm{\theta})$ and the observations $\bm{y}^{E}(\bm{x})$ in the equation:

where $\bm{\epsilon}$ is the observation error that is usually assumed to be independent and identically distributed as $\mathcal{N}(0,\sigma_{exp}^{2})$. It should be noted that this assumption may not always hold in reality, especially in time-dependent problems (Wang et al. (2018b)). $\bm{\delta}(\bm{x})$ is the model discrepancy term, which is caused by incomplete or inaccurate physics employed in the model. Here, the parametric uncertainty is derived from the $\bm{\theta}$ parameter, and other forms of uncertainties are incorporated in the model discrepancy term.

Following the model updating equation, the posterior PDF of the calibration parameters can be found using Bayes’ rule:

where $p(\bm{\theta})$ is the prior distribution of the calibration parameters, and $p(\bm{y}^{E},\bm{y}^{M}|\bm{\theta})$ is the likelihood function. From equation 1, we know that $\bm{\epsilon}=\bm{y}^{E}(\bm{x})-\bm{y}^{M}(\bm{x},\bm{\theta})-\bm{\delta}(\bm{x})$ follows a multivariate normal distribution. So the posterior can be written as:

The covariance matrix $\bm{\Sigma}_{t}$ is defined as:

where $\bm{\Sigma}_{exp}$ is the experimental uncertainty, $\bm{\Sigma}_{\delta}$ is the model uncertainty, and $\bm{\Sigma}_{code}$ is the code uncertainty, introduced by using surrogate models to replace the full model. $p(\bm{\theta})$ is the prior distribution of calibration parameters and can be treated as a non-informative uniform distribution over a certain range.

It can be very challenging to estimate the model discrepancy $\bm{\delta}$. In cases where model discrepancy is negligible, we can simply ignore the $\bm{\delta}$ term. However, ignoring the model discrepancy when it does exist can cause over-fitting issue of the calibration parameters, meaning that the IUQ process will favor a $\bm{\theta}$ distribution that yields the best agreement between the measurement data and the model simulation, rather than the “true” value. A method called the Modular Bayesian Approach (MBA) was developed to tackle the model discrepancy problem in IUQ (Wu et al. (2018b) Wu et al. (2018c) Wang et al. (2018a)).

One question we need to answer at the very beginning of IUQ is the identification of input parameters requiring study. Sensitivity Analysis (SA) examines how uncertain inputs contribute to variations in model QoIs. There are typically two classes of SA: local SA and global SA (Saltelli (2002)). Local SA addresses the local impact of input variations, whereas global SA considers the entire variation range of inputs and can also provide insights into how interactions between input parameters influence QoIs. Generally, local SA is computationally less expensive than global SA; thus, in the IUQ framework, we use local SA for initial parameter selection and global SA for a more accurate study to determine the final list of inputs to be investigated.

Sobol indices method provides a straightforward measure of sensitivity in arbitrarily complex computational models. It is a variance-based method where the variance of the output is decomposed as a sum of contributions of each input variable and their combinations. This decomposition of variance is referred to as ANOVA (ANalysis Of VAriance). Sobol indices method has been successfully applied to many UQ related applications (Aly et al. (2019) Wang et al. (2017) Wu et al. (2017b)) and the details of the method will not be repeated in this paper. The calculation of Sobol indices in this paper will follow the samping-based approach outlined by Saltelli et al. (2010). This approach requires $N(d+2)$ number of samples. $N$ is usually a value of several thousand for sufficient accuracy and $d$ is the dimension of the input, and GP surrogate model is used here due to a large number of code runs required.

Surrogate models, also known as metamodels or emulators, play a critical role in Bayesian calibration of computer models, particularly when dealing with computationally demanding simulations such as nuclear reactor analyses. These surrogate models are constructed using a limited number of full model runs at specifically chosen input parameter values, combined with a learning algorithm. Many types of surrogate models have been utilized in the Bayesian calibration process, Different types of surrogate models have different characteristics, so they may fit in different engineering scenarios. Factors to consider when selecting an appropriate surrogate model include input dimensionality, the non-linearity of the input-output relationship, and whether the problem is time-dependent or steady-state.

In this work, we primarily employ two types of surrogate models: GP and PR. GP has been widely used as surrogate models since the seminal work of Bayesian calibration by Kennedy and O’Hagan (2001) utilizing GP modeling. One notably characteristic of GP is that it not only provides a mean estimator at any untried input location, but also the corresponding variance estimator, which can serve as the code uncertainty term ($\bm{\Sigma}_{code}$) in Equation 4. In other surrogate models, since this term is not available, we will need to treat the covariance term as a calibration parameter. Details of the GP surrogate model will not be repeated in this paper, interested readers may refer to Rasmussen (2004) for more details on GP and Wu et al. (2018b) Wang et al. (2019a) for similar IUQ applications using GP.

In addition to the GP model, as a more convenient but less accurate alternative, Polynomial Regression model is also used in the current work. The incentive of using PR derives from the requirement of gradient information in some advanced and efficient MCMC algorithms. As a parametric regression model, it is straightforward to calculate the gradient at any given location. The PR model with degree of $d$ for two variables $x_{1},x_{2}$ has the following form:

where all d^th order polynomials and interaction terms are considered. The parameter $\bm{w}=[w_{1},w_{2},...]$ can be easily estimated by Ordinary Least Squares method, similar to other linear regression models.

The accuracy of the surrogate models needs to be assessed before use. We are particularly interested in the predictive accuracy at untried points, which can be achieved by quantifying the predictive error using an additional set of validation data. Various metrics, such as the Mean Squared Error, Leave-One-Out-Error, Coefficient of Determination, and the Mean Absolute Error, can be employed to measure the discrepancy between the predicted values and the actual values. By comparing these metrics, we can ascertain the effectiveness of the surrogate model in capturing the underlying patterns and relationships in the data, ensuring its reliability for subsequent analysis and applications.

In the traditional Bayesian calibration setting, we have assumed that the observations are independent of each other so that the joint likelihood function can easily be formulated as the product of each individual likelihood. However, in many situations, such independence may not hold. Hierarchical model is also called a multi-level model or mixed-effect model. A simple illustration of the structure of the hierarchical model is shown in Figure 2. Observations $\bm{y}$ are from different clusters determined by (1) cluster parameter $c_{m}$ which might be different among clusters, and (2) parameter $\bm{\theta}$ shared across clusters. If the probability distribution of $c_{m}$ can be parameterized by $\bm{\Sigma}_{c}$, $\bm{\Sigma}_{c}$ can be seen as the cluster-specific parameter, or per-group parameter, to be estimated.

Consider the simple example illustrated in Figure 2, where the observations are assumed to follow a normal distribution. In this scenario, the parameters to be estimated are the mean and variance. Let’s assume that all clusters share the same variance, denoted as $\sigma^{2}$, but have different means $\mu_{i}$ for each cluster. So the $\sigma^{2}$ is a shared parameter ($\theta$ in Figure 2) and $\mu_{i}$ is per-group parameter ($c_{i}$ in Figure 2). Now we would like to specify the distribution over the cluster-specific $\mu_{i}$. We can assume the distributions are also normal, or other task-specification distributions, without loss of generality. If normal distributions are assumed, two parameters (global mean $\mu$ and global standard deviation $\sigma_{y}$, correspond to the $\Sigma_{c}$ in Figure 2) would be required. So we can describe the parameters $\mu_{i}$ and observations $y_{ij}$ as:

Now we can re-arrange the structure of the above hierarchical model to better illustrate the process of the two preceding equations: the per-group parameters are generated from shared parameters, and the observations are generated from per-group parameters. This relationship is shown in Figure 3.

In summary, the application of hierarchical models is highly beneficial when there’s an observable existence of group-level characteristics, and the data’s structure can be effectively modeled or assumed. When dealing with hierarchically structured data, utilizing a non-hierarchical model may not be appropriate for two reasons: (1) a model containing insufficient parameters may not successfully fit the dataset, and (2) a model with an abundance of parameters may cause over-fitting of the dataset. In contrast, hierarchical models, having a substantial number of parameters for fitting the data accurately and a population distribution to model the inter-dependencies of the parameters, can avoid the issue of over-fitting. This is achieved as the parameters are “constrained” by the population distribution, preventing them from exactly replicating the data.

In this section we describe how to estimate parameters in a hierarchical Bayesian model structure. The hierarchical model structure in Figure 2 can be plotted into a more succinct and formal graphical representation, as shown in Figure 4. In the figure, $p(\bm{\theta})$ is the prior distribution of $\bm{\theta}$, $i$ is the cluster index, $\bm{y}$ is the observation, and there are $N$ observations in total. The per-group parameter $b_{i}$ is governed by its distribution $\bm{\Sigma_{b}}$, which has a prior distribution $p(\bm{\Sigma_{b}})$.

It should be noted that within this model, our primary goal is to learn about $\bm{\theta}$ and $\bm{\Sigma_{b}}$, as opposed to the individual group-specific parameters $b_{i}$, but $b_{i}$ also has to be estimated in order to estimate $\bm{\Sigma_{b}}$. So we need to marginalize over the $b_{i}$ parameters after we have obtained the joint posterior distributions. There are generally two approaches to estimate the posteriors, we can either use Maximum Likelihood Estimation (MLE) to obtain best-fit point estimations or use sampling to get full posterior distributions for these parameters.

In the MLE method, the likelihood of $\theta$ and $\Sigma_{b}$, marginalized over $b$, can be expressed as:

In case full posterior distributions are desired, we can integrate the prior distribution and apply Bayes’ rule. The prior distribution can be decomposed into the following form:

which is based on the fact that $\bm{\Sigma_{b}}$ influences observations only through $b_{i}$. So the marginalized joint posterior distribution can be expressed as:

Hierarchical models provide a flexible framework for modeling the complex interactions but come with higher computational requirements in inference. The parameter space to be estimated can be high-dimensional: we need to estimate cluster-specific parameter $b_{i}$, which may contain up to hundreds. Traditional random-walk based MCMC algorithms such as Metropolis-Hastings (MH) or advanced MH algorithms tailored for relatively high dimensional problems still suffer from serious convergence issues because that the random behavior of the proposal function in very inefficient in high-dimension domains. Hamiltonian Monte Carlo (HMC) and No-U-Turn Sampler (NUTS), which will be introduced in the following section, would be useful in such high-dimensional conditions.

Estimating the posterior in Equation 3 requires numerical sampling method as analytical solution is not available. Statistical tools that can be used here include MLE, Maximum A Posteriori (MAP), Variational Inference (VI), Laplace Approximation (LA), and various MCMC methods. MLE, MAP, LA all give point estimates and VI, LA are both approximate solutions rather than exact solutions. MCMC, although more time consuming, can estimate the exact posteriors for arbitrary distributions. It is a class of algorithms for sampling from an arbitrary probability distribution, for example, the probability distribution of equation 3, where the density function needs to be known only up to a normalizing constant.

The MH algorithm is one of the most widely used MCMC methods for sampling from a target distribution. In this approach, a family of potential transitions from one state to the next is determined using a proposal distribution. However, when the number of parameters to be estimated increases to 20 or 30, the random walk-based methods, such as the traditional MH algorithm, suffer from a very low acceptance rate. To address these limitations, the HMC method has been developed as an alternative to the conventional MH algorithm. HMC is an MCMC technique that avoids the random walk behavior and the sensitivity to correlated parameters by utilizing first-order gradient information to inform the sequence of steps taken during the sampling process (Andrieu and Thoms (2008) Hoffman and Gelman (2014)). By employing an auxiliary variable scheme, HMC can transform the problem of sampling from a target distribution into simulating Hamiltonian dynamics, effectively suppressing the random walk behavior. As a result, HMC offers a more efficient approach to sampling from high-dimensional distributions and provides a smoother transition between states, thus enhancing the performance of MCMC methods in scenarios with a large number of parameters to be estimated.

In HMC, an auxiliary momentum variable $r_{d}$ is introduced for each model model parameter $\theta_{d}$. In most cases, these momentum variables are drawn independently from standard normal distribution, yielding the joint density $p(\theta,r)\propto\exp(L(\theta)-\frac{1}{2}r\cdot r)$, where $L(\cdot)$ is the logarithm of the joint density of the variables of interest $\theta$. We can interpret this model in physical terms as a fictitious Hamiltonian system where $\theta$ is a particle’s position in D-dimensional space, $r_{d}$ is the momentum of that particle in the d-th dimension, $L$ is a position-dependent negative potential energy function, $\frac{1}{2}r\cdot r$ is the kinetic energy of the particle, and $log(p(\theta,r))$ is the negative energy of the particle. The evolution of the Hamiltonian dynamics of this system can be simulated via the “leapfrog” integrator (Hoffman and Gelman (2014)), which will update the $r$ and $\theta$ accordingly. The details of the algorithm is shown in Algorithm 1.

Suppose we are drawing $M$ samples using the HMC algorithm. For each sample, we first re-sample the momentum variables from $N(0,I)$, which can be interpreted as a Gibbs sampling update. Then, the leapfrog function is used $l$ times to update the variables $\theta$ and $r$, and a proposal position-momentum pair of $\tilde{\theta}$ and $\tilde{r}$ is generated. The proposal can be accepted or rejected according to the Metropolis algorithm. The leapfrog functions plays an important role in improving the efficiency of sampling. By taking $l$ leapfrog steps per sample, the proposal generated for $\theta$ has a high probability of being accepted. However, as we can see from the definition of leapfrog function, derivative of the joint density function of $\theta$ is required. The performance of the HMC also depends on hand-chosen values of $\epsilon$ and $l$. If the two parameters are not properly selected, the Markov chain may be slow to move between regions of high and low densities.

NUTS is an extension of HMC that eliminates the need to manually tune these parameters in HMC. The idea is to use a recursive algorithm to build a set of likely candidate points that spans the target distribution, and stop automatically when it starts to retrace its steps. A detailed description about NUTS can be found in Hoffman and Gelman (2014). In this work, the NUTS is implemented with an open-source package PyMC3.8 (Salvatier et al. (2016)). A significant advantage of NUTS is that it can be used without any hand-tuning, making it suitable for a wide range of engineering applications. However, as can be seen in the algorithm, the gradient of the log-likelihood with respect to the sampling parameter is required. Therefore, the surrogate model should also be capable of providing the gradient of outputs if NUTS is to be employed.

The international OECD/NRC BWR full-size fine-mesh bundle test (BFBT) benchmark (Neykov et al. (2006)) was created to encourage advancement in sub-channel analysis of two-phase flow in rod bundles, which has great relevance to the nuclear reactor safety evaluation. The BFBT test program includes experiments performed under steady-state and transient conditions, with measurements for single and two-phase pressure losses, void fraction, and critical power. The benchmark has been widely used in the validation process of various computational approaches (Wang et al. (2019b), Borowiec et al. (2017) Wu et al. (2018c)).

The void fractions are measured at four elevations along the test section, three of them are measured by X-ray densitometer, and the upper one is measured by a X-ray CT scanner. The void fraction data measured from the lower to upper location are referred to as ‘Void Fraction 1’ to ‘Void Fraction 4’, respectively. BFBT benchmark includes various assembly types and in each assembly type, experiments were conducted at various boundary conditions and void fraction data were measured.

The BFBT benchmark consists of 5 assembly types in all. The Assembly-4 is chosen in the current work. It consists of 86 experimental cases under various boundary conditions (flow rate, inlet temperature, power, and inlet pressure). The TRACE model is constructed according to the experimental geometry and boundary conditions. A comparison between the simulated void fraction results and experimental measurements are shown in Figure 9. The experimental data on ‘Void Fraction 1’, ‘Void Fraction 2’, and ‘Void Fraction 3’ by densitometers are corrected according to Glück (2008). We can see that TRACE model results agree well with experimental data and no obvious model discrepancy exists.

TRACE provides the option to adjust 36 physical model parameters (PMPs) by a multiplicative factor, however, not all of those parameters will be active in the BFBT bundle assembly model because some parameters involve phenomena that do not occur in BFBT benchmark, e.g. reflood. A simple perturbation method is firstly used to perturb each parameter in the range of $(0,5)$ while fixing other parameters. 50 uniform samples in that range are used to test the effect of this parameter on the simulated void fraction data. The resulting output variance is calculated for each parameter. The results show that most of the variances are 0 or very close to 0. Finally, eight parameters with variances larger than $10^{-3}$ are selected and shown in Table 1.

The Sobol indices method is subsequently employed for an further screening, and the first and total Sobol indices are shown in Figure 10. The different colors represent the corresponding indices for a certain measurement location. We can see that four out of eight have more significant influences on the model outputs, so the four parameters P1008, P1012, P1022, and P1028 are selected in the sensitivity analysis step and will be treated as uncertain inputs.

In this section, we provide the specifics of surrogate models that were introduced earlier. The task is to create a surrogate model with 4 inputs and 4 outputs. In order to determine the number of samples required for an accurate surrogate model, a convergence study is needed. Below we show a convergence study to find the sufficient number of Latin Hypercube Sampling (LHS) samples required by surrogates. Two types of surrogates are considered, GP and PR with a degree of 2. An additional 50 samples are taken as a validation set. Figure 11 displays the Mean Absolute Error (MAE). The graph indicates that the out-of-sample errors plateau once the sample number surpasses 100. Therefore, in this study, we choose 100 LHS samples to build the surrogate models. For the results presented below, PR model is used because they can provide convenient gradient information.

In order to determine the range of the input parameter for sampling, an iterative sampling procedure is used in this work. In the posterior sampling phase, it’s possible that the posterior distribution of the calibration parameter $\bm{\theta}$ may exceed the specified prior range. This occurrence is particularly common in the context of certain PMPs, such as the interfacial drag coefficient, where deviations from their nominal value can be large. In these instances, the posterior distribution may be artificially truncated by the upper bound of the initial range.

To circumvent this issue, we have employed an iterative resampling procedure in Wang (2020). If the posterior distribution is truncated by the upper bound, this suggests that the true posterior distribution may exceed the current range. In such a scenario, it is necessary to extend the initial upper bound of this parameter to a higher value. Doing so ensures that the prior range can accommodate the final posterior distribution. This process involves performing new LHS of the original code, the construction of a surrogate model, and the application of MCMC sampling. As the calibration parameters selected in this paper are all multiplicative factors, the lower bound remains 0.

In TRACE, dimensionless multipliers are usually used because it provides a straightforward way to compare those coefficients with the nominal value of $1.0$, and there is no need to calculate their absolute values. Traditional calibration methods implicitly assume that $\bm{\theta}$ is a global variable for all experimental cases $\bm{y}^{E}$. This assumption greatly simplifies the problem being considered, however, it may be problematic in the context of closure models in TH codes. The closure models are derived from Separate Effects Tests and are strongly dependent on boundary conditions and flow regimes, thus the “true” parameter $\bm{\theta}$ may have different distributions (from run to run over the physical experiments), resulting in different posteriors (including the uncertainty information as well as the best-fitting value information of the parameter) sampled by MCMC algorithms.

Ignoring this fact and sampling a single global variable might lead to unrealistic posteriors. An accurate way to treat this problem is to partition the experimental data and simulation outputs according to each correlation equation in a certain range of boundary conditions. However, the work will be tedious and is not meaningful in BEPU applications, because simple input distributions are desired rather than many input distributions dependent on other variables. As a middle ground between the two treatments, we propose to use hierarchical models to account for both similarities and differences among calibration parameters over all considered physical experiments.

The hierarchical model used in this work allows for the calibration parameters $\bm{\theta}_{i}$ to vary for each experimental case $i$, reflecting the different boundary conditions. At the same time, we think that those parameters (multipliers) should not display substantial disparities and they can follow a shared distribution, for example, normal distribution. This treatment permits variability for input parameters, which is more realistic and can help reduce the over-fitting issue.

In the hierarchical model proposed for this application, we assume that each experimental case forms a ‘group‘ or ‘cluster‘ since they have different boundary conditions but they have the same assembly type and similar underlying physics. Thus, we can view the 86 experimental cases as forming 86 distinct clusters. Furthermore, each group includes four measurements, leading to the model structure shown in Figure 12.

In Figure 12, $\theta_{i},i=1,2,3,4$ represent our selected four PMPs, ‘P1008’, ‘P1012’, ’P1022’, ’P1028’, respectively. M is the number of groups and N is the total number of measurement points. In this case $N=4M$ because there are four measurements in each group. $\sigma$ is measurement error. For each calibration parameter $\theta_{i}$, since they can be different across groups, so we suppose they are from a common normal distribution:

This approach accounts for possible variations of $\theta_{i}$ across different experiments that could arise due to potential errors or discrepancies. It is worth noting that some experimental data may exhibit a significant discrepancy for reasons that remain unknown or unaccounted for. Such data points could disproportionately influence the likelihood function, thereby making the posterior distribution highly sensitive to these points. We can view these points as “outliers”. By employing a hierarchical framework, we can focus on the distribution of $\bm{\theta}$, making the model robust against outliers. Adding suitable prior information to $\bm{\mu}$ and $\bm{\sigma}$ allows us to efficiently perform Bayesian inference for these parameters using the NUTS algorithm.