- Email: [email protected]
Estimation of probability Density Functions for model input parameters using inverse uncertainty quantification with bias terms
Annals of Nuclear Energy 133 (2019) 1–8
Contents lists available at ScienceDirect
Annals of Nuclear Energy journal homepage: www.elsevier.com/locate/anucene
Estimation of probability Density Functions for model input parameters using inverse uncertainty quantification with bias terms Rabie A. Abu Saleem a,⇑, Tomasz Kozlowski b a b
Nuclear Engineering Department, Jordan University of Science and Technology, P.O. Box 3030, Irbid 22110, Jordan Department of Nuclear, Plasma and Radiological Engineering, University of Illinois at Urbana-Champaign, 216 Talbot Laboratory, 104 S. Wright St., Urbana, IL 61801, USA
a r t i c l e
i n f o
Article history: Received 14 November 2018 Received in revised form 2 May 2019 Accepted 4 May 2019
Keywords: Inverse uncertainty quantification Maximum Likelihood Estimation Maximum A Posterior estimation Bias terms BFBT
a b s t r a c t The documentation of most nuclear thermal-hydraulics codes does not provide sufficient information on uncertainty of physical models (e.g. interfacial heat transfer coefficients). These models were derived based on experimental data and implemented as empirical correlations in the computational code. The uncertainty quantification for the relevant output quantity (e.g. Peak Cladding Temperature) requires estimation of the Probability Density Functions (PDFs) of the code inputs, such as physical models. In this paper, we investigate the effect of boundary conditions (outlet pressure, inlet liquid temperature, and inlet flow rate) on the uncertainty of two physical models (the interfacial friction coefficient and the wall to liquid friction coefficient). The boundary conditions effect was accounted for by adding a bias term to the mathematical framework of two existing methods for Inverse Uncertainty Quantification (IUQ): the Maximum Likelihood Estimation (MLE) method and the Maximum A Posterior (MAP) method. The two methods were demonstrated using the BFBT benchmark, experimental data was compared to code predictions of the RSTART thermal-hydraulics code for two different cases: without and with bias term. The results show an evident improvement in code prediction when the bias term is used. Finally, a validation set of experimental data was used to investigate the possibility of data overfitting, and the proposed methodology showed absence of overfitting when bias terms are used. Ó 2019 Elsevier Ltd. All rights reserved.
1. Introduction Uncertainty Quantification (UQ) plays a crucial role in the verification and validation of nuclear codes. One of the main goals of UQ is to develop rigorous methods to properly characterize and quantify the variability in model input parameters and output results of a simulation code (Shrestha and Kozlowski, 2015). The use of Best Estimate Plus Uncertainty (BEPU) approach for thermal-hydraulics transient calculations has been promoted by the U.S. NRC (Nuclear Regulatory Committee) over the conservative approach for the purpose of understanding the performance of Emergency Core Cooling Systems (ECCSs) during a scenario of a reactor transient (D’Auria et al., 2006). In the conservative approach, the output response of a code is calculated using extreme values of input parameters. In other words, the values used in such codes are believed to produce the most unfavorable results. The reason behind this approach is to model the physical phenomena such that it always predicts the worst-case scenario. On the other hand, the BEPU approach is meant to predict the ⇑ Corresponding author. E-mail addresses: [email protected] (R.A. Abu Saleem), [email protected] (T. Kozlowski). https://doi.org/10.1016/j.anucene.2019.05.005 0306-4549/Ó 2019 Elsevier Ltd. All rights reserved.
results as realistically as possible by predicting the mean of experimental data and considering the effects of all the important variables. With this approach, some form of uncertainty is provided to account for omitting the variables that are not possible or practical to consider for a given phenomenon (Hu and Kozlowski, 2016). BEPU approach enables more precise calculation of safety margins and aids in regulating the licensing and operation of nuclear reactors with less cost compared to the conservative approach. There are two types of uncertainty quantification; forward Uncertainty Quantification (UQ) and Inverse Uncertainty Quantification (IUQ). Forward uncertainty quantification aims at quantifying the overall uncertainty in the results (output) of a code by propagating the uncertainties of the input model parameters of that code. Inverse uncertainty quantification, on the other hand, entails the process of determining the values and ranges of the input code parameters of a code by using experimental data along with code prediction data (Kennedy and O’Hagan, 2001). In other words, it determines Probability Density Functions (PDFs) of the code’s input model parameters. These PDFs can then be used as the start point for forward UQ analysis leading to more precise safety margins. For forward UQ problems, Monte Carlo methods, perturbation methods, and spectral expansion methods are some of the most commonly used (Lee and Chen, 2009; Robert and
2
R.A. Abu Saleem, T. Kozlowski / Annals of Nuclear Energy 133 (2019) 1–8
Casella, 2013; Mooney, 1997). Common methods for inverse problems include Maximum-Likelihood Estimation (MLE), Maximum A Posteriori (MAP), and Markov Chain Monte Carlo (MCMC) (Shrestha and Kozlowski, 2015; Gelman et al., 2014; McLachlan and Krishnan, 2007; Gilks, 2005). Many researchers have investigated using IUQ in verifying and validating results of different nuclear codes. In some studies, an Expectation-Maximization (EM) algorithm was applied within a mathematical framework based on MLE and MAP methods (Shrestha and Kozlowski, 2015; D’Auriaa and Mazzantini, 2009). These methods were constrained with some limitations based on the assumptions used to derive them, some of these constraints are: assuming a linear correlation between the code output and the input parameter and the need for local sensitivity analysis (Wu and Kozlowski, 2017). Some other researchers established a mathematical framework for estimating and quantifying the uncertainties in input parameters and code results concurrently, and they were able to verify their method for time-dependent nonlinear systems of large scale (Cacuci et al., 2010; Badea et al., 2012; Petruzzi et al., 2010; Cacuci and Arslan, 2014). Furthermore, this predictive method was extended to develop an innovative mathematical procedure for predictive modeling of coupled multi-physics systems (Cacuci, 2014), and was applied to different benchmarks to demonstrate its ability to produce best estimate responses and parameter values with reduced uncertainties (Cacuci and Badea, 2014; Latten and Cacuci, 2014). Other studies explored the implementation of inference methods with metamodels rather than system codes. One study investigated using Bayesian inference methods for the IUQ problem with surrogate models derived based on polynomial chaos expansion (Wu and Kozlowski, 2017). This method was applied for a point reactor kinetics equation coupled with a lumped parameter thermal-hydraulics equation. It demonstrated high efficiency and accuracy in quantifying the uncertainties of three important input parameters: external reactivity, Doppler reactivity coefficient, and coolant temperature coefficient. In another study, Gaussian process metamodel was developed as a surrogate to a TRACE model based on the FEBA benchmark, then the probabilistic framework based on Bayesian calibration was used to quantify uncertainties of physical model parameters (Wicaksono et al., 2018). Usually, the uncertainties of the boundary conditions, initial conditions, and geometry are determined by measurement uncertainty or manufacturing tolerances. In this paper, statistical parameters for the boundary conditions are estimated based on local sensitivity coefficients and their effect on the Probability Density Functions (PDFs) of input parameters is accounted for by adding a bias term to the mathematical formulation of the problem. The paper is outlined as follows: Section 2 contains the theory and the mathematical derivation for MLE and MAP algorithms with the addition of a bias term. Section 3 contains a brief description of the BFBT benchmark and the method implementation. Section 4 contains the set of obtained results with discussion and finally, Section 5 contains conclusions and suggestions for future work. 2. Theory and mathematical derivation
! h , then the likelihood function representing the probability of ! observing n independent output variables given h is calculated as follows:
! ! L h jy ¼ f y1 ; y2 ; ; yn j h
n ! ! ! Y ! f yi j h ¼ f y1 j h :f y2 j h :f yn j h ¼
To find the set of statistical parameters of input variables (e.g. ! physical model) that make the experimental data Y most probable, the logarithm of the likelihood function is maximized:
h ! io log L h jy ¼0 n Q ! @ log f yi j h ¼0 @hi i¼1 n h i P ! @ log f yi j h ¼0 @hi @ @hi
2.1. Maximum likelihood estimate (MLE) ! If f yi j h is the PDF representing the probability of observing the ith set of output variables yi given a set of statistical parameters
n
ð2Þ
i¼1
! The code output Y X can be assumed to be a linear combina! tion of a set of input parameters X with a corresponding sensitivity ! matrix Ax and a set of boundary conditions B with a corresponding sensitivity matrix Ab :
! ! ! ! Y ¼ Y 0 þ Ax X þ Ab B
ð3Þ
In Eq. (3), the sensitivity matrices Ax and Ab consist of the local sensitivity coefficients of the input parameters (e.g. physical models) and the boundary conditions, respectively. Since the linearity assumption in Eq. (3) is a fundamental part of the following derivation, it is important to mention that this assumption is valid for parameter values close to the nominal value, but becomes non-linear far from the nominal value (Hu, 2015). Assuming Gaussian distributions for the input parameters and ! ! the boundary conditions with mean vectors l x and l b and covari ! ance matrices Rx and Rb , respectively, it follows that f yi j h is a Gaussian distribution:
2 0
6 @ 6 6 6e ! f yi j h ¼ 6 6 6 6 4
1 3 !! y l y;i y;i A7 2 7 7 7 7 pffiffiffiffiffiffiffiffiffiffi 7 2p jRy;i j 7 7 5
!! y l y;i
T
R1
ð4Þ
with statistical parameters:
!
!
l y;i ¼ Y 0;i þ Ax;i ! l x þ Ab;i ! l b Ry;i ¼ Ax;i Rx ATx;i þ Ab;i Rb ATb;i þ 2Ax;i Rx;b ATb;i
ð5Þ
In Eq. (5), ly;i and Ry;i are the mean and the covariance matrix of ! f yi j h , respectively. If we assume the boundary conditions and
the input parameters to be independent (Rx;b ¼ 0Þ, we end up with:
! Two algorithms were considered in this paper, the Maximum Likelihood Estimate (MLE) algorithm and the Maximum A Posterior (MAP) algorithm. The following two sections contain a detailed derivation of the two methods.
ð1Þ
i¼1
!
l y;i ¼ Y 0;i þ Ax;i ! l x þ Ab;i ! l b Ry;i ¼ Ax;i Rx ATx;i þ Ab;i Rb ATb;i
ð6Þ
By substituting Eq. (6) into Eq. (4), the likelihood function is: 2 0 6 @ 6 ! 6 6e f yi j h ¼ 6 6 6 6 4
1 ! 1 3 ! ! T ! ! !! ! Ax;i Rx AT þAb;i Rb AT y Y 0;i Ab;i l b Ax;i l x y Y 0;i Ab;i l b Ax;i l x b;i x;i A7 2 7 7 7 7 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 7 T T 7 2p jAx;i Rx Ax;i þ Ab;i Rb Ab;i j 7 5 ð7Þ
R.A. Abu Saleem, T. Kozlowski / Annals of Nuclear Energy 133 (2019) 1–8
After that, a Maximization-Expectation (ME) algorithm is implemented by iteratively searching for optimal values for both the mean and the variance of the input parameters that would maximize the likelihood function (McLachlan and Krishnan, 2007). Detailed mathematical formulation for the ME algorithm is presented in the appendix of this paper. 2.2. Maximum A Posteriori estimate (MAP) With MAP, maximization is performed for the posterior distribution rather than the likelihood function. Based on Bayesian theorem, the posterior distribution for the statistical parameters of the physical models given a set of observable data ðyÞ can be calculated as follows:
! ! ! ! ! L h jy p h p h jy ¼ R ! ! ! ¼ K ðyÞp h L h jy L h jy p h d h
ð8Þ
where K ðyÞ is integration constant. Now, the log of the posterior is maximized as follows:
h ! io log p h jy ¼0 h !i P h !i n @ log ½ K ð y Þ þ log p h þ log f yi j h ¼0 @hi @ @hi
n
ð9Þ
i¼1
Because K ðyÞ is constant, Eq. (9) simplifies into: n h !i X h !i @ @ log p h þ log f yi j h ¼0 @hi @hi i¼1
ð10Þ
The second part of Eq. (10) is the same as that of the MLE method. As for the first part, a prior distribution is needed. The posterior is assumed to follow a normal distribution with unknown mean and standard deviation. In this case, the conjugate prior distribution is given by a normal inverse Gamma distribution:
p l; r ¼ p ljr p r ¼ N b; r 2
2
2
2
=k C1 ðr; kÞ
Explicitly:
pffiffiffi r kk
p l; r2 ¼ pffiffiffiffiffiffiffi rð2rþ2þ1Þ e 2pCðr Þ
ðkðlbÞ2 þ2kÞ 2r 2
ð11Þ
The prior knowledge of the statistical parameters is dictated by the assignment of the hyper parameters (r; k; b, and k) of the distribution in Eq. (11). Taking the logarithm of Eq. (11) and derivative ! with respect to the mean vector l yields:
h !i @ ! ! ¼ K R1 lb x ! log p h @l
ð12Þ
Eq. (12) can be used now for implementing the maximization expectation algorithm. Starting with the maximization step:
h !i @ @ þ ! ! log p h @ lx @ lx 2 0 ! ! 1 3 ! ! ! ! T ! ! y Y 0;i Ab;i l b Ax;i l x Ry 1 y Y 0;i Ab;i l b Ax;i l x 6 @ A7 2 7 6 7 6 n X 6e 7 7 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi log 6 7 6 T T 7 6 i¼1 jA p R A þ A R A j 2 x x;i b;i b x;i b;i 7 6 5 4 ¼0 ð13Þ This can be written in matrix form as:
! !
! M þ M prior l x ¼ B þ B prior
3
ð14Þ
! where M and B are defined in Eq. (A.4) (see appendix), and Mprior ! and B prior are given by:
Mprior ¼ K R1 x ! ! B prior ¼ K R1 x b
ð15Þ
The expectation step for MAP is the same as for the MLE method. The mathematical formulation derived in this section was applied to determine the probability density functions of two physical models (the interfacial friction coefficient and the wall to liquid friction coefficient) in the thermal-hydraulics code RSTART based on data from the BFBT benchmark. Methodology implementation is described in detail in Section 3 of this paper. 3. BFBT benchmark and methodology implementation The OECD/NEA BWR Full-size Fine-mesh Bundle Test (BFBT) benchmark is one of the most valuable databases for the modeling of Boiling Water Reactor (BWR) thermal-hydraulics. This publicly available benchmark contains sub-channel void fraction measurements in a full-scale BWR fuel assembly and has been used in several studies to perform IUQ of physical model parameters for different thermal-hydraulics codes (Hu and Kozlowski, 2016; Cacuci and Arslan, 2014; Neykov, 2006). In this study, IUQ was performed based on BFBT data to estimate the PDFs of the input parameters of the thermal-hydraulics code RSTART developed by Hu and Kozlowski (2018a,b). A brief description of the BFBT benchmark and methodology implementation is provided in the following sections. 3.1. The BFBT benchmark BFBT is a benchmark related to the thermal-hydraulics subchannel analysis in BWR rod bundles with two-phase flow. The full-scale BWR fuel assembly is simulated by an electrically heated rod bundle (Neykov, 2006). Void fraction and critical power measurements were performed with typical BWR reactor power and boundary conditions. With the specifications of the BFBT benchmark, code predictions of void fractions and critical power can be compared to full-scale experimental data in a prototypical BWR rod bundle. The BFBT benchmark contains different data sets of void fraction and critical power measurements. These data sets correspond to different inlet and outlet boundary conditions (pressure, inlet temperature, flow rate, and power). The ranges of the experimental conditions are shown in Table 1. For void fraction, there were 392 data sets of steady-state measurements and 2 data sets of transient measurements. Five different types of bundle assemblies were tested in the BFBT benchmark (Assemblies 0–4) corresponding to different combinations of geometries and axial power distributions. In this study, we used 86 data sets of steady-state void fraction measurements corresponding to assembly 4 with uniform axial power profile. There were two types of void distribution measurement systems used in BFBT: an X-ray CT scanner and an X-ray densitometer. Three X-ray densitometers (DEN #3, DEN #2, and DEN #1) were located at three axial elevations within the heated section of the bundle, whereas the X-ray CT scanner was located 50 mm above the heated section (see Fig. 1). The steady-state cross-sectional averaged void fraction measurements from the three X-ray densitometers are used in this study for comparison with the RSTART predictions.
4
R.A. Abu Saleem, T. Kozlowski / Annals of Nuclear Energy 133 (2019) 1–8 Table 1 Variations of experimental conditions for the BFBT benchmark. Parameter
Variation
Power Pressure Flow rate Inlet temperature Inlet subcooling Exit quality Exit void fraction
0.62–7.3 (MW) 3.9–8.7 (MPa) 10–70 (t/h) 238–292 (°C) 50–56 (kJ/kg) 8–25 (%) 45–90 (%)
3.2. Methodology implementation The RSTART model was composed of a single pipe with 24 cells equal in length with inlet and outlet boundary conditions corresponding to the 86 data sets of the BFBT benchmark. No approximations were made to the geometry in the RSTART model and the uncertainty due to geometry was assumed negligible. Data related to the three X-ray densitometers (DEN #3, DEN #2, and DEN #1) was divided into two groups. The first group of data constitutes experimental measurements and code prediction for DEN #2. This data set was used only for MLE and MAP implementation. The second group of data that constitutes experimental measurements and code predictions for DEN #1 and DEN #3. This data set was used only for methodology validation. It was entirely anonymous to the MLE and MAP algorithms throughout the estimation process and was used later to validate the results and examine the possibility of result overfitting. The first step towards IUQ implementation is to perform a sensitivity analysis and calculate the local sensitivity coefficients for the parameters and boundary conditions of interest. RSTART code determines the local sensitivity coefficients for 14 input parameters and boundary conditions based on discrete adjoint methods (Hu and Kozlowski, 2018c–e). From a simple investigation of the values of the local sensitivity coefficients, the interfacial friction
Fig. 2. Methodology implementation to the BFBT benchmark.
coefficient (f i ) and the wall to liquid friction coefficient (f wl ) exhibited two of the largest average values and, therefore, were considered for this study as the most important input parameters. A flow chart depicting the methodology implementation to the BFBT benchmark is shown in Fig. 2. 4. Results and discussion Two input parameters (physical models) uncertainty were estimated in this study based on their high sensitivity coefficients, the interfacial friction coefficient (f i ) and the wall to liquid friction coefficient (f wl ). To study the effect of the bias term (boundary conditions), two sets of results were obtained from the thermal-hydraulics code. The first set was obtained using MLE and MAP algorithms without
Table 2 Statistical parameters using MLE. Parameter
Perturbation mean
Variance
Interfacial friction coefficient (f i ) Wall to liquid friction coefficient (f wl )
0.88 0.30
0.94 0.80
Table 3 Statistical parameters using MAP. Parameter
Perturbation mean
Variance
Interfacial friction coefficient (f i ) Wall to liquid friction coefficient (f wl )
0.86 0.30
0.94 0.80
Table 4 Statistical parameters for boundary conditions obtained using MLE.
Fig. 1. Void fraction measurement at four axial elevations denoted by: DEN #3, DEN #2, DEN #1 and CT (Neykov, 2006).
Boundary condition
Perturbation mean
Variance
Outlet pressure (pout ) Inlet liquid temperature (Tl;in ) Inlet liquid velocity (ul;in )
5.8 E-05 0.003 0.17
3.6 E-07 1.9 E-07 0.003
5
R.A. Abu Saleem, T. Kozlowski / Annals of Nuclear Energy 133 (2019) 1–8 Table 5 Updated statistical parameters using MLE with bias term. Parameter
Perturbation mean
Variance
Interfacial friction coefficient (f i ) Wall to liquid friction coefficient (f wl )
0.30 0.26
0.004 0.09
Table 6 Updated statistical parameters using MAP with bias term. Parameter
Perturbation mean
Variance
Interfacial friction coefficient (f i ) Wall to liquid friction coefficient (f wl )
0.05 0.10
0.006 0.07
Table 7 Average absolute difference between experimental data and code prediction for void fraction. Method used
Without bias term
With bias term
Nominal code values MLE MAP
0.030 0.020 0.020
0.014 0.013
accounting for the bias term. The second set was obtained using MLE and MAP algorithms with bias terms taken into account. For the first set of results, the statistical parameters of the perturbation for the two input parameters were obtained using MLE
and MAP algorithms without accounting for the bias term. Tables 2 and 3 contain the estimated statistical parameters using MLE and MAP, respectively. For the second set of results, statistical parameters for the boundary conditions were needed. To obtain reasonable estimates for such statistical parameters, the MLE algorithm was applied to the three most influential boundary conditions, namely, the outlet pressure (pout ), the inlet liquid temperature (T l;in ) and the inlet liquid velocity (ul;in ). Results are shown in Table 4. The statistical parameters of the boundary conditions in Table 4 were then used in a new calculation for both MLE and MAP algorithms to get updated estimates for the statistical parameters of the two input parameters, the interfacial friction coefficient (f i ) and the wall to liquid friction coefficient (f wl ). Results of the updated estimates using MLE and MAP algorithms are shown in Tables 5 and 6, respectively. To investigate the effect of adding the bias term, the average absolute difference between code prediction of void fraction and experimental data was calculated for each case, results of this comparison are reported in Table 7. It can be seen that even without bias term MLE and MAP algorithms led to a noticeable decrease in the average absolute difference. However, results from MLE and MAP algorithms with added bias term led to further reduction in the average absolute difference. A comparison between code prediction and experimental data with and without the bias term for MLE and MAP algorithms is shown in Fig. 3.
Fig. 3. Comparison between code prediction of void fraction and experimental data using (a) MLE algorithm (b) MAP algorithm.
Table 8 Average absolute difference between code prediction of void fraction and experimental results for validation data. Method used
Validation data set 1 Without bias
Nominal code values MLE MAP
0.049 0.037 0.037
Validation data set 2 With bias
Without bias
With bias
0.015 0.016
0.021 0.017 0.017
0.013 0.013
6
R.A. Abu Saleem, T. Kozlowski / Annals of Nuclear Energy 133 (2019) 1–8
Fig. 4. MLE and MAP results for validation set 1.
Fig. 5. MLE and MAP results for validation set 2.
The above results were calculated with the same data as used for the MLE and MAP algorithm, which gives a possibility that the statistical parameters for the input (physical models) and bias (boundary conditions) were over fitted to this particular data set.
To investigate the possibility of over fitting, the average absolute difference between code prediction and experimental data was calculated for each case using the new (validation) data. This data was not used for parameter estimation using MLE and MAP, i.e. void
7
R.A. Abu Saleem, T. Kozlowski / Annals of Nuclear Energy 133 (2019) 1–8
fraction data corresponding to DEN #1 and DEN #3. Results of this comparison are reported in Table 8. It can be seen from Table 8 that although the validation data was anonymous to the MLE and MAP algorithms, it still led to a better code prediction results. Here, the effect of the bias term in reducing the average absolute error is evident. A comparison between code prediction and experimental data for validation data sets 1 and 2 is shown in Figs. 4 and 5, respectively. 5. Summary and conclusions Proper quantification of the uncertainties in the input parameters (e.g. physical models) of a simulation code is an essential part of forward Uncertainty Quantification (UQ). Such quantification can be systematically achieved through an Inverse Uncertainty Quantification (IUQ) process, where the code results and experimental measurements are used to estimate uncertainty of the code input parameters. There are several factors affecting the estimation results of IUQ. Some of these are geometry uncertainties, measurement errors, and uncertainties in boundary conditions. In this paper, a mathematical framework with additional bias term was developed based on the Maximum Likelihood Estimate (MLE) and Maximum A Posterior (MAP) algorithms to investigate the effect of uncertainties in boundary conditions on the input uncertainty estimation. The methodology was tested with the well-known BFBT benchmark and a new thermal-hydraulics code RSTART. The sensitivity analysis and calculation of local sensitivity coefficients were obtained from the code’s adjoint results. Experimental data for void fraction distribution from three X-ray densitometers (DEN #3, DEN #2, and DEN #1) were used. The IUQ problem with and without bias term was implemented using the data from DEN #2. Data from the other two densitometers was used for validation and investigation of overfitting. By calculating the average absolute difference between experimental data and code prediction, the results showed a noticeable reduction in the average error for the case with the bias term. It was also shown that the IUQ results led to improved results when used for the validation data, reducing the possibility of overfitting.
where:
M¼ ! B ¼
n P i¼1 n P i¼1
ATx;i R1 y Ax;i ! ! ! ATx;i R1 y Y 0;i ATx;i R1 y y Ab;i l b
ðA:4Þ
After the completion of the maximization step, an expectation step is implemented to obtain an estimate of the covariance matrix based on an updated value for the mean vector. Using the Linear Minimum Mean Square Error (LMMSE) estimate:
j 2
r
new
2 ! 2 ! b j lj jY b j jY X ¼ E Xj X i þE i h !i bj X b j lj jY i þ 2E X j X
ðA:5Þ
b j is the LMMSE of X j conditioned on ! X Y i and defined by:
! !1 ! ! b j ¼ lj þ Cov X j ; ! Y i Cov Y i ; Y i Y i l y;i X
ðA:6Þ
It has the following properties: ! ! !1 2 ! ! j j b j jY E Xj X Cov Y i ;X j Cov X j ; Y i Cov Y i ; Y i i ¼ Cov X ;X ! ! !1 ! 2 ! 2 ! jY i ¼ Cov X j ; Y i Cov Y i ; Y i Y i l y;i h !i bj X b j lj jY E Xj X i ¼0
E
X j lj
ðA:7Þ
Substituting Eq. (14) into Eq. (12), and considering N number of observations, the updated covariance matrix is given by: N n 1 X ! T 1 ! RTxy;i R1 y i l y;i y;i Rxy;i þ Rxy;i Ry;i N i¼1 T ! ! ðA:8Þ y i l y;i RTxy;i R1 y;i
Rx;new ¼ Rx;old þ
where:
! !
Rxy;i ¼ Cov X ; Y i ¼ Ax;i Rx
ðA:9Þ
Appendix
Appendix A. Supplementary data
The Maximization-Expectation algorithm
Supplementary data to this article can be found online at https://doi.org/10.1016/j.anucene.2019.05.005.
The algorithm starts with a maximization step where the likelihood function is maximized as follows: 8 > > > > > > > >
2 0
6 @ 6 6 6e @ 6 log ! > 6 @ lx> 6 i¼1 > > 6 > > 4 > > :
1 3 9 !! !! ! ! ! ! > > y Y 0;i Ab;i l b Ax;i l x y Y 0;i Ab;i l b Ax;i l x Ax;i Rx AT þAb;i Rb AT x;i b;i > A 7> 2 7> > > 7> 7= 7 rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 7> 7> > 2p Ax;i Rx ATx;i þ Ab;i Rb ATb;i 7> > 5> > > ; T
1
¼0 ðA:1Þ
This leads to: n X
!
ATx;i
X T ! ! ! Ry Ax;i ! lx ¼ Ax;i Ry 1 y Y 0;i ATx;i Ry 1 Ab;i l b
i¼1
n
1
i¼1
ðA:2Þ or, in a matrix form:
! Mlx ¼ B !
ðA:3Þ
References Badea, M.C., Cacuci, D.G., Badea, A.F., 2012. Best-estimate predictions and model calibration for reactor thermal hydraulics. Nucl. Sci. Eng. 172 (1), 1–19. Cacuci, D.G., 2014. Predictive modeling of coupled multi-physics systems: I. theory. Ann. Nucl. Energy 70, 266–278. Cacuci, D.G., Arslan, E., 2014. Reducing uncertainties via predictive modeling: FLICA4 calibration using BFBT benchmarks. Nucl. Sci. Eng. 176 (3), 339–349. Cacuci, D.G., Badea, M.C., 2014. Predictive modeling of coupled multi-physics systems: II. Illustrative application to reactor physics. Ann. Nucl. Energy 70, 279–291. Cacuci, D.G., Ionescu-Bujor, M., 2010. Sensitivity and uncertainty analysis, data assimilation, and predictive best-estimate model calibration. In: Cacuci, D.G. (Ed.), Handbook of Nuclear Engineering. Springer New York/Berlin, pp. 1913– 2051 (Chapter 17 in vol. 3, ISBN: 978-0-387-98150-5). D’Auria, F., Bousbia-Salah, A., Petruzzi, A., Nevo, A.D., 2006. State of the art in using best estimate calculation tools in nuclear technology. Nucl. Eng. Technol. 38 (1), 11–32. D’Auriaa, F., Mazzantini, O., 2009. The Best-Estimate Plus Uncertainty (BEPU) challenge in the licensing of current generation of reactors. International conference on opportunities and challenges for water cooled reactors in the 21 century. Vienna, Austria. Gelman, A., Carlin, J.B., Stern, H.S., Dunson, D.B., Vehtari, A., Rubin, D.B., 2014. Bayesian Data Analysis. Chapman & Hall/CRC, Boca Raton, FL, USA.
8
R.A. Abu Saleem, T. Kozlowski / Annals of Nuclear Energy 133 (2019) 1–8
Gilks, W.R., 2005. Markov Chain Monte Carlo. John Wiley & Sons Ltd.. Hu, G., and Kozlowski, T., 2018. Application of discrete adjoint method to sensitivity and uncertainty analysis in steady-state two-phase flow simulations. arXiv preprint, arXiv:1805.01451. Hu, G., and Kozlowski, T., 2018. Assessment of continuous and discrete adjoint method for sensitivity analysis in two-phase flow simulations. arXiv preprint, arXiv:1805.08083. Hu, G., Kozlowski, T., 2016. Inverse uncertainty quantification of trace physical model parameters using BFBT benchmark data. Ann. Nucl. Energy 96, 197–203. Hu, G., Kozlowski, T., 2018c. Application of continuous adjoint method to steadystate two-phase flow simulation. Ann. Nucl. Energy 117, 202–212. Hu, G., Kozlowski, T., 2018b. Application of implicit Roe-type scheme and Jacobianfree Newton-Krylov method to two-phase problems. Ann. Nucl. Energy 119, 180–190. Hu, G., Kozlowski, T., 2018a. A Roe-type numerical solver for the two-phase twofluid six-equation model with realistic equation of state. Nucl. Eng. Des. 326, 354–370. Hu, G., 2015. Inverse Uncertainty Quantification of TRACE Physical Model Parameters using Bayesian Analysis (thesis). Kennedy, M.C., O’Hagan, A., 2001. Bayesian calibration of computer models. J. R. Stat. Soc. B 63 (3), 425–464. Latten, C., Cacuci, D.G., 2014. Predictive modeling of coupled systems: uncertainty reduction using multiple reactor physics benchmarks. Nucl. Sci. Eng. 178 (2).
Lee, S.H., Chen, W., 2009. A comparative study of uncertainty propagation methods for black-box-type problems. Struct. Multidiscip. Opt. 37 (3), 239–253. McLachlan, G.J., Krishnan, T., 2007. The EM Algorithm and Extensions. John Wiley & Sons. Mooney, C.Z., 1997. Monte Carlo Simulation. Sage Publications. Neykov, B., 2006. NUPEC BWR Full-size Fine-mesh Bundle Test (BFBT) Benchmark: Specifications. Nuclear Energy Agency Organization for Economic Cooperation and Development. Petruzzi, A., Cacuci, D.G., D’Auria, F., 2010. Best-estimate model calibration and prediction through experimental data assimilation-II: application to a blowdown benchmark experiment. Nucl. Sci. Eng. 165 (1), 45–100. Robert, C., Casella, G., 2013. Monte Carlo Statistical Methods. Springer Science & Business Media. Shrestha, R., Kozlowski, T., 2015. Inverse uncertainty quantification of input model parameters for thermal-hydraulics simulations using expectation– maximization under Bayesian framework. J. Appl. Stat., 1–16 Wicaksono, D.C., Pautz, A., Zerkak, O,. 2018. Bayesian Uncertainty Quantification of Physical Models in Thermal-Hydraulics System Codes (thesis). Wu, X., Kozlowski, T., 2017. Inverse uncertainty quantification of reactor simulations under the Bayesian framework using surrogate models constructed by polynomial chaos expansion. Nucl. Eng. Des. 313, 29–52.
Contents lists available at ScienceDirect
Annals of Nuclear Energy journal homepage: www.elsevier.com/locate/anucene
Estimation of probability Density Functions for model input parameters using inverse uncertainty quantification with bias terms Rabie A. Abu Saleem a,⇑, Tomasz Kozlowski b a b
Nuclear Engineering Department, Jordan University of Science and Technology, P.O. Box 3030, Irbid 22110, Jordan Department of Nuclear, Plasma and Radiological Engineering, University of Illinois at Urbana-Champaign, 216 Talbot Laboratory, 104 S. Wright St., Urbana, IL 61801, USA
a r t i c l e
i n f o
Article history: Received 14 November 2018 Received in revised form 2 May 2019 Accepted 4 May 2019
Keywords: Inverse uncertainty quantification Maximum Likelihood Estimation Maximum A Posterior estimation Bias terms BFBT
a b s t r a c t The documentation of most nuclear thermal-hydraulics codes does not provide sufficient information on uncertainty of physical models (e.g. interfacial heat transfer coefficients). These models were derived based on experimental data and implemented as empirical correlations in the computational code. The uncertainty quantification for the relevant output quantity (e.g. Peak Cladding Temperature) requires estimation of the Probability Density Functions (PDFs) of the code inputs, such as physical models. In this paper, we investigate the effect of boundary conditions (outlet pressure, inlet liquid temperature, and inlet flow rate) on the uncertainty of two physical models (the interfacial friction coefficient and the wall to liquid friction coefficient). The boundary conditions effect was accounted for by adding a bias term to the mathematical framework of two existing methods for Inverse Uncertainty Quantification (IUQ): the Maximum Likelihood Estimation (MLE) method and the Maximum A Posterior (MAP) method. The two methods were demonstrated using the BFBT benchmark, experimental data was compared to code predictions of the RSTART thermal-hydraulics code for two different cases: without and with bias term. The results show an evident improvement in code prediction when the bias term is used. Finally, a validation set of experimental data was used to investigate the possibility of data overfitting, and the proposed methodology showed absence of overfitting when bias terms are used. Ó 2019 Elsevier Ltd. All rights reserved.
1. Introduction Uncertainty Quantification (UQ) plays a crucial role in the verification and validation of nuclear codes. One of the main goals of UQ is to develop rigorous methods to properly characterize and quantify the variability in model input parameters and output results of a simulation code (Shrestha and Kozlowski, 2015). The use of Best Estimate Plus Uncertainty (BEPU) approach for thermal-hydraulics transient calculations has been promoted by the U.S. NRC (Nuclear Regulatory Committee) over the conservative approach for the purpose of understanding the performance of Emergency Core Cooling Systems (ECCSs) during a scenario of a reactor transient (D’Auria et al., 2006). In the conservative approach, the output response of a code is calculated using extreme values of input parameters. In other words, the values used in such codes are believed to produce the most unfavorable results. The reason behind this approach is to model the physical phenomena such that it always predicts the worst-case scenario. On the other hand, the BEPU approach is meant to predict the ⇑ Corresponding author. E-mail addresses: [email protected] (R.A. Abu Saleem), [email protected] (T. Kozlowski). https://doi.org/10.1016/j.anucene.2019.05.005 0306-4549/Ó 2019 Elsevier Ltd. All rights reserved.
results as realistically as possible by predicting the mean of experimental data and considering the effects of all the important variables. With this approach, some form of uncertainty is provided to account for omitting the variables that are not possible or practical to consider for a given phenomenon (Hu and Kozlowski, 2016). BEPU approach enables more precise calculation of safety margins and aids in regulating the licensing and operation of nuclear reactors with less cost compared to the conservative approach. There are two types of uncertainty quantification; forward Uncertainty Quantification (UQ) and Inverse Uncertainty Quantification (IUQ). Forward uncertainty quantification aims at quantifying the overall uncertainty in the results (output) of a code by propagating the uncertainties of the input model parameters of that code. Inverse uncertainty quantification, on the other hand, entails the process of determining the values and ranges of the input code parameters of a code by using experimental data along with code prediction data (Kennedy and O’Hagan, 2001). In other words, it determines Probability Density Functions (PDFs) of the code’s input model parameters. These PDFs can then be used as the start point for forward UQ analysis leading to more precise safety margins. For forward UQ problems, Monte Carlo methods, perturbation methods, and spectral expansion methods are some of the most commonly used (Lee and Chen, 2009; Robert and
2
R.A. Abu Saleem, T. Kozlowski / Annals of Nuclear Energy 133 (2019) 1–8
Casella, 2013; Mooney, 1997). Common methods for inverse problems include Maximum-Likelihood Estimation (MLE), Maximum A Posteriori (MAP), and Markov Chain Monte Carlo (MCMC) (Shrestha and Kozlowski, 2015; Gelman et al., 2014; McLachlan and Krishnan, 2007; Gilks, 2005). Many researchers have investigated using IUQ in verifying and validating results of different nuclear codes. In some studies, an Expectation-Maximization (EM) algorithm was applied within a mathematical framework based on MLE and MAP methods (Shrestha and Kozlowski, 2015; D’Auriaa and Mazzantini, 2009). These methods were constrained with some limitations based on the assumptions used to derive them, some of these constraints are: assuming a linear correlation between the code output and the input parameter and the need for local sensitivity analysis (Wu and Kozlowski, 2017). Some other researchers established a mathematical framework for estimating and quantifying the uncertainties in input parameters and code results concurrently, and they were able to verify their method for time-dependent nonlinear systems of large scale (Cacuci et al., 2010; Badea et al., 2012; Petruzzi et al., 2010; Cacuci and Arslan, 2014). Furthermore, this predictive method was extended to develop an innovative mathematical procedure for predictive modeling of coupled multi-physics systems (Cacuci, 2014), and was applied to different benchmarks to demonstrate its ability to produce best estimate responses and parameter values with reduced uncertainties (Cacuci and Badea, 2014; Latten and Cacuci, 2014). Other studies explored the implementation of inference methods with metamodels rather than system codes. One study investigated using Bayesian inference methods for the IUQ problem with surrogate models derived based on polynomial chaos expansion (Wu and Kozlowski, 2017). This method was applied for a point reactor kinetics equation coupled with a lumped parameter thermal-hydraulics equation. It demonstrated high efficiency and accuracy in quantifying the uncertainties of three important input parameters: external reactivity, Doppler reactivity coefficient, and coolant temperature coefficient. In another study, Gaussian process metamodel was developed as a surrogate to a TRACE model based on the FEBA benchmark, then the probabilistic framework based on Bayesian calibration was used to quantify uncertainties of physical model parameters (Wicaksono et al., 2018). Usually, the uncertainties of the boundary conditions, initial conditions, and geometry are determined by measurement uncertainty or manufacturing tolerances. In this paper, statistical parameters for the boundary conditions are estimated based on local sensitivity coefficients and their effect on the Probability Density Functions (PDFs) of input parameters is accounted for by adding a bias term to the mathematical formulation of the problem. The paper is outlined as follows: Section 2 contains the theory and the mathematical derivation for MLE and MAP algorithms with the addition of a bias term. Section 3 contains a brief description of the BFBT benchmark and the method implementation. Section 4 contains the set of obtained results with discussion and finally, Section 5 contains conclusions and suggestions for future work. 2. Theory and mathematical derivation
! h , then the likelihood function representing the probability of ! observing n independent output variables given h is calculated as follows:
! ! L h jy ¼ f y1 ; y2 ; ; yn j h
n ! ! ! Y ! f yi j h ¼ f y1 j h :f y2 j h :f yn j h ¼
To find the set of statistical parameters of input variables (e.g. ! physical model) that make the experimental data Y most probable, the logarithm of the likelihood function is maximized:
h ! io log L h jy ¼0 n Q ! @ log f yi j h ¼0 @hi i¼1 n h i P ! @ log f yi j h ¼0 @hi @ @hi
2.1. Maximum likelihood estimate (MLE) ! If f yi j h is the PDF representing the probability of observing the ith set of output variables yi given a set of statistical parameters
n
ð2Þ
i¼1
! The code output Y X can be assumed to be a linear combina! tion of a set of input parameters X with a corresponding sensitivity ! matrix Ax and a set of boundary conditions B with a corresponding sensitivity matrix Ab :
! ! ! ! Y ¼ Y 0 þ Ax X þ Ab B
ð3Þ
In Eq. (3), the sensitivity matrices Ax and Ab consist of the local sensitivity coefficients of the input parameters (e.g. physical models) and the boundary conditions, respectively. Since the linearity assumption in Eq. (3) is a fundamental part of the following derivation, it is important to mention that this assumption is valid for parameter values close to the nominal value, but becomes non-linear far from the nominal value (Hu, 2015). Assuming Gaussian distributions for the input parameters and ! ! the boundary conditions with mean vectors l x and l b and covari ! ance matrices Rx and Rb , respectively, it follows that f yi j h is a Gaussian distribution:
2 0
6 @ 6 6 6e ! f yi j h ¼ 6 6 6 6 4
1 3 !! y l y;i y;i A7 2 7 7 7 7 pffiffiffiffiffiffiffiffiffiffi 7 2p jRy;i j 7 7 5
!! y l y;i
T
R1
ð4Þ
with statistical parameters:
!
!
l y;i ¼ Y 0;i þ Ax;i ! l x þ Ab;i ! l b Ry;i ¼ Ax;i Rx ATx;i þ Ab;i Rb ATb;i þ 2Ax;i Rx;b ATb;i
ð5Þ
In Eq. (5), ly;i and Ry;i are the mean and the covariance matrix of ! f yi j h , respectively. If we assume the boundary conditions and
the input parameters to be independent (Rx;b ¼ 0Þ, we end up with:
! Two algorithms were considered in this paper, the Maximum Likelihood Estimate (MLE) algorithm and the Maximum A Posterior (MAP) algorithm. The following two sections contain a detailed derivation of the two methods.
ð1Þ
i¼1
!
l y;i ¼ Y 0;i þ Ax;i ! l x þ Ab;i ! l b Ry;i ¼ Ax;i Rx ATx;i þ Ab;i Rb ATb;i
ð6Þ
By substituting Eq. (6) into Eq. (4), the likelihood function is: 2 0 6 @ 6 ! 6 6e f yi j h ¼ 6 6 6 6 4
1 ! 1 3 ! ! T ! ! !! ! Ax;i Rx AT þAb;i Rb AT y Y 0;i Ab;i l b Ax;i l x y Y 0;i Ab;i l b Ax;i l x b;i x;i A7 2 7 7 7 7 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 7 T T 7 2p jAx;i Rx Ax;i þ Ab;i Rb Ab;i j 7 5 ð7Þ
R.A. Abu Saleem, T. Kozlowski / Annals of Nuclear Energy 133 (2019) 1–8
After that, a Maximization-Expectation (ME) algorithm is implemented by iteratively searching for optimal values for both the mean and the variance of the input parameters that would maximize the likelihood function (McLachlan and Krishnan, 2007). Detailed mathematical formulation for the ME algorithm is presented in the appendix of this paper. 2.2. Maximum A Posteriori estimate (MAP) With MAP, maximization is performed for the posterior distribution rather than the likelihood function. Based on Bayesian theorem, the posterior distribution for the statistical parameters of the physical models given a set of observable data ðyÞ can be calculated as follows:
! ! ! ! ! L h jy p h p h jy ¼ R ! ! ! ¼ K ðyÞp h L h jy L h jy p h d h
ð8Þ
where K ðyÞ is integration constant. Now, the log of the posterior is maximized as follows:
h ! io log p h jy ¼0 h !i P h !i n @ log ½ K ð y Þ þ log p h þ log f yi j h ¼0 @hi @ @hi
n
ð9Þ
i¼1
Because K ðyÞ is constant, Eq. (9) simplifies into: n h !i X h !i @ @ log p h þ log f yi j h ¼0 @hi @hi i¼1
ð10Þ
The second part of Eq. (10) is the same as that of the MLE method. As for the first part, a prior distribution is needed. The posterior is assumed to follow a normal distribution with unknown mean and standard deviation. In this case, the conjugate prior distribution is given by a normal inverse Gamma distribution:
p l; r ¼ p ljr p r ¼ N b; r 2
2
2
2
=k C1 ðr; kÞ
Explicitly:
pffiffiffi r kk
p l; r2 ¼ pffiffiffiffiffiffiffi rð2rþ2þ1Þ e 2pCðr Þ
ðkðlbÞ2 þ2kÞ 2r 2
ð11Þ
The prior knowledge of the statistical parameters is dictated by the assignment of the hyper parameters (r; k; b, and k) of the distribution in Eq. (11). Taking the logarithm of Eq. (11) and derivative ! with respect to the mean vector l yields:
h !i @ ! ! ¼ K R1 lb x ! log p h @l
ð12Þ
Eq. (12) can be used now for implementing the maximization expectation algorithm. Starting with the maximization step:
h !i @ @ þ ! ! log p h @ lx @ lx 2 0 ! ! 1 3 ! ! ! ! T ! ! y Y 0;i Ab;i l b Ax;i l x Ry 1 y Y 0;i Ab;i l b Ax;i l x 6 @ A7 2 7 6 7 6 n X 6e 7 7 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi log 6 7 6 T T 7 6 i¼1 jA p R A þ A R A j 2 x x;i b;i b x;i b;i 7 6 5 4 ¼0 ð13Þ This can be written in matrix form as:
! !
! M þ M prior l x ¼ B þ B prior
3
ð14Þ
! where M and B are defined in Eq. (A.4) (see appendix), and Mprior ! and B prior are given by:
Mprior ¼ K R1 x ! ! B prior ¼ K R1 x b
ð15Þ
The expectation step for MAP is the same as for the MLE method. The mathematical formulation derived in this section was applied to determine the probability density functions of two physical models (the interfacial friction coefficient and the wall to liquid friction coefficient) in the thermal-hydraulics code RSTART based on data from the BFBT benchmark. Methodology implementation is described in detail in Section 3 of this paper. 3. BFBT benchmark and methodology implementation The OECD/NEA BWR Full-size Fine-mesh Bundle Test (BFBT) benchmark is one of the most valuable databases for the modeling of Boiling Water Reactor (BWR) thermal-hydraulics. This publicly available benchmark contains sub-channel void fraction measurements in a full-scale BWR fuel assembly and has been used in several studies to perform IUQ of physical model parameters for different thermal-hydraulics codes (Hu and Kozlowski, 2016; Cacuci and Arslan, 2014; Neykov, 2006). In this study, IUQ was performed based on BFBT data to estimate the PDFs of the input parameters of the thermal-hydraulics code RSTART developed by Hu and Kozlowski (2018a,b). A brief description of the BFBT benchmark and methodology implementation is provided in the following sections. 3.1. The BFBT benchmark BFBT is a benchmark related to the thermal-hydraulics subchannel analysis in BWR rod bundles with two-phase flow. The full-scale BWR fuel assembly is simulated by an electrically heated rod bundle (Neykov, 2006). Void fraction and critical power measurements were performed with typical BWR reactor power and boundary conditions. With the specifications of the BFBT benchmark, code predictions of void fractions and critical power can be compared to full-scale experimental data in a prototypical BWR rod bundle. The BFBT benchmark contains different data sets of void fraction and critical power measurements. These data sets correspond to different inlet and outlet boundary conditions (pressure, inlet temperature, flow rate, and power). The ranges of the experimental conditions are shown in Table 1. For void fraction, there were 392 data sets of steady-state measurements and 2 data sets of transient measurements. Five different types of bundle assemblies were tested in the BFBT benchmark (Assemblies 0–4) corresponding to different combinations of geometries and axial power distributions. In this study, we used 86 data sets of steady-state void fraction measurements corresponding to assembly 4 with uniform axial power profile. There were two types of void distribution measurement systems used in BFBT: an X-ray CT scanner and an X-ray densitometer. Three X-ray densitometers (DEN #3, DEN #2, and DEN #1) were located at three axial elevations within the heated section of the bundle, whereas the X-ray CT scanner was located 50 mm above the heated section (see Fig. 1). The steady-state cross-sectional averaged void fraction measurements from the three X-ray densitometers are used in this study for comparison with the RSTART predictions.
4
R.A. Abu Saleem, T. Kozlowski / Annals of Nuclear Energy 133 (2019) 1–8 Table 1 Variations of experimental conditions for the BFBT benchmark. Parameter
Variation
Power Pressure Flow rate Inlet temperature Inlet subcooling Exit quality Exit void fraction
0.62–7.3 (MW) 3.9–8.7 (MPa) 10–70 (t/h) 238–292 (°C) 50–56 (kJ/kg) 8–25 (%) 45–90 (%)
3.2. Methodology implementation The RSTART model was composed of a single pipe with 24 cells equal in length with inlet and outlet boundary conditions corresponding to the 86 data sets of the BFBT benchmark. No approximations were made to the geometry in the RSTART model and the uncertainty due to geometry was assumed negligible. Data related to the three X-ray densitometers (DEN #3, DEN #2, and DEN #1) was divided into two groups. The first group of data constitutes experimental measurements and code prediction for DEN #2. This data set was used only for MLE and MAP implementation. The second group of data that constitutes experimental measurements and code predictions for DEN #1 and DEN #3. This data set was used only for methodology validation. It was entirely anonymous to the MLE and MAP algorithms throughout the estimation process and was used later to validate the results and examine the possibility of result overfitting. The first step towards IUQ implementation is to perform a sensitivity analysis and calculate the local sensitivity coefficients for the parameters and boundary conditions of interest. RSTART code determines the local sensitivity coefficients for 14 input parameters and boundary conditions based on discrete adjoint methods (Hu and Kozlowski, 2018c–e). From a simple investigation of the values of the local sensitivity coefficients, the interfacial friction
Fig. 2. Methodology implementation to the BFBT benchmark.
coefficient (f i ) and the wall to liquid friction coefficient (f wl ) exhibited two of the largest average values and, therefore, were considered for this study as the most important input parameters. A flow chart depicting the methodology implementation to the BFBT benchmark is shown in Fig. 2. 4. Results and discussion Two input parameters (physical models) uncertainty were estimated in this study based on their high sensitivity coefficients, the interfacial friction coefficient (f i ) and the wall to liquid friction coefficient (f wl ). To study the effect of the bias term (boundary conditions), two sets of results were obtained from the thermal-hydraulics code. The first set was obtained using MLE and MAP algorithms without
Table 2 Statistical parameters using MLE. Parameter
Perturbation mean
Variance
Interfacial friction coefficient (f i ) Wall to liquid friction coefficient (f wl )
0.88 0.30
0.94 0.80
Table 3 Statistical parameters using MAP. Parameter
Perturbation mean
Variance
Interfacial friction coefficient (f i ) Wall to liquid friction coefficient (f wl )
0.86 0.30
0.94 0.80
Table 4 Statistical parameters for boundary conditions obtained using MLE.
Fig. 1. Void fraction measurement at four axial elevations denoted by: DEN #3, DEN #2, DEN #1 and CT (Neykov, 2006).
Boundary condition
Perturbation mean
Variance
Outlet pressure (pout ) Inlet liquid temperature (Tl;in ) Inlet liquid velocity (ul;in )
5.8 E-05 0.003 0.17
3.6 E-07 1.9 E-07 0.003
5
R.A. Abu Saleem, T. Kozlowski / Annals of Nuclear Energy 133 (2019) 1–8 Table 5 Updated statistical parameters using MLE with bias term. Parameter
Perturbation mean
Variance
Interfacial friction coefficient (f i ) Wall to liquid friction coefficient (f wl )
0.30 0.26
0.004 0.09
Table 6 Updated statistical parameters using MAP with bias term. Parameter
Perturbation mean
Variance
Interfacial friction coefficient (f i ) Wall to liquid friction coefficient (f wl )
0.05 0.10
0.006 0.07
Table 7 Average absolute difference between experimental data and code prediction for void fraction. Method used
Without bias term
With bias term
Nominal code values MLE MAP
0.030 0.020 0.020
0.014 0.013
accounting for the bias term. The second set was obtained using MLE and MAP algorithms with bias terms taken into account. For the first set of results, the statistical parameters of the perturbation for the two input parameters were obtained using MLE
and MAP algorithms without accounting for the bias term. Tables 2 and 3 contain the estimated statistical parameters using MLE and MAP, respectively. For the second set of results, statistical parameters for the boundary conditions were needed. To obtain reasonable estimates for such statistical parameters, the MLE algorithm was applied to the three most influential boundary conditions, namely, the outlet pressure (pout ), the inlet liquid temperature (T l;in ) and the inlet liquid velocity (ul;in ). Results are shown in Table 4. The statistical parameters of the boundary conditions in Table 4 were then used in a new calculation for both MLE and MAP algorithms to get updated estimates for the statistical parameters of the two input parameters, the interfacial friction coefficient (f i ) and the wall to liquid friction coefficient (f wl ). Results of the updated estimates using MLE and MAP algorithms are shown in Tables 5 and 6, respectively. To investigate the effect of adding the bias term, the average absolute difference between code prediction of void fraction and experimental data was calculated for each case, results of this comparison are reported in Table 7. It can be seen that even without bias term MLE and MAP algorithms led to a noticeable decrease in the average absolute difference. However, results from MLE and MAP algorithms with added bias term led to further reduction in the average absolute difference. A comparison between code prediction and experimental data with and without the bias term for MLE and MAP algorithms is shown in Fig. 3.
Fig. 3. Comparison between code prediction of void fraction and experimental data using (a) MLE algorithm (b) MAP algorithm.
Table 8 Average absolute difference between code prediction of void fraction and experimental results for validation data. Method used
Validation data set 1 Without bias
Nominal code values MLE MAP
0.049 0.037 0.037
Validation data set 2 With bias
Without bias
With bias
0.015 0.016
0.021 0.017 0.017
0.013 0.013
6
R.A. Abu Saleem, T. Kozlowski / Annals of Nuclear Energy 133 (2019) 1–8
Fig. 4. MLE and MAP results for validation set 1.
Fig. 5. MLE and MAP results for validation set 2.
The above results were calculated with the same data as used for the MLE and MAP algorithm, which gives a possibility that the statistical parameters for the input (physical models) and bias (boundary conditions) were over fitted to this particular data set.
To investigate the possibility of over fitting, the average absolute difference between code prediction and experimental data was calculated for each case using the new (validation) data. This data was not used for parameter estimation using MLE and MAP, i.e. void
7
R.A. Abu Saleem, T. Kozlowski / Annals of Nuclear Energy 133 (2019) 1–8
fraction data corresponding to DEN #1 and DEN #3. Results of this comparison are reported in Table 8. It can be seen from Table 8 that although the validation data was anonymous to the MLE and MAP algorithms, it still led to a better code prediction results. Here, the effect of the bias term in reducing the average absolute error is evident. A comparison between code prediction and experimental data for validation data sets 1 and 2 is shown in Figs. 4 and 5, respectively. 5. Summary and conclusions Proper quantification of the uncertainties in the input parameters (e.g. physical models) of a simulation code is an essential part of forward Uncertainty Quantification (UQ). Such quantification can be systematically achieved through an Inverse Uncertainty Quantification (IUQ) process, where the code results and experimental measurements are used to estimate uncertainty of the code input parameters. There are several factors affecting the estimation results of IUQ. Some of these are geometry uncertainties, measurement errors, and uncertainties in boundary conditions. In this paper, a mathematical framework with additional bias term was developed based on the Maximum Likelihood Estimate (MLE) and Maximum A Posterior (MAP) algorithms to investigate the effect of uncertainties in boundary conditions on the input uncertainty estimation. The methodology was tested with the well-known BFBT benchmark and a new thermal-hydraulics code RSTART. The sensitivity analysis and calculation of local sensitivity coefficients were obtained from the code’s adjoint results. Experimental data for void fraction distribution from three X-ray densitometers (DEN #3, DEN #2, and DEN #1) were used. The IUQ problem with and without bias term was implemented using the data from DEN #2. Data from the other two densitometers was used for validation and investigation of overfitting. By calculating the average absolute difference between experimental data and code prediction, the results showed a noticeable reduction in the average error for the case with the bias term. It was also shown that the IUQ results led to improved results when used for the validation data, reducing the possibility of overfitting.
where:
M¼ ! B ¼
n P i¼1 n P i¼1
ATx;i R1 y Ax;i ! ! ! ATx;i R1 y Y 0;i ATx;i R1 y y Ab;i l b
ðA:4Þ
After the completion of the maximization step, an expectation step is implemented to obtain an estimate of the covariance matrix based on an updated value for the mean vector. Using the Linear Minimum Mean Square Error (LMMSE) estimate:
j 2
r
new
2 ! 2 ! b j lj jY b j jY X ¼ E Xj X i þE i h !i bj X b j lj jY i þ 2E X j X
ðA:5Þ
b j is the LMMSE of X j conditioned on ! X Y i and defined by:
! !1 ! ! b j ¼ lj þ Cov X j ; ! Y i Cov Y i ; Y i Y i l y;i X
ðA:6Þ
It has the following properties: ! ! !1 2 ! ! j j b j jY E Xj X Cov Y i ;X j Cov X j ; Y i Cov Y i ; Y i i ¼ Cov X ;X ! ! !1 ! 2 ! 2 ! jY i ¼ Cov X j ; Y i Cov Y i ; Y i Y i l y;i h !i bj X b j lj jY E Xj X i ¼0
E
X j lj
ðA:7Þ
Substituting Eq. (14) into Eq. (12), and considering N number of observations, the updated covariance matrix is given by: N n 1 X ! T 1 ! RTxy;i R1 y i l y;i y;i Rxy;i þ Rxy;i Ry;i N i¼1 T ! ! ðA:8Þ y i l y;i RTxy;i R1 y;i
Rx;new ¼ Rx;old þ
where:
! !
Rxy;i ¼ Cov X ; Y i ¼ Ax;i Rx
ðA:9Þ
Appendix
Appendix A. Supplementary data
The Maximization-Expectation algorithm
Supplementary data to this article can be found online at https://doi.org/10.1016/j.anucene.2019.05.005.
The algorithm starts with a maximization step where the likelihood function is maximized as follows: 8 > > > > > > > >
2 0
6 @ 6 6 6e @ 6 log ! > 6 @ lx> 6 i¼1 > > 6 > > 4 > > :
1 3 9 !! !! ! ! ! ! > > y Y 0;i Ab;i l b Ax;i l x y Y 0;i Ab;i l b Ax;i l x Ax;i Rx AT þAb;i Rb AT x;i b;i > A 7> 2 7> > > 7> 7= 7 rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 7> 7> > 2p Ax;i Rx ATx;i þ Ab;i Rb ATb;i 7> > 5> > > ; T
1
¼0 ðA:1Þ
This leads to: n X
!
ATx;i
X T ! ! ! Ry Ax;i ! lx ¼ Ax;i Ry 1 y Y 0;i ATx;i Ry 1 Ab;i l b
i¼1
n
1
i¼1
ðA:2Þ or, in a matrix form:
! Mlx ¼ B !
ðA:3Þ
References Badea, M.C., Cacuci, D.G., Badea, A.F., 2012. Best-estimate predictions and model calibration for reactor thermal hydraulics. Nucl. Sci. Eng. 172 (1), 1–19. Cacuci, D.G., 2014. Predictive modeling of coupled multi-physics systems: I. theory. Ann. Nucl. Energy 70, 266–278. Cacuci, D.G., Arslan, E., 2014. Reducing uncertainties via predictive modeling: FLICA4 calibration using BFBT benchmarks. Nucl. Sci. Eng. 176 (3), 339–349. Cacuci, D.G., Badea, M.C., 2014. Predictive modeling of coupled multi-physics systems: II. Illustrative application to reactor physics. Ann. Nucl. Energy 70, 279–291. Cacuci, D.G., Ionescu-Bujor, M., 2010. Sensitivity and uncertainty analysis, data assimilation, and predictive best-estimate model calibration. In: Cacuci, D.G. (Ed.), Handbook of Nuclear Engineering. Springer New York/Berlin, pp. 1913– 2051 (Chapter 17 in vol. 3, ISBN: 978-0-387-98150-5). D’Auria, F., Bousbia-Salah, A., Petruzzi, A., Nevo, A.D., 2006. State of the art in using best estimate calculation tools in nuclear technology. Nucl. Eng. Technol. 38 (1), 11–32. D’Auriaa, F., Mazzantini, O., 2009. The Best-Estimate Plus Uncertainty (BEPU) challenge in the licensing of current generation of reactors. International conference on opportunities and challenges for water cooled reactors in the 21 century. Vienna, Austria. Gelman, A., Carlin, J.B., Stern, H.S., Dunson, D.B., Vehtari, A., Rubin, D.B., 2014. Bayesian Data Analysis. Chapman & Hall/CRC, Boca Raton, FL, USA.
8
R.A. Abu Saleem, T. Kozlowski / Annals of Nuclear Energy 133 (2019) 1–8
Gilks, W.R., 2005. Markov Chain Monte Carlo. John Wiley & Sons Ltd.. Hu, G., and Kozlowski, T., 2018. Application of discrete adjoint method to sensitivity and uncertainty analysis in steady-state two-phase flow simulations. arXiv preprint, arXiv:1805.01451. Hu, G., and Kozlowski, T., 2018. Assessment of continuous and discrete adjoint method for sensitivity analysis in two-phase flow simulations. arXiv preprint, arXiv:1805.08083. Hu, G., Kozlowski, T., 2016. Inverse uncertainty quantification of trace physical model parameters using BFBT benchmark data. Ann. Nucl. Energy 96, 197–203. Hu, G., Kozlowski, T., 2018c. Application of continuous adjoint method to steadystate two-phase flow simulation. Ann. Nucl. Energy 117, 202–212. Hu, G., Kozlowski, T., 2018b. Application of implicit Roe-type scheme and Jacobianfree Newton-Krylov method to two-phase problems. Ann. Nucl. Energy 119, 180–190. Hu, G., Kozlowski, T., 2018a. A Roe-type numerical solver for the two-phase twofluid six-equation model with realistic equation of state. Nucl. Eng. Des. 326, 354–370. Hu, G., 2015. Inverse Uncertainty Quantification of TRACE Physical Model Parameters using Bayesian Analysis (thesis). Kennedy, M.C., O’Hagan, A., 2001. Bayesian calibration of computer models. J. R. Stat. Soc. B 63 (3), 425–464. Latten, C., Cacuci, D.G., 2014. Predictive modeling of coupled systems: uncertainty reduction using multiple reactor physics benchmarks. Nucl. Sci. Eng. 178 (2).
Lee, S.H., Chen, W., 2009. A comparative study of uncertainty propagation methods for black-box-type problems. Struct. Multidiscip. Opt. 37 (3), 239–253. McLachlan, G.J., Krishnan, T., 2007. The EM Algorithm and Extensions. John Wiley & Sons. Mooney, C.Z., 1997. Monte Carlo Simulation. Sage Publications. Neykov, B., 2006. NUPEC BWR Full-size Fine-mesh Bundle Test (BFBT) Benchmark: Specifications. Nuclear Energy Agency Organization for Economic Cooperation and Development. Petruzzi, A., Cacuci, D.G., D’Auria, F., 2010. Best-estimate model calibration and prediction through experimental data assimilation-II: application to a blowdown benchmark experiment. Nucl. Sci. Eng. 165 (1), 45–100. Robert, C., Casella, G., 2013. Monte Carlo Statistical Methods. Springer Science & Business Media. Shrestha, R., Kozlowski, T., 2015. Inverse uncertainty quantification of input model parameters for thermal-hydraulics simulations using expectation– maximization under Bayesian framework. J. Appl. Stat., 1–16 Wicaksono, D.C., Pautz, A., Zerkak, O,. 2018. Bayesian Uncertainty Quantification of Physical Models in Thermal-Hydraulics System Codes (thesis). Wu, X., Kozlowski, T., 2017. Inverse uncertainty quantification of reactor simulations under the Bayesian framework using surrogate models constructed by polynomial chaos expansion. Nucl. Eng. Des. 313, 29–52.