Optimal Uncertainty Quantification of a risk measurement from a thermal-hydraulic code using Canonical Moments

TL;DR

This paper introduces a method to robustly quantify risk in nuclear safety computer simulations by optimizing the maximum quantile over distributions constrained by moments, using canonical moments for efficient computation.

Contribution

It develops a novel approach employing canonical moments to optimize risk quantification under moment constraints, enhancing robustness in uncertainty analysis.

Findings

Effective optimization of risk measures under moment constraints.

Use of canonical moments simplifies complex constrained problems.

Improved robustness in nuclear safety risk assessment.

Abstract

We study an industrial computer code related to nuclear safety. A major topic of interest is to assess the uncertainties tainting the results of a computer simulation. In this work we gain robustness on the quantification of a risk measurement by accounting for all sources of uncertainties tainting the inputs of a computer code. To that extent, we evaluate the maximum quantile over a class of distributions defined only by constraints on their moments. Two options are available when dealing with such complex optimization problems: one can either optimize under constraints; or preferably, one should reformulate the objective function. We identify a well suited parameterization to compute the optimal quantile based on the theory of canonical moments. It allows an effective, free of constraints, optimization.

Figures

Tables

Table 1. Table 1: Initial distribution of the 4 inputs of the hydraulic model.

Variable	Description	Distribution
$Q$	annual maximum flow rate	$Gumbel\left(mode=1013,scale=558\right)$
$K_s$	Manning-Strickler coefficient	$𝒩\left(\overline{x}=30,\sigma =7.5\right)$
$Z_v$	Depth measure of the river downstream	$𝒰(49,51)$
$Z_m$	Depth measure of the river upstream	$𝒰(54,55)$

Table 2. Table 2: Corresponding moment constraints of the 4 inputs of the hydraulic model.

Variable	Bounds	Mean	Second order moment	Third order moment
$Q$	$[160,3580]$	$1320.42$	$2.1632\times {10}^6$	$4.18\times {10}^9$
$K_s$	$[12.55,47.45]$	$30$	$949$	$31422$
$Z_v$	$[49,51]$	$50$	$2500$	$125050$
$Z_m$	$[54,55]$	$54.5$	$2970$	$161892$

Table 3. Table 3: Corresponding moment constraints of the 9 most influential inputs of the CATHARE model [ 2 ] .

Variable	Bounds	Mean	Second moment	Initial distribution (truncated on bounds)
$n^{\circ }10$	$[0.1,10]$	$1.33$	$3.02$	$LogNormal(0,0.76)$
$n^{\circ }22$	$[0,12.8]$	$6.4$	$45.39$	$Normal(6.4,4.27)$
$n^{\circ }25$	$[11.1,16.57]$	$13.83$	$192.22$	$Normal(13.79$
$n^{\circ }2$	$[-44.9,63.5]$	$9.3$	$1065$	$Uniform(-44.9,63.5)$
$n^{\circ }12$	$[0.1,10]$	$1.33$	$3.02$	$LogNormal(0,0.76)$
$n^{\circ }9$	$[0.1,10]$	$1.33$	$3.02$	$LogNormal(0,0.76)$
$n^{\circ }14$	$[0.235,3.45]$	$0.99$	$1.19$	$LogNormal(-0.1,0.45)$
$n^{\circ }15$	$[0.1,3]$	$0.64$	$0.55$	$LogNormal(-0.6,0.57)$
$n^{\circ }13$	$[0.1,10]$	$1.33$	$3.02$	$LogNormal(0,0.76)$

Table 4. Table 4: Results for the 95 % percent \% -quantile estimates.

	Empirical	Plug-in	Full-Gp	OUQ
Mean	746.80	735.83	741.46	788
90$\%$-CI	$[736.7,747.41]$		$[738.76,744.17]$

Equations

$$\mathcal{A}$$

$${\mathcal{A}}_{\Delta }$$

$${\overline{Q}}_{\mathcal{A}}(p)$$

$${\overline{Q}}_{\mathcal{A}}(p)=\inf \left\{h\in \mathbb{R}|\underset{\mu \in \mathcal{A}}{\inf }{F}_{\mu }(h)\ge p\right\}.$$

$$\mathcal{A}:=\left\{(G,\mu)\ \begin{array}[]{|l}G:\mathcal{X}\rightarrow\mathcal{Y}\text{, is a real valued measurable function},\\ \mu=\mu_{1}\otimes\dots\otimes\mu_{p}\in\bigotimes_{i=1}^{d}\mathcal{P}(\mathcal{X}_{i})\ ,\\ \text{for each G, and for some measurable functions }\\ \varphi_{l}:\mathcal{X}\rightarrow\mathbb{R}\text{ and }\varphi_{i}^{(j)}:\mathcal{X}_{i}\rightarrow\mathbb{R}\ ,\\ \qquad\qquad\bullet\ \mathbb{E}_{\mu}[\varphi_{l}]\leq 0\text{ for }l=1,\dots,N_{0}\ ,\\ \qquad\qquad\bullet\ \mathbb{E}_{\mu_{i}}[\varphi_{i}^{(j)}]\leq 0\text{ for }j=1,\dots,N_{i}\text{ and }i=1,\dots,d\end{array}\right\}$$

$${\mathcal{A}}_{\Delta }:=\left\{(G,\mu )\in \mathcal{A}|{\mu }_i\in {\Delta }_{N_0+N_i}({\mathcal{X}}_i)\right\}.$$

$$\underset{(G,\mu )\in \mathcal{A}}{\sup }{\mathbb{E}}_{\mu }[q(X,G(X))]=\underset{(G,\mu )\in {\mathcal{A}}_{\Delta }}{\sup }{\mathbb{E}}_{\mu }[q(X,G(X))].$$

$$\displaystyle\inf_{\mu\in\mathcal{A}}F_{\mu}(h)\$$

$$c_{n+1}^{+}=\max \left\{c\in \mathbb{R}:(c_1,\dots ,c_n,c)\in M_{n+1}\right\},$$

$$c_{n+1}^{-}=\min \left\{c\in \mathbb{R}:(c_1,\dots ,c_n,c)\in M_{n+1}\right\},$$

$$p_n=p_n(\mathbf{c})=\frac{c_n-c_n^{-}}{c_n^{+}-c_n^{-}}.$$

$${c^{\prime }}_i^{(j)}=\frac{1}{(u_i-l_i{)}^j}\sum _{k=0}^j\left(\genfrac{}{}{0ex}{}{j}{k}\right)(-l_i{)}^{j-k}{c}_i^{(k)},\text{ for }j=1,\dots ,N_i.$$

$$p_1=c_1$$

$$S(z)=S(z,\mu )={\int }_a^b\frac{d\mu (x)}{z-x},(z\in \mathbb{C}\backslash \{\text{supp}(\mu )\}).$$

$$S(z)={\int }_a^b\frac{d\mu (x)}{z-x}=\sum _{i=1}^n\frac{{\omega }_i}{z-x_i},$$

$$S(z)=\frac{Q^{(n-1)}(z)}{P_{\ast }^{(n)}(z)},$$

$${\omega }_i=\frac{Q^{(n-1)}(x_i)}{\frac{d}{dx}{P}_{\ast }^{(n)}(x){|}_{x=x_i}}.$$

$$S(z)$$

$$\displaystyle=\ \begin{array}[]{c}\lx@intercol\hfil 1\hfil\lx@intercol\vrule\lx@intercol\\ \hline\cr\vrule\lx@intercol\hfil z-a-\zeta_{1}(b-a)\hfil\lx@intercol\end{array}\ \ -\ \ \begin{array}[]{c}\lx@intercol\hfil\zeta_{1}\zeta_{2}(b-a)^{2}\hfil\lx@intercol\vrule\lx@intercol\\ \hline\cr\vrule\lx@intercol\hfil z-a-(\zeta_{2}+\zeta_{3})(b-a)\hfil\lx@intercol\end{array}\$$

$$\displaystyle-\ \ \begin{array}[]{c}\lx@intercol\hfil\zeta_{3}\zeta_{4}(b-a)^{2}\hfil\lx@intercol\vrule\lx@intercol\\ \hline\cr\vrule\lx@intercol\hfil z-a-(\zeta_{4}+\zeta_{5})(b-a)\hfil\lx@intercol\end{array}\ \ -\ \dots$$

$$P_{\ast }^{(k+1)}(x)=(x-a-(b-a)({\zeta }_{2k}+{\zeta }_{2k+1}))P_{\ast }^{(k)}(x)-(b-a{)}^2{\zeta }_{2k-1}{\zeta }_{2k}{P}_{\ast }^{(k-1)}(x),$$

$${\mathcal{A}}_{\Delta }=\left\{\mu =\sum _{i=1}^{n+1}{\omega }_i{\delta }_{x_i}\in \mathcal{P}([a,b])|\text{mathbbm}{E}_{\mu }(x^j)=c_j,j=1,\dots ,n\right\}.$$

$$\Gamma =\left\{(p_{n+1},\dots ,p_{2n+1})\in [0,1{]}^{n+1}|p_i\in \{0,1\}\Rightarrow p_k=0,k>i\right\}.$$

$$\left\{\begin{array}[]{lclcll}\omega_{i}^{(1)}&+&\dots&+&\omega_{i}^{(N_{i}+1)}&=1\\ \omega_{i}^{(1)}x_{i}^{(1)}&+&\dots&+&\omega_{i}^{(N_{i}+1)}x_{i}^{(N_{i}+1)}&=c_{i}^{(1)}\\ \lx@intercol\hfil\vdots\hfil\lx@intercol&&&&\lx@intercol\hfil\vdots\hfil\lx@intercol&\lx@intercol\hfil\vdots\hfil\lx@intercol\\ \omega_{i}^{(1)}{x_{i}^{(1)}}^{N_{i}}&+&\dots&+&\omega_{i}^{(N_{i}+1)}{x_{i}^{(N_{i}+1)}}^{N_{i}}&=c_{i}^{(N_{i})}\end{array}\right.$$

$${\mathbb{P}}_{\mu }(G(X)\le h)$$

$$\mathcal{A}=\left\{\mu =\otimes {\mu }_i\in \underset{i=1}{\overset{d}{⨂}}\mathcal{P}([l_i,u_i])|{\alpha }_i^{(j)}\le {\mathbb{E}}_{{\mu }_i}[x^j]\le {\beta }_i^{(j)},1\le j\le N_i\right\}.$$

$$H={\left(\frac{Q}{300K_s\sqrt{\frac{Z_m-Z_v}{5000}}}\right)}^{3/5}.$$

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

$$\text{for all }h\ge b;\underset{\mu \in \mathcal{A}}{\inf }{F}_{\mu }(h)\ge p,$$

$$\iff$$

$$\Rightarrow$$

$$\text{for all }h\ge a;\text{ for all }\mu \in \mathcal{A};F_{\mu }(h)\ge p,$$

Full text

††volume: Volume x, Issue x, 2019

1 Introduction

Computer codes are increasingly used to measure safety margins in nuclear accident management analysis instead of conservative procedures [1]. 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 physic and the model input parameters. The Best Estimate Plus Uncertainty (BEPU) methods [2] were developed in safety analyses, especially for the large break loss of coolant accident (see [3], [4]). Its principles rely mainly on a probabilistic modeling of the model input uncertainties, on Monte Carlo sampling for running the thermal-hydraulic computer code on sets of input, and on the application of statistical tools to infer high quantiles of the scalar output variables of interest [5] .

This takes place in a more general setting, known as Uncertainty Quantification (UQ) methods [6]. Quantitative assessment of the uncertainties tainting the results of computer simulations is a major topic of interest in both industrial and scientific communities. In the context of nuclear safety, the computer models are expensive to run. Uncertainty propagation, risk measurement such as high quantile inference, or system robustness analysis become a difficult task using such models. In order to circumvent this problem, a widely accepted method consists in replacing the cpu time expensive numerical simulations by inexpensive mathematical functions called metamodels. This metamodel is build from a set of computer code simulations that must be as representative as possible of the code in the variation domain of its uncertain inputs. Generally, space-filling designs of experiments are created with a given budget of $n$ code evaluations, that provide a full coverage of the input space [7]. In the presence of a high number of input parameters, screening strategies are then performed in order to identify the Primary Influential Inputs (PII) on the model output variability and rank them by decreasing influence. From the learning sample, a metamodel is therefore built to fit the simulator output, considering only the PII as the explanatory inputs, while the remaining inputs remain fixed to a default value. Among all the metamodel-based solutions (polynomials, splines, neural networks, etc.), Gaussian process metamodeling, also known as Kriging [8], has been very attractive. It makes the assumption that the response is a realization of a Gaussian process, conditioned on code observations. This approach provides the basis for statistical inference. In that, we dispose of simple analytical formulations of the predictor and the mean squared error of the predictions. The metamodel is then validated before being used.

Once the predictive metamodel has been built, it can be used to perform uncertainty propagation and in particular, estimate probabilities or, as here, quantiles (see [9], [10]). This measure of risk will be designated from now on as our quantity of interest. The most trivial approach to estimate a quantile with a Gp metamodel, called plug-in approach, is to apply the quantile definition to the predictor of the metamodel. As the expectation of the Gp mean is a deterministic function of the input, this provides a deterministic expression of the quantile but no confidence intervals. Moreover for high quantiles, this methods tends to substantially underestimate the true quantile [9]. To assess this problem, [11] has proposed to take into account the covariance structure of the Gp metamodel. The quantile definition is therefore applied to the global Gp metamodel and yields a random variable, whose expectation can be considered as the quantile estimator and its variance an indicator of the accuracy of its prediction. This full-Gp approach leads to confidence intervals. In practice, the estimation of a quantile with the full Gp approach is based on stochastic simulations (conditional on the learning sample) of the Gp metamodel.

This overall methodology yields the estimation of the $p$-quantile of the model output. In nuclear safety, as in other engineering domains, methods of conservative computation of quantiles have been largely studied [12]. Though the above construction work largely increases the robustness of the metamodel, the evaluation of the quantile remains tainted by the uncertainty of the input distributions. For simplicity, the inputs probability densities are usually chosen in parametric families (uniform, normal, log-normal, etc. See for instance Table 3). The distribution’s parameters are itself set with the available information coming from data and/or an expert opinion. This information is often reduced to an input’s mean value or a variance. Nevertheless, the distribution choice differs inevitably from reality. This uncertainty on the input probability densities is propagated to the quantile, hence, different choices of distributions will lead to different quantile values, thus different safety margins.

In this work, we propose to gain robustness on the quantification of this measure of risk. We aim to account for the uncertainty on the input distributions by evaluating the maximum quantile over a class of probability measures $\mathcal{A}$. In this optimization problem, the set $\mathcal{A}$ must be large enough to effectively represent our uncertainty on the inputs, but not too large in order to keep the estimation of the quantile representative of the physical phenomena. For example, the maximum quantile over the very large class $\mathcal{A}=\left\{\text{all distributions}\right\}$, proposed in [13], will certainly be too conservative to remain physically meaningful. Several articles which discuss possible choices of classes of distributions can be found in the literature of Bayesian robustness (see [14]). [15], consider a class of measures specified by a type of upper and lower envelope on their density. [16] study the class of unimodal distributions. In more recent work, [17] propose to optimize the measure of risk over a class of distributions specified by constraints on their generalized moments. Their work (and ours) is called Optimal Uncertainty Quantification (OUQ), because given a set of assumptions and information, there exist optimal bounds on uncertainties. In practical engineering cases, the available information on an input distribution is often reduced to the knowledge of its mean and/or variance. This is why in this paper, we are interested in a specific case of the framework introduced by [17]. We consider the class of measures known by some of their classical moments, which we refer to as the moment class:

[TABLE]

where $\mathcal{P}([l_{i},u_{i}])$ denotes the set of scalar probability measure on the interval $[l_{i},u_{i}]$. This set traduces simply that each random input is restricted to a given range and has some of its moments fixed. The tensorial product of measure sets traduces the mutual independence of the $d$ inputs. The choice of range enforcement is not a very strong hypothesis, as input variables represent physical parameters which are rarely unbounded.

The solution to our optimization problem is numerically computed thanks to the OUQ reduction theorem ([17], [18]). This theorem states that the measure corresponding to the extremal CDF (can be extended to quantile through Proposition 2.1), is located on the extreme points of the distribution set. In the context of the moment class, the extreme distributions are located on the $d$-fold product of finite convex combinations of Dirac masses:

[TABLE]

To be more specific it holds that when $N$ pieces of information are available on the moments of a scalar measure $\mu$, it is enough to pretend that the measure is supported on at most $N+1$ points. This powerful theorem gives the basis for practical optimization of our optimal quantity of interest. In this matter, Semi-Definite-Programming [19] has been already already explored by [20] and [21], but the deterministic solver used rapidly reaches its limitation as the dimension of the problem increases. One can also find in the literature a Python toolbox developed by [22] called Mystic framework that fully integrates the OUQ framework. However, it was built as a generic tool for generalized moment problems and the enforcement of the moment constraints is not optimal. By restricting the work to our moment class, we propose an original and practical approach based on the theory of canonical moments [23]. Canonical moments of a measure can be seen as the relative position of its moment sequence in the moment space. It is inherent to the measure and therefore presents many interesting properties. It allows to explore very efficiently the optimization space $\mathcal{A}_{\Delta}$, were the maximum quantile is to be found. Hence, we rewrite the optimization problem on the highly constrained set $\mathcal{A}_{\Delta}$ into a simplified and constraints free optimization problem.

The paper proceeds as follows. Section 2 describes the OUQ framework and the OUQ reduction theorem. In Section 3, we then describe step by step the algorithm calculating our quantity of interest with the canonical moments parameterization. We present in Section 4, an extended algorithm to deal with inequality constraints on the moments. Section 5 and 6 are dedicated to the application of our algorithm to a toy example, and to the peak cladding temperature for the IBLOCA application presented in the introduction. Section 7 gives some conclusions and perspectives.

2 OUQ principles

2.1 Duality transformation

In this work, we consider the quantile of the output of a computer code $G:\mathbb{R}^{d}\rightarrow\mathbb{R}$, seen as a black box function. As we said, in order to gain robustness on the risk measurement, our goal is to find the maximum quantile over the moment class $\mathcal{A}$ described in Equation (1). The objective value writes:

[TABLE]

The objective value written as a quantile (2.1) is not very convenient to work with. In order to applied the OUQ reduction Theorem 2.1 [17], one must optimize an affine functional of the measure. In particular, it is necessary to optimize a probability instead of a quantile. The following result, illustrated in Figure 1, can be interpreted as a duality transformation of our optimization problem (2.1), into the optimization of a probability of failure (p.o.f). The proof is postponed to A.

Proposition 2.1**.**

The following duality result holds

[TABLE]

Proposition 2.1 reads like this: the optimal quantile over a class of distributions is equal to the quantile of the CDF lower envelope. Our problem is therefore to evaluate the lowest probability of failure $\inf_{\mu\in\mathcal{A}}\mathbb{P}_{\mu}(G(X)\leq h)$ for a given threshold $h$.

2.2 Reduction Theorem

Under the form of Proposition 2.1, the OUQ reduction theorem applies (see [17], [18]). It states that the optimal solution of the p.o.f optimization is a product of discrete measures. A general form of the theorem reads as follows:

Theorem 2.1** (OUQ reduction [17, p.37]).**

Suppose that $\mathcal{X}:=\mathcal{X}_{1}\times\dots\times\mathcal{X}_{d}$ is a product of Radon spaces. Let

[TABLE]

Let $\Delta_{n}(\mathcal{X})$ be the set of all discrete measure supported on at most $n+1$ points of $\mathcal{X}$, and

[TABLE]

Let $q$ be a measurable real function on $\mathcal{X}\times\mathcal{Y}$. Then

[TABLE]

This theorem derives from the work of [18], who has shown that the extreme measures of a moment class $\{\mu\in\mathcal{P}(\mathcal{X})\ |\ \mathbb{E}_{\mu}[\varphi_{1}]\leq 0,\allowbreak\dots,\mathbb{E}_{\mu}[\varphi_{n}]\allowbreak\leq 0\}$ are the discrete measures that are supported on at most $n+1$ points. The strength of Theorem 2.1 is that it extends the result to a tensorial product of moment sets. The proof relies on a recursive argument using Winkler’s classification on every set $\mathcal{X}_{i}$. A remarkable fact is that, as long as the quantity to be optimized is an affine function of the underlying measure $\mu$, this theorem remains true whatever the function $G$ and the quantity of interest $q$ are. Applying Theorem 2.1 to the optimization of the probability of failure, it is rewritten as:

[TABLE]

3 Parameterization using canonical moments

The optimization problem in Equation (4) shows that the weights and positions of the input distributions provide a natural parameterization for the computation of the p.o.f. However, in order to compute the lowest p.o.f, one must be able to explore the whole set of admissible measures $\mathcal{A}_{\Delta}$. Two ways to handle the problem appear. The first one consists in optimizing the objective value $F_{\mu}(h)$ under constraints, that is $\mu\in\mathcal{A}_{\Delta}$: this is the approach taken by [22] with the Mystic Framework. The second option, always favored when possible, consists in reformulating the objective function. This requires to identify a new parameterization adapted to the problem. Here, canonical moments [23] provide a surprisingly well tailored reparameterization.

The work on canonical moments was first introduced by [24]. His main contribution covered the original study of the geometric aspect of general moment space [25], [26]. In a number of further papers, Skibinsky proves numerous other interesting properties of the canonical moments. [23] have shown the intrinsic relation between a measure $\mu$ and its canonical moments. They highlight the interest of canonical moments in many areas of statistics, probability and analysis such as problem of design of experiments, or the Hausdorff moment problem [27]. In the following, we describe step by step the algorithm used to transform the optimization problem of Equation 4 under the canonical moments parameterization.

3.1 Step 1. From classical moments to canonical moments

We enforce some moments on the input distributions of the code $G$. In this section, we present how to transform these classical moment constraints, into canonical moments constraints.

We define the moment space $M:=M(a,b)=\{\mathbf{c}(\mu)\ |\ \mu\in\mathcal{P}([a,b])\}$ where $\mathbf{c}(\mu)$ denote the sequence of all moments of some measure $\mu$. The $n$th moment space $M_{n}$ is defined by projecting $M$ onto its first $n$ coordinates, $M_{n}=\{\mathbf{c}_{n}(\mu)=(c_{1},\dots,c_{n})\;|\;\mu\in\mathcal{P}([a,b])\}$. $M_{2}$ is depicted in Figure 2. We first define the extreme values,

[TABLE]

which represent the maximum and minimum values of the $(n+1)$th moment that a measure can have, when its moments up to order $n$ equal to $\mathbf{c}_{n}$ . The $n$th canonical moment is then defined recursively as

[TABLE]

Note that the canonical moments are defined up to the degree $N=N(\mathbf{c})=\min\left\{\right.n\in\mathbb{N}\ |\ \mathbf{c}_{n}\in\partial M_{n}\left.\right\}$, and $p_{N}$ is either [math] or $1$. Indeed, we know from [23, Theorem 1.2.5] that $\mathbf{c}_{n}\in\partial M_{n}$ implies that the underlying $\mu$ is uniquely determined, so that, $c_{n}^{+}=c_{n}^{-}$. We also introduce the quantity $\zeta_{n}=(1-p_{n-1})p_{n}$ that will be of some importance in the following. The very nice properties of canonical moments are that they belong to $[0,1]$ and are invariant by any affine transformation of the support of the underlying measures. Hence, we may restrict ourselves to the case $a=0$, $b=1$.

Therefore, for every $i=1,\dots,d$, the support of the measure $\mu_{i}$ is transformed into $[0,1]$ using the affine transformation $y=l_{i}+(u_{i}-l_{i})x$. The sequences of moments of the corresponding measures are written $\mathbf{c}^{\prime}_{i}=({c^{\prime}}_{i}^{(1)},\dots,{c^{\prime}}_{i}^{(N_{i})})$ where ${c^{\prime}}_{i}^{(j)}$ reads

[TABLE]

Given a sequence of moment constraints $(c_{i}^{(j)})_{1\leq j\leq N_{i}}$ enforced to the $i$th input, it is then possible to calculate the corresponding sequence of canonical moments $\mathbf{p}_{i}=(p_{i}^{(1)},\dots,p_{i}^{(N_{i})})$. [23, p. 29] propose a recursive algorithm named Q-D algorithm that allows this computation. It drastically fastens the computational time compared to the raw formula that consists of computing Hankel determinants [23, p. 32]. In practical applications, we wish to enforce low order of moments, typically order 2 (see for instance Table 1). In this case we dispose of the simple analytical formulas

[TABLE]

One can easily see that enforcing $N_{i}$ moments or $N_{i}$ canonical moments to the $i$th input is equivalent. Indeed, Equations (5) and (6) can be inverted.

3.2 Step 2. From canonical moments to support points

From a given sequence of canonical moments, one wishes to reconstruct the support of a discrete measure. We introduce the Stieltjes Transform, which connects canonical moments of a measure to its support. The Stieltjes transform [23] of a scalar measure $\mu$ is defined as

[TABLE]

The transform $S(z,\mu)$ is an analytic function of $z$ in $\mathbb{C}\backslash\text{supp}(\mu)$. If $\mu$ has a finite support then

[TABLE]

where the support points of the scalar measure $\mu$ are distinct and denoted by $x_{1},\allowbreak\dots,x_{n}$, with corresponding weights $\omega_{1},\dots,\omega_{n}$. Alternatively, the weights are given by $\omega_{i}=\lim_{z\rightarrow x_{i}}(z-x_{i})S(z)$. We can rewrite the transform as a ratio of two polynomials with no common zeros. The zeros of the denominator being the support of $\mu$.

[TABLE]

where $P_{*}^{(n)}(z)=\prod_{i=1}^{n}(z-x_{i})$ and

[TABLE]

The Stieljes transform can also be written as a continuous fraction, some basic definitions and properties of continuous fraction are postponed to B.

Theorem 3.1** ([23, Theorem 3.3.1]).**

Let $\mu$ be a probability measure on the interval $[a,b]$ and $z\in\mathbb{C}\backslash[a,b]$, then the Stieltjes transform of $\mu$ has the continued fraction expansion (see Appendix B for notation):

[TABLE]

Where we recall that $\zeta_{n}:=p_{n-1}(1-p_{n})$.

Theorem 3.1 states that the Stieltjes transform can be computed when one knows the canonicals moments. It immediately follows from Equation (7), Theorem 3.1 and Lemma B.1 that we have the following recursive formula for $P_{\ast}^{(n)}$

[TABLE]

where $P_{\ast}^{(-1)}=0$, $P_{\ast}^{(0)}=1$. The support of $\mu$ thus consists of the roots of $P_{\ast}^{(n)}$. This obviously leads to the following theorem.

Theorem 3.2** ([23, Theorem 3.6.1]).**

Let $\mu$ denote a measure on the interval $[a,b]$ supported on $n$ points with canonical moments $p_{1},p_{2},\dots$. Then, the support of $\mu$ is the set of $\{x:P_{\ast}^{(n)}(x)=0\}$ defined by Equation (8).

In the following we consider a fixed sequence of moments $\mathbf{c}_{n}=(c_{1},\dots,\allowbreak c_{n})\allowbreak\in M_{n}$, let $\mu$ be a measure supported on at most $n+1$ points, such that its moments up to order $n$ coincide to $\mathbf{p}_{n}=(p_{1},\dots,p_{n})$ the corresponding sequence of canonical moments related to $\mathbf{c}_{n}$, as described in Section 3.1. Corollary 3.3 is the moment version of Theorem 3.2. The only difficulty compared to Theorem 3.2 is that one tries to generate admissible measures supported on at most $n+1$ Dirac masses. Given a measure supported on strictly less than $n+1$ points, the question is therefore to know whether it makes sense to evaluate the $n+1$ roots of $P_{\ast}^{(n+1)}$. A limit argument is used for the proof.

Corollary 3.3**.**

Consider a sequence of moments $\mathbf{c}_{n}=(c_{1},\allowbreak\dots,\allowbreak c_{n})\in M_{n}$, and the set of measures

[TABLE]

We define

[TABLE]

Then there exists a bijection between $\mathcal{A}_{\Delta}$ and $\Gamma$.

Proof 3.4**.**

Without loss of generality we can always assume $a=0$ and $b=1$ as the problem is invariant using affine transformation. We first consider the case where $\text{card}(\text{supp}(\mu))$ is exactly $n+1$. From Theorem 3.2, the polynomial $P_{\ast}^{(+1)}$ is well defined with $n+1$ distinct roots corresponding to the support of $\mu$. Notices that this implies that $(p_{1},\dots,p_{2n-1})$ belongs to $]0,1[^{2n-1}$ and that $p_{2n},p_{2n+1}$ or $p_{2n+2}$ belong to $\{0,1\}$.

Now, the functions $g(x,z)=1/(z-x)$ are equicontinuous for $z$ in any compact region which has a positive distance from $[0,1]$. The Stieljes transform is a finite sum of equicontinuous functions and therefore also equicontinuous. Thus if a measure $\mu$ converges weakly to $\mu^{*}$, the convergence must be uniform in any compact set with positive distance from $[0,1]$ (see [28]). It is then always possible to restrict ourselves to measures of cardinal $m<n+1$, by letting $p_{k}$ converge to [math] or $1$ for $2m-2\leq k\leq 2m$. Note that by doing so the polynomials $P_{\ast}^{(m)}$ and $P_{\ast}^{(n+1)}$ will have the same roots. But, $P_{\ast}^{(n+1)}$ and $Q^{(n-1)}$ will have some others roots of multiplicity strictly equal (see Equation (7) and (8)). The corresponding weights of these roots are vanishing, so that the measures extracted from $P_{\ast}^{(m)}$ and $P_{\ast}^{(n+1)}$ are the same.

Remark 3.5**.**

From a computational point of view, as the proof relies on a limit argument, we can always generate $p_{k}\in]0,1[\ ,\ \text{ for }n+1\leq k\leq 2n+1$. This prevents the condition $p_{k}\in\{0,1\}\Rightarrow p_{j}=0$ for $j>k$.

We use Section 3.1 to transform the $N_{i}$ constraints on the moments of the $i$th input into $N_{i}$ canonical moment constraints. However, the construction of the polynomial $P_{i\ast}^{(N_{i}+1)}$ requires the sequence $(p_{1}^{(i)},\dots,p_{2N_{i}+1}^{(i)})$. The $N_{i}$ first canonical moments of this sequence are known by the constraints, while the canonical moments $(p_{k})_{N_{i}+1\leq k\leq 2N_{i}+1}\in\Gamma$ constitute $N_{i}+1$ free parameters, in equal number to the cardinal of $\mu_{i}$. The computation of $\Gamma$ is very simple, it is basically done by random generation of $N_{i}+1$ numbers in $]0,1[$, yet it allows to generate the support of all the measures in $\mathcal{A}_{\Delta}$. This provides a very nice parameterization of the problem that takes naturally into account the constraints.

We also provide a geometrical interpretation of this results. Consider the support of point $(x_{1},\dots,x_{n+1})$ of a discrete measure in $\mathcal{A}_{\Delta}$. The set of support point for all measure in $\mathcal{A}_{\Delta}$ is a manifold in $\mathbb{R}^{n+1}$. Corollary 3.3 described the structure of the manifold and states it is an algebraic variety. Indeed, it is the zero locus of the set of polynomial $P_{\ast}^{(n)}$ built from a sequence of canonical moments in $(p_{1},\dots,p_{n})+\Gamma$.

3.3 Step 3. From support points to weights

From the positions of a discrete measure $\mu$ in $\mathcal{A}_{\Delta}$ generated in Section 3.2. We easily recover the associated weights. Indeed, we enforce $N_{i}$ constraints on the moments of a scalar measure $\mu_{i}=\sum_{j=1}^{N_{i}+1}\omega_{j}^{(i)}\delta_{x_{j}^{(i)}}$, supported by at most $N_{i}+1$ points according to Theorem 2.1. A noticeable fact is that as soon as the $N_{i}+1$ support points of the distribution are set, the corresponding weights are uniquely determined. Indeed, the $N_{i}$ constraints lead to $N_{i}$ equations, and one last equation derives from the measure mass equals to 1. For each $1\leq i\leq d$, the following $N_{i}+1$ linear equations holds

[TABLE]

The determinant of the previous system is a Vandermonde matrix. Hence, the system is invertible as long as the $(x_{i}^{(j)})_{j}$ are distinct.

3.4 Step 4. Computation of the objective function

Thanks to Sections 3.1, 3.2, and 3.3, we can compute the positions $(x_{i}^{(j)})_{1\leq j\leq N_{i}+1}$ and the weights $(\omega_{i}^{(j)})_{1\leq j\leq N_{i}+1}$ of the $i$th input of some $\mu\in\mathcal{A}_{\Delta}$. We can therefore compute the following probability of failure (in Equation (4)):

[TABLE]

We recall that the positions and consequently the weights, were determined using Corollary 3.3 from a sequence of canonical moments $(p_{k})_{N_{i}+1\leq k\leq 2N_{i}+1}\allowbreak\in\Gamma$. So that, the exploration of $\mathcal{A}_{\Delta}$ is parameterized with canonical moments in $\Gamma$. No constraints need to be enforced, as a discrete measure generated from the canonical moments naturally satisfies the moment constraints. The p.o.f is then optimized globally using a differential evolution solver [29]. Algorithm 1 summarizes step 1 to step 4 in order to compute the lowest probability of failure (4). The main cost of the algorithm arises from the high number of metamodel calls for $G$, evaluated on a $d$-dimensional grid of size $\prod_{i=1}^{d}(N_{i}+1)$.

4 Modified algorithm for inequality constraints

In the following, we consider inequality constraints for the moments. The optimization set reads

[TABLE]

One can notice that $\alpha_{i}^{(j)}\leq\mathbb{E}_{\mu_{i}}[x^{j}]\leq\beta_{i}^{(j)}$ is equivalent to enforcing two constraints, thus drastically increasing the dimension of the problem. However, it is possible to restrict ourselves to one constraint. Considering the convex function $\varphi_{i}^{(j)}:x\mapsto(x^{j}-\alpha_{i}^{(j)})(x^{j}-\beta_{i}^{(j)})$, Jensen’s inequality states that $\varphi_{i}^{(j)}(\mathbb{E}_{\mu_{i}}(x))\leq\mathbb{E}_{\mu_{i}}(\varphi_{i}^{(j)}(x))$. Therefore, the sole constraint $\mathbb{E}(\varphi_{i}^{(j)}(x))\leq 0$ ensures $\alpha_{i}^{(j)}\leq\mathbb{E}_{\mu_{i}}[x^{j}]\leq\beta_{i}^{(j)}$. Without loss of generality we still consider measures $\mu_{i}$ that are convex combinations of $N_{i}+1$ Dirac masses, for $i=1,\dots,d$.

We now propose a modified version of Algorithm 1 to solve the problem with inequality constraints. For $i=1,\dots,d$, we denote the moments lower bounds $\bm{\alpha}_{i}=(\alpha_{i}^{(1)},\dots,\alpha_{i}^{(N_{i})})$ and the moments upper bounds $\bm{\beta}_{i}=(\beta_{i}^{(1)},\dots,\beta_{i}^{(N_{i})})$. We use Equation (6) to calculate the corresponding moment sequence $\bm{\alpha}^{\prime}_{i}$ and $\bm{\beta}^{\prime}_{i}$ after affine transformation to $[0,1]$.

The p.o.f of algorithm 2 has $d+2\times\sum_{i=1}^{d}N_{i}$ arguments. The new parameters are actually the first ${(N_{i})}_{|i=1,\dots,d}$th moments of the inputs that were previously fixed. A new step in the algorithm is needed to calculate the canonical moments up to degree $N_{i}$ for $i=1,\dots,d$. This ensures that the constraints are satisfied while the canonical moments from degree $N_{i}+1$ up to degree $2N_{i}+1$ can vary between $]0,1[$ in order to generate all possible measures. The increase of the dimension does not affect the computational times neither the complexity. Indeed, the main cost still arises from the large number of evaluation of the code $G$, that remains equal to $\prod_{i=1}^{d}(N_{i}+1)$. Once again this new p.o.f function can be optimized using any global solver.

5 Numerical tests on a toy example

5.1 Presentation of the hydraulic model

In the following, we address a simplified hydraulic model [30]. This code calculates the water height $H$ of a river subject to a flood event. It takes four inputs whose initial joint distribution is detailed in Table 1. It is always possible to calculate the plug-in quantiles for those particular distributions. The result is given in Figure 4, which present the initial CDF. However, as we desire to evaluate the robust quantile over a class of measures, we present in Table 2 the corresponding moment constraints that the variables must satisfy. The constraints are calculated based on the initial distributions, while the bounds are chosen in order to match the initial distributions most representative values.

The height of the river $H$ is calculated through the analytical model

[TABLE]

We are interested in the flood probability $\sup_{\mu\in\mathcal{A}}\mathbb{P}(H\geq h)$.

5.2 Maximum constraints order influence

We compare the influence of the constraint order on the optimum. The initial distributions and the constraints enforced are available in Table 2. The value of the constraints correspond to the moments of the initial distributions. Figure 4 shows how the size of the optimization space $\mathcal{A}$ decreases by adding new constraints. A differential evolution solver was used to perform the optimization. The initial CDF was computed with a Monte Carlo algorithm. One can observe that enforcing only one constraint on the mean will give a robust quantile significantly larger than the one of the initial distribution. On the other hand, adding three constraints on every inputs reduces quite drastically the space so that the optimal quantile found are close to the one of the initial CDF. A good compromise is to enforce two constraints on the first two moments. This is equivalent to enforcing the mean and the variance on the input distributions, which is a natural way of proceeding. Using this rule, the moment space is not to large so that the worst case quantile remains physically viable.

5.3 Comparison with the Mystic framework

We highlight the interest of the canonical moments parameterization by comparing its performances with the Mystic framework [22]. Mystic is a Python toolbox suitable for OUQ. In Figure 5 one can see the comparison between Mystic and our algorithm. Both computations were realized with an identical solver, and computational times were similar ($\approx$30 min). We enforced one constraint on the mean of each input (see Table 2). The performance of the Mystic framework is outperformed by our algorithm. Indeed, the generation of the weights and support points of the input distributions is not optimized in the Mystic framework. Hence, an intermediary transformation of the measure is needed in order to respect the constraints. During this transformation, the support points can be send out of bounds so that the measure is no more admissible. Many population vectors are rejected, which reduces the overall performance of the algorithm. Meanwhile, our algorithm warrants the exploration of the whole admissible set of measure without any vector rejection.

6 Application to the use-case

6.1 Presentation of the use-case

Our use-case consists in thermal-hydraulic computer experiments, typically used in support of regulatory work and nuclear power plant design and operation. The numerical model is based on code CATHARE 2 (V2.5$\_$3mod3.1) which simulates the time evolution of physical quantities during a thermal hydraulic transient. The simulated accidental transient is an Intermediate Break Loss Of Coolant Accident (IBLOCA) with a break on the cold leg and no safety injection on the broken leg. One run of the computer code takes around 20 minutes on an ordinary computer. Nevertheless, in regards of both physical phenomena and dimensions of the system, we could consider our case under study as a simplified model with respect to a realistic modeling of a reactor. In this use-case, $d=27$ scalar inputs variables of CATHARE code are uncertain and defined by their probability density function. They correspond to various physical parameters, for instance: interfacial friction, critical flow rates, heat transfer coefficients, etc. They are all considered mutually independent. The output variable $Y$ is a single scalar which is the maximal peak cladding temperature (PCT) during the accident transient, see an example in Figure 7.

The number $n$ of simulations chosen for the design of experiments is a compromise between the CPU time required for each simulation and the number of input parameters. For uncertainty propagation and metamodel building purpose, it is a common rule to chose $n$ at least 10 times the dimension $d$ of the input vector [31]. Here $n=1000$ simulations were performed using a space filling Latin Hypercube Sample (LHS) in dimension 27, thus providing a nice coverage of the high-dimensional input space [7].

A screening based on the Hilbert-Schmidt Independence Criterion (HSIC) dependence measure [32], was performed on the $n=1000$ learning simulations. The hypothesis “$\mathcal{H}_{0}^{(k)}$: the input $X_{k}$ and the output $Y$ are independent" was rejected for 9 inputs with a significance level $\alpha=0.1$. Those inputs, designated as PII, are given by Table 3. The screening based on the HSIC takes into account the whole output variability. However, in some pathological case a low influential input for global variability may have an important impact on a $p$-quantile. A perspective to solve this issue, is to consider the conditional and target HSIC [33, p.23], this will be the aim of a future work.

The input distributions described in Table 3 derive from the CIRCE method [34, 35], where in a Bayesian setting the posterior distributions are computed with respect to an experimental database. The moment presented in Table 3 correspond to the moments of the the calibrated distributions. Those are the moments enforced to the measures in the moment class. This way, we can compare the quantiles obtained through the plug-in and full-Gp approaches with the maximal quantile of the OUQ framework. Indeed, with these constraints the initial distribution is also an element of the moment class.

The Gaussian process (Gp) is then build on the PII reduced space, conditioned from the available $n=1000$ learning simulations. We usually consider in computer experiments an anisotropic (stationary) covariance, and the covariance kernel is here chosen as a Matérn 5/2. The metamodel accuracy is evaluated using the predictivity coefficient $Q^{2}$ [36]:

[TABLE]

where $(x_{i})_{1\leq i\leq n_{test}}$ is a test sample, $(y_{i})_{1\leq i\leq n_{test}}$ are the corresponding observed outputs and $(\widehat{y}_{i})_{1\leq i\leq n_{test}}$ are the metamodel predictions. We use a Leave one out strategy in order to perform the validation on the learning sample, we obtain $Q^{2}=0.92$. As already discussed in the introduction, once the predictive metamodel has been built, it can be used to perform uncertainty propagation and in particular, estimate quantiles.

6.2 Results on the use-case

Two constraints were enforced on the first two moments of each inputs as depicted in Table 3. Indeed, informations on the input distributions are often reduced to the knowledge of the means and variances. Moreover, in regards to the toy example, only one constraints on the mean defines a moment class too large, which leads to a worst case scenario not physically viable. At last, more constraints increases the number of support points of each discrete measure on the extreme points, thus increases the computational cost.

We successfully applied the methodology on the 9 dimensional restricted Gp metamodel of the CATHARE code. However, the computation was one day long for each threshold. We restricted the computation of the CDF to a small specific area of interest (high quantile 0.5-0.99) and we parallelized the task so that the computation did not exceed one week. One can compare, in Figure 8, the results of the computation realized with the Mystic framework and our algorithm. Once again, the computations were performed with identical solvers. It confirms the difficulty for the Mystic framework to explore the whole space of admissible measures. Indeed, the maximal quantile obtained with the Mystic framework is lower than the result of our optimization. This proves the efficiency of the canonical moments parameterization to solve this optimization problem.

One can compare the estimation of the 95$\%$-quantile of the peak cladding temperature for the IBLOCA application in Table 4. A 90$\%$-confidence interval for the empirical quantile estimator was constructed with a bootstrap method. The plug-in and full-Gp approaches were defined in the introduction and correspond to a classical estimation of the quantile using respectively the predictor of the metamodel and the full Gaussian process [2]. The OUQ method corresponds to the maximal quantile, when the input distributions are only defined by two of their moments (see Table 3). It is optimal in a sense that it minimize the uncertainty on the input distribution considering the available information. With this information, industrials are able to quantify the worst impact of the inputs uncertainty on the measure of risk, and adapt their choice of safety margins.

The overall optimization can also be realized via Algorithm 2, by relaxing the equality constraints into inequality constraints. This way, it is possible to quantify the sensitivity of the maximal quantile to the moment values.

7 Summary

Metamodels are widely used in industry to perform uncertainty propagation, in particular to evaluate measures of risk such as high quantiles. In this work, we successfully increased the robustness of the quantile evaluation by removing the main sources of uncertainties tainting the inputs of the computer code. We evaluated the maximum measure of risk over a class of distribution. We focus on set of measures only known by some of their moments, and adapted the theory of canonical moments into an improved methodology for solving OUQ problems. Our objective function has been parameterized with the canonical moments, which allows the natural integration of the constraints. The optimization can therefore be performed free of constraints, thus drastically increasing its efficiency. The restriction to moment constraints suits most of practical engineering cases. We also provide an algorithm to deal with inequality constraints on the moments, if an uncertainty lies in their values. Our algorithm shows very good performances and great adaptability to any constraints order. However, the optimization is subject to the curse of dimension and should be kept under 10 input parameters.

The joint distribution of the optimum is a discrete measure. One can criticize that it hardly corresponds to a physical, real world, interpretation. In order to address this issue, we will search for new optimization sets whose extreme points are not discrete measure. The unimodal class found in the literature of robust Bayesian analysis or the $\varepsilon$-contamination class, might be of some interest in this situation. New measures of risk will also be explored, for instance, superquantiles [37], and Bayesian estimates associated to a given utility or loss function [38], which are of particular industrial interest.

Appendix A Proof of duality proposition 2.1

Proof A.1**.**

we denote by \displaystyle a=\sup_{\mu\in\mathcal{A}}\bigg{[}\inf\left\{h\in\mathbb{R}\ ;\ F_{\mu}(h)\geq p\right\}\bigg{]}\ and $\displaystyle\ b=\inf\left\{h\in\mathbb{R}\;|\;\inf_{\mu\in\mathcal{A}}F_{\mu}(h)\geq p\right\}$. In order to prove $a=b$, we proceed in two step. First step, we have

[TABLE]

so that $b\geq a$. Second step, because $a$ is the sup of the quantiles,

[TABLE]

so that

[TABLE]

and $b\leq a$.

Appendix B Basic properties of continuous fraction

Lemma B.1**.**

A finite continued fraction is an expression of the form

[TABLE]

The quantities $A_{n}$ and $B_{n}$ are called the $n$th partial numerator and denominator. There are basic recursive relations for the quantities $A_{n}$ and $B_{n}$ given by

[TABLE]

for $n\geq 1$ with initial conditions

[TABLE]