Advanced methodology for uncertainty propagation in computer experiments with large number of inputs

TL;DR

This paper introduces an advanced statistical methodology combining Gaussian process metamodels and screening techniques to efficiently propagate uncertainties and perform sensitivity analysis in high-dimensional nuclear safety simulations, reducing computational costs.

Contribution

It proposes a novel approach that effectively handles high-dimensional inputs by combining screening, group-wise Gaussian process modeling, and residual modeling, improving accuracy and efficiency in uncertainty quantification.

Findings

Accurately estimates Sobol' sensitivity indices with few simulations

Provides precise high quantile estimates compared to empirical methods

Outperforms simple Gaussian process models in high-dimensional cases

Abstract

In the framework of the estimation of safety margins in nuclear accident analysis, a quantitative assessment of the uncertainties tainting the results of computer simulations is essential. Accurate uncertainty propagation (estimation of high probabilities or quantiles) and quantitative sensitivity analysis may call for several thousand of code simulations. Complex computer codes, as the ones used in thermal-hydraulic accident scenario simulations, are often too cpu-time expensive to be directly used to perform these studies. A solution consists in replacing the computer model by a cpu inexpensive mathematical function, called a metamodel, built from a reduced number of code simulations. However, in case of high dimensional experiments (with typically several tens of inputs), the metamodel building process remains difficult. To face this limitation, we propose a methodology which combines several advanced statistical tools: initial space-filling design, screening to identify the non-influential inputs, Gaussian process (Gp) metamodel building with the group of influential inputs as explanatory variables. The residual effect of the group of non-influential inputs is captured by another Gp metamodel. Then, the resulting joint Gp metamodel is used to accurately estimate Sobol' sensitivity indices and high quantiles (here $95\%$-quantile).The efficiency of the methodology to deal with a large number of inputs and reduce the calculation budget is illustrated on a thermal-hydraulic calculation case simulating with the CATHARE2 code a Loss Of Coolant Accident scenario in a Pressurized Water Reactor. A predictive Gp metamodel is built with only a few hundred of code simulations and allows the calculation of the Sobol' sensitivity indices. This Gp also provides a more accurate estimation of the 95%-quantile and associated confidence interval than the empirical approach, at equal calculation budget. Moreover, on this test case, the joint Gp approach outperforms the simple Gp.

Tables

Table 1. Table I : List of the 27 27 27 uncertain input parameters and associated physical models in CATHARE2 code.

Type of inputs	Inputs	pdf aaaU, LU, N and LN respectively stands for Uniform, Log-Uniform, Normal and Log-Normal probability distributions.	Physical models
Heat transfer	$X_1$	N	Departure from nucleate boiling
in the core	$X_2$	U	Minimum film stable temperature
	$X_3$	LN	HTCbbbHeat Transfer Coefficient. for steam convection
	$X_4$	LN	Wall-fluid HTC
	$X_5$	N	HTC for film boiling
Heat transfer in the steam	$X_6$	LU	HTC forced wall-steam convection
generators (SG) U-tube	$X_7$	N	Liquid-interface HTC for film condensation
Wall-steam friction in core	$X_8$	LU
Interfacial friction	$X_9$	LN	SG outlet plena and crossover legs together
	$X_{10}$	LN	Hot legs (horizontal part)
	$X_{11}$	LN	Bend of the hot legs
	$X_{12}$	LN	SG inlet plena
	$X_{13}$	LN	Downcomer
	$X_{14}$	LN	Core
	$X_{15}$	LN	Upper plenum
	$X_{16}$	LN	Lower plenum
	$X_{17}$	LN	Upper head
Condensation	$X_{18}$	LN	Downcomer
	$X_{19}$	U	Cold leg (intact)
	$X_{20}$	U	Cold leg (broken)
	$X_{27}$	U	Jet
Break flow	$X_{21}$	LN	Flashing (undersaturated)
	$X_{22}$	N	Wall-liquid friction (undersaturated)
	$X_{23}$	N	Flashing delay (undersaturated)
	$X_{24}$	LN	Flashing (saturated)
	$X_{25}$	N	Wall-liquid friction (saturated)
	$X_{26}$	LN	Global interfacial friction (saturated)

Table 2. Table II : HSIC-based sensitivity indices R HSIC 2 subscript superscript 𝑅 2 HSIC R^{2}_{\text{HSIC}} for the PII (Primary Influential Inputs), identified by independence test (Step 2).

PII	X10	X2	X12	X22	X15	X13	X9	X5	X14	X26	X27
$R_{\text{HSIC}}^2$	0.39	0.04	0.02	0.02	0.01	0.01	0.01	0.01	0.01	0.01	0.01

Table 3. Table III : Evolution of G ​ p m 𝐺 subscript 𝑝 𝑚 Gp_{m} metamodel predictivity during the sequential process building, for each new additional PII.

PII	X10	X2	X12	X22	X15	X13	X9	X5	X14	X26	X27
$Q^2$	0.60	0.64	0.70	0.79	0.81	0.83	0.85	0.85	0.87	0.87	0.87

Table 4. Table IV : First and second Sobol’ indices of PII (in %), estimated with G ​ p m 𝐺 subscript 𝑝 𝑚 Gp_{m} of the joint Gp metamodel.

PII	X10	X2	X12	X22	X15	X13	X9	X5	X14	X26	X27
$1^{\text{st}}$-order Sobol’ indices	59	3	8	8	2	1	2	0	2	0	0

Table 5. Table V : Results for the 95 % percent 95 95\% -quantile estimates of the PCT (in   ∘ ​ C superscript   C \text{ }^{\circ}\text{C} ) with its 90 % percent 90 90\% -confidence interval (CI) when available.

	ReferencefffThe reference method uses the full sample of size $n=1100$.	EmpiricalgggEmpirical and Gp-based methods are applied from the learning sample of size $n=500$.	Plug-ing,hhhThe plug-in and full-Gp estimators are based on either the homoscedastic or heteroscedastic Gp metamodel.		Full-Gpg,h (conditional simulations)
			homo-Gp	hetero-Gp	homo-Gp	hetero-Gp
	${\hat{q}}_{95}^{\text{ref}}$	${\hat{q}}_{95}^{\text{emp}}$	${\hat{q}}_{95}^{G{p}_{m,1}\text{-pi}}$	${\hat{q}}_{95}^{G{p}_{m,2}\text{-pi}}$	${\hat{q}}_{95}^{G{p}_{m,1}\text{-full}}$	${\hat{q}}_{95}^{G{p}_{m,2}\text{-full}}$
Mean	$742.28$	$746.80$	$736.26$	$735.83$	$747.11$	$741.46$
$90\%$-CI	—	$[736.7;747.41]$	—	—	$[742.93;751.32]$	$[738.76;744.17]$

Equations

$$Y=g\left(X_1,\dots ,X_d\right)$$

$$Y_m({\mathbf{X}}_{\mathbf{exp}})=\mathbb{E}(Y|{\mathbf{X}}_{\mathbf{exp}})$$

$$Y_d({\mathbf{X}}_{\mathbf{exp}})=\text{mbox}Var(Y|{\mathbf{X}}_{\mathbf{exp}})=\mathbb{E}\left[(Y-Y_m({\mathbf{X}}_{\mathbf{exp}}){)}^2|{\mathbf{X}}_{\mathbf{exp}}\right].$$

$$Q^2=1-\frac{{\sum }_{i=1}^{n_{\text{test}}}{\left(y^{(i)}-{\hat{y}}^{(i)}\right)}^2}{{\sum }_{i=1}^{n_{\text{test}}}{\left(y^{(i)}-\frac{1}{n_{\text{test}}}{\sum }_{i=1}^{n_{\text{test}}}{y}^{(i)}\right)}^2}$$

$${\hat{q}}_{95}^{\text{ref}}=742.28{}^{\circ }\text{C}.$$

$$⟨f_{X_k},C_{X_k,Y}{g}_Y{⟩}_{{\mathcal{F}}_k}=\text{Cov}\left(f_{X_k},g_Y\right).$$

$$\text{HSIC}(X_k,Y{)}_{{\mathcal{F}}_k,\mathcal{G}}=\Vert C_k{\Vert }_{HS}^2.$$

$$R_{\text{HSIC},k}^2=\frac{\text{HSIC}(X_k,Y)}{\sqrt{\text{HSIC}(X_k,X_k)\text{HSIC}(Y,Y)}}.$$

$$g(X)$$

$$S_u=\frac{\text{mbox}Var\left[g_u(X_u)\right]}{\text{mbox}Var\left[g(X)\right]}\text{ and }{S}_u^T=\sum _{v\supset u}{S}_v.$$

$$\text{mbox}Var[Y({\mathbf{X}}_{\mathbf{exp}},X_{\epsilon })]=\text{mbox}{Var}_{{\mathbf{X}}_{\mathbf{exp}}}\left[Y_m({\mathbf{X}}_{\mathbf{exp}})\right]+{\mathbb{E}}_{{\mathbf{X}}_{\mathbf{exp}}}\left[Y_d({\mathbf{X}}_{\mathbf{exp}})\right]$$

$$\text{mbox}Var(Y)=V_{\epsilon }(Y)+\sum _{u\subset {\mathbf{X}}_{\mathbf{exp}}}[V_u(Y)+V_{u\epsilon }(Y)]$$

$$\text{mbox}Var(Y_m)=\sum _{u\subset {\mathbf{X}}_{\mathbf{exp}}}{V}_u(Y_m).$$

$$S_u=\frac{V_u(Y_m)}{\text{mbox}Var(Y)}\text{for any}u\subset {\mathbf{X}}_{\mathbf{exp}}.$$

$$S_{\epsilon }^T=\frac{V_{\epsilon }(Y)+{\sum }_{u\subset {\mathbf{X}}_{\mathbf{exp}}}{V}_{u\epsilon }(Y)}{\text{mbox}Var(Y)}=\frac{{\mathbb{E}}_{{\mathbf{X}}_{\mathbf{exp}}}[Y_d({\mathbf{X}}_{\mathbf{exp}})]}{\text{mbox}Var(Y)}.$$

Full text

Keywords — Gaussian process, Metamodel, Quantile, Screening, Sobol’ indices

I.   Introduction

Best-estimate computer codes are increasingly used to estimate safety margins in nuclear accident management analysis instead of conservative procedures [1, 2]. In this context, it is essential to evaluate the accuracy of the numerical model results, whose uncertainties come mainly from the lack of knowledge of the underlying physics and the model input parameters. The so-called Best Estimate Plus Uncertainty (BEPU) methods were then developed and introduced in safety analyses, especially for thermal-hydraulic issues, and even more precisely for the large-break loss-of-coolant accident [3, 4, 5]. Its main principles rely on a probabilistic modeling of the model input uncertainties, on Monte Carlo sampling for running the thermal-hydraulic computer code on sets of inputs, and on the application of statistical tools (based for example on order statistics as the Wilks’ formula) to infer high quantiles of the output variables of interest [6, 7, 8, 9, 10].

Not restricted to the nuclear engineering domain, the BEPU approach is more largely known as the uncertainty analysis framework [11]. Quantitative assessment of the uncertainties tainting the results of computer simulations is indeed a major topic of interest in both industrial and scientific communities. One of the key issues in such studies is to get information about the output when the numerical simulations are expensive to run. For example, one often faces up with cpu time consuming numerical models and, in such cases, uncertainty propagation, sensitivity analysis, optimization processing and system robustness analysis become difficult tasks using such models. In order to circumvent this problem, a widely accepted method consists in replacing cpu-time expensive computer models by cpu inexpensive mathematical functions (called metamodels) based, e.g., on polynomials, neural networks, or Gaussian processes [12]. This metamodel is built from a set of computer code simulations and must be as representative as possible of the code in the variation domain of its uncertain parameters, with good prediction capabilities. The use of metamodels has been extensively applied in engineering issues as it provides a multi-objective tool [13]: once estimated, the metamodel can be used to perform global sensitivity analysis, as well as uncertainty propagation, optimization, or calibration studies. In BEPU-kind analyses, several works [14, 15, 16] have introduced the use of metamodels and shown how this technique can help estimate quantiles or probability of failure in thermal-hydraulic calculations.

However, the metamodeling technique is known to be relevant when simulated phenomena are related to a small number of random input variables (see [17] for example). In case of high dimensional numerical experiments (with typically several tens of inputs), and depending on the complexity of the underlying numerical model, the metamodel building process remains difficult, or even impracticable. For example, the Gaussian process (Gp) model [18] which has shown strong capabilities to solve practical problems, has some caveats when dealing with high dimensional problems. The main difficulty relies on the estimation of Gp hyperparameters. Manipulating pre-defined or well-adapted Gp kernels (as in [19, 20]) is a current research way, while several authors propose to couple the estimation procedure with variable selection techniques [21, 22, 23].

In this paper, following the latter technique, we propose a rigorous and robust method for building a Gp metamodel with a high-dimensional vector of inputs. Moreover, concerning the practical use of this metamodel, our final goal is twofold: to be able to interpret the relationships between the model inputs and outputs with a quantitative sensitivity analysis, and to have a reliable high-level quantile estimation method which does not require additional runs of the computer code.

In what follows, the system under study is generically denoted

[TABLE]

where $g(\cdot)$ is the numerical model (also called the computer code), whose output $Y$ and input parameters $X_{1},\ldots,X_{d}$ belong to some measurable spaces $\mathcal{Y}$ and $\mathcal{X}_{1},\ldots,\mathcal{X}_{d}$ respectively. $\mathbf{X}=\left(X_{1},\ldots,X_{d}\right)$ is the input vector and we suppose that $\mathcal{X}=\prod_{k=1}^{d}\mathcal{X}_{k}\subset\mathbb{R}^{d}$ and $\mathcal{Y}\subset\mathbb{R}$. For a given value of the vector of inputs $\mathbf{x}=\left(x_{1},\ldots,x_{d}\right)\in\mathbb{R}^{d}$, a simulation run of the code yields an observed value $y=g(\mathbf{x})$.

To meet our aforementioned objectives, we propose a sequential methodology which combines several relevant statistical techniques. Our approach consists in four steps:

Step 1: Design of experiments. Knowing the variation domain of the input variables, a design of $n$ numerical experiments is firstly performed and yields $n$ model output values. The obtained sample of inputs/outputs constitutes the learning sample. The goal is here to explore, the most efficiently and with a reasonable simulation budget, the domain of uncertain parameters $\mathcal{X}$ and get as much information as possible about the behavior of the simulator output $Y$. For this, we use a space-filling design (SFD) of experiments, providing a full coverage of the high-dimensional input space [12, 23]. 2. 2.

Step 2: Preliminary screening. From the learning sample, a screening technique is performed in order to identify the Primary Influential Inputs (PII) on the model output variability and rank them by decreasing influence. To achieve it, we use dependence measures with associated statistical tests. These measures have several advantages: they can be directly estimated on a SFD, the sensitivity indices that they provide are interpretable quantitatively and their efficiency for screening purpose has been recently illustrated by [24, 25, 26]. 3. 3.

Step 3: Building and validation of joint metamodel. From the learning sample, a metamodel is built to fit the simulator output $Y$. For this, we propose to use a joint Gp metamodel [27], by considering only the PII as the explanatory inputs while the other inputs (screened as non significantly influential) are considered as a global stochastic (i.e. unknown) input. Moreover, we use a sequential process to build the joint Gp where the ranked PII are successively included as explanatory inputs in the metamodel (ranking from Step 2). At each iteration, a first Gp model [22], only depending on the current set of explanatory inputs, is built to approximate the mean component. The residual effect of the other inputs is captured using a second Gp model, also function of the explanatory inputs, which approximates the variance component. The accuracy and prediction capabilities of the joint metamodel are controlled on a test sample or by cross-validation. 4. 4.

Step 4: Use of the metamodel for sensitivity analysis and uncertainty propagation. A quantitative sensitivity analysis (Step 4A) and an uncertainty propagation (Step 4B) are performed using the joint metamodel instead of the computer model, leading to a large gain of computation time [28, 27]. In this work, we are particularly interested in estimating variance-based sensitivity indices (namely Sobol’ indices) and the $95\%$-quantile of model output.

For ease of reading and understanding of the methodology, Figure 1 gives a general workflow of the articulation of the main steps. For each step, the dedicated sections, the main notations that will be properly introduced later and the key equations are also referenced, in order to provide a guideline for the reader. The paper is organized as follows. In the following section, the nuclear test case which constitutes the guideline application of the paper is described. It consists in a thermal-hydraulic calculation case which simulates an accidental scenario in a nuclear Pressurized Water Reactor. Then, each step of the above statistical methodology is detailed in a dedicated section and illustrated as the same time on the use-case. A last section concludes the work.

II.   Thermal-hydraulic test case

Our use-case consists in thermal-hydraulic computer experiments, typically used in support of regulatory work and nuclear power plant design and operation. Indeed, some safety analysis considers the so-called “Loss Of Coolant Accident” which takes into account a double-ended guillotine break with a specific size piping rupture. The test case under study is a simplified one, as regards both physical phenomena and dimensions of the system, with respect to a realistic modeling of the reactor. The numerical model is based on code CATHARE2 (V2.5_3mod3.1) which simulates the time evolution of physical quantities during a thermal-hydraulic transient. It models a test carried out on the Japanese mock-up “Large Scale Test Facility” (LSTF) in the framework of the OECD/ROSA-2 project, and which is representative of an Intermediate Break Loss Of Coolant Accident (IBLOCA) [29]. This mock-up represents a reduced scale Westinghouse PWR ($1/1$ ratio in height and $1/48$ in volume), with two loops instead of the four loops on the actual reactor and an electric powered heating core ($10$ MWe), see Figure 2. It operates at the same pressure and temperature values as the reference PWR. The simulated accidental transient is an IBLOCA with a break on the cold leg and no safety injection on the broken leg. The test under study reproduces a PWR $17\%$ (of cold leg cross-sectional area) cold leg IBLOCA transient with total failure of the auxiliary feedwater, single failure of diesel generators and three systems only available in the intact loop (high pressure injection, accumulator and low pressure injection).

CATHARE2 is used to simulate this integral effect test (see [30] for the full details of the modeling). During an IBLOCA, the reactor coolant system minimum mass inventory and the Peak Cladding Temperature (PCT) are obtained shortly after the beginning of the accumulators’ injection. Figure 3 shows the CATHARE2 prediction and the experimental values of the maximal cladding temperature (also called maximal heater rod temperature) obtained during the test. The conclusion of [30], which also presents other results, is that the CATHARE2 modeling of the LSTF allows to reproduce the global trends of the different physical phenomena during the transient of the experimental test. In our study, the main output variable of interest will be a single scalar which is the PCT during the accident transient (as an example, see the peak in Figure 3). This quantity is derived from the physical outputs provided by CATHARE2 code.

The input parameters of the CATHARE2 code correspond to various system parameters as boundary conditions, some critical flow rates, interfacial friction coefficients, condensation coefficients, heat transfer coefficients, etc. In our study, only uncertainties related to physical parameters are considered and no uncertainty on scenario variables (initial state of the reactor before the transient) is taken into account. All uncertain physical models identified in a IBLOCA transient of a nuclear power plant are supposed to apply to the LSTF, except phenomena related to fuel behavior because of the fuel absence in the LTSF. A physical model uncertainty consists in an additive or multiplicative coefficient associated to a physical model. Finally, $d=27$ scalar input variables of CATHARE2 code are considered uncertain and statistically independent of each other. They are then defined by their marginal probability density function (uniform, log-uniform, normal or log-normal). Table I gives more details about these uncertain inputs and their probability density functions (pdf). The nature of these uncertainties appears to be epistemic since they come from a lack of knowledge on the true value of these parameters.

Our objective with this use-case is to provide a good metamodel for sensitivity analysis, uncertainty propagation and, more generally, safety studies. Indeed, the cpu-time cost of one simulation is too important to directly perform all the statistical analysis which are required in a safety study and for which many simulations are needed. To overcome this limitation, an accurate metamodel, built from a reduced number of direct code calculations, will make it possible to develop a more complete and robust safety demonstration.

III.   Step 1: Design of experiments

This initial step of sampling is to define a design of $n$ experiments for the inputs and performing the corresponding runs with the numerical model $g$. The obtained sample of inputs/outputs will constitute the learning sample on which the metamodel will then be fitted. The objective is therefore to investigate, most efficiently and with few simulations, the whole variation domain of the uncertain parameters in order to build a predictive metamodel which approximates as accurately as possible the output $Y$.

For this, we propose to use a space-filling design (SFD) of a $n$ of experiments, this kind of design providing a full coverage of the high-dimensional input space [12]. Among SFD types, a Latin Hypercube Sample (LHS, [31])) with optimal space-filling and good projection properties [23] would be well adapted. In particular, [12, 32] have shown the importance of ensuring good low-order sub-projection properties. Maximum projection designs [33] or low-centered $L^{2}$ discrepancy LHS [34] are then particularly well-suited.

Mathematically, the experimental design corresponds to a $n$-size sample $\left\{\mathbf{x}^{(1)},\ldots,\mathbf{x}^{(n)}\right\}$ which is performed on the model (or code) $g$. This yields $n$ model output values denoted $\left\{y^{(1)},\ldots,y^{(n)}\right\}$ with $y^{(i)}=g(\mathbf{x}^{(i)})$. The obtained learning sample is denoted $\left(X_{s},Y_{s}\right)$ with $X_{s}=\left[{\mathbf{x}^{(1)}}^{T},\ldots,{\mathbf{x}^{(n)}}^{T}\right]^{T}$ and $Y_{s}=\left[y^{(1)},\ldots,y^{(n)}\right]^{T}$. Then, the goal is to build an approximating metamodel of $g$ from the $n$-sample $\left(X_{s},Y_{s}\right)$.

The number $n$ of simulations is a compromise between the CPU time required for each simulation and the number of input parameters. For uncertainty propagation and metamodel-building purpose, some rules of thumb propose to choose $n$ at least as large as $10$ times the dimension $d$ of the input vector [35, 22].

To build the metamodel for the IBLOCA test case, $n=500$ CATHARE2 simulations are performed following a space-filling LHS built in dimension $d=27$. The histogram of the obtained values for the output of interest, namely the PCT, is given by Figure 4 (temperature is in *∘*C). A kernel density estimator [36] of the data is also added on the plot to provide an estimator of the probability density function. A bimodality seems to be present in the histogram. It underlines the existence of bifurcation or threshold effects in the code, probably caused by a phenomenon of counter current flow limitation between the bend of hot legs and the steam generator inlet plena.

Remark III.1

Note that the input values are sampled following their prior distributions defined on their variation ranges. Indeed, as we are not ensured to be able to build a sufficiently accurate metamodel, we prefer sample the inputs following the probabilistic distributions in order to have at least a probabilized sample of the uncertain output, on which statistical characteristics could be estimated. Moreover, as explained in the next section, dependence measures can be directly estimated on this sample, providing first usable results of sensitivity analysis.

Finally, the bimodality that has been observed on the PCT distribution strengthens the use of advanced sensitivity indices (i.e. more general than linear ones or variance-based ones) in our subsequent analysis.

IV.   Step 2: Preliminary screening based on dependence measure

From the learning sample, a screening technique is performed in order to identify the primary influential inputs (PII) on the variability of model output. It has been recently shown that screening based on dependence measures [24, 25] or on derivative-based global sensitivity measures [37, 38] are very efficient methods which can be directly applied on a SFD. Moreover, beyond the screening job, these sensitivity indices can be quantitatively interpreted and used to order the PII by decreasing influence, paving the way for a sequential building of metamodel. In the considered IBLOCA test case, the adjoint model is not available and the derivatives of the model output are therefore not computed because of their costs. The screening step will then be based only on dependence measures, more precisely on Hilbert-Schmidt independence criterion (HSIC) indices, directly estimated from the learning sample.

The dependence measures for screening purpose have been proposed by [24] and [25]. These sensitivity indices are not the classical ones based on the decomposition of output variance (see [39] for a global review). They consider higher order information about the output behavior in order to provide more detailed information. Among them, the Hilbert-Schmidt independence criterion (HSIC) introduced by [40] builds upon kernel-based approaches for detecting dependence, and more particularly on cross-covariance operators in reproducing kernel Hilbert spaces, see Appendix A for mathematical details.

From the estimated $R^{2}_{\text{HSIC}}$ [24], independence tests can be performed for a screening purpose. The objective is to separate the inputs into two sub-groups, the significant ones and the non-significant ones. For a given input $X_{k}$, statistical HSIC-based tests aims at testing the null hypothesis “$\mathcal{H}_{0}^{(k)}$: $X_{k}$ and $Y$ are independent”, against its alternative “$\mathcal{H}_{1}^{(k)}$: $X_{k}$ and $Y$ are dependent”. The significance levelcccThe significance level of a statistical hypothesis test is the probability of rejecting the null hypothesis $\mathcal{H}_{0}$ when it is true. of these tests is hereinafter noted $\alpha$. Several HSIC-based statistical tests are available: asymptotic versions based on an approximation with a Gamma law (for large sample size), spectral extensions and permutation-based versions for non-asymptotic case (case of small sample size). All these tests are described and compared in [25]; a guidance to use them for a screening purpose is also proposed.

So, at the end of this HSIC-based screening step, the inputs are clustered into two subgroups, PII and non-influential inputs, and the PII are ordered by decreasing $R^{2}_{\text{HSIC}}$.This order will be used for the sequential metamodel building in step 3.

On the IBLOCA test case, from the learning sample of $n=500$ simulations, $R^{2}_{\text{HSIC}}$ dependence measures are estimated and bootstrap independence tests with $\alpha=0.1$ are performed. The independence hypothesis $\mathcal{H}_{0}$ is rejected for eleven inputs, which are now designated as PII (Primary Influential Inputs). These inputs are given by Table II with their estimated $R^{2}_{\text{HSIC}}$. Ordering them by decreasing $R^{2}_{\text{HSIC}}$ reveals:

• •

the large influence of the interfacial friction coefficient in the horizontal part of the hot legs (X10),

• •

followed by the minimum stable film temperature in the core X2, the interfacial friction coefficient in the SG inlet plena X12 and the wall to liquid friction (in under-saturated break flow conditions) in the break line X22,

• •

followed by seven parameters with a lower influence: the interfacial friction coefficients in the upper plenum X15, the downcomer X13, the core X14 and the SG outlet plena and crossover legs together X9, the heat transfer coefficient in the core for film boiling X5, the interfacial friction coefficient of the saturated break flow X26 and the condensation coefficient in the jet during the injection X27.

These results clearly underline the predominance influence of the uncertainties on various interfacial friction coefficients.

From the learning sample, some scatterplots of the PCT with respect to some well-chosen inputs (the three most influential ones: X2, X10 and X12) are displayed in Figure 5. An additional local regression using weighted linear least squares and a first degree polynomial model (moving average filter) is added on each scatterplot to extract a possible tendency. One can observe that larger values of the interfacial friction coefficient in the horizontal part of the hot legs (X10) lead to larger values of the PCT. This can be explained by the increase of vapor which brings the liquid in the horizontal part of hot legs, leading to a reduction of the liquid water return from the rising part of the U-tubes of the SG to the core (through the hot branches and the upper plenum). Since the amount of liquid water available to the core cooling is reduced, higher PCT are observed. In addition, we notice a threshold effect concerning this input: beyond a value of $2$, the water non-return effect seems to have been reached, and X10 no longer appears to be influential. We also note that the minimum stable film temperature in the core (X2) shows a trend: the more it increases, the lower the PCT. This is explained by the fact that in the film boiling regime in the core (i.e. when the rods are isolated from the liquid by a film of vapor), X2 represents (with a decrease in heat flux) the temperature from which the thermal transfer returns to the nucleate boiling regime. Thus, the larger X2, the faster the re-wetting of the rods, the faster the cladding temperature excursion is stopped, and thus the lower the PCT.

The same kind of physical analysis can be made for other PII by looking at their individual scatterplot. Finally, it is important to note that the estimated HSIC and the results of significant tests are relatively stable when the learning sample size varies from $n=300$ to $n=500$. Only two or three selected variables with a very low HSIC ($R^{2}_{\text{HSIC}}$ around 0.01) can differ. This confirms the robustness, with respect to the sample size, of the estimated HSIC and the results of the associated significance tests. Their relevance for qualitative ranking and screening purpose is emphasized.

In the next steps, only the eleven PII are considered as explanatory variables in the joint metamodel and will be successively included in the building process. The other sixteen variables will be joined in a so-called uncontrollable parameter.

V.   Step 3: Joint Gp metamodel with sequential building process

Among all the metamodel-based solutions (polynomials, splines, neural networks, etc.), we focus our attention on the Gaussian process (Gp) regression, which extends the kriging principles of geostatistics to computer experiments by considering the correlation between two responses of a computer code depending on the distance between input variables. The Gp-based metamodel presents some real advantages compared to other metamodels: exact interpolation property, simple analytical formulations of the predictor, availability of the mean squared error of the predictions and the proved capabilities for modeling of numerical simulators (see [41], [18] or [22]). The reader can refer to [42] for a detailed review on Gp metamodel.

However, for its application to complex industrial problems, developing a robust implementation methodology is required. Indeed, it often implies the estimation of several hyperparameters involved in the covariance function of the Gp (e.g. usual case of anisotropic stationary covariance function). Therefore, some difficulties can arise from the parameter estimation procedure (instability, high number of hyperparameters, see [22] for example). To tackle this issue, we propose a progressive estimation procedure based on the result of the previous screening step and using a joint Gp approach [27]. The interest of the previous screening step becomes twofold. First, the input space, on which each component of the joint Gp is built, can be reduced to the PII space (only the PII are explanatory inputs of Gp). Secondly, the joint Gp is built with a sequential process where the ranked PII are successively included as explanatory inputs in the metamodel. It is expected that these two uses of screening results could significantly make the joint Gp building easier and more efficient.

V.A.   Sequential Gp-building process based on successive inclusion of PII as explanatory variables

At the end of the screening step, the PII are ordered by decreasing influence (decreasing $R^{2}_{\text{HSIC}}$). They are successively included as explanatory inputs in the Gp metamodel while the other inputs (the remaining PII and the other non-PII inputs) are joined in a single macro-parameter which is considered as an uncontrollable parameter (i.e. a stochastic parameter, notion detailed in section V.B). Thus, at the $j^{th}$ iteration, a joint Gp metamodel is built with, as explanatory inputs, the $j$ first ordered PII. The definition and building procedure of a joint Gp is fully described in [27] and summarized in the next subsection.

However, building a Gp or a joint Gp involves to perform a numerical optimization in order to estimate all the parameters of the metamodel (covariance hyperparameters and variance parameter). As we usually consider in computer experiments anisotropic (stationary) covariance, the number of hyperparameters linearly increases with the number of inputs. In order to improve the robustness of the optimization process and deal with a large number of inputs, the estimated hyperparameters obtained at the $(j-1)^{th}$ iteration are used, as starting points for the optimization algorithm. This procedure is repeated until the inclusion of all the PII. Note that this sequential estimation process is directly adapted from the one proposed by [22].

V.B.   Joint Gp metamodeling

We propose to use a joint Gp metamodeling to handle the group of non-PII inputs and capture their residual effect. In the framework of stochastic computer codes, [28] proposed to model the mean and dispersion of the code output by two interlinked Generalized Linear Models (GLM), called “joint GLM”. This approach has been extended by [27] to several nonparametric models and the best results on several tests are obtained with two interlinked Gp models, called “joint Gp”. In this case, the stochastic input is considered as an uncontrollable parameter denoted $\mathbf{X}_{\varepsilon}$ (i.e. governed by a seed variable).

We extend this approach to a group of non-explanatory variables. More precisely, the input variables $\mathbf{X}=(X_{1},\ldots,X_{d})$ are divided in two subgroups: the explanatory ones denoted $\mathbf{X_{exp}}$ and the others denoted $\mathbf{X}_{\varepsilon}$. The output is thus defined by $Y=g(\mathbf{X_{exp}},\mathbf{X}_{\varepsilon})$ and the metamodeling process will now focus on fitting the random variable $Y|\mathbf{X_{exp}}$ddd$Y|\mathbf{X_{exp}}$ (i.e. $Y$ knowing $X_{exp}$) is a random variable as its value depends on the uncontrollable random variable $\mathbf{X}_{\varepsilon}$.. Under this hypothesis, the joint metamodeling approach yields building two metamodels, one for the mean $Y_{m}$ and another for the dispersion component $Y_{d}$:

[TABLE]

[TABLE]

where $\mathbb{E}[Z]$ is the usual notation for the expected (i.e. mean) value of a random variable $Z$.

To fit these mean and dispersion components, we propose to use the methodology proposed by [27]. To summarize, it consists in the following steps. First, an initial Gp denoted $Gp_{m,1}$ is estimated for the mean component with homoscedastic nugget effecteeeBorrowed from geostatistics to refer to an unexpected nugget of gold found in a mining process, a constant nugget effect assumes an additive white noise effect, whose variance constitutes the nugget parameter. Most often, this variance is assumed to be constant, independent from the inputs (here $X_{exp}$), and the nugget effect is called homoscedastic. When this variance depends on the value of $x$ (i.e. is a function of $X$), the nugget effect is called heteroscedastic.. A nugget effect is required to relax the interpolation property of the Gp metamodel. Indeed, this property, which would yield zero residuals for the whole learning sample, is no longer desirable as the output $Y|\mathbf{X_{exp}}$ is stochastic. Then, a second Gp, denoted $Gp_{v,1}$, is built for the dispersion component with, here also, an homoscedastic nugget effect. $Gp_{v,1}$ is fitted on the squared residuals from the predictor of $Gp_{m,1}$. Its predictor is considered as an estimator of the dispersion component. The predictor of $Gp_{v,1}$ provides an estimation of the dispersion at each point, which is considered as the value of the heteroscedastic nugget effect. The homoscedastic hypothesis is so removed and a new Gp, denoted $Gp_{m,2}$, is fitted on data, with the estimated heteroscedastic nugget. This nugget, as a function of $\mathbf{X_{exp}}$, accounts for a potential interaction between $\mathbf{X_{exp}}$ and the uncontrollable parameter $\mathbf{X}_{\varepsilon}$. Finally, the Gp on the dispersion component is updated from $Gp_{m,2}$ following the same methodology as for $Gp_{v,1}$.

Remark V.1

Note that some parametric choices are made for all the Gp metamodels: a constant trend and a Matérn stationary anisotropic covariance are chosen. All the hyperparameters (covariance parameters) and the nugget effect (when homoscedastic hypothesis is done) are estimated by maximum likelihood optimization process.

V.C.   Assessment of metamodel accuracy

To evaluate the accuracy of a metamodel, we use the predictivity coefficient $Q^{2}$:

[TABLE]

where $(x^{(i)})_{1\leq i\leq n_{\text{test}}}$ is a test sample, $(y^{(i)})_{1\leq i\leq n_{\text{test}}}$ are the corresponding observed outputs and $(\hat{y}^{(i)})_{1\leq i\leq n_{\text{test}}}$ are the metamodel predictions. $Q^{2}$ corresponds to the coefficient of determination in prediction and can be computed on a test sample independent from the learning sample or by cross-validation on the learning sample. The closer to one the $Q^{2}$, the better the accuracy of the metamodel. In the framework of joint Gp-modeling, $Q^{2}$ criterion will be used to assess the accuracy of the mean part of the joint Gp, namely $Gp_{m,\bullet}$, whether $Gp_{m,\bullet}$ is a homoscedastic ($Gp_{m,1}$) or heteroscedastic Gp ($Gp_{m,2}$). This quantitative information can be supplemented with a plot of predicted against observed values ($\hat{y}^{(i)}$ with respect to $y^{(i)}$) or a quantile-quantile plot.

To evaluate the quality of the dispersion part of a joint metamodel (denoted $Gp_{v,\bullet}$), we use the graphical tool introduced in [27] to assess the accuracy of the confidence intervals predicted by a Gp metamodel. For a given Gp metamodel, It consists in evaluating the proportions of observations that lie within the $\alpha$-theoretical confidence intervals which are built with the mean squared error of Gp predictions (the whole Gp structure is used and not only the conditional mean). These proportions (i.e. the observed confidence intervals) can be visualized against the $\alpha$-theoretical confidence intervals, for different values of $\alpha$.

V.D.   Application on IBLOCA test case

The joint Gp metamodel is built from the learning sample of $n=500$: the eleven PII identified at the end of the the screening step are considered as the explanatory variables while the sixteen others are considered as the uncontrollable parameter. Gps on mean and dispersion components are built using the sequential building process described in section V.A where PII ordered by decreasing $R^{2}_{\text{HSIC}}$ are successively included in Gp. $Q^{2}$ coefficient of mean component $Gp_{m}$ is computed by cross validation at each iteration of the sequential building process. The results which are given by Table III show an increasing predictivity until its stabilization around 0.87, which illustrates the robustness of the Gp building process. The first four PII make the major contribution yielding a $Q^{2}$ around 0.8, the four following ones yield minor improvements (increase of 0.02 on average for each input) while the three last PII does not improve the Gp predictivity. Note that, in this application, these results remain unchanged whether we consider homoscedastic Gp ($Gp_{m,1}$) or heteroscedastic Gp ($Gp_{m,2}$), the heteroscedastic nugget effect not significantly improving the Gp predictor for the mean component. Thus, only $13\%$ of the output variability remains here not explained by $Gp_{m}$, this includes both the inaccuracy of the $Gp_{m}$ (part of $Y_{m}$ not fitted by Gp) and the total effect of the uncontrollable parameter, i.e. the group of non-selected inputs.

For this exercise, 600 remaining CATHARE2 simulations, different from the learning sample, are also available. As they are not used to build the Gp metamodel, this set of simulations will constitutes a test sample. The $Q^{2}$ computed on this test sample is $Q^{2}=0.90$ for both $Gp_{m,1}$ and $Gp_{m,2}$, which is consistent with the estimations by cross-validation.

Now, to assess the quality of dispersion part of the joint Gp, the predicted confidence intervals are compared with the theoretical ones (cf. Section V.C). Figure 6 gives the results obtained by $Gp_{m,1}$ and $Gp_{m,2}$ on the learning sample (by cross-validation) and the test sample, since available here. It clearly illustrates both the interest of considering a heteroscedatic nugget effect and the efficiency of using a joint Gp model to fit and predict this nugget. It can be seen that the joint Gp yields the most accurate confidence intervals in prediction, especially for the test sample. Indeed, all its points are close to the theoretical $y=x$ line (a deviation is only observed for the learning sample for the highest $\alpha$), while the homoscedastic Gp tends to give too large confidence intervals. Thus, in this case, the heteroscedasticity hypothesis is justified and, consequently, the proposed joint Gp model is clearly more competitive than the simple one.

VI.   Step 4A: Variance-based sensitivity analysis

Sensitivity Analysis methods allow to answer the question “How do the input parameters variations contribute, qualitatively or quantitatively, to the variation of the output?” [43]. These tools can detect non-significant input parameters in a screening context, determine the most significant ones, measure their respective contributions to the output or identify an interaction between several inputs which impacts strongly the model output. From this, engineers can guide the characterization of the model by reducing the output uncertainty: for instance, they can calibrate the most influential inputs and fix the non-influential ones to nominal values. Many surveys on sensitivity analysis exist in the literature, such as [44], [45] or [46]. Sensitivity analysis can be divided into two sub-domains: the local sensitivity analysis and the global sensitivity analysis. The first one studies the effects of small input perturbations around nominal values on the model output [47] while the second one considers the impact of the input uncertainty on the output over the whole variation domain of uncertain inputs [43]. We focus here on one of the most widely used global sensitivity indices, namely Sobol’ indices which are based on output variance decomposition.

A classical approach in global sensitivity analysis is to compute the first-order and total Sobol’ indices which are based on the output variance decomposition [48, 49], see Appendix B for mathematical details. Sobol’ indices are widely used in global sensitivity analysis because they are easy to interpret and directly usable in a dimension reduction approach. However, their estimation (based on Monte-Carlo methods for example) requires a large number of model evaluations, which is intractable for time expensive computer codes. An classical alternative solution consists in using a metamodel to compute these indices. Note that a connection can be made between the estimation error of Sobol’ indices when using a metamodel and the predictivity coefficient $Q^{2}$ of the metamodel. Indeed, when the $Q^{2}$ is estimated on a probabilized sample of the inputs (in other words when it is computed under the probability distribution of the inputs), it provides an estimation of the part of variance unexplained by the metamodel. This can be kept in mind when interpreting the Sobol’ indices estimated with the metamodel.

VI.A.   Sobol’ indices with a joint Gp metamodel

In the case where a joint Gp metamodel is used to take into account an uncontrollable input $X_{\varepsilon}$, [50] and [27] have shown how to deduce Sobol’ indices from this joint metamodel, see Appendix C for mathematical details.

Therefore, from a joint Gp, it is only possible to estimate Sobol’ indices of any subset of $\mathbf{X_{exp}}$ (equation (12)) and the total Sobol’ index of $X_{\varepsilon}$ (equation (13)). The latter is interpreted as the total sensitivity index of the uncontrollable process. The individual index of $X_{\varepsilon}$ or any interaction index involving $X_{\varepsilon}$ are not directly reachable from joint Gp; their contributions in $S_{\varepsilon}^{T}$ can not be distinguished. This constitutes a limitation of this approach. However, the potential interactions between $X_{\varepsilon}$ and inputs of $\mathbf{X_{exp}}$ could be pointed out, considering all the primary and total effects of all the other parameters. The sensitivity analysis of $Y_{d}$ can also be a relevant indicator: if a subset $X_{u}$ of $\mathbf{X_{exp}}$ is not influential on $Y_{d}$, we can deduce that $S_{u\varepsilon}$ is equal to zero. Note that in practice, $\mbox{Var}(Y)$ which appears in both equations (12) and (13) can be estimated directly from the learning sample (empirical estimator) or from the fitted joint Gp, using equation (9).

VI.B.   Results on IBLOCA test case

From the joint Gp built in section V.D, Sobol’ indices of PII are estimated from $Gp_{m}$ by using equation (12), $\mbox{Var}(Y)$ being estimated with $Gp_{m}$ and $Gp_{d}$ using equation (9). For this, intensive Monte Carlo methods are used (see e.g. the pick-and-freeze estimator of [51]): tens of thousands simulations of Gp are done. Remind that predictions of Gp are very inexpensive (few seconds for several thousand simulations), especially with respect to the thermal-hydraulic simulator. The obtained first Sobol’ indices of PII are given by Table IV and represent $85\%$ of the total variance of the output. Qualitative results of HSIC indices are confirmed and refined: X10 remains the major influential input with $59\%$ of explained variance, followed to a lesser extend by X12 and X22 with for each of them $8\%$ of variance. The partial total Sobol’ indices involving only PII and derived from $Gp_{m}$ show that additional $4\%$ of variance is due to interaction between X10, X12 and X22. The related second order Sobol’ indices are estimated and the significant ones are given in Table IV. The other PII have negligible influence. In short, the set of PII explain a total $89\%$ of the output variance, of which $79\%$ is only due to X10, X12 and X22.

For the physical interpretation, these results confirm those revealed in Section IV, with a rigorous quantification of inputs’ importance: the interfacial friction coefficient in the horizontal part of the hot legs (X10) is the main contributor to the uncertainty of the PCT. Moreover, some results have not been revealed by the HSIC-based screening analysis of Table II. At present, Sobol’ indices clearly indicate that the interfacial friction coefficient in the SG inlet plena X12 and the wall to liquid friction (in under-saturated break flow conditions) in the break line X22 are more influential than the minimum stable film temperature in the core X2. X22 has a significant influence on the PCT because its low values lead to higher break flow rates, resulting in a loss of the primary water mass inventory at the higher break, and thus a more significant core uncovering (then higher PCT). For X12, its higher values lead to a greater supply (by the vapor) of liquid possibly stored in the water plena to the rest of the primary loops (then lower PCT). Table IV also shows that there are some small interaction effects (possibly antagonistic) between X10 and X12, as well as between X10 and X22. Let us remark that deepening this question (which is outside the scope of this paper) would be possible via plotting, from the Gp metamodel, the conditional expectations of the PCT as a function of the interaction variables.

At present, from $Gp_{d}$ and using equation (13), the total effect of the group of the sixteen non-PII inputs (i.e. $X_{\varepsilon}$) is estimated to $9.7\%$. This includes its effect alone and in interaction with the PII. To further investigate these interactions, Sobol indices of $Y_{d}$ are estimated and HSIC-based statistical dependence tests are applied on $Y_{d}$. They reveal that only X10, X14, X2, X22 and X4 potentially interact with the group of non-PII inputs. At this stage of our analysis, this result cannot be physically interpreted.

VII.   Step 4B: Quantile estimation

As already said in the introduction, once a predictive Gp metamodel has been built, it can be used to perform uncertainty propagation and in particular, estimate probabilities or, as here, quantiles.

VII.A.   Gp-based quantile estimation

The most trivial and intuitive approach to estimate a quantile with a Gp metamodel is to apply the quantile definition to the predictor of the metamodel. This direct approach yields a so called plug-in estimator. More precisely, with a Gp metamodel, this approach consists in using only the predictor of the Gp metamodel (i.e. the expectation conditional to the learning sample) in order to estimate the quantile. As the expectation of the Gp mean is a deterministic function of the inputs, this provides a deterministic expression of the quantile but no confidence intervals are available. Moreover, for high (resp. low) quantiles, this methods tends to substantially underestimate (resp. overestimate) the true quantile because the metamodel predictor is usually constructed by smoothing the computer model output values (see an analysis of this phenomenon in [14]).

To overcome this problem, [52] has proposed to take advantage of the probabilistic-type Gp metamodel by using its entire structure: not only the mean of the conditional Gp metamodel but also its covariance structure are taken into account. In this full-Gp based approach also called modular Bayesian approach, the quantile definition is therefore applied to the global Gp metamodel and leads to a random variable. The expectation of this random variable can be then considered as a quantile estimator. Its variance and, more generally, all its distribution can then be used as an indicator of the accuracy of the quantile estimate. Confidence intervals can be deduced, which constitutes a great advantage of this full-Gp based approach. Moreover, the efficiency of this approach for high quantile (of the order of $95\%$) has been illustrated by [53].

In practice, the estimation of quantile with the full-Gp approach is based on stochastic simulations (conditionally to the learning sample) of the Gp metamodel, by using the standard method of conditional simulations [54]. We recall just that, to generate conditional Gaussian simulations on a given interval, a discretization of the interval is first considered, then the vector of the conditional expectation and the matrix of the conditional covariance on the discretized interval are required. Note that [55] uses this approach with Gp conditional simulations for the estimation of Sobol’ indices and their associated confidence intervals.

In this paper, from our joint Gp model, we will compare the full-Gp approach applied to either the homoscedastic Gp or the heteroscedastic Gp and will proceed as follows:

• •

For the homoscedastic Gp, the standard technique of conditional simulations is applied to $Gp_{m,1}$ (built to estimate $Y_{m}$, with a constant nugget effect).

• •

For the heteroscedastic Gp, we propose a new technique to simulate the conditional Gp trajectories:

The heteroscedastic Gp built for $Y_{m}$ (namely $Gp_{m,2}$) provides the conditional expectation vector and a preliminary conditional covariance matrix. 2. 2.

The Gp built for $Y_{d}$ (namely $Gp_{d,2}$) is predicted and provides the heteroscedastic nugget effect which is added to the diagonal of the conditional covariance matrix of the previous step. 3. 3.

Conditional simulations are then done using the standard method [54].

VII.B.   Results on IBLOCA test case

In this section, we focus on the estimation of the $95\%$-quantile of the PCT for the IBLOCA application. From the learning sample of size $n=500$ and the joint Gp, we compare here the following approaches to estimate the PCT quantile:

• •

The classical empirical quantile estimator, denoted $\hat{q}_{95}^{\text{emp}}$. A bootstrap method (see for example [56]) makes it possible to obtain in addition a $90\%$-confidence interval for this empirical quantile.

• •

The plug-in estimators from the homoscedastic or the heteroscedastic Gp, denoted $\hat{q}_{95}^{Gp_{m,1}\text{-pi}}$ and $\hat{q}_{95}^{Gp_{m,2}\text{-pi}}$. No confidence intervals are obtained using this estimation method.

• •

The full-Gp estimators from the homoscedastic or the heteroscedastic Gp, denoted $\hat{q}_{95}^{Gp_{m,1}\text{-full}}$ and $\hat{q}_{95}^{Gp_{m,2}\text{-full}}$ respectively. As explained in the previous section, confidence intervals can be deduced with this full-Gp approach.

Table V synthesizes all the results obtained for the PCT $95\%$-quantile with these different approaches given above. In addition, $90\%$-confidence intervals are given when they are available.

As explained in Section V.D, we also have 600 other CATHARE2 simulations, for a total of 1100 PCT simulations (learning sample plus test sample). A reference value of the $95\%$-quantile is obtained from this full sample with the classical empirical quantile estimator:

[TABLE]

The empirical estimator based on the learning sample is imprecise but its bootstrap-based confidence interval is consistent with the reference value. As expected, the plug-in approaches provide poor estimations of the quantile, which is here strongly underestimated. The full-Gp approach based on the homoscedastic assumption overestimates the quantile; this is consistent with the analysis of Gp confidence intervals in Figure 6 (too large confidence intervals provided by homoscedastic Gp). Finally, the best result is obtained with the conditional simulation method based on the heteroscedastic Gp metamodel, built with the joint Gp method. This approach yields a more accurate prediction than the usual homoscedastic Gp and outperforms the empirical estimator in terms of confidence interval. Once again, the heteroscedasticity hypothesis is clearly relevant in this case.

VIII.   Conclusion

In the framework of the estimation of safety margins in nuclear accident analysis, it is essential to quantitatively assess the uncertainties tainting the results of Best-estimate codes. In this context, this paper has been focused on an advanced statistical methodology for Best Estimate Plus Uncertainty (BEPU) analysis, illustrated by a high dimensional thermal-hydraulic test case simulating accidental scenario in a Pressurized Water Reactor (IBLOCA test case). Some statistical analyses such as the estimation of high-level quantiles or quantitative sensitivity analysis (e.g., estimation of Sobol’ indices based on variance decomposition) may call in practice for several thousand of code simulations. Complex computer codes, as the ones used in thermal-hydraulic accident scenario simulations, are often too cpu-time expensive to be directly used to perform these studies.

To cope with this limitation, we propose a methodology mainly based on a predictive joint Gp metamodel, built with an efficient sequential algorithm. First, an initial screening step based on advanced dependence measures and associated statistical tests enabled to identify a group of significant inputs, allowing a reduction of the dimension. The efforts of optimization when fitting the metamodel can then be concentrated on the main influential inputs. The robustness of metamodeling is thus increased. Moreover, thanks to the joint metamodel approach, the non-selected inputs are not completely removed: the residual uncertainty due to dimension reduction is integrated in the metamodel and the global influence of non-selected inputs is so controlled. From this joint Gp metamodel, accurate uncertainty propagation and quantitative sensitivity analysis, not feasible with the numerical model because of its computational cost, become accessible. Thus, the uncertainties of model inputs are propagated inside the joint Gp to estimate Sobol’ indices, failure probabilities and/or quantiles, without additional code simulations.

Thus, on the IBLOCA application, a predictive Gp metamodel is built with only a few hundred of code simulations (500 code simulations for 27 uncertain inputs). From this joint Gp, a quantitative sensitivity analysis based on variance decomposition is performed without additional code simulation: Sobol’ indices are computed and reveal that the output is mainly explained by four uncertain inputs. One of them, namely the interfacial friction coefficient in the hot legs, is strongly influential with around $60\%$ of output variance explained, the three others being of minor influence. The quite less influence of all the other inputs is also confirmed. Note that a direct and accurate computation of Sobol’ indices with the thermal-hydraulic code would have required tens of thousands of code simulations.

The physical interpretation of the results obtained with the screening step and the variance-based sensitivity analysis step are known to be useful for modelers and engineers. This study has demonstrated this once again, in the particular case of an IBLOCA safety analysis, by revealing some physical effects on the PCT of influential inputs which cannot be understood without a global statistical approach (e.g. the threshold effect due to the interfacial friction coefficient in the horizontal part of the hot legs). [57] has also underlined the importance of sensitivity analysis in a validation methodology in order to identify relevant physical model uncertainties on which safety engineers must focus their efforts. Counter-intuitive effects are also present during an IBLOCA transient and only a global understanding of physical phenomena can help. As a perspective of the present work, extending the screening and sensitivity analysis tools to a joint analysis of several relevant output variables of interest (as the PCT time, the primary pressure and the core differential pressure) would be essential.

In the IBLOCA test case, we are particularly interested by the estimation of the $95\%$-quantile of the model output temperature. In nuclear safety, as in other engineering domains, methods of conservative computation of quantiles [6, 58] have been largely studied. We have shown in the present work how to use and simulate the joint Gp metamodel to reach the same objective: the uncertainty of the influential inputs are directly and accurately propagated through the mean component of the joint metamodel while a confidence bound is derived from the dispersion component in order to take into account the residual uncertainty of the other inputs. Results on the IBLOCA test case show that joint Gp provides a more accurate estimation of the $95\%$-quantile than the empirical approach, at equal calculation budget. Besides, this estimation is very close from the reference value obtained with 600 additional code simulations. Furthermore, the interest of the heteroscedastic approach in joint Gp is also emphasized: the estimated quantile and the associated confidence interval are much better than those of the homoscedastic approach.

As a future application of modern statistical methods on IBLOCA safety issues, one should mention the use of Gp metamodels to identify penalizing thermal-hydraulic transients, with respect to some particular scenario inputs. As in the present work, strong difficulties are raised by the cpu time cost of the computer code and the large number of uncertain (and uncontrollable) inputs.

Acknowledgments

We are grateful to Henri Geiser and Thibault Delage who performed the computations of the CATHARE2 code.

Appendix A Appendix: HSIC dependence measures

If we consider two reproducing kernel Hilbert spaces $\mathcal{F}_{k}$ and $\mathcal{G}$ of functions $\mathcal{X}_{k}\rightarrow\mathbb{R}$ and $\mathcal{Y}\rightarrow\mathbb{R}$ respectively, the crossed-covariance $C_{X_{k},Y}$ operator associated to the joint probabilistic distribution of $\left(X_{k},Y\right)$ is the linear operator defined for every $f_{X_{k}}\in\mathcal{F}_{k}$ and $g_{Y}\in\mathcal{G}$ by:

[TABLE]

$C_{X_{k},Y}$ generalizes the covariance matrix by representing higher order correlations between $X_{k}$ and $Y$ through nonlinear kernels. The HSIC criterion is then defined by the Hilbert-Schmidt norm of the cross-covariance operator:

[TABLE]

From this, [24] introduces a normalized version of the HSIC which provides a sensitivity index of $X_{k}$, lying in $[0,1]$:

[TABLE]

The closer to one the $R^{2}_{\text{HSIC},k}$, the stronger the dependence between $X_{k}$ and $Y$. In practice, [40] propose a Monte Carlo estimator of $\text{HSIC}(X_{k},Y)$ and a plug-in estimator can be deduced for $R^{2}_{\text{HSIC},k}$. Note that Gaussian kernel functions with empirical estimations of the variance parameter are used in our application (see [40] for more details).

Appendix B Appendix: Sobol’ indices

If $X=(X_{1},\ldots,X_{d})$ with independence between the variables $X_{1},\ldots,X_{d}$ and if $\mathbb{E}[g^{2}(X)]<+\infty$, we can apply the Hoeffding decomposition to the random variable $g(X)$ [59]:

[TABLE]

where $g_{\emptyset}=\mathbb{E}[g(X)]$, $g_{i}(X_{i})=\mathbb{E}[g(X)|X_{i}]-g_{\emptyset}$ and $g_{u}(X_{u})=\mathbb{E}[g(X)|X_{u}]-\sum_{v\subsetneq u}g_{v}(X_{v})$, with $X_{u}=(X_{i})_{i\in u}$, for all $u\subset\{1,\ldots,d\}$. All the $2^{d}$ terms in (8) have zero mean and are mutually uncorrelated with each other. This decomposition is unique and leads to the Sobol’ indices. These are the elements of the $g(X)$ variance decomposition according to the different groups of input parameter interactions in (8). More precisely, for each $u\subset\{1,\ldots,d\}$, the first-order and total Sobol’ indices of $X_{u}$ are defined by:

[TABLE]

$S_{u}$ represents the part of the output variance explained by $X_{u}$, independently from the other inputs, and $S_{u}^{T}$ is the part of the output variance explained by $X_{u}$ considered separately and in interaction with the other input parameters.

In practice, we are usually interested in the first-order sensitivity indices $S_{1},\ldots,S_{d}$, the total ones $S_{1}^{T},\ldots,S_{d}^{T}$ and sometimes in the second-order ones $S_{ij}$, $1\leq i<j\leq d$. The model $g$ is devoid of interactions if $\sum_{i=1}^{d}S_{i}\approx 1$.

Appendix C Appendix: Sobol’ indices from a joint metamodel

By adopting the same notations as in Section V.B where $\mathbf{X_{exp}}$ is the vector of all the controllable (or explanatory) inputs and $X_{\varepsilon}$ denotes the uncontrollable input (group of non explanatory inputs), we can show that the variance of the output variable $Y(\mathbf{X_{exp}},X_{\varepsilon})$ can be rewritten and deduced from the two components $Y_{m}$ and $Y_{d}$ (equations (2) and (3)):

[TABLE]

where $\mathbb{E}_{X}$ (resp. $\mbox{Var}_{X}$) denotes the mean (resp. variance) operator with respect to the probability density function of $X$. Furthermore, the variance of $Y$ is the sum of the contributions of both all the explanatory inputs in $\mathbf{X_{exp}}$ and $X_{\varepsilon}$:

[TABLE]

where $V_{\varepsilon}(Y)=\mbox{Var}_{X_{\varepsilon}}[\mathbb{E}_{\mathbf{X_{exp}}}(Y|X_{\varepsilon})]$, $V_{u}(Y)=\mbox{Var}_{\mathbf{X_{u}}}[\mathbb{E}_{\mathbf{X_{-u}}}(Y|X_{u})]-\sum_{v\subsetneq u}V_{v}(Y)$, $V_{u\varepsilon}(Y)=\mbox{Var}_{X_{u}X_{\varepsilon}}[\mathbb{E}_{\mathbf{X_{exp\;-u}}}(Y|X_{u}X_{\varepsilon})]-V_{u}(Y)-V_{\varepsilon}(Y)$.

Variance of the mean component $Y_{m}(\mathbf{X_{exp}})$ denoted hereafter $Y_{m}$ can be also decomposed:

[TABLE]

As $V_{u}(Y_{m})=\mbox{Var}_{\mathbf{X_{u}}}\mathbb{E}_{\mathbf{X_{exp\,-u}}}[\mathbb{E}_{X_{\varepsilon}}(Y|\mathbf{X_{exp}})|X_{u}]=V_{u}(Y)$, Sobol’ indices according to any subset of input variables of $\mathbf{X_{exp}}$ can be derived and estimated from $Y_{m}$:

[TABLE]

Similarly, the total sensitivity index of $X_{\varepsilon}$ is given by:

[TABLE]

Note that, as $Y_{d}(\mathbf{X_{exp}})$ is a positive random variable, the positivity of $S_{\varepsilon}^{T}$ is guaranteed.