Emmanuel Vazquez

Julien Bect

A sequential Bayesian algorithm to estimate a probability of failure

Keywords: Sequential Bayesian planning, Design of experiments, Kriging, Probability of failure

This paper deals with the problem of estimating the probability of failure of a system, in the challenging case where only an expensive-to-simulate model is available. In this context, the budget for simulations is usually severely limited and therefore classical Monte Carlo methods ought to be avoided. We present a new strategy to address this problem, in the framework of sequential Bayesian planning. The method uses kriging to compute an approximation of the probability of failure, and selects the next simulation to be conducted so as to reduce the mean square error of estimation. By way of illustration, we estimate the probability of failure of a control strategy in the presence of uncertainty about the parameters of the plant.

INTRODUCTION

Uncertainty on the parameters of a system have to be taken into account during its design to insure a desired level of reliability. In particular, it is important to estimate the probability of the system to work under abnormal or dangerous operating conditions due to unknown random dispersions of its parameters. The probability of failure of a system is usually expressed as the probability of the excursion set of a function of performance above a fixed threshold. More precisely, let f be a measurable real function defined over a probability space (X, X , P X ) and let u ∈ R be a threshold. The problem to be considered in this paper is the estimation of the volume under P X of the excursion set A(f ) := {x ∈ X : f (x) ≥ u} (1) of the function f above the level u. In the context of robust design, the volume α(f ) := P X (A(f )) can be viewed as the probability of failure of a system. Then, X plays the role of an input or factor space, the probability P X models the uncertainty on the inputs or the factors, and f is some function of performance of the outputs of the system, whose evaluation may involve complex and timeconsuming simulations or fabrications of prototypes. When f is expensive to evaluate, the estimation of α(f ) must be carried out with a restricted number of evaluations of f , generally excluding the estimation of the probability of excursion by a Monte Carlo approach. Indeed, consider the empirical estimator

α m (f ) := 1 m m i=1 1 {f (Xi)≥u} , (2)
where the X i s are independent random variables with distribution P X . According to the strong law of large numbers, the estimator α m (f ) converges to α(f ) almost surely when m increases. Moreover, it is an unbiased estimator of α(f

), i.e. E(α m (f )) = α(f ). Its mean square error is E (α m (f ) -α(f )) 2 = 1 m α(f ) 1 -α(f ) .
If the probability of failure α(f ) is small, then the standard deviation of α m (f ) is approximately (α(f )/m) 1/2 . To achieve a given standard deviation δα(f ) thus requires approximately 1/(δ 2 α(f )) evaluations, which can be prohibitively high if δ is small. Of course, many importance sampling methods have been proposed to improve over the basic Monte Carlo convergence rate -for instance, methods based on cross-entropy [START_REF] Rubinstein | The cross-entropy method for combinatorial and continuous optimization[END_REF], or subset sampling [START_REF] Au | Estimation of small failure probabilities in high dimensions by subset simulation[END_REF]. They will not be considered here for the sake of brevity and because the required number of function evaluations is still very high.

Methods to estimate a probability of failure that do not require a large number of evaluations of f are usually based on an approximation of the frontier ∂A(f ) of A(f ) by a simple geometrical shape. Compared to Monte Carlo methods, this class of methods does achieve a large reduction in the required number of system evaluations, but the accuracy of the estimator depends on the actual shape of ∂A(f ) and its resemblance to the approximating shape. The popular first-order reliability method (FORM) uses an hyperplane as an approximation of ∂A(f ) and the second-order reliability method (SORM) uses a paraboloid (see, e.g., [START_REF] Bjerager | On computational methods for structural reliability analysis[END_REF]. Nevertheless, ∂A(f ) may depart significantly from such elementary shapes. In fact, this class of methods does not provide statistically consistent estimators of the probability of failure.

This paper presents a new approach. The method proposed, which is outlined in Section 2, is in essence a Bayesian sequential search algorithm similar in spirit to Bayesian algorithms for global optimization [START_REF] Zilinskas | A review of statistical models for global optimization[END_REF][START_REF] Mockus | Application of Bayesian approach to numerical methods of global and stochastic optimization[END_REF] such as the Expected Improvement algorithm [START_REF] Schonlau | Global optimization with nonparametric function fitting[END_REF] and the IAGO algorithm Villemonteix et al. (to appear.). Here, we transpose the main ideas and propose an algorithm to choose sequentially the evaluations of f to produce a consistent and rapidly converging estimation of α(f ). This algorithm is detailed in Section 3. Section 4 will illustrate the potential evaluation savings of the methodology on a simple but archetypal example, in which the probability of failure of a control strategy of an uncertain plant is estimated.

SEQUENTIAL BAYESIAN SEARCH ALGORITHM

Sequential Bayesian optimal design

Let (G n ) be any given sequence of functions,

G n : (X× R) n → R + , such that G n ((x 1 , f (x 1 )), . . . , (x n , f (x n ))
) yields an approximation of α(f ). Our goal is to choose locations X n = X n (f ) sequentially in such a way that the approximations

n := G n (X 1 , f (X 1
)), . . . , (X n , f (X n )) will converge as fast as possible to the true value α(f ). By doing so, it is hoped that a reasonably good approximation of α(f ) will be obtained using a limited number of evaluations. Obviously, in such a procedure, the n th location X n = X n (f ) is only allowed to depend on the unknown function f through the past evaluations f (X i ), 1 ≤ i < n. For the sake of clarity, the explicit dependence on f will be dropped in the rest of the paper.

We formulate the problem in the framework of sequential Bayesian optimal design, where the information from the past evaluations can be combined with prior information on the system in order to take the "best" sequence of decisions (in a sense to be defined later). Let F = R X be the space of all real functions X → R and denote by U the cylinder σ-algebra on F . Prior information on the unkown function f is taken into account in the Bayesian approach under the form of a probability measure P on (F, U). This amounts to considering f as a random function, the distribution of which is given by the prior probability P before evaluations are made, and by the posterior probability P { • |F n } after n evaluations, where F n is the σ-algebra generated by f (X 1 ), . . . , f (X n ). For each n ≥ 0, the location X n+1 is assumed to be chosen measurably from the past evaluations (in other words, X n+1 is an F n -measurable random variable). It must be pointed out that despite the semantic shift we are still dealing with deterministic algorithms: randomness only comes from the fact that we are now considering f as a random function.

Given a cost function C : R + × R + → R and a budget of N evaluations of f , a strategy X = (X 1 , . . . , X N ) is said to be optimal if it minimizes the Bayes risk E (C (α(f ), α N )) associated to the prior P. It is well known [START_REF] Berger | Statistical Decision Theory and Bayesian Analysis[END_REF] that, at least in principle, the optimal strategy can be found using a form of dynamic programming. In practice, finding the optimal strategy is usually not feasible, especially in problems like the one at hand where there is a continuum of possible states and decisions at each stage n. In this paper we handle the difficulty with the following suboptimal greedy strategy (sometimes called a one-step look-ahead procedure) :

X n = argmin xn∈X E C (α(f ), α n ) | F n-1 , (3)
for all n ∈ {1, . . . , N }, where the predictor α n is computed using the past evaluations f (X i ), 1 ≤ i < n and the uncertain outcome f (x n ) of an evaluation at x n .

The reader is referred to [START_REF] Chaloner | Bayesian experimental design: A review[END_REF] for a survey on both sequential and non-sequential Bayesian experimental design (from a statistical point of view) and to [START_REF] Pronzato | Optimal experimental design and some related control problems[END_REF] for a survey on the connections between the fields of design of experiments and automatic control.

Gaussian process priors and kriging

The proposed method will now be specified following the Bayesian approach explained in Section 2.1, using a Gaussian process prior as the probability P. The idea of modelling an unknown deterministic function f by a stochastic process has originally been introduced in the field of geostatistics (see [START_REF] Chilès | Geostatistics: Modeling Spatial Uncertainty[END_REF], and the references therein), and also plays a central role in the design and analysis of "computer experiments" (see, e.g., [START_REF] Sacks | Design and analysis of computer experiments[END_REF][START_REF] Currin | Bayesian prediction of deterministic functions, with applications to the design and analysis of computer experiments[END_REF][START_REF] Welch | Screening, predicting and computer experiments[END_REF].

The distribution P of a Gaussian process f is uniquely determined by its mean m(x) := E(f (x)) and its covariance function k(x, y) := E ((f (x)m(x))(f (y)m(y))).

Restricting our framework to Gaussian processes makes it possible to deal with the conditional distributions P { • |F n } that appear in equation (3). Another benefit of using such a prior is that the posterior mean

f n (x) := E f (x) | F n , (4
) also called the kriging predictor of f (x), can easily be computed using linear algebra (see, e.g., [START_REF] Chilès | Geostatistics: Modeling Spatial Uncertainty[END_REF] and provides a cheap surrogate model for the true function f . Computing the probability of failure α(f n ) of the surrogate model is a natural choice for the approximations G n of subsection 2.1, which leads us to set α n = α(f n ). Another important benefit of using Gaussian process priors is that the posterior variance σ 2 n (x) is also easily computed, which allows to provide prediction intervals.

Practically, the probability α(f n ) can be approximated by a Monte Carlo estimator such as α m (f n ), which requires m evaluations of f n . Since the evaluation of the surrogate model is much faster than the evaluation of the true model, m can be chosen large enough to make α m (f n )α(f n ) negligible with respect to α(f n )α(f ). Thus, we shall only focus on the convergence of α(f n ) to α(f ) and assume that α(f n ) can be computed with the required accuracy in a reasonable time.

Finally, choosing the quadratic cost function C(α, α ′ ) = (α-α ′ ) 2 , the proposed strategy is to sequentially minimize the mean square error of approximation:

Program 1 X n = argmin xn∈X Υ n (x n ) := E (α(f ) -α( fn )) 2 | F n-1 .
Here and subsequently, we denote by fn (x) the kriging predictor of f (x) based on f (X 1 ), . . . , f (X n-1 ) and f (x n ). This strategy can be seen as an instance of the so-called stepwise uncertainty reduction principle [START_REF] Geman | An active testing model for tracking roads in satellite images[END_REF]Villemonteix et al., to appear.). The criterion Υ n (x n ) provides a measure of the expected approximation error if a new evaluation is performed at x n . It depends of course on the past evaluations f (X 1 ), . . . , f (X n-1 ), since it is defined as a conditional expectation with respect to F n , but this depency is kept implicit for simplicity of notation. It must be noted that the actual result f (x n ) of the additional evalution at x n is integrated out by the conditional expectation with respect to F n-1 .

TOWARD A NUMERICALLY TRACTABLE SEQUENTIAL ALGORITHM

This section presents a numerically tractable approximation of Program 1, proceeding in two steps. In the first step an upper bound Υ ′ n on the criterion Υ n is derived. This upper bound is simpler than Υ n but still not in closed form. In the second step we provide a computationally efficient procedure to approximate Υ ′ n , the result of which is our final criterion Υ ′′ n .

Upper bound on the criterion Υ n

Using the chain rule for conditional expectations, observe that Υ n (x n ) can also be written as

Υ n (x n ) = E Γ n (x n , f (x n )) 2 F n-1 , (5)
where the function Γ n is defined by

Γ n (x n , z n ) := E n α(f ) -α( fn ) 2 1/2 .
For abbreviation, the notation E n has been introduced to denote the conditional expectation operator 5) is simply a one-dimensional integral with respect to the value of f (x n ). Unfortunately, Γ n (x n , z n ) has no analytical expression. A numerical approximation could in principle be obtained using a Monte Carlo approach, but would require a very large number of conditional simulations (see Chilès and Delfiner, 1999, Chap. 7) of the random function f . Instead, we choose to replace Γ n (x n , z n ) by an upperbound which will be easier to compute and minimize. Remark that

E n := E{ • F n-1 , f (x n ) = z n }. Note that Equation (
α(f ) -α( fn ) = X 1 f (x)>u -1 fn(x)>u P X (dx).
An upper bound is then readily provided by the generalized Minkowski inequality (see, e.g., Vestrup, 2003, Section 10.7):

Γ n (x n , z n ) ≤ Γ ′ n (x n , z n ) := X υ n (x; x n , z n ) 1 2 P X (dx) , with υ n (x; x n , z n ) := E n (1 f (x)>u -1 fn(x)>u ) 2 .
Replacing Γ n by Γ ′ n in Program 1 yields the following modified strategy:

Program 2 X n = argmin xn∈X Υ ′ n (x n ) := E Γ ′ n (x n , f (x n )) 2 F n-1 .
This strategy is indeed simpler than the previous one, since υ n (x; x n , f (x n )) is given by

υ n (x; x n , f (x n )) = Ψ u -fn (x) σn (x) ,
where Ψ is the complementary cumulative distribution function of the normal distribution and σn (x) the standard deviation of fn (x) conditionally to F n-1 and f (x n ).

An informal interpretation of Program 2 is that X n minimizes, conditionally to the past evaluations, the expected posterior error of prediction of 1 f (x)>u by 1 fn(x)>u , averaged over X (using the distribution P X ).

Approximation of the upper bound

In the following paragraphs, we show that an approximate version of Program 2 can be solved with an acceptable computational complexity.

A first approximation is required for the computation of Γ ′ n . Assuming that the distribution P X is easy to sample from, we propose to use a simple Monte Carlo approximation for the integral with respect to P X :

Γ ′ n (x n , z n ) ≈ Γ ′′ n (x n , z n ) := 1 m m i=1 υ n (Y i ; x n , z n ) 1/2 ,
where (Y i ) i≥1 is a sequence of independent and identically distributed variables with common distribution P X .

Another approximation is required for the computation of

Υ ′ n (x n ) ≈ E Γ ′′ n (x n , f (x n )) 2 F n-1 . (6)
Using a tagged partition (∆ j n , z j n ), 1 ≤ j ≤ Q of R, we replace the random variable f (x n ) in the right-hand side of ( 6) by the quantized variable ξ n := j z j n 1 f (xn)∈∆ j n . The proposed algorithm is thus Program 3

X n = argmin xn∈X Υ ′′ n (x n ) := E Γ ′′ n (x n , ξ n ) 2 F n-1 .
Practically, the design criterion is computed as

Υ ′′ n (x n ) = Q j=1 Γ ′′ (x n , z j n ) P f (x n ) ∈ ∆ j n F n-1 , using the fact that f (x n ) is normally distributed condition- naly to F n-1 , with mean f n-1 (x n ) and variance σ 2 n-1 (x n ).

NUMERICAL EXAMPLES

Illustration

This section provides a compendious one-dimensional illustration of the algorithm. We wish to estimate α(f ), where f (x) is a given function defined over R, which is endowed with the probability distribution P X = N (0, σ 2 ), σ = 0.4 (see Figure 1). After a few iterations, the unknown function f has been sampled so that the probability of excursion P({f : f (x) > u} | F n ) is either 0 or 1 in the region where the density of P X is high (as shown on Figure 1). This example illustrates the effectiveness of the proposed algorithm. Note that in practice, a parametrized covariance has to be chosen for specifying P and its parameters should be estimated from the data, using, for instance, a maximum likelihood approach (e.g. [START_REF] Stein | Interpolation of Spatial Data: Some Theory for Kriging[END_REF]. for the Monte Carlo integration and Q = 20 for the tagged partition; kriging approximation f n (thick line); 95% confidence intervals computed from the kriging variance (dashed lines). Middle: probability of excursion (solid line); probability density of P X (dotted line). Bottom: graph of Υ ′′ (x i ), i = 1, . . . , m = 800, the minimum of which indicates where the next evaluation of f should be done (i.e., at approximately 0.75).

Probability of failure of a regulated system

A typical example for which our method is particularly relevant is the estimation of a probability of failure of a system described by partial differential equations with complex boundary conditions and a few uncertain parameters.

We choose, however, to consider a much simpler illustrative problem, for two reasons. First, it is possible, within the space allowed, to give enough details to allow the reader to use it and to compare the performance of our approach with other methods not considered here. Second, nothing is lost by considering such an example, as the methodology would be strictly the same for a more expensive-to-simulate model.

We consider a simple continuous-time dynamical system, corresponding to a two-compartment model (Figure 2).

x 1

x 2 q 1 q 2

x 3 u Fig. 2. A two-compartment model.

The state variables q 1 and q 2 represent an amount of material in each compartment, the positive real numbers x 1 , x 2 , x 3 correspond to exchange rates and u is the control input of the system. The state vector is governed by the ordinary differential equation q1

q2 = -x 1 -x 3 x 2 x 1 -x 2 q 1 q 2 + u 1 0 .
We assume that x 3 is known and equal to 0.35. However, we suppose that the parameters x 1 and x 2 are uncertain. We model this uncertainty by a probability distribution P X centered on (x 1 , x 2 ) = (0.6, 0.15) (see Figure 3).

Our objective is to regulate the output variable q 2 , by means of a proportional-integral controller, i.e. a linear feedback corrector whose transfer function can be written as G(s) = K a + s s .

The value of K and a of the corrector were chosen to insure satisfactory phase and gain margins for (x 1 , x 2 ) = (0.6, 0.15). We obtained K = 0.93 and a = 0.23. To measure the performance of the system, we consider the settling time f (x 1 , x 2 ) of the step response of systemthe response is considered to be settled at time s when |q 2 (t)lim t→∞ q 2 (t)| becomes smaller than a fraction equal to 0.02 of its final value, ∀t > s. For instance, we have f (0.6, 0.15) = 10.95. We would like to estimate the probability α(f ) = P X (A(f )) of the excursion set A(f ) = {(x 1 , x 2 ) ∈ R 2 : f (x 1 , x 2 ) ≥ u = 12.5} , i.e. the probability that the system has a settling time greater than 12.5, given the probability distribution P X on the uncertain parameters. The result of our sequential algorithm is presented on Figure 3. After n = 40 iterations, the unknown settling-time function f has been sampled so that the probability of excursion P({f : f (x 1 , x 2 ) > u} | F n ) is very close to either 0 or 1 in the region where the density of P X is high.

Fig. 1 .

 1 Fig. 1. Top: threshold u = 1 (horizontal solid line);function f (thin line); n=10 evaluations (squares) as obtained by the proposed algorithm using m = 800 for the Monte Carlo integration and Q = 20 for the tagged partition; kriging approximation f n (thick line); 95% confidence intervals computed from the kriging variance (dashed lines). Middle: probability of excursion (solid line); probability density of P X (dotted line). Bottom: graph of Υ ′′ (x i ), i = 1, . . . , m = 800, the minimum of which indicates where the next evaluation of f should be done (i.e., at approximately 0.75).

Fig. 3 .

 3 Fig. 3. Top: gray-level plot of the probability density of P X along the parameters x 1 and x 2 of the twocompartment model. The bright areas correspond to high values of the density. Bottom: gray-level plot of the probability of excursion conditioned on the evaluations of the settling time, whose locations are indicated by square dots. The solid line indicates the frontier ∂A of the failure domain. One can see that the probability of excursion is very close to either 0 or 1 in the region where the density of P X is high.

ACKNOWLEDGEMENTS

This work was partially supported by the French Agence Nationale de la Recherche in the context of the opuscoprin project (ref. anr-07-cis7-010 and anr-07-tlog-015).