Vincent Chabridon
email: vincent.chabridon@onera.fr

Mathieu Balesdent

Jean-Marc Bourinet

Jérôme Morio

Nicolas Gayton

Reliability-based sensitivity estimators of rare event probability in the presence of distribution parameter uncertainty

Keywords: Distribution parameter uncertainty, Rare event simulation, Adaptive importance sampling, Reliability sensitivity analysis, Score functions

This paper aims at presenting sensitivity estimators of a rare event probability in the context of uncertain distribution parameters (which are often not known precisely or poorly estimated due to limited data). Since the distribution parameters are also affected by uncertainties, a possible solution consists in considering a second probabilistic uncertainty level. Then, by propagating this bi-level uncertainty, the failure probability becomes a random variable and one can use the mean estimator of the distribution of the failure probabilities (i.e. the "predictive failure probability", PFP) as a new measure of safety. In this paper, the use of an augmented framework (composed of both basic variables and their probability distribution parameters) coupled with an Adaptive Importance Sampling strategy is proposed to get an efficient estimation strategy of the PFP. Consequently, double-loop procedure is avoided and the computational cost is decreased. Thus, sensitivity estimators of the PFP are derived with respect to some deterministic hyper-parameters parametrizing a priori modeling choice. Two cases are treated: either the uncertain distribution parameters follow an unbounded probability law, or a bounded one. The method efficiency is assessed on two different academic test-cases and a real space system computer code (launch vehicle stage fallback zone estimation).

Introduction

Reliability analysis and sensitivity analysis are two major steps in uncertainty quantification of complex systems. For a large variety of applications, assessing the reliability of complex engineering systems (such as aerospace ones) implies, first, to build a dedicated computer code whose aim is to mimick the behavior of the real system. This code can be high-fidelity and consequently, costly-to-evaluate. In uncertainty quantification, it is often considered as an input-output black-box.

Then, one needs to track down and quantify the uncertainties affecting the basic input variables (i.e. physical variables) or those arising in the model itself. Finally, one can propagate these uncertainties through the simulation code and estimate, with some dedicated methods, a so-called failure probability associated to an unsafe and undesired state of the system [START_REF] Morio | Estimation of Rare Event Probabilities in Complex Aerospace and Other Systems: A Practical Approach[END_REF]. In a context of highly safe systems (e.g., systems implying potential risks in terms of human security, environmental impact and/or huge financial loss), the low failure probability requires a huge computational cost to be estimated by crude Monte Carlo (CMC) [START_REF] Rubinstein | Simulation and the Monte Carlo Method[END_REF] which can make these calculations intractable, especially for time-demanding simulation models.

For these reasons, several methods are available to handle this problem of rareness: approximation methods of the failure region [START_REF] Lemaire | Structural Reliability[END_REF], simulation methods based on Monte Carlo simulations or on quasi-random sampling [START_REF] Morio | A survey of rare event simulation methods for static input-output models[END_REF], and, finally, surrogatebased methods [START_REF] Sudret | Meta-models for structural reliability and uncertainty quantification[END_REF].

After the uncertainty propagation phase, it is often relevant to examinate how sensitive some output quantities are with respect to (w.r.t.) the variability affecting input quantities. One the one hand, one can apply sensitivity analysis of model output (SAMO) whose aim is to study how the model output of the computer code varies regarding the input uncertainty [START_REF] Saltelli | Sensitivity Analysis in Practice: A Guide to Assessing Scientific Models[END_REF][START_REF] Helton | Survey of sampling-based methods for uncertainty and sensitivity analysis[END_REF][START_REF] Storlie | Implementation and evaluation of nonparametric regression procedures for sensitivity analysis of computationally demanding models[END_REF][START_REF] Iooss | A Review on Global Sensitivity Analysis Methods[END_REF][START_REF] Wei | Variable importance analysis: a comprehensive review[END_REF][START_REF] Borgonovo | Sensitivity analysis: a review of recent advances[END_REF]. On the other hand, reliability-based sensitivity analysis (RbSA) aims at quantifying the impact of the variability affecting any input quantity on the estimated measure of safety. In this second type of sensitivity analysis, the quantity of interest (QoI) is no longer the model output but a reliability measure. Differences appearing between SAMO and RbSA are due, not only to the different QoIs under consideration, but also to the potential numerical issues appearing in the rare event probability estimation procedure (definition of the limit-state function, isoprobabilistic transformation, estimation strategy, global estimation cost) which can make the overall numerical procedure more complicated compared to SAMO. In the reliability framework, several methods have been developed to quantify the sensitivity of the failure probability w.r.t. deterministic distribution parameters of the basic variables. On the one hand, several global sensitivity methods have been developed: for instance, variance-decomposition-based methods (through the use of Sobol' indices) in [START_REF] Morio | Influence of input PDF parameters of a model on a failure probability estimation[END_REF], methods based on density perturbation in [START_REF] Liu | Relative Entropy Based Method for Probabilistic Sensitivity Analysis in Engineering Design[END_REF][START_REF] Lemaître | Density modification-based reliability sensitivity analysis[END_REF] or, more recently, quantile-oriented sensitivity estimators based on constrast functions in [START_REF] Browne | Estimate of quantile-oriented sensitivity indices[END_REF]. On the other hand, local sensitivity methods (often based on partial derivatives of the QoI) are also available. However, their formulation and implementation procedure mostly depend on the method used to estimate the failure probability.

One can cite the so-called "sensitivity factors" [START_REF] Bjerager | Parametric Sensitivity in First Order Reliability Theory[END_REF][START_REF] Melchers | A fast approximate method for parameter sensitivity estimation in Monte Carlo structural reliability[END_REF][START_REF] Bourinet | FORM Sensitivities to Distribution Parameters with the Nataf Transformation[END_REF] and "omission factors" [START_REF] Madsen | Omission sensitivity factors[END_REF] which can be seen as by-products of applying the First-Order Reliability Method (FORM). Other measures specifically dedicated for sampling-based methods are available: for instance, one can cite [START_REF] Wu | Computational Methods for Efficient Structural Reliability and Reliability Sensitivity Analysis[END_REF][START_REF] Wu | Adapative Importance Sampling (AIS)-Based System Reliability Sensitivity Analysis Method[END_REF] for Adaptive Importance Sampling (AIS), [START_REF] Lu | Reliability sensitivity method by line sampling[END_REF] for Line Sampling, [START_REF] Song | Reliability sensitivity by method of moments[END_REF] for the Method of Moments and [START_REF] Song | Subset simulation for structural reliability sensitivity analysis[END_REF] for Subset Simulations. Among these local measures for RbSA, some of them rely on the mathematical concept of score function (SF) [START_REF] Rubinstein | Simulation and the Monte Carlo Method[END_REF] (also known as kernel function [START_REF] Millwater | Universal properties of kernel functions for probabilistic sensitivity analysis[END_REF]). From a general point of view, the SF method relies on using the derivatives of the probability density function (pdf) in an integral equation (e.g., an expected value equation) and enables to estimate simultaneously (i.e. using the same samples) the performance function and all its derivatives w.r.t. distribution parameters. In the reliability context, this approach has been used, considering deterministic distribution parameters, in [START_REF] Wu | Computational Methods for Efficient Structural Reliability and Reliability Sensitivity Analysis[END_REF][START_REF] Rahman | Stochastic sensitivity analysis by dimensional decomposition and score functions[END_REF][START_REF] Millwater | Development of a localized probabilistic sensitivity method to determine random variable regional importance[END_REF]. Finally, in the two papers [START_REF] Millwater | Probabilistic Sensitivity Analysis With Respect to Bounds of Truncated Distributions[END_REF][START_REF] Lee | Sampling-Based Stochastic Sensitivity Analysis Using Score Functions for RBDO Problems With Correlated Random Variables[END_REF], the problem of estimating the sensitivities of a failure probability regarding deterministic distribution parameters which are bounds of truncated distributions (e.g., uniform or truncated Gaussian) is described.

To sum up, the previous methods appear to be dedicated to problems where distribution parameters (e.g., mean or standard deviation of a Gaussian distribution) are assumed to be deterministic (they are set to a certain value). However, in the case where distribution parameters are affected by uncertainties (due to lack of knowledge or because of limited data), one needs to consider both uncertainty sources (i.e. the one affecting the basic variables and the one affecting the distribution parameters) in the probability estimation, especially for "high-consequence" systems, as stressed in [START_REF] Helton | Uncertainty and sensitivity analysis in the presence of stochastic and subjective uncertainty[END_REF][START_REF] Helton | Probability of loss of assured safety in systems with multiple time-dependent failure modes: Representations with aleatory and epistemic uncertainty[END_REF][START_REF] Helton | Conceptual structure and computational organization of the 2008 performance assessment for the proposed high-level radioactive waste repository at Yucca Mountain[END_REF]. Bayesian framework has been used and recommended in [START_REF] Der Kiureghian | Aleatory or epistemic? Does it matter?[END_REF][START_REF] Pasanisi | Some useful features of the Bayesian setting while dealing with uncertainties in industrial practice[END_REF][START_REF] Pasanisi | Estimation of a quantity of interest in uncertainty analysis: Some help from Bayesian decision theory[END_REF] for assessing reliability under distribution parameter uncertainty. The effect of parameter uncertainty in reliability assessment has also been studied within some extra-probabilistic frameworks, also known as "imprecise probabilities" [START_REF] Beer | Imprecise probabilities in engineering structures[END_REF][START_REF] Schöbi | Surrogate models for uncertainty quantification in the context of imprecise probability modelling[END_REF]. Several papers [START_REF] Helton | An exploration of alternative approaches to the representation of uncertainty in model predictions[END_REF][START_REF] Helton | Sensitivity analysis in conjunction with evidence theory representations of epistemic uncertainty[END_REF][START_REF] Helton | A sampling-based computational strategy for the representation of epistemic uncertainty in model predictions with evidence theory[END_REF][START_REF] Helton | Representation of analysis results involving aleatory and epistemic uncertainty[END_REF][START_REF] Helton | Quantification of margins and uncertainties: Alternative representations of epistemic uncertainty[END_REF] propose notional examples illustrating and comparing the use of several different frameworks (probabilistic and extra-probabilistic such as evidence theory and possibility theory) in the characterization of epistemic uncertainty affecting distributions of input aleatory uncertainty.. Among these contributions, one can cite for instance [START_REF] Qiu | Probabilistic interval reliability of structural systems[END_REF] in which reliability assessment under interval-based analysis of distribution parameters is treated, [START_REF] Schöbi | Structural reliability analysis for p-boxes using multi-level meta-models[END_REF] in which the propagation of parameter uncertainty through probability-boxes is achieved, [START_REF] Li | On upper and lower bounds for quantity of interest in problems subject to epistemic uncertainty[END_REF] for the derivation of probability bounds and [START_REF] Pedroni | Hierarchical propagation of probabilistic and non-probabilistic uncertainty in the parameters of a risk model[END_REF] for integrating either probabilistic and non-probabilistic variables in an hybrid framework. If the topic of reliability under parameter uncertainty is now studied with probabilistic and non-probabilistic approaches, it appears that the works can hardly be compared and contrasted because of two main characteristics: first, the initial assumptions are different and the information sources for modeling input uncertainties are treated in different ways; second, the output QoI of the reliability analysis is also very different quantitatively (single vs. multiple measures of safety) and qualitatively (point measure, bounds, full distributions, truncated distributions, etc.). Consequently, in this work, the Bayesian view is adopted [START_REF] Pasanisi | Estimation of a quantity of interest in uncertainty analysis: Some help from Bayesian decision theory[END_REF] and the following measure of safety incorporating this bi-level uncertainty can be considered: the predictive failure probability, as described in [START_REF] Chabridon | Evaluation of failure probability under parameter epistemic uncertainty: application to aerospace system reliability assessment[END_REF]. Thus, dedicated RbSA methods have to be developed to examinate sensitivities of this measure of safety w.r.t. the deterministic distribution hyper-parameters which parametrize the a priori distribution on the parameters of the basic variables. In this paper, local sensitivity estimators of the predictive failure probability w.r.t. these hyper-parameters are derived. In a context of low failure probability estimation, since the calculations are performed in the joint space of stochastic basic variables and stochastic distribution parameters (i.e. the "augmented space"), an augmented formulation of Adaptive Importance Sampling can be used so as to estimate jointly the predictive failure probability and the sensitivities at a reduced cost compared to CMC. This paper is organized as follows. Section 2 adresses the problem of reliability assessment under distribution parameter uncertainty and aims at introducing the formal concepts and notations. Section 3 describes both the derived sensitivity estimators and their implementation within an augmented Adaptive Importance Sampling strategy. Section 4 illustrates the benefits of such a methodology on different test-cases (two academic test-cases and one black-box computer code issued from aerospace research) and a synthesis gathering the key aspects and issues of the proposed approach is provided at the end of this section. Finally, a conclusion gathering the most important outcomes of this paper is given in Section 5.

Reliability assessment under probability distribution law parameter uncertainty via an augmented sampling

Usual formulation of a time-invariant reliability problem

Let us consider a complex system represented by a static input-output black-box simulation model M(•) whose numerical evaluation can be achieved pointwise. Thus, it can be seen as a real-valued function defined on D X ⊆ d such that M : D X → . In this paper, for the sake of simplicity, the model output is assumed to be reduced to a scalar y

= M(x) ∈ D Y ⊆ .
Moreover, the model is assumed to be deterministic. Assessing reliability to such a system traditionally relies on setting a probabilistic framework in which one assumes that the input basic variables (i.e. variables that are directly observable or which can be directly linked to empirical data [START_REF] Der Kiureghian | Aleatory or epistemic? Does it matter?[END_REF]) can be modeled by a set of random variables X 1 , X 2 , . . . , X d , where d is the input dimension. Under a perfect state of knowledge about the distribution family and the distribution parameters, and provided a perfect knowledge about the dependence between input variables (which can be set up through the copula framework [START_REF] Lebrun | Do Rosenblatt and Nataf isoprobabilistic transformations really differ?[END_REF]), a parametrized probabilistic model can be assumed for the d-dimensional random input vector X through the cumulative distribution function (cdf) F X (•) and corresponding joint pdf f X (•). By propagating the uncertainties through the model via the equation Y = M(X), the ouput is also a random variable which can be characterized by the distribution F Y (•). For reliability assessment, the behavior of the system can be characterized by a so-called limit-state function (lsf) (note that the common "limit-state" analysis performed in reliability is also known as "margin analysis" in nuclear engineering [START_REF] Helton | Quantification of margins and uncertainties: Alternative representations of epistemic uncertainty[END_REF]) g : D X → . A common formulation for the lsf can be g(X) def = y th -Y = y th -M(X) with y th ∈ a given scalar threshold beyond which the system falls into a failure state. The aim of g(•) is to split the space of realizations of the random vector X into a failure domain F x = {x ∈ D X : g(x) ≤ 0} and a safe domain S x = {x ∈ D X : g(x) > 0}. Hence, the failure probability p f is given by:

p f = [g(X) ≤ 0] = F x f X (x)dx = D X F x (x) f X (x)dx = f X F x (X) (1)
where dx = dx 1 . . . dx d and F x (•), the indicator function of the failure domain, is defined by:

F x (x) = 1 if x ∈ F x
and F x (x) = 0 otherwise. Estimating this probability can be achieved by a large variety of methods [START_REF] Lemaire | Structural Reliability[END_REF]. Among these techniques, some require to reformulate the reliability problem in the standard normal space where all the basic variables are independent and standard-Gaussian distributed (e.g., FORM, SORM, Line Sampling or Subset Simulations in the original version by [START_REF] Au | Estimation of small failure probabilities in high dimensions by subset simulation[END_REF]). This reformulation is achieved by means of a transformation which can be, in general, either the Nataf one [START_REF] Nataf | Détermination des distributions dont les marges sont données[END_REF] or the Rosenblatt one [START_REF] Rosenblatt | Remarks on a Multivariate Transformation[END_REF], depending on the available probabilistic information (i.e. marginals and potentially the copula) [START_REF] Lebrun | Do Rosenblatt and Nataf isoprobabilistic transformations really differ?[END_REF].

Bayesian framework for reliability assessment under parameter uncertainty

Among information available to construct a parametrized probabilistic model for the input basic variables, one often has only access to limited data, possibly unadapted literature-based recommendations and finally subjective expert opinions [START_REF] Ditlevsen | Structural Reliability Methods[END_REF]. Thus, an imperfect state of knowledge [START_REF] Ditlevsen | Generalized Second Moment Reliability Index[END_REF][START_REF] Der Kiureghian | Measures of Structural Safety Under Imperfect States of Knowledge[END_REF] may lead to a misestimation of the failure probability and lead to dramatic consequences in terms of risk mitigation. Statistical uncertainty arises in the estimation procedure of the probability distribution parameters when one can only deal with unsufficient measures or data. In some cases, it may also happen that neither data nor expert judgment are available, which imposes to the engineer yet to make a choice for the values of parameters. This is often encountered in the field of complex systems for which data acquisition is difficult.

To take distribution parameter uncertainty into account in the reliability assessment, the previous framework set up in subsection 2.1 needs to be reinterpreted.

In this paper, two levels of uncertainty are considered: the first one represents the variability in the basic input variables and thus affects the input random vector X when the second one represents the lack of knowledge affecting the distribution parameters Θ. To do so, we consider the following Bayesian hierarchical model [START_REF] Gelman | Bayesian Data Analysis[END_REF]:

X ∼ f X|Θ (x|θ) : D X ⊆ d → + (first layer) (2a) Θ ∼ f Θ|ξ (θ|ξ) : D Θ ⊆ k → + (second layer) (2b) ξ = (ξ 1 , ξ 2 , . . . , ξ q ) ⊤ ∈ D ξ ⊆ q (third layer). (2c)
In this hierarchical representation, one can distinguish three layers of inputs:

• the first layer is constituted by the random vector X gathering the stochastic basic variables. Based on prior knowledge, a probability distribution can be assumed through the choice of a parametric model. This random vector can be possibly of great dimension and may involve a complex dependence structure [START_REF] Lebrun | Do Rosenblatt and Nataf isoprobabilistic transformations really differ?[END_REF];

• the second layer is constituted by uncertain and deterministic (i.e. supposed to be known accurately enough) distribution parameters. Adding such a layer is coherent with the Bayesian point of view of modeling either "uncertain" (in the sense of stochastic) or "unknown but fixed" parameters [START_REF] Gelman | Bayesian Data Analysis[END_REF]. If one can theoretically consider a possible dependence structure for the vector Θ, it is often, from a pragmatic point of view (if we do not care about the deterministic parameters in it), of rather small dimension compared to the basic random vector and structured such that quadrature schemes or quasi-random sampling can be easily used to sample over the space D Θ [START_REF] Chabridon | Evaluation of failure probability under parameter epistemic uncertainty: application to aerospace system reliability assessment[END_REF]. Other techniques, such as the association of a rank correlation procedure developed in [START_REF] Iman | Rank correlation plots for use with correlated input variables[END_REF] and quasi-random sampling have been successfully used in literature [START_REF] Helton | Latin hypercube sampling and the propagation of uncertainty in analyses of complex systems[END_REF]. In the present paper, the prior distribution is mostly assumed to be derived from expert judgment or from a limited set of data. Thus, despite the fact that Bayes' theorem is not used as an updating procedure (but could be, if more data or a better characterization was available), one can consider that this prior distribution characterizes epistemic uncertainty affecting the distribution parameters. As a remark, no deep investigation about the type of prior distributions (neither in terms of the choice of a parametric family, nor in terms of informativeness of the prior) is claimed in this study. This topic is beyond the scope of the paper and the interested reader is invited to refer to [START_REF] Gelman | Bayesian Data Analysis[END_REF] and [START_REF] Roos | Sensitivity Analysis for Bayesian Hierarchical Models[END_REF] for any further information;

• the third layer is composed of deterministic hyper-parameters gathered in ξ (which can be some moments or bounds) of prior distributions of stochastic parameters.

Following the paradigm proposed in [START_REF] Der Kiureghian | Aleatory or epistemic? Does it matter?[END_REF], the distinction between aleatory and epistemic uncertainties in a given set of modeling, is achieved here by considering that the probabilistic modeling of X could represent aleatory uncertainty whereas the probabilistic prior modeling of Θ could represent epistemic uncertainty. As an example, let us consider a problem involving some aleatory uncertainty defined as a single random variable X following a normal distribution f X |Θ (first layer).

One of the distribution parameters (e.g., the mean value denoted by Θ) can also be affected by epistemic uncertainty (second layer). Thus one can add a probability distribution f Θ|ξ on Θ indexed by its own distribution hyper-parameters ξ (third layer).

Under this bi-level uncertainty (because of the first two layers), the failure probability p f as written in Eq. (1) only represents a realization of the random variable P f for a given realization θ of the distribution parameters. One can express this conditional failure probability P f (θ) such that:

P f (θ) = [g(X) ≤ 0 | Θ = θ] = D X F x (x) f X|Θ (x|θ)dx = f X|Θ F x (X) | Θ = θ . (3)
Then, considering the mathematical expectation of the previous random variable w.r.t. the prior density f Θ|ξ allows to get the so-called "predictive failure probability" P f , as introduced in [START_REF] Der Kiureghian | Measures of Structural Safety Under Imperfect States of Knowledge[END_REF], which is a measure of reliability taking into account the effect of the lack of knowledge about distribution parameters:

P f (ξ) def = f Θ|ξ [P f (Θ)] = f Θ|ξ f X|Θ F x (X) | Θ = θ = D Θ P f (θ) f Θ|ξ (θ|ξ)dθ. (4)
This predictive failure probability is thus the mean of all the conditional failure probabilities over the prior distribution on uncertain parameters. This quantity can be estimated by two different approaches [START_REF] Chabridon | Evaluation of failure probability under parameter epistemic uncertainty: application to aerospace system reliability assessment[END_REF]. This indicator allows to have a single aggregated measure of uncertainty gathered in both Θ and X. Following the discussion in [START_REF] Chabridon | Evaluation of failure probability under parameter epistemic uncertainty: application to aerospace system reliability assessment[END_REF], to have disagreggated effects, one should reconstruct the full distribution of P f (θ) and use a more conservative measure (e.g., a quantile) than the mean estimator. Other formalisms (e.g., p-boxes, p-bounds) could provide such an information [START_REF] Schöbi | Structural reliability analysis for p-boxes using multi-level meta-models[END_REF]. However, the estimation of such indicators is not the scope of this paper.

Nested vs. augmented approaches for estimating the predictive failure probability

On the one hand, the double expectation in Eq. ( 4) can be estimated via a double-loop estimation procedure over both integration domains. This strategy is called here "nested reliability approach" (NRA). The outer sampling loop over D Θ can be achieved, depending on the structure of the vector Θ, either by deterministic sampling methods such as quadrature schemes [START_REF] Davis | Methods of Numerical Integration[END_REF], or by CMC or any other advanced methods for design of experiments (DOE) [START_REF] Pronzato | Design of computer experiments: space filling and beyond[END_REF]. For the inner loop (i.e.

the nested loop), the conditional failure probability can be estimated by any approximation or simulation method such as those cited in subsection 2.1. Various applications of NRA for rare event probability estimation and safety assessment under parameter uncertainty can be found in [START_REF] Helton | Conceptual structure and computational organization of the 2008 performance assessment for the proposed high-level radioactive waste repository at Yucca Mountain[END_REF][START_REF] Limbourg | Accelerated uncertainty propagation in two-level probabilistic studies under monotony[END_REF][START_REF] Balesdent | Rare Event Probability Estimation in the Presence of Epistemic Uncertainty on Input Probability Distribution Parameters[END_REF].

On the other hand, a different strategy can be adopted by considering an "augmented" random vector

Z def = (X, Θ) ⊤ defined on D Z = D X × D Θ (where × is the Cartesian product) with joint pdf f Z|ξ (z|ξ) def = f (X,Θ)|ξ (x, θ|ξ) = f X|Θ (x|θ) f Θ|ξ (θ|ξ)
such that the expression in Eq. ( 4) can be rewritten as follows:

P f (ξ) = D Θ D X F x (x) f X|Θ (x|θ) f Θ|ξ (θ|ξ)dxdθ = D Z F z (z) f Z|ξ (z|ξ)dz = f Z|ξ F z (Z)|ξ (5)
where

F z = {z ∈ D Z : g(z) ≤ 0}.
This augmented formulation numerically implies to estimate the expected value in Eq (5). This is called here the "augmented reliability approach" (ARA) or in a wider Bayesian context the "predictive approach" since it relies on the use of predictive distributions [START_REF] Geisser | Aspects of the Predictive and Estimative Approaches in the Determination of Probabilities[END_REF][START_REF] Christensen | Bayesian Point Estimation Using the Predictive Distribution[END_REF]. From a practical point of view, this approach has been used in the context of FORM analysis in [START_REF] Hong | Evaluation of the Probability of Failure with Uncertain Distribution Parameters[END_REF][START_REF] Most | Efficient structural reliability methods considering incomplete knowledge of random variable distributions[END_REF].

Recently, in [START_REF] Chabridon | Evaluation of failure probability under parameter epistemic uncertainty: application to aerospace system reliability assessment[END_REF], a numerical comparison between NRA and ARA has been achieved and the coupling between both approaches and advanced simulation-based method such as the Subset Simulations has been investigated. Possibilities offered by ARA in terms of simulation cost reduction, estimation accuracy and robustness w.r.t. several numerical challenges concerning real aerospace test-cases have been illustrated in this work. The interested reader may refer to [START_REF] Chabridon | Evaluation of failure probability under parameter epistemic uncertainty: application to aerospace system reliability assessment[END_REF] and [START_REF] Pasanisi | Estimation of a quantity of interest in uncertainty analysis: Some help from Bayesian decision theory[END_REF] for any further details.

Among technical steps involved in ARA, the transformation between the physical space and the standard normal space is an issue that needs to be adressed by adapting the traditional Rosenblatt transformation to an "augmented" version of it, i.e. by applying the transformation on the vector Z. One needs to apply first the transformation to the stochastic distribution parameters and then to the stochastic basic variables, conditioned by the distribution parameters [START_REF] Chabridon | Evaluation of failure probability under parameter epistemic uncertainty: application to aerospace system reliability assessment[END_REF]. In the following sections, the use of the prefix "ARA" in front of any name implies that the sampling strategy is set in the augmented framework.

Sensitivity analysis of predictive failure probability with respect to distribution hyper-parameters

The use of local sensitivities in this study is motivated mainly by two reasons. The first one corresponds to the way the problem is set as explained above (i.e. one wants to measure the sensitivity w.r.t. a local choice of ξ). Finally, the second reason, which is a key constraint in this work, remains the limited allowable extra simulation budget one can afford to get sensitivities while the rare event probability estimation can be very expensive too (without any consideration here of any use of a metamodel).

Derivation of sensitivity estimators in the augmented framework

The gradient of the predictive failure probability P f w.r.t. the vector of the hyper-parameters ξ is defined as follows:

∇ P f (ξ) = ∂ P f (ξ) ∂ ξ j , j = 1, . . . , q ⊤ . ( 6
)
Depending on the nature of the hyper-parameter ξ j , two cases are considered:

• Case #1: ξ j is an hyper-parameter of a prior distribution with an unbounded support;

• Case #2: ξ j is an hyper-parameter of a prior distribution with a bounded/truncated support.

Derivations for Case #1

The partial derivative of the predictive failure probability w.r.t. the j-th component of ξ is given by:

∂ P f (ξ) ∂ ξ j = ∂ ∂ ξ j D Θ P f (θ) f Θ|ξ (θ|ξ)dθ = D Θ P f (θ) ∂ f Θ|ξ (θ|ξ) ∂ ξ j dθ. ( 7
)
Note that, in the previous derivations, the differential and integral operators are switched due to Lebesgue's dominated convergence theorem [START_REF] Jacod | Probability Essentials. Universitext[END_REF]. Following the idea given in [START_REF] Rubinstein | The score function approach for sensitivity analysis of computer simulation models[END_REF], one can use the so-called importance sampling trick so as to get an expectation w.r.t. the same probability measure as the one used for the failure probability estimation. It thus comes:

∂ P f (ξ) ∂ ξ j = D Θ P f (θ) ∂ ∂ ξ j f Θ|ξ (θ|ξ) f Θ|ξ (θ|ξ) f Θ|ξ (θ|ξ)dθ (8a) = D Θ D X F x (x) f X|Θ (x|θ)dx ∂ ln f Θ|ξ (θ|ξ) ∂ ξ j f Θ|ξ (θ|ξ)dθ (8b) = D Θ D X F x (x) s j (θ, ξ) f X|Θ (x|θ)dx f Θ|ξ (θ|ξ)dθ (8c) = f Z|ξ F z (Z) s j (Θ, ξ) (8d)
where s j (θ, ξ)

def = ∂ ln f Θ|ξ (θ|ξ) ∂ ξ j
is called the "score function" (SF). As recalled in Introduction, the SF approach has been widely used in the RbSA literature under single-level uncertainty [START_REF] Wu | Computational Methods for Efficient Structural Reliability and Reliability Sensitivity Analysis[END_REF][START_REF] Rahman | Stochastic sensitivity analysis by dimensional decomposition and score functions[END_REF][START_REF] Millwater | Development of a localized probabilistic sensitivity method to determine random variable regional importance[END_REF][START_REF] Millwater | Probabilistic Sensitivity Analysis With Respect to Bounds of Truncated Distributions[END_REF][START_REF] Lee | Sampling-Based Stochastic Sensitivity Analysis Using Score Functions for RBDO Problems With Correlated Random Variables[END_REF]. It is used here, in the context of bi-level uncertainty, for three main reasons: first, it enables to provide the targeted local sensitivity in Eq. ( 7) as a simple derivation (assuming rather simple cases for prior distribution); then, it fits naturally to an importance-sampling-based estimation framework; and finally, it provides an efficient sampling-based estimator which allows to avoid finite difference schemes.

One should notice that, to avoid any confusion, in the above equations and in the rest of the paper, the vector Θ is explicitely written instead of Z = (X, Θ) ⊤ since the dependence w.r.t. ξ is through Θ. An example of a SF associated to an unbounded normal prior for an uncertain distribution parameter Θ j is given in Table 1. Then, considering a sample {Z (i) } N i=1 of N independent and identically distributed (i.i.d.) copies of the augmented vector Z, one can derive the following Monte Carlo estimator:

∂ P f (ξ) ∂ ξ j ≈ MC 1 N N i=1 F z (Z (i) ) s j (Θ (i) , ξ). ( 9
)
As a remark, if the probability is estimated with ARA, its gradient in Eq. ( 9) may be estimated at a reduced cost. However, due to central limit theorem, one can show that the predictive failure probability and its sensitivities may have different convergence rates. Thus, the variance associated to the two asymptotic distributions will necessarily differ.

Normal (Case #1) 1 ξ 2 2π exp -1 2 θ j -ξ 1 ξ 2 2 ξ 1 = µ Θ j , ξ 2 = σ Θ j 1 ξ 2 θ j -ξ 1 ξ 2 1 ξ 2 θ j -ξ 1 ξ 2 2 -1 Uniform (Case #2) 1 ξ 2 -ξ 1 [ξ 1 , ξ 2 ] 1 ξ 2 -ξ 1 -1 ξ 2 -ξ 1 3.1.2. Derivations for Case #2
Here, at least, one of the basic variables X j follows a parametrized distribution whose one parameter Θ j follows a bounded or truncated distribution denoted by f Θ j |ξ (θ j |ξ), with ξ = (ξ 1 , ξ 2 , . . . , ξ q ) ⊤ the vector of the hyper-parameters.

Here, ξ j ∈ ξ could be either a bound or a moment. In the next derivations, one considers that the distribution parameters

Θ i , i = 1, . . . , k are independent, which leads to f Θ|ξ (θ|ξ) = k i=1 f Θ i |ξ (θ i |ξ).
Before deriving the sensitivities, one should notice that, in this case, the support D Θ is a function of ξ j as the support of Θ j is either bounded or truncated (note that, without any loss of generality, the bounds are denoted by a(ξ j ) and b(ξ j ) in Eq. ( 10b)). In the following, the notation Θ -j (respectively θ -j ) denotes the vector without the j-th component Θ j (repectively θ j ) which depends on the hyper-parameter ξ j . Thus, the derivations start as follows:

∂ P f (ξ) ∂ ξ j = ∂ ∂ ξ j D Θ (ξ j ) P f (θ) f Θ|ξ (θ|ξ)dθ (10a) = ∂ ∂ ξ j k-1 b(ξ j ) a(ξ j ) P f (θ) f Θ j |ξ (θ j |ξ) f Θ -j |ξ (θ -j |ξ)dθ j dθ -j (10b) = k-1 ∂ ∂ ξ j b(ξ j ) a(ξ j ) P f (θ) f Θ j |ξ (θ j |ξ)dθ j f Θ -j |ξ (θ -j |ξ)dθ -j (10c) = k-1 ∂ ∂ ξ j I P f (θ), ξ f Θ -j |ξ (θ -j |ξ)dθ -j (10d)
where I P f (θ), ξ is an integral whose bounds (denoted by a and b) depend on the parameter ξ j . Indeed, using the Leibniz integral rule for differentiation of a definite integral whose limits are functions of the differential variables, one gets:

∂ ∂ ξ j I P f (θ), ξ = ∂ ∂ ξ j b(ξ j ) a(ξ j ) P f (θ) f Θ j |ξ (θ j |ξ)dθ j (11a) = b(ξ j ) a(ξ j ) P f (θ) ∂ f Θ j |ξ (θ j |ξ) ∂ ξ j dθ j + P f θ -j , θ j = b(ξ j ) f Θ j |ξ (b(ξ j )|ξ) ∂ b(ξ j ) ∂ ξ j -P f θ -j , θ j = a(ξ j ) f Θ j |ξ (a(ξ j )|ξ) ∂ a(ξ j ) ∂ ξ j ( 11b
)
where P f θ -j , θ j = •(ξ j ) represents the failure probability estimated with θ j fixed to one of the integration bounds (i.e.

•(ξ j ) = a(ξ j ) or b(ξ j )).

To illustrate the previous reasoning, one can apply these derivations to a test-case involving a distribution parameter Θ j following a continuous uniform distribution such that Θ j ∼ U ([a, b]). In this example, ξ j can be either a bound (a or b) or even a moment (e.g., µ Θ j = a+b 2 or σ Θ j = b-a 12 ). Let us first assume that ξ j = a:

∂ ∂ a I P f (θ), ξ = b a P f (θ) ∂ f Θ j |ξ (θ j |ξ) ∂ a dθ j + 0 -P f (θ -j , θ j = a) f Θ j |ξ (a|ξ) × 1 (12a) = b a P f (θ) ∂ f Θ j |ξ (θ j |ξ) ∂ a dθ j - 1 b -a P f (θ -j , θ j = a). (12b)
Hence, using the same trick as in paragraph 3.1.1, one gets:

∂ P f (ξ) ∂ a = k P f (θ) ∂ ∂ a f Θ j |ξ (θ j |ξ) f Θ j |ξ (θ j |ξ) f Θ j |ξ (θ j |ξ) f Θ -j |ξ (θ -j |ξ)dθ j dθ -j - 1 b -a k-1 P f (θ -j , θ j = a) f Θ -j |ξ (θ -j |ξ)dθ -j (13a) = k d F x (x) f X|Θ (x|θ)dx ∂ ln f Θ j |ξ (θ j |ξ) ∂ a f Θ|ξ (θ|ξ)dθ - 1 b -a k-1 d F x (x) f X|Θ (x|θ -j , θ j = a)dx f Θ -j |ξ (θ -j |ξ)dθ -j (13b) = k d F x (x) s j (θ j , a) f X|Θ (x|θ)dx f Θ|ξ (θ|ξ)dθ - 1 b -a k-1 d F x (x) f X|Θ (x|θ -j , θ j = a)dx f Θ -j |ξ (θ -j |ξ)dθ -j (13c) = f Z|ξ F z (Z) s j (Θ j , ξ) - 1 b -a f Z|Θ j =a,ξ F z (Z -j )|Θ j = a (13d)
where Z -j stands for (X, Θ -j ) ⊤ , which means that the integration is achieved over a (n aug -1)-dimensional hypersurface with

n aug = k × d.
Using formulas given in Table 1 and the linearity of expectation, one obtains the following two sensitivities:

∂ P f (ξ) ∂ a = 1 b -a f Z|ξ F z (Z) -f Z|Θ j =a,ξ F z (Z -j )|Θ j = a = 1 b -a P f (ξ) -P a f,aux ; (14a) ∂ P f (ξ) ∂ b = - 1 b -a f Z|ξ F z (Z) -f Z|Θ j =b,ξ F z (Z -j )|Θ j = b = 1 b -a P b f,aux -P f (ξ) . ( 14b
)
where P a f,aux and P b f,aux are two "auxiliary" failure probabilities which have to be estimated over the (n aug -1)-dimensional hypersurface. These derivations are consistent with those provided in [START_REF] Millwater | Probabilistic Sensitivity Analysis With Respect to Bounds of Truncated Distributions[END_REF] and [START_REF] Lee | Sampling-Based Stochastic Sensitivity Analysis Using Score Functions for RBDO Problems With Correlated Random Variables[END_REF] in the context of RbSA under single-level uncertainty (i.e. deterministic distribution parameters and traditional failure probability as a QoI). As for the interpretation, the sensitivity estimators obtained in Eqs. (14a) and (14b) are close to those derived in [START_REF] Millwater | Probabilistic Sensitivity Analysis With Respect to Bounds of Truncated Distributions[END_REF] and similarly involve flux integrals, here P a f,aux and P b f,aux , over a (n aug -1)-dimensional space. As a final remark, one should highlight the fact that for uniform prior distributions, one can calculate the sensitivities of the predictive failure probability w.r.t. the moments

µ Θ j = a+b 2 or σ Θ j = b-a
12 by combining the previous sensitivities w.r.t. the bounds such that:

∂ P f (ξ) ∂ µ Θ j = ∂ P f (ξ) ∂ a ∂ a(µ Θ j ) ∂ µ Θ j + ∂ P f (ξ) ∂ b ∂ b(µ Θ j ) ∂ µ Θ j = 2 ∂ P f (ξ) ∂ a + ∂ P f (ξ) ∂ b ; (15a) ∂ P f (ξ) ∂ σ Θ j = ∂ P f (ξ) ∂ a ∂ a(σ Θ j ) ∂ σ Θ j + ∂ P f (ξ) ∂ b ∂ b(σ Θ j ) ∂ σ Θ j = 12 ∂ P f (ξ) ∂ b - ∂ P f (ξ) ∂ a . ( 15b
)
As a remark, one could notice that for Case #2, the computational cost (i.e. estimating two probabilities) is similar to firstly, FDM is an approximation method to compute the gradient; secondly, the type of finite difference scheme (forward, backward or centered) may influence the results; thirdly, choosing an optimal perturbation step can be problematic. The proposed method allows to overcome these difficulties by providing (assuming the SF is available for the prior distribution)

an exact formulation to get the gradient of the predictive failure probability and estimate it independently of any choice for perturbation step.

Up to now, the previous framework allows to estimate jointly, within the same sampling phase and with limited extra computational effort, a predictive failure probability and its derivatives w.r.t. a priori deterministic hyper-parameters.

However, even with the ARA strategy, ARA/CMC is not able to handle rare event probability estimation regarding real engineering system safety assessment. The idea of the next subsection is to propose a dedicated numerical methodology combining efficient sampling strategy in the augmented space and the above sensitivity estimators.

Proposed methodology (ARA/AIS) for reliability-based sensitivity analysis in a context of rare event probability estimation

Estimating a rare event probability with CMC can be cumbersome and can even become intractable for costly-to-evaluate computer codes. Importance Sampling (IS) is now a well-known variance-reduction technique [START_REF] Rubinstein | Simulation and the Monte Carlo Method[END_REF]. The idea is to use a so-called "auxiliary density" h(•) to generate samples such that, if this density is the optimal one, one gets a zero variance of the IS estimator of the rare event probability. To introduce it, one can start from the observation that the following equality holds:

D Z F z (z) f Z|ξ (z|ξ)dz = D Z F z (z) f Z|ξ (z|ξ) h(z) h(z)dz = D Z F z (z)w(z)h(z)dz (16)
where w(z)

def = f Z|ξ (z|ξ) h(z)
is called the likelihood ratio [START_REF] Rubinstein | Simulation and the Monte Carlo Method[END_REF]. This weight is introduced in the probability estimator and takes into account the change in the pdf to generate samples. Thus, considering {Z (i) } N i=1 as an i.i.d. sample drawn according to f Z|ξ (z|ξ), the IS estimators for both the probability and its sensitivities can be derived such that:

P f (ξ) ≈ IS 1 N N i=1 F z (Z (i) ) w(Z (i) ) (17a) ∂ P f (ξ) ∂ ξ j ≈ IS 1 N N i=1 F z (Z (i) ) w(Z (i) ) s j (Θ (i) , ξ). ( 17b
)
The estimator P f of P f given in the right hand side in Eq. (17a) is unbiased (i.e. h P f = P f ) and its variance P f reduces to zero as the density h(•) equals the optimal auxiliary density h * (•) which is given by:

h * (z) = F z (z) f Z|ξ (z|ξ) P f . ( 18
)
Since this quantity depends on the predictive probability P f one would like to estimate, this intricate problem can be solved by using Adaptive Importance Sampling (AIS) techniques [START_REF] Tokdar | Importance sampling: a review[END_REF]. These techniques aim at, using different adaptive strategies, to sequentially approximate the optimal auxiliary density.

In this paper, we propose to adapt two existing AIS techniques, namely the Nonparametric Adaptive Importance Sampling (NAIS) [START_REF] Morio | Non-parametric adaptive importance sampling for the probability estimation of a launcher impact position[END_REF] and the Adaptive Importance Sampling by Cross-Entropy (AIS-CE) [START_REF] Homem-De Mello | Estimation of rare event probabilities using cross-entropy[END_REF] to estimate, within the augmented framework, both the predictive failure probability and its sensitivities w.r.t. deterministic hyper-parameters at a reduced cost compared to CMC. Thus, one presents in this paper two different methods which are called respectively "ARA/NAIS" and "ARA/AIS-CE". Two generic algorithms are given in Algorithm 1 (for ARA/NAIS) and Algorithm 2 (for ARA/AIS-CE).

Since these algorithms are already existing but adapted to the augmented framework, a complete description of the details is not given in this paper for the sake of conciseness. The interested reader can refer to [START_REF] Morio | Estimation of Rare Event Probabilities in Complex Aerospace and Other Systems: A Practical Approach[END_REF] for any further detail about practical implementation of the original version of these AIS algorithms. However, one should insist on the fact that these two rare event estimation algorithms rely on some assumptions and imply some parameters (e.g., the choice of a given kernel K(•) for ARA/NAIS, the choice of an initial parametric family for the auxiliary pdf h λ (•) with λ ∈ Λ for ARA/AIS-CE, the choice of the ρ-quantiles for both methods) whose tuning and performance optimization is not treated in the core discussion of the present paper (interested reader may refer to [START_REF] Morio | Estimation of Rare Event Probabilities in Complex Aerospace and Other Systems: A Practical Approach[END_REF] for more information about it). Note that the algorithms' parameters used in this paper are given as footnotes of, respectively Algorithms 1 and 2 .

The global methodology is able to handle both Case #1 and Case #2 as detailed in subsection 3.1. However, in Case #2, estimating the predictive failure probability is not enough to get the sensitivities w.r.t. the bounds. An estimation of the auxiliary failure probability, as shown in Eqs. (14a) and (14b) for the uniform case, is required. To do so, one needs to apply a second time the ARA/AIS algorithm (with either NAIS or AIS-CE as the core algorithm) to estimate this quantity and allow an accurate estimation of the sensitivity.

As an illustration of the ARA/AIS sampling strategy, one can consider a simple test case involving two basic random variables, similar to a so-called "Resistance -Sollicitation" problem (similar to the one presented further in numerical applications in subsection 4.2). Assuming these variables are Gaussian, one can imagine that, due to limited data, their distribution parameters are affected by epistemic uncertainty: for instance, the standard deviation of the resistance variable and the mean value of the sollicitation one are not perfectly known. Some prior distributions are considered to model the a priori knowledge about these parameters. The propagation of this bi-level uncertainty (to get an estimate of the predictive failure probability and its sensitivities) using the ARA/AIS strategy is presented only with the ARA/NAIS method in Figure 1 for the sake of conciseness.

In Figure 1a, one can see two sets of N = 10 3 realizations in the original space (denoted as "X-space") of the vector

Θ = (Θ 1 , Θ 2 ) ⊤ .
The first set of samples corresponds to the first iteration of the ARA/NAIS algorithm. The second set corresponds to the final iteration of the algorithm. By comparing them, one can notice the modification between the initial sampling density and the final one. Such a modification is also noticeable in Figure 1b which represents the same samples in the standard normal space (denoted as "U-space"). One can see the shrinkage of the initial standard normal density to the optimal one. Based on these realizations, one can observe the corresponding realizations of the basic variables gathered in X. Figure 1c and Figure 1d shows the corresponding samples plotted respectively in X-space and in U-space. Again, from the first iteration to the final one, convergence of the density towards the optimal one is noticeable. One can also highlight that considering a second level of uncertainty does affect the realizations of the basic variables by changing the shape of the distribution of X.

In the following section, the numerical efficiency of the ARA/AIS method is demonstrated on two academic test-cases and on a real aerospace test-case.

Application examples

The following numerical applications have been implemented in Matlab ® and performed using a rare event simulation toolbox developed at ONERA -The French Aerospace Lab.

Figure 1: Illustration of the ARA/NAIS method on a Resistance -Sollicitation test-case (similar to the one described in Example #1 (cf. 4.2)). In this example, two different sets of samples, drawn at iterations #1 and #4 of the ARA/NAIS method (Algorithm 1) are presented. One can see the evolution of the samples reflecting the adaptive evolution of the augmented sampling density towards a near optimal one.

Methodology and comparison metrics

The overall methodology developed in this paper, as summed up in Table 2, relies on four blocks detailed below:

• firstly, the ARA/AIS method has been implemented and tested on three different test-cases (two academic and one industrial) of increasing difficulty (see the first column and the footnotes below Table 2 for the specifications of each test-case);

• secondly, reference results for the estimation of the predictive failure probability and its sensitivities are obtained using a CMC with large sample size performed by ARA (see the second column denoted ARA/CMC, see [START_REF] Chabridon | Evaluation of failure probability under parameter epistemic uncertainty: application to aerospace system reliability assessment[END_REF] for details about the coupling between ARA and sampling methods);

• thirdly, as shown in the third column of Table 2, the two different methods, namely ARA/AIS-CE and ARA/NAIS are applied on the three test-cases for the sake of comparison, but also to illustrate the modular aspect of the methodology (the black squares stand for the performed calculations). Finally, after validation of the method, the impact of the increasing rareness of the failure event (regarding a limited simulation budget available) is studied as an extension of the last two test-cases (see the four black stars ). This last extension is called "rare event context" in the rest of the paper. One should notice that, for these specific cases, the reference calculation by ARA/CMC is considered as computationally "intractable" (see the corresponding crosses ×).

Algorithm 1 -Generic algorithm for ARA/AIS (with an ARA/NAIS plug-in box in this example).

Start

Define: pdf f Z|ξ , budget N , model M(•), threshold y th , quantile ρ ∈ [0, 1] a , kernel K(•) b Set: k = 1 and h 0 = f Z|ξ Generate: i.i.d. samples z (i) 1 of {Z (i) 1 } N i=1 ∼ h 0 Evaluate: Y (1) 1 = M(z (1) 1 ), . . . , Y (N ) 1 = M(z (N ) 1 ) Compute: empirical ρ-quantile γ 1 of the samples {Y (i) 1 } N i=1 While γ k < y th do Estimate: l k = 1 kN k j=1 N i=1 {M(z (i) j )≥γ k } (Z (i) j ) f Z|ξ (Z (i) j ) h j-1 (Z (i) j )
and set w j (Z (i)

j ) = {M(z (i) j )≥γ k } (Z (i) j ) f Z|ξ (Z (i) j ) h j-1 (Z (i) j ) Update c : h k+1 (z) = 1 kN l k det(B 1/2 k ) k j=1 N i=1 w j (Z (i) j )K B 1/2 k (z -Z (i) j ) Set: k = k + 1 Generate: i.i.d. samples z (i) k of {Z (i) k } N i=1 ∼ h k Evaluate: Y (1) k = M(z (1) k ), . . . , Y (N ) k = M(z (N ) k ) Compute: empirical ρ-quantile γ k of the samples {Y (i) k } N i=1
Estimate:

P f = 1 N N i=1 {M(z (i) k )>y th } (Z (i) k ) f Z|ξ (Z (i) k ) h k (Z (i) k ) ∂ P f /∂ ξ j = 1 N N i=1 {M(z (i) k )>y th } (Z (i) k ) f Z|ξ (Z (i) k ) h k (Z (i) k ) s j (Θ (i) k , ξ) End
a In this paper: ρ = 0.9. b In this paper: the Gaussian kernel is used. c B k is a symmetric positive definite bandwidth matrix optimized with asymptotic integrated square error (AMISE) criterion, see [START_REF] Morio | Estimation of Rare Event Probabilities in Complex Aerospace and Other Systems: A Practical Approach[END_REF]. Finally, for a comparison in terms of numerical efficiency of the method w.r.t. CMC estimation, the following standard coefficient ν ARA/AIS is used:

ν ARA/AIS = N ARA/CMC sim N ARA/AIS sim ( 19
)
where N ARA/CMC sim is the required number of CMC samples to reach the same coefficient of variation δ on the probability Algorithm 2 -Generic ARA/AIS-CE plug-in box.

Define: . . . , parametric family of pdfs h λ (•) with λ ∈ Λ a . . .

Set: k = 1 Generate: i.i.d. samples z (i) of {Z (i) } N i=1 ∼ h λ 0 . . . While γ k < y th do
Optimize b : the parameters of the auxiliary pdf family

λ k = argmax λ∈Λ 1 N N i=1 {M(z (i) j )≥γ k } (Z (i) j ) f Z|ξ (Z (i) ) h λ k-1 (Z (i) ) ln[h λ (Z (i) )] Set: k = k + 1 Generate: i.i.d. samples z (i) of {Z (i) } N i=1 ∼ h λ k-1 Evaluate: Y (1) = M(z (1) ), . . . , Y (N ) = M(z (N ) ) Compute: empirical ρ-quantile γ k of the samples {Y (i) } N i=1 . . .
a In this paper: the Gaussian parametric family is used.

b For particular density families (e.g. Gaussian), the pdf optimal parameters have analytical formulas, see [START_REF] Morio | Estimation of Rare Event Probabilities in Complex Aerospace and Other Systems: A Practical Approach[END_REF]. In this paper: both mean and standard deviation are optimized.

estimate for both methods. Thus, it leads to the following expression:

ν ARA/AIS = 1 -P ARA/AIS f N ARA/AIS sim × P ARA/AIS f × δ 2 . ( 20
)
A value of ν ARA/AIS > 1 indicates that the method ARA/AIS is more efficient than CMC for the given test-case. In other words, ν ARA/AIS indicates the quantity by which we can divide the initial CMC simulation budget for a same level of accuracy.

Example #1: a resistance-sollicitation toy-case

Description. The goal of this first academic test-case is to validate the method regarding the estimation accuracy of both the predictive failure probability and its sensitivities w.r.t. the hyper-parameters. Table 3 gives the input data. The lsf

g 1 (X) = R -S = X 1 -X 2 is
linear and involves two independent Gaussian random variables. One assumes that both the standard deviation of the first variable and the mean of the second one are affected by epistemic uncertainty. Thus, two prior distributions (one unbounded and one bounded) are assumed for the uncertain parameters (e.g., based either on limited data, literature-based information or expert opinion). As a remark, one could argue that the choice of a normal prior for Θ 1 may be inappropriate regarding physical constraints, conjugacy and informativeness. In this paper, common priors (i.e.

normal and uniform) are set for the sake of illustration to characterize epistemic uncertainty and for a certain convenience in sampling, without loss of generality. For a pure Bayesian approach (involving possible updating), the reader should refer to [START_REF] Gelman | Bayesian Data Analysis[END_REF][START_REF] Hamada | Bayesian Reliability[END_REF][START_REF] Straub | Bayesian Updating with Structural Reliability Methods[END_REF].

Variable Distribution Parameter #1

Parameter #2

X 1 = R Normal µ X 1 = 7 σ X 1 uncertain X 2 = S Normal µ X 2 uncertain a σ X 2 = 1 Θ 1 = σ X 1 Normal ξ 1 = µ σ X 1 = 0.7 ξ 2 = σ σ X 1 = 0.07 Θ 2 = µ X 2 Uniform ξ 3 = a µ X 2 = 1.5 ξ 4 = b µ X 2 = 2.5
a For fixed values σ X 1 = 0.7 and µ X 2 = 2, p f = 2.50 × 10 -5 . Results. Table 4 gathers numerical values for probability estimates, sensitivities and efficiencies. For comparison purpose, a CMC with N x,θ = 10 8 samples is performed. The coefficient of variation (cv) for any estimate is calculated using a hundred replicates of each algorithm. However, for the sake of conciseness, only the results associated to the highest numerical efficiency ν ARA/AIS are presented in Table 4 for the two proposed methods (here, N x,θ = 10 4 samples/step). As a remark, one can see that both the predictive failure probability and its sensitivities are well estimated with both methods compared to ARA/CMC. Concerning the estimation of the predictive failure probability, one can see that P f = 3.72 × 10 -5 which shows that taking a second uncertainty level into account implies a slight increase compared to the failure probability under single-level uncertainty p f = 2.50 × 10 -5 (see below

Example #2: a two d.o.f. primary/secondary damped oscillator

Description. This nonlinear oscillator is based on a two-degree-of-freedom primary-secondary system, as shown in Figure 2, excited by a white noise base acceleration [START_REF] Der Kiureghian | Efficient Algorithm for Second-Order Reliability Analysis[END_REF]. The system's behavior is characterized by the masses m p and m s , the spring stiffnesses k p and k s , the natural frequencies ω p = (k p /m p ) 1/2 and ω s = (k s /m s ) 1/2 and the damping ratios ζ p and ζ s , where the subscripts p and s respectively refer to the primary and secondary oscillators. As for F s , it denotes the force capacity of the secondary spring. Thus, the reliability of the system can be evaluated using the following lsf:

g 2 (X) = F s -3k s πS 0 4ζ s ω 3 s   ζ a ζ s ζ p ζ s 4ζ 2 a + r 2 + γζ 2 a ζ p ω 3 p + ζ s ω 3 s ω p 4ζω 4 a   = F s -F acc . (21)
This equation defines the potential failure event {F s ≤ F acc } occuring if the force induced by the white noise base acceleration overcomes the force capacity in the secondary spring. Table 5 gives the input data (the set of parameters here for the basic variables are mean values and coefficients of variation). The lsf g 2 (•) is highly nonlinear and involves eight independent lognormal random variables. One assumes that both the mean of the second mass is not precisely known due to measure uncertainty and the mean of the force capacity of the secondary spring is also affected by epistemic uncertainty. Thus, two prior distributions (Gaussian and uniform) are assumed for the uncertain parameters based on expert judgment.

Variable a

Distribution Parameter #1 Parameter #2

X 1 = m p (kg) Lognormal µ X 1 = 1.5 δ X 1 = 10% X 2 = m s (kg) Lognormal µ X 2 uncertain b δ X 2 = 10% X 3 = k p (N.m -1 ) Lognormal µ X 3 = 1 δ X 3 = 20% X 4 = k s (N.m -1 ) Lognormal µ X 4 = 0.01 δ X 4 = 20% X 5 = ζ p (1) Lognormal µ X 5 = 0.05 δ X 5 = 40% X 6 = ζ s (1) Lognormal µ X 6 = 0.02 δ X 6 = 50% X 7 = F s (N) Lognormal µ X 7 uncertain δ X 7 = 10% X 8 = S 0 (m.s -2 ) Lognormal µ X 8 = 100 δ X 8 = 10% Θ 1 = µ X 7 (N) Normal ξ 1 = µ µ X 7 = 21.5 / 27.5 ( ) c ξ 2 = σ µ X 7 = 2.15 / 2.75 ( ) Θ 2 = µ X 2 (kg) Uniform ξ 3 = a µ X 2 = 0.008 ξ 4 = b µ X 2 = 0.012
a The basic variables are independent.

b For fixed values µ X 2 = 0.01 and µ X 7 = 21.5 / 27.5, p f = 4.78 × 10 -5 / 3.78 × 10 -7 . c The second value is for the rare event case, cf. Table 2.

Results. Numerical results summarized in Table 6 show that, for a moderately rare failure event, both ARA/AIS-CE and ARA/NAIS give accurate results in the predictive failure probability estimation compared to ARA/CMC. Again, the coefficient of variation (cv) for any estimate is calculated using a hundred replicates of each algorithm. It first reveals the variations between the failure probability under single-level uncertainty p f and the predictive one P f (here, it increases from 4.78×10 -5 to 2.35 × 10 -4 ) due to the bi-level uncertainty. In this case, considering uncertainty on a distribution parameter makes the system less safe, which can be an important indicator for design or re-design purposes. Concerning sensitivities, most of them are correctly estimated, except for ∂ P f /∂ ξ 2 for which one can observe a tiny relative bias between the proposed methods and reference results. In terms of comparison, once again, since the problem involves two uncertain distribution parameters (Θ 1 unbounded and Θ 2 bounded), one can use Eqs. (15a) and (15b) to get sensitivities w.r.t. moments of Θ 2 instead of its bounds. Finally, in this case, the predictive failure probability seems to be slightly more sensitivive to the a priori choice modelling the uncertain mean of the secondary mass. Thus, the lack of knowledge about the mean value of the mass plays a key role in terms of system safety. As for the convergence of the results, Figure 3a compares the estimated sensitivities (to avoid any redundancy, only the ARA/NAIS plots are presented) to the reference results obtained by ARA/CMC. One can notice the convergence w.r.t. the increasing number of samples per step and a low variability of the two last sensitivities as mentioned previously. Finally, for a moderate rareness of the failure event, the efficiency of the method is promising:

ν ARA/AIS is equals to 14 for ARA/AIS-CE and 17 for ARA/NAIS, meaning the equivalent ARA/CMC simulation budget can be still divided while ensuring a given target accuracy in the estimation.

In the rare event context ( ) (see Table 7), reference results are supposed to be intractable. Again, a hundred replicates were used to get samples' statistics. One can first observe a slight increase for the predictive failure probability (around 7 × 10 -6 ) compared to the failure probability under single-level uncertainty (equals to 3.78 × 10 -7 as given below Table 5).

With the proposed method, one can notice that the estimated values show relatively low coefficient of variation. Comparing sensitivities leads to notice that the rareness of the failure event (i.e. between Table 6 and Table 7) does not impact the relative order in terms of influence. The moments of Θ 2 are still the most influent hyper-parameters of the predictive failure probability. The method ensures a higher efficiency than as the rareness of the probability increases (ν ARA/AIS from 231 for ARA/AIS-CE to 276 for ARA/NAIS). Finally, the global convergence is observed on Figure 3b (for ARA/NAIS only).

Example #3: application to a launch vehicle stage fallback zone estimation

Description. The role of a launch vehicle is to carry a payload (e.g., a satellite) from the Earth's surface to a given orbit.

A traditional expendable space launcher is composed of multiple stages, equipped with their propulsion systems. An example of such a launch vehicle is given in Figure 4. During the flight, uncertainties can affect several variables in multiple disciplines (e.g., on the dynamics perturbations or stage combustion). For instance, focusing on the optimal trajectory assessesment leads to consider the separation point (see "stage separation" in Figure 4) as a key point in terms of uncertainty analysis. Dynamic perturbations, varying unburned propellant left mass, error measurements due to sensors, are all types of uncertainties that are of prime importance to estimate and predict the fallback zone into the ocean. A misestimation can have dramatic consequences in terms of launcher safety and in terms of human security and environmental impact. To assess reliability of such a system, the black-box model M(•) considered here is a trajectory simulation code of the dynamic fallback phase of a generic launcher first stage [START_REF] Chabridon | Reliability-based sensitivity analysis of aerospace systems under distribution parameter uncertainty using an augmented approach[END_REF]. As input, the state vector gathers six independent random variables representing some perturbations of the initial conditions and launcher stage characteristics at separation point:

• X 1 : stage altitude perturbation at separation (∆h (m));

• X 2 : velocity perturbation at separation (∆v (m.s -1 ));

• X 3 : flight path angle perturbation at separation (∆γ (rad));

• X 4 : azimuth angle perturbation at separation (∆ψ (rad));

• X 5 : propellant mass perturbation at separation (∆m (kg));

• X 6 : drag force error perturbation (∆C d dimensionless).

The input probabilistic model is given in Table 8. The code output is the distance D code between the theoretical fallback position into the ocean and the estimated one. The lsf g 3 (•) can be written as follows: [START_REF] Lu | Reliability sensitivity method by line sampling[END_REF] for which the rareness of the failure event depends on the safety threshold distance d safe . To get a moderate target failure probability, one can set the safety distance to d safe = 15, 000 meters while to increase the rareness, one can set it to 20, 000 meters ( ).

g 3 (X) = d safe -M(X) = d safe -D code

Variable a

Distribution Parameter #1 Parameter #2 In addition to this first level of uncertainty, one assumes that epistemic uncertainty is affecting two mean values, respectively µ X 2 and µ X 3 , i.e. the mean values of the perturbations affecting the velocity at separation and the perturbations affecting the flight path angle at separation as shown in Table 8. These physical quantities are difficult to measure and to control in real conditions. For the sake of demonstration, the numerical values (distribution parameters and hyper-parameters) used in this example are hypothetic. The analysis results for the proposed method are averaged over a hundred replications of the algorithm (this code is a simplified version of a real trajectory code with a reasonable computational cost).

X 1 = ∆h (m) Normal µ X 1 = 0 σ X 1 = 1650 X 2 = ∆v (m.s -1 ) Normal µ X 2 uncertain b σ X 2 = 3.7 X 3 = ∆γ (rad) Normal µ X 3 uncertain σ X 3 = 0.001 X 4 = ∆ψ (rad) Normal µ X 4 = 0 σ X 4 = 0.0018 X 5 = ∆m (kg) Normal µ X 5 = 0 σ X 5 = 70
X 6 = ∆C d (1) Normal µ X 6 = 0 σ X 6 = 0.1 Θ 2 = µ X 2 (m.s -1 ) Normal ξ 1 = µ µ X 2 = 0 ξ 2 = σ µ X 2 = 3.7 Θ 3 = µ X 3 (rad) Normal ξ 3 = µ µ X 3 = 0 ξ 4 = σ µ X 3 = 0.
Results. Numerical results gathered in Table 9 show that both methods manage to correctly estimate the predictive failure probability. As a first remark, this predictive failure probability is slightly greater than the failure probability under singlelevel uncertainty p f recalled below Table 8. This shows again that, in this case, the distribution parameter makes the system be less safe. Then, one can highlight the fact that ARA/NAIS gets closer results to reference ones (obtained by ARA/CMC) than ARA/AIS-CE. The problem here is inherent to the cross-entropy method as it shows some difficulty to converge with a correct accuracy when the problem involves multiple failure regions. However, in [START_REF] Kurtz | Cross-entropy-based adaptive importance sampling using Gaussian mixture[END_REF], a cross-entropy-based adaptive importance sampling algorithm using Gaussian mixture has been proposed and seems to manage to compensate this issue.

In terms of comparison, the predictive failure probability seems to be more sensitivive to the hyper-parameters ξ 4 and ξ 3 which are respectively the standard deviation and the mean of Θ 3 = µ X 3 . Here, the lack of knowledge affecting the mean value of the flight path angle perturbation really plays a key role on the final predictive failure probability. This can be a relevant information for refining the a priori probabilistic model for Θ 3 (especially in terms of variance reduction) and set up an investigation policy about the possible reduction of epistemic (statistical) uncertainty affecting Θ 3 . Concerning the efficiencies, while ARA/AIS-CE is inefficient in this specific case where AIS-CE is not accurate enough, ARA/NAIS manages to provide a precise estimation for both the probability and the sensitivities. The convergence plot in Figure 5a illustrates these results.

In the rare event context ( , as given in Table 10, one can still see that ARA/NAIS provides better results than ARA/AIS-CE, even if this one still manages to get relevant orders of magnitude for both the probability and the sensitivities. However, the efficiency of ARA/AIS-CE is annealed by the poor accuracy of the estimation while ARA/NAIS outperforms ARA/CMC by allowing to reduce the simulation budget by 207. Similar comparisons can be drawn to the previous case regarding the relative influence of the hyper-parameters. However, one can still notice that increasing the rareness of the failure event decreased, in proportion, the relative influence of ξ 4 . Finally, the global convergence of ARA/NAIS sensitivity estimation is represented in Figure 5b.

Synthesis about numerical results and discussion

The aim of this subsection is to give a synthesis for the interested reader to get the main advantages and drawbacks of the proposed approach. According to the numerical results, one can sum up the following characteristics:

• ARA/AIS-CE and ARA/NAIS lead to similar results for predictive failure probability and reliability-based sensitivities estimation to the reference approach, namely ARA/CMC. They both enable to reduce the variance of estimation and limit the number of calls to the black-box model;

• ARA/AIS-CE and ARA/NAIS, coupled to the derived sensitivity estimators, enable to estimate sensitivities w.r.t. distribution parameters appearing in both unbounded and bounded priors;

• ARA/AIS-CE suffers from the possible multimodality of the optimal auxiliary density (i.e. problem with multiple failure regions), which is a known issue of this method [START_REF] Kurtz | Cross-entropy-based adaptive importance sampling using Gaussian mixture[END_REF];

• ARA/NAIS is dedicated to problems whose input dimensionality (including the stochastic distribution parameters) is a few dozen while ARA/AIS-CE can handle higher input dimensions. In terms of interpretation of the sensitivities derived in this paper, one should remember that local sensitivities are not quantities that can be easily ranked. Using local sensitivities can only give a qualitative result which should be considered as an indication for the user to know whether he/she should get more information about the lack of knowledge affecting one or more distribution parameters. It cannot be used as a variable importance measure for factor fixing purposes [START_REF] Wei | Variable importance analysis: a comprehensive review[END_REF]. However, they can be useful to highlight some features and underlying behaviors associated, either to the input probabilistic modeling, or to the model behavior at failure. To do so, one should compare sensitivities of the same nature, for instance sensitivities to a mean together, sensitivities to a standard deviation together, sensitivities to bounds together. Another possible solution is to build what is called an elasticity which is a normalized sensitivity [START_REF] Lemaire | Structural Reliability[END_REF]. In reliability literature, several authors advocated different types of elasticities. For instance, Table 11 gathers three different formulas and the corresponding literature sources for computing elasticities (adapted here to the bi-level uncertainty framework). From the authors knowledge, there is no consensus about the elasticity formulation. For the sake of generality, the sensitivity indicator provided in this paper do not take into account any normalization factor. All the derivative-based formulas (and associated numerical results) described in this paper can be combined with the formulas given in Table 11 to obtain elasticities. ξ j = µ Θ j or ξ j = σ Θ j [START_REF] Wu | Computational Methods for Efficient Structural Reliability and Reliability Sensitivity Analysis[END_REF] (f2) e j = ∂ P f ∂ ξ j × ξ j ξ j = µ Θ j or ξ j = σ Θ j [START_REF] Millwater | Probabilistic Sensitivity-Based Ranking of Damage Tolerance Analysis Elements[END_REF] (f3) e j =

∂ P f ∂ ξ j × ξ j P f ξ j = µ Θ j or ξ j = σ Θ j [3]
It is important to notice that the choice of prior distribution is a crucial step in practice and can be lead by various techniques (expert judgment, data analysis, experimental results). For the sake of conciseness, in this paper, the sensitivity estimator formula (and application cases) are explicitely given using only normal and uniform probability distributions as priors. When facing more complex priors (e.g., truncated normal distribution to avoid unrealistic physical data), these formula can be adapted. Finally, for some priors, score functions may not be available [START_REF] Rahman | Stochastic sensitivity analysis by dimensional decomposition and score functions[END_REF][START_REF] Millwater | Probabilistic Sensitivity Analysis With Respect to Bounds of Truncated Distributions[END_REF][START_REF] Lee | Sampling-Based Stochastic Sensitivity Analysis Using Score Functions for RBDO Problems With Correlated Random Variables[END_REF].

Conclusion -perspectives

In this article, new reliability-based sensitivity estimators are proposed using the so-called "augmented reliability approach". The proposed derivative-based local sensitivity estimators of predictive failure probability, with respect to deterministic distribution hyper-parameters, are derived for two cases: firstly, when all the stochastic distribution parameters follow some unbounded prior probability distributions; secondly, when at least one distribution parameter follows a bounded prior.

Thus, this method allows to get sensitivities with either none (in the unbounded case) or a moderate extra computational effort (in the bounded case). To enhance the efficiency of the method, these estimators are derived using an adaptive importance sampling scheme, either using a parametric algorithm (ARA/AIS-CE) or a nonparametric one (ARA/NAIS). Then, three numerical applications on both academic and industrial test-cases (launch vehicle stage fallback zone estimation) are considered, the final ones beeing representative of the complex simulation codes used in aerospace engineering. The comparison with a reference method (ARA/CMC) demonstrates the convergence and the performance of the proposed method.

Finally, this study shows the benefits of using an ARA/AIS strategy when the failure event becomes very rare, especially for complex models.

A first enhancement track could be to study the possibility of enhancing ARA/AIS-CE strategy for reliability and sensitivity assessment under bi-level uncertainty by adapting recent works on AIS-CE [START_REF] Kurtz | Cross-entropy-based adaptive importance sampling using Gaussian mixture[END_REF][START_REF] Wang | Cross-entropy-based adaptive importance sampling using von Mises-Fisher mixture for high dimensional reliability analysis[END_REF] to handle both multiple failure regions and high-dimensionality. In order to further reduce the computational cost of these estimations, another possible enhancement track could be to couple the ARA/AIS strategy to a metamodel such as proposed in [START_REF] Dubourg | Metamodel-based importance sampling for reliability sensitivity analysis[END_REF] to possibly extend it to the bi-level uncertainty problem. However, metamodels also bring their own contribution to the overall uncertaintiy by adding some "model uncertainty". This type of uncertainty has to be handled and should not cover the uncertainty introduced by considering a prior on the distribution parameters. Some recent metamodel-based strategies also consider a so-called "augmented framework" by considering that the metamodel can handle all kinds of epistemic uncertainties in its own definition [START_REF] Walter | Moving particles: A parallel optimal multilevel splitting method with application in quantiles estimation and meta-model based algorithms[END_REF]. These possible tracks are left for future work.

  X|Θ samples in U-space 1e3 samples at iteration #1 1e3 samples at iteration #4

Figure 2 :

 2 Figure 2: Two-degree-of-freedom damped oscillator with primary and secondary systems.

  (a) Estimated sensitivities (error bars) vs. reference results obtained by ARA/CMC (dashed lines). Pf w.r.t. hyper-parameters∂ Pf /∂ξ1 ∂ Pf /∂ξ2 ∂ Pf /∂ξ3 ∂ Pf /∂ξ4(b) Estimated sensitivities (error bars) in a context of rare event ( ).

Figure 3 :

 3 Figure 3: Convergence plots obtained by ARA/NAIS for Example #2.

001 a

 001 The basic variables are independent. b For fixed values µ X 2 = 0 and µ X 3 = 0: for d safe = 15 (km), p f = 1.36 × 10 -4 while for d safe = 20 (km), p f = 2.31 × 10 -7 .

Figure 4 :

 4 Figure 4: Illustration scheme of a first stage fallback phase (adapted from an infographic by Jon Ross (http://www.zlsadesign.com)).

( b )

 b Estimated sensitivities (error bars) in a context of rare event ( ).

Figure 5 :

 5 Figure 5: Convergence plots obtained by ARA/NAIS for Example #3.

Figure 6

 6 Figure6summarizes the evolution of the numerical efficiency ν ARA/NAIS (only for the ARA/NAIS method here) for the two last examples (Example #2 and Example #3) in their two configurations -initial problem ( ) and increased rareness of the failure event ( ) -as a function of the number of samples per step (from 1e3 to 1e4). If the efficiencies for the initial problems reach values around ν ARA/NAIS ≈ 10 to 15 (meaning that one can divide the simulation budget by 15), the efficiencies in a rare event context reach high values such as ν ARA/NAIS ≈ 50 to 270 which indicates that the proposed method can efficiently handle problems with very low predictive failure probabilities whose estimation by ARA/CMC is definitely intractable. As a remark, one can notice that the increasing curves as a function of the number of samples only represent the underlying tradeoff of the proposed method between accuracy (coefficient of variation decreasing as the number of samples increases) and the global computational cost.

Figure 6 :

 6 Figure 6: Efficiency.

Table 1 :

 1 Score functions for normal (Case #1) and uniform (Case #2) prior distributions on an uncertain parameter Θ j .

Distribution

pdf f Θ j |ξ (θ j |ξ 1 , ξ 2 ) Hyper-parameters Score function s 1 (θ j , ξ 1 ) Score function s 2 (θ j , ξ 2 )

Table 2 :

 2 Overall methodology.

	Test-case	Reference	ARA/AIS
		ARA/CMC ARA/AIS-CE ARA/NAIS
	Example #1: Resistance -Sollicitation a	(cf. 4.2)
	Example #2: Nonlinear oscillator b	(cf. 4.3)	×
	Example #3: Launch vehicle fallback zone estimation c (cf. 4.4)	×
	a 2 basic variables, 2 uncertain parameters (1 unbounded & 1 bounded), g(•) linear.
	b 8 basic variables, 2 uncertain parameters (1 unbounded & 1 bounded), g(•) nonlinear.
	c 6 basic variables, 2 uncertain parameters (2 unbounded), g(•) nonlinear.

Table 3 :

 3 Input probabilistic model for Example #1.

Table 4 :

 4 Results for Example #1.

		ARA/CMC	ARA/AIS-CE	ARA/NAIS
		(N x,θ = 10 8 samples)	(N x,θ = 10 4 samples/step) (N x,θ = 10 4 samples/step)
		Estimate	cv	Estimate	cv	Estimate	cv
	P f	3.72 × 10 -5 (1.60 %) 3.73 × 10 -5 (4.39 %)	3.71 × 10 -5 (3.62 %)
	∂ P f /∂ ξ 1	2.77 × 10 -4 (3.50 %) 2.81 × 10 -4 (14.2 %)	2.75 × 10 -4 (13.7 %)
	ν ARA/AIS	-	-	103	-	158	-

∂ P f /∂ ξ 2 1.37 × 10 -4 (10.8 %) 1.45 × 10 -4 (62.9 %) 1.35 × 10 -4 (59.1 %) ∂ P f /∂ ξ 3 3.30 × 10 -5 (1.91 %) 3.32 × 10 -5 (4.91 %) 3.30 × 10 -5 (4.06 %) ∂ P f /∂ ξ 4 8.91 × 10 -5 (1.41 %) 8.87 × 10 -5 (5.99 %) 8.89 × 10 -5 (4.74 %) ∂ P f /∂ µ Θ 2 a 2.44 × 10 -4 -2.44 × 10 -4 -2.44 × 10 -4 -∂ P f /∂ σ Θ 2 1.94 × 10 -4 -1.92 × 10 -4 -1.94 × 10 -4a cf. Eqs. (15a) and (15b) for the uniform case.

Table 3 )

 3 . From these results, one can notice that the sensitivities w.r.t. bounds (i.e. ∂ P f /∂ ξ 3 and ∂ P f /∂ ξ 4 ) show a reduced coefficient of variation compared to those for the unbounded distribution (i.e. ∂ P f /∂ ξ 1 and ∂ P f /∂ ξ 2 ). This could be due to the fact that they are estimated as a difference of two probabilites estimated both by the method which lead to a reduced variance. In terms of comparison, since the problem involves two uncertain distribution parameters (Θ 1 unbounded and Θ 2 bounded), one can use Eqs. (15a) and (15b) to get sensitivities w.r.t. moments of Θ 2 instead of its bounds. Finally, in this case, the predictive failure probability seems to be

	S(t)	x p	x s
	k p	k s
		m p	m s
	ζ p	ζ s

slightly more sensitivive to the a priori choice of the mean values of both distribution parameters than to the choice of the standard deviations. As for the ν ARA/AIS coefficients, they are very high compared unity for both methods. This implies that, for the same level of accuracy, one can reduce the CMC simulation budget by 103 times if one uses ARA/AIS-CE and by 158 times if one uses ARA/NAIS. This gain can be of practical importance for applications involving rare event probabilities.

Table 5 :

 5 Input probabilistic model for Example #2.

Table 6 :

 6 Results for Example #2.

			ARA/CMC	ARA/AIS-CE	ARA/NAIS
			(N x,θ = 10 8 samples)	(N x,θ = 10 4 samples/step) (N x,θ = 10 4 samples/step)
			Estimate	cv	Estimate	cv	Estimate	cv
	P f		2.35 × 10 -4	(0.603 %) 2.36 × 10 -4	(4.69 %)	2.38 × 10 -4	(4.57 %)
	∂ P f /∂ ξ 1		-1.61 × 10 -4 (0.661 %) -1.60 × 10 -4 (4.80 %)	-1.61 × 10 -4 (5.32 %)
	∂ P f /∂ ξ 2		9.04 × 10 -3	(0.602 %) 1.08 × 10 -2	(4.69 %)	1.09 × 10 -2	(4.57 %)
	∂ P f /∂ ξ 3		5.57 × 10 -2	(0.650 %) 5.58 × 10 -2	(4.89 %)	5.63 × 10 -2	(4.81 %)
	∂ P f /∂ ξ 4		1.36 × 10 -1	(0.525 %) 1.35 × 10 -1	(6.45 %)	1.34 × 10 -1	(6.22 %)
	∂ P f /∂ µ Θ 2	a	3.82 × 10 -1	-	3.82 × 10 -1	-	3.82 × 10 -1	-
	∂ P f /∂ σ Θ 2		2.77 × 10 -1	-	2.75 × 10 -1	-	2.71 × 10 -1	-
	ν ARA/AIS		-	-	14	-	17	-

a cf. Eqs. (15a) and (15b) for the uniform case.

Table 7 :

 7 Results for Example #2 considering the influence of the failure event rareness.

					ARA/AIS-CE ( )	ARA/NAIS ( )
					(N x,θ = 10 4 samples/step) (N x,θ = 10 4 samples/step)
					Estimate		cv	Estimate	cv
		P f			6.71 × 10 -6	(6.25 %)	6.82 × 10 -6	(5.97 %)
		∂ P f /∂ ξ 1		-4.95 × 10 -6 (6.58 %)	-4.98 × 10 -6 (5.56 %)
		∂ P f /∂ ξ 2		2.41 × 10 -4	(6.25 %)	2.45 × 10 -4	(5.97 %)
		∂ P f /∂ ξ 3		1.62 × 10 -3	(6.51 %)	1.65 × 10 -3	(6.14 %)
		∂ P f /∂ ξ 4		4.63 × 10 -3	(6.74 %)	4.59 × 10 -3	(9.01 %)
		∂ P f /∂ µ Θ 2	a	1.25 × 10 -2	-	1.25 × 10 -2	-
		∂ P f /∂ σ Θ 2		1.04 × 10 -2	-	1.02 × 10 -2	-
		ν ARA/AIS			231		-	276	-
		a cf. Eqs. (15a) and (15b) for the uniform case.
		0.2	Sensitivities of Pf w.r.t. hyper-parameters
		0.18					∂ Pf /∂ξ1 ∂ Pf /∂ξ2
		0.16					∂ Pf /∂ξ3 ∂ Pf /∂ξ4
		0.14
	Sensitivities	0.08 0.1 0.12
		0.06
		0.04
		0.02
		0
		1e3		2e3	3e3	5e3	7e3	1e4
				Number of samples per step

Table 8 :

 8 Input probabilistic model for Example #3.

Table 9 :

 9 Results for Example #3.

		ARA/CMC	ARA/AIS-CE	ARA/NAIS
		(N x,θ = 10 6 samples)	(N x,θ = 10 4 samples/threshold) (N x,θ = 10 4 samples/threshold)
		Estimate	cv	Estimate	cv	Estimate	cv
	P f	4.40 × 10 -3	(1.38 %) 4.41 × 10 -3	(10.3 %)	4.40 × 10 -3	(2.08 %)
	∂ P f /∂ ξ 1 -9.13 × 10 -4 (3.44 %) -8.68 × 10 -4 (27.7 %)	-9.12 × 10 -4 (5.90 %)
	∂ P f /∂ ξ 2 2.95 × 10 -3	(2.32 %) 3.02 × 10 -3	(14.8 %)	2.95 × 10 -3	(3.22 %)
	∂ P f /∂ ξ 3 -2.31	(3.88 %) -2.29	(25.7 %)	-2.30	(5.82 %)
	∂ P f /∂ ξ 4 6.43	(2.18 %) 6.26	(14.6 %)	6.41	(3.77 %)
	ν ARA/AIS	-	-	0.5	-	13	-

Table 10 :

 10 Results for Example #3 considering the influence of the failure event rareness.

						ARA/AIS-CE ( )		ARA/NAIS ( )
					(N x,θ = 10 4 samples/step) (N x,θ = 10 4 samples/step)
					Estimate		cv	Estimate	cv
			P f		1.00 × 10 -4	(29.7 %)	1.19 × 10 -4	(2.85 %)
			∂ P f /∂ ξ 4 2.51 × 10 -1	(25.2 %)	3.10 × 10 -1	(4.24 %)
			ν ARA/AIS	2			-	207	-
		8	Sensitivities of Pf w.r.t. hyper-parameters				0.4	Sensitivities of Pf w.r.t. hyper-parameters
		6								0.3
							∂ Pf /∂ξ1
							∂ Pf /∂ξ2						∂ Pf /∂ξ1
	Sensitivities	2 4					∂ Pf /∂ξ3 ∂ Pf /∂ξ4		Sensitivities	0.2 0.1			∂ Pf /∂ξ2 ∂ Pf /∂ξ3 ∂ Pf /∂ξ4
		0								0
		-2								-0.1
		-4								-0.2
		1e3	2e3	3e3	5e3	7e3	1e4			1e3	2e3	3e3	5e3	7e3	1e4
			Number of samples per step					Number of samples per step
	(a) Estimated sensitivities (error bars) vs. reference results obtained
	by ARA/CMC (dashed lines).

∂ P f /∂ ξ 1 -4.65 × 10 -5 (46.9 %) -3.66 × 10 -5 (7.00 %)

∂ P f /∂

ξ 2 1.21 × 10 -4 (40.3 %) 1.41 × 10 -4 (3.62 %) ∂ P f /∂ ξ 3 -1.19 × 10 -1 (34.4 %) -9.18 × 10 -2 (7.67 %)

Table 11 :

 11 Different formulas for elasticities.

	Formula	Hyper-parameters	Source
	(f1) e j =	∂ P f ∂ ξ j ×	σ Θ j P f

the one required by applying a finite difference method (FDM). However, using FDM can be difficult for several reasons: 206

Acknowledgements

The authors would primarily like to thank the editor and the anonymous reviewers for the constructive remarks and valuable suggestions. The first author is enrolled in a PhD program co-funded by ONERA -The French Aerospace Lab and SIGMA Clermont. Their financial supports are gratefully acknowledged.


Eminent domains. Boris Pétric. Le Monde diplomatique (English edition). Le Monde diplomatique. English edition. July, 2014

China has this year become the world's biggest purchaser of Bordeaux wines. It consumes as yet just over a litre of wine per person annually, compared with 47 litres in France and 37 in Italy (1), so there's plenty of room for growth, for foreign exporters and local producers. The area devoted to vineyards in China has doubled over the past 15 years to 600,000 hectares (2), while Europe is losing vineyards -it now has 3.5m hectares, 800,000 of them in France.

Although grapes have been grown in China for 2,000 years, volume winemaking only started in the 1980s. Until 1990, wine was produced solely in Hebei Province (around Beijing), Shandong Province and the Xinjiang Uyghur Autonomous Region, and was controlled by a few major state enterprises -Changyu, Dragon Seal, Great Wall and Suntime -which still dominate the sector. In the early 2000s, China developed highly successful joint ventures with foreign firms, and agreements were signed with multinationals including Miguel Torres, Domecq, Pernod Ricard and Castel. But while China has opened up its economy, it still holds that agricultural land is inalienable property, and access to it only possible through long-term leases contracted with the government or with local partners.

Joint ventures provided China with knowhow, funds and technology, while limiting imports. The French government invested 2m in establishing a vineyard in Hebei Province, a project later abandoned, although while it ran it formed partnerships, exchanged knowledge about types of wine, and trained the first Chinese wine specialists, such as Li Demei, who is now a leading consultant oenologist. French companies that invested had mixed results. Pernod Ricard withdrew from a complicated partnership, but has kept a presence in the Ningxia Hui Autonomous Region with its Helan Mountains domain. Castel still has a partnership in Shandong Province.

Since 2000, wine growing has been integrated into a government plan to develop western China, which has fallen behind in the economic boom. The Xinjiang, Ningxia and Inner Mongolia autonomous regions and Shanxi Province have become priority investment regions. Local authorities hope that winegrowing will help limit the rural exodus and create jobs. New vineyards, some very small, have sprung up with the ambition to produce quality wines. Irrigation, especially near the Yellow River, makes viticulture possible in desert regions with harsh, freezing winters. (The vines are buried during the winter to protect them from the cold.) Ecologically, the vines help to anchor the sand and therefore reduce the sandstorms that hit China's great cities.

Vineyards in resort complexes

For the past decade, China's nouveaux riches have hired wine and agriculture specialists, and French architects, to set up French-style vineyards -they want to gain AOC status. The appellation d'origine controlée guarantees the origin of a wine and ensures standards in production and alcohol content. Domaines Barons de Rothschild (DBR) has established vineyards in the Penglai region (Shandong Province) in a partnership with China International Trust and Investment Corporation. LVMH Moët Hennessy Louis Vuitton is in Yunnan Province in the south, where it has registered the brand name Shangri-La in partnership with Vats Group.

Wine tourism venues, like the replica of Saint-Emilion in Dalian, Liaoning Province, cater to the urban middle class and provide investment opportunities for promoters. In these vast resort complexes, the vineyards are secondary to château-style hotels, golf courses and tennis courts. We were told: "Planting a vineyard gives you the right to have a property lease on land. Businessmen often don't care about the vineyard. With the anticipated urban growth, the vines can be pulled up and they will still have the land rights, which is a big advantage."

For the government, wine symbolises that China is now part of a globalised economy. It wants Chinese wines to have cultural dimensions -wine must be a national brand in the worldwide battle of national brands. A grower in Ningxia said: "I want my wine to become the pride of China. We are becoming capable of producing great wines and we are going to prove it to the whole world." Investment in wine is justified as contributing to China's collective ambition, which leaves little room for individual dreams, although young Chinese oenologists are making visible careers, such as Emma Gao, chief winemaker at the Silver Heights domain and Zhang Jing of the Jiabelan domain, both in Ningxia, distinguished in Revue du vin de France and Decanter .

There is a major obstacle -lack of suitable land. Consultants are searching China's massive landmass for ideal locations because, as expert Gérard Colin said, "the three main winegrowing regions, Hebei, Shandong and Xinjiang, have no future. The problem in Xinjiang is that burying the vines every winter is very expensive in an area with a shortage of manpower. I won't even mention the problems of irrigation and water shortages. In Hebei, the cities are encroaching at an impressive pace and farmers prefer to leave the land and work on building sites, where they earn more money." In Shandong, according to Bruno Paumard, "very high heat and humidity trigger many diseases and there's a limit to what you can do with chemicals... we have yet to find the ideal place to make wine." Despite that, in 2012 China became the fifth largest producer in the world and should rapidly make it to the top three [START_REF]France, Italy, Spain and the US still dominate the market[END_REF].

Agrifood groups buy abroad

China's major agrifood groups are acquiring vineyards abroad. Trading giant Cofco owns Château de Viaud in Bordeaux as well as the 800 hectares of the Bodega Bisquertt in Chile (which has a trading agreement with China) and is trying to acquire vineyards in Australia. Chinese companies are also buying up foreign brands to distribute them directly at home. Wang Quang owns Bellefont-Belcier, a major Saint-Emilion grand cru classé. Chinese purchasers have bought 50 domains in the Bordeaux area, in some cases to protect their wealth from political reversals in China. Although there are more than 7,000 vineyards in Bordeaux, those acquisitions caused an outcry in France. Whether Chinese companies abroad or foreign companies in China, the industry is rapidly becoming transnational (4).

In China, the increase in production is primarily to satisfy growing domestic demand -more than 80% of which is for Chinese wine (5). There is a difference between what people drink and what they buy, because wine has become part of the ubiquitous gift-giving culture. Most wines tasted in cellars cannot be bought on the market. A cellar master told us: "We have no problem selling our wines because my boss has connections with the local authorities and we give them a large amount of our wine. Then the party bosses oblige companies to buy our production, so our bottles end up as gifts at political banquets. This is not a supply and demand market... our entire economy functions like that. Getting into the real consumer market is another matter." The interdependence of companies creates mutual obligations to purchase and beverages are no exception to the dominant logic of guanxi (relationship networks). Very often the bottles end up as status symbols in the homes of local cadres, never to be opened.

Wine culture has been encouraged as a social practice. The government has run media campaigns demonstrating the health benefits of wine in the hope of reducing the consumption of strong drinks such as baiju, a grain-based spirit. Television reports stigmatise Party members seen drunk at banquets where the baiju flows, and the government would prefer land devoted to its production to be used for growing food. Yet social drinking is a prerequisite among the political and business classes for forming relationships of trust. "If you don't drink it will be taken as a lack of respect for your business partner," said a French wine trader. "You have to get completely drunk and share a moment of intimacy before you sign a contract. Things are going well for me, but my liver is paying the price." "We don't drink wine, it's too expensive. We prefer baiju, " said a peasant working in a vineyard in Hebei Province. Baiju is still popular and China's 600 million working-class men and women will always prefer it, or beer, to wine. But among the wealthy, wine confers social distinction, and foreign wines, especially French, are prestigious. A young privileged, urban elite are adopting different social practices in a relatively new mixed-gender context. They distrust domestic food and drink, because of China's food safety scares, and prefer foreign produce.

The new government's anti-corruption campaign has had an effect since last year on the gift-giving economy, slowing imports of Bordeaux. Nevertheless, wine is becoming accepted by Chinese society. This year the Chinese became the world's leading consumers of red wine, ahead of the French, with more than 1.8bn bottles sold (6).

Notes :

Boris Pétric is an anthropologist and works for the Centre National de la Recherche Scientifique (CNRS).


Torbjørn Cunis

Comments on the computation of multiple Lyapunov-like functions for switched hybrid systems

Introduction

Switched hybrid systems representing a wide range of physical and cyber-physical systems, the analysis of hybrid stability is of particular interest yet the more difficult. The use of multiple Lyapunov functions here is well known [START_REF] Michael | Multiple Lyapunov functions and other analysis tools for switched and hybrid systems[END_REF]. In the work of [Zheng et al . 2], the authors aim to iteratively compute multiple Lyapunov-like functions to estimate the domain of attraction. Notwithstanding, the approach bears formal and conceptual shortcomings, which we shall present in the following. In consequence, as we argue, the proposed algorithm might not yield a correct estimate. Notation. The state-space in mode i be D i ; the Lyapunov-function candidates V i : R n → R; and the estimated domain of attraction D. The set of sum of squares polynomials is denoted as Σ [x] and the zero-polynomial as 0(x) ≡ 0.

Corrections

Increase of the estimate

In order to properly enlarge the estimate from iteration k to k + 1, the authors require [2, Eq. 3]

D k D k+1 , (1)
where

D k = i (D k,i ∩ D i ), and therefore check [2, cf. Eq. 8] D k,i D k+1,i (2)
for all i. However, (2) does not imply (1) in general, which can be easily seen:

Counterexample. Let D k+1,i = D k,i ∪ {x i } for all i with x i ∈ D k,i and suppose x i ∈ D i ; then, i (D k+1,i ∩ D i ) = i (D k,i ∩ D i ) 1 although D k,i D k+1,i .
Even we assumed x i ∈ D i -D k,i , this would not guarantee x i ∈ D k as the statespaces D i are not pairwise disjunct, thus disproving the implication. Whether or not the implication might still hold for the assumptions taken, has not been proven here. The iteration might therefore still yield a larger estimate. The restriction to D k+1,i = {x |V k+1,i (x) ≤ 1 } as partial estimate further unnecessarily neglects possibly stable level sets of V k+1,i (•) larger than 1. It is finally questionable whether the iteration will actually approximate the "true" domain of attraction, since a set D which is by an arbitrary measure larger than D k does not necessarily imply D k ⊂ D .

Continuity of the Lyapunov-functions

Multiple Lyapunov-function candidates must be continuous along the switching surfaces, here SS i,j ⊂ S i,j = r {x |h i,j,r x = 0 } for i = j. The constraint [2, Eq. 6]

S i,j ⊆ {x |V i (x) ≥ V j (x) } (3)
is then supposed to be encoded as [2, Eq. 11]

-s i,j (x) r h i,j,r x -(V j (x) -V i (x)) ∈ Σ [x] (4)
and s i,j (•) ∈ Σ [x] for all i = j. Yet, (4) actually implies [3, Lemma 2]

x r h i,j,r x ≥ 0 ⊆ {x |V i (x) -V j (x) ≥ 0 } , (5)
which is a much stronger constraint then (3).

Continuity along the boundary

Finally, the authors require [2, Eq. 7]

SS i,j ∩ {x |V i (x) = 1 } ⊆ {x |V j (x) = 1 } , (6)
which they encode as [2, Eq. 12]

-s 1,i,j (x)

r h i,j,r x -s 2,i,j (x) (1 -V i (x)) ± s 3,i,j (x) (1 -V j (x)) ∈ Σ [x] (7)
and

s 1,i,j (•) , s 2,i,j (•) , s 3,i,j (•) ∈ Σ [x].
Eq. ( 7), we are afraid to say, does not imply (6) at all. We recall the primal lemma [3, Lemma 2]: Lemma 1. Let g 0 , g 1 , . . . , g m be polynomials in x; if there are ς

1 , . . . , ς m ∈ Σ [x] such that g 0 (x) - m i=1 ς i (x) g i (x) ∈ Σ [x]
, then {x |g 1 (x) , . . . , g m (x) ≥ 0 } ⊆ {x |g 0 (x) ≥ 0 }.

In the case of (7), we have g 0 = 0(x) and, subsequently, {x |g 0 (x) ≥ 0 } = R n . Indeed, s 1,i,j = s 2,i,j = s 3,i,j = 0(x) satisfy (7) for any h i,j,r , V i (•), and V j (•). That is, the constraint in (7) is a tautology.

Conclusions

The estimation of the domain of attraction of a switched hybrid system does not work as the authors [2] suppose. With (1) not proven, an increase of D k+1 is possibly not enforced. The effectively imposed constraint of ( 5) is much stricter than suggested. In any case, as (7) is void, the output is not necessarily a domain of attraction in the sense of the author's first lemma.