{'fraction_non_alphanumeric': 0.0403020099989797, 'fraction_numerical': 0.012549739822467095, 'mean_word_length': 4.56726308813784, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
CAMERA: A Method for Cost-aware, Adaptive, Multifidelity, Efficient Reliability Analysis

Estimating probability of failure in aerospace systems is a critical requirement for flight certification and qualification. Failure probability estimation involves resolving tails of probability distribution, and Monte Carlo sampling methods are intractable when expensive high-fidelity simulations have to be queried. We propose a method to use models of multiple fidelities that trade accuracy for computational efficiency. Specifically, we propose the use of multifidelity Gaussian process models to efficiently fuse models at multiple fidelity, thereby offering a cheap surrogate model that emulates the original model at all fidelities. Furthermore, we propose a novel sequential \emph{acquisition function}-based experiment design framework that can automatically select samples from appropriate fidelity models to make predictions about quantities of interest in the highest fidelity. We use our proposed approach in an importance sampling setting and demonstrate our method on the failure level set estimation and probability estimation on synthetic test functions as well as two real-world applications, namely, the reliability analysis of a gas turbine engine blade using a finite element method and a transonic aerodynamic wing test case using Reynolds-averaged Navier--Stokes equations. We demonstrate that our method predicts the failure boundary and probability more accurately and computationally efficiently while using varying fidelity models compared with using just a single expensive high-fidelity model.


Introduction
Flight certification by analysis is a new paradigm in which airworthiness certificates for aircraft and handling qualities compliance for aircraft subsystems can be obtained via calculations with quantifiable accuracy. 1 This opens up possibilities for leveraging physics-based high-fidelity models in aircraft design, certification, and qualification. However, the high computational costs associated with high-fidelity models necessitates the development of advanced mathematical methods that judiciously spend computational resources on problems where such models may have to be queried many (≥ O(10 2 ) − O(10 6 )) times.
One such problem we are interested in is reliability analysis [36]. Specifically, we want to estimate probabilities of failure in aerospace subsystems, whose performance are governed by models that can be queried by perturbing some control parameter. Rigorous quantification of the probability of subsystem failure is a critical requirement for flight certification. Failure probability estimation (FPE) involves accurately resolving failure boundaries (contours or level sets) in the design space and tails of probability distributions. The popular Monte Carlo (MC) [15] sampling-based methods are intractable when expensive high-fidelity simulations have to be queried. Furthermore, the high dimensionality of the control parameters exacerbates the tractability issue. In this work we focus on a method that circumvents direct MC and instead seeks to develop a surrogate model of the mapping between the input parameters and the quantity of interest. Our hope is that a surrogate model can be constructed with far fewer queries to the expensive high-fidelity model, following which we use conventional variance-reducing sampling-based techniques for FPE.
Computational models typically include one or more customizable fidelity parameters that allow trading model accuracy for a gain in computational speedup. For example, the exact Navier-Stokes equations simulating fluid flows require numerical mesh resolutions on the order of the fine Kolmogorov scale, which shares an inverse nonlinear power relationship to the flow Reynolds number [58], resulting in a large number of degrees of freedom that in turn makes the computational cost of evaluating such models intractable. However, filtering out the smaller scales, leading to simplified models such as Reynolds-averaged Navier-Stokes and large eddy simulations, significantly reduces computational cost for a (typically acceptable) loss of accuracy. In general, physics-based models governing conservation laws can be computationally simplified by using coarse-grained mesh resolutions and/or fewer solver iterations for convergence (in the iterative solver) to achieve a similar effect. We exploit the availability of models with tunable fidelities and propose multifidelity methods [57,21,23,56] that can judiciously utilize the cost-accuracy trade-off in models to estimate rare-event probabilities in aerospace design.
Surrogate models for conservation laws can be efficiently constructed via reducedorder models [42,41,43,46] for state variables and general regression/interpolation techniques [44,45,37,38] for scalar quantities of interest (QoI). In this work we are interested in a setting where a data-driven surrogate model for the QoI is adaptively constructed; that is, starting with an initial seed sample set (observations of the QoI), the surrogate model is sequentially updated with new samples that are optimally chosen to improve the prediction of the the failure boundaries in the design space. Furthermore, we account for the fact that our model of the subsystem has tunable fidelity parameters that trade predictive accuracy for computational cost. Therefore, our framework develops a surrogate model that is cost-aware while learning the relationships between the predictions between any two given fidelities. Additionally, the adaptive model construction chooses the appropriate control parameter and fidelity level to query at each step of the model-building process. For these reasons and since we are interested in reliability analysis with scalar QoI in this work, we use a Gaussian process (GP) [40] regression model to learn the mapping between the QoI and the augmented input-fidelity space.
The use of surrogate models for FPE is not new. Li and Xiu [26] show that a hybrid approach of using both surrogate models and high-fidelity function evaluations can lead to computationally cheap failure probability estimates. Moreover, they show that if the chosen surrogate model is of known accuracy, their proposed approach provides a consistent failure probability estimate (i.e., approaching the true function in the limit). Extending this work, Li et al. [27] show that the same hybrid approach can be used with importance sampling and a cross-entropy method to estimate biasing distributions. However, neither [26] nor [27] addresses the choice of high-fidelity model queries necessary to construct the surrogate model (that is, they do not actively train their surrogate model); furthermore they do not consider the availability of multiple-fidelity models. Peherstorfer et al. [32] demonstrate the idea of using a surrogate model to construct biasing distributions for importance sampling (IS), but that does not include details on adaptively constructing the surrogate model.
In the context of using GP regression as surrogate models, a large body of work exists, of which we briefly review what we believe are most relevant to our work. Seminal work on using adaptive GP models for contour finding and FPE was performed by Ranjan et al. [39], who use Jones et al. [19]'s expected improvement framework to propose an improvement function that seeks to place points close (within a userspecified distance) to the target level set. They show that their acquisition function efficiently balances exploration and exploitation, and they prove that their approach asymptotically predicts the contour exactly. Bichon et al. [5] then modify the formulation of [39]; both [39] and [5] have essentially the same meaning. As we will show later, our framework uses the acquisition function from [5] (and hence indirectly [39])-which we consider single-fidelity acquisition functions-to develop multifidelity extensions. Indeed, our acquisition function can extend those beyond [5] and [39] to the multifidelity setting.
Other noteworthy works that propose an acquisition function for adaptive GP construction are as follows. Bect et al. [4] propose a lookahead acquisition function using GP models, called stepwise uncertainty reduction. Echard et al. [10] use a sequential sampling strategy that greedily (myopically) seeks points with high probability of being in the failure domain. Picheny et al. [34] propose a weighted integrated mean-squared prediction error criterion, where the weight function is designed to be very high in regions where the GP mean is close to the target level set with high certainty, or when the GP uncertainty is high. Gotovos et al. [14] use GP to construct an acquisition function that picks the smallest distance to the target level set from one of the confidence bounds. They leverage the theoretical results of GP-upper confidence bound (UCB) [53] to derive a bound on the number of samples required to achieve a certain accuracy in resolving the level set. Dubourg et al. [9] use an adaptive Kriging surrogate model along with importance sampling to construct a quasi-optimal biasing distribution. Similar to [26], they propose a hybrid approach for FPE that combines the use of both the true limit state function and its (cheap) surrogate prediction. Wang et al. [60] use GP to approximate the forward model. Then they use Bayesian experimental design to sequentially add points toward the contour estimation. They use Simultaneous perturbation stochastic approximation (SPSA) [52] to appropriately optimize their stochastic acquisition function. More recently, Marques et al. [30] develop an entropy-based acquisition function that seeks to reduce the entropy of the failure contour. Furthermore, their method (CLoVER) allows for multiple information sources with varying costs of querying, and their acquisition function allows them to choose the information source. Their method, however, involves an integral in the input-dimension space, which can become intractable in higher dimensions. Recently, Cole et al. [7] improve upon [30] by developing a closed-form entropy expression for the contour that is (relatively) easier to optimize. While their method does not naturally include accounting for multiple information sources, they show that they outperform CLoVER on several problems.
Our work differs from the aforementioned work that use GP models in the following way. We propose the use of multifidelity GP [40] models to efficiently fuse models at multiple fidelity, where the fidelity space can be continuous or discrete. Even though [30] allows for multifidelity models, their approach depends on a Kennedy-O'Hagan [22]-like framework (see also [35]), where they learn a linear correction factor between each fidelity level and the highest-fidelity model. This has two drawbacks: (i) when the fidelity space is continuous, the approach can become intractable, and (ii) even for discrete fidelity spaces, they have to query each fidelity at a regular grid to learn the discrepancies. Our approach does not face either of these limitations. The method proposed in [24] involves the fusion of unbiased estimators from multiple fidelity models; however, the authors do not consider the choice of input and fidelities to query for an adaptive construction of the surrogate model, as we do. An adaptive construction of the surrogate model potentially results in a more judicious choice of input-fidelity pairs to query and thus a more parsimonious use of the computational resources. Furthermore, as stated previously, we propose a novel acquisition function that can provide a multifidelity extension to any single-fidelity acquisition function 2 , for example, [39,5]. Our acquisition function optimization allows for sequentially selecting samples from appropriate fidelity models to make predictions about QoI in the highest fidelity. We use our proposed approach in an importance sampling setting, similar to [32], where biasing distributions are constructed with the cheap multifidelity problems.
The remainder of the article is organized as follows. In Section 2 we provide the background and preliminaries including details about FPE and multifidelity GP models. In Section 3 we provide details of our method and discuss key theoretical properties in Section 4. In Section 5 we show results of our method demonstrated on the synthetic test functions and a real-world stress analysis of a gas turbine blade. We provide concluding remarks in Section 6.

Background and Preliminaries
Let x ∈ X ⊂ R d be the uncertain input to the function f : X → R. We assume that evaluations of f (x) are computationally expensive. We define the limit state function g : R → R that operates on f (x); and we define failure to occur, without loss of generality, when g(f (x)) < 0. The failure level set is then defined as For example, for a given failure threshold a ∈ R, g( Remark 1. We restrict our analysis to affine limit state functions of the form g(x) = ρf (x) − a, ρ ∈ R, a ∈ R, in order to preserve the Gaussianity of their sample paths, when we place a GP prior on f . This is not a restrictive assumption, however, because for nonaffine limit state functions we simply place a GP prior on g(x) instead of f (x).

Failure probability estimation
Consider the probability space (S, F, P) = (X , B(X ), P), where S denotes the sample space that in our setting is X ⊂ R d ; B(X ) denotes the Borel sets in X ; and P denotes a probability measure. We model the uncertain parameter x as a random variable in R d and assume that x is distributed according to the distribution P x , which has density q x . The failure hypersurface defined by C divides the parameter space into a failure set F = {x : g(f (x)) ≤ 0} and safe set F := X \F. We are interested in estimating the probability of failure p F that is defined by where I {x∈F } denotes the indicator function that takes the value 1 when x ∈ F and 0 otherwise. MC methods are often used to estimate p F ; the idea is to draw N independent and identically distributed (iid) samples of x, {x i , i = 1, . . . , N }, x i ∼ P x and approximate p F as where x i is the ith realization of x. The MC estimator is unbiased, as shown below.
The variance of the MC estimator is given by Using the approach described by (3) requires a large number of samples in order to obtain a failure probability estimate of desired accuracy. For instance, consider p F = 10 −3 and a root mean-squared error in the estimate of 10 −4 (and therefore a variance of 10 −8 ). Then the number of samples necessary to achieve this accuracy via MC is ≈ 10 5 . This could be prohibitively expensive in the reliability analysis of complex engineered systems with expensive forward models. A common means to (partially) circumvent this shortcoming is to employ an auxiliary distribution that concentrates the samples around the failure region [20,54]; this approach is referred to as importance sampling. Equation (2) can be equivalently written as where q x is the density of the biasing distribution P x . The key idea is that samples drawn from the biasing distribution are more likely to be from the set F. From (4), the IS estimate of p F can be written as where w(x i ) = qx(x i ) q x (x i ) are the importance weights and x i ∼ P x . The estimator in (5) is also unbiased, provided that q x (x) > 0 whenever As mentioned previously, we will construct the biasing distribution via a surrogate model in order to keep the computational costs tractable. As we will show later, similar to the construction of our surrogate model, our biasing distribution is also adaptively improved. We now introduce the surrogate modeling technique via GP regression.

Multifidelity Gaussian process regression
Crucially, we assume that f (x) is approximated by modelsf (x, s), where x are common inputs to all the models and s ∈ S ⊂ R is a tunable fidelity parameter for each model. For the sake of simplicity, in this work we set s to be scalar valued and without loss of generality S = [0, 1], where s = 1 and s = 0 represent the model at the highest and lowest fidelity, respectively; therefore, Additionally, we assume there is known a cost function c(s) : S → R, which is monotonic increasing in s and models the computational cost of queryingf at a specific fidelity. Note that we assume c(·) is independent of x for the sake of simplicity, but our framework applies to more general c(x, s) as well. Overall, we are interested in computing (4) by queryingf , while keeping the overall cost lower than computing (4) by querying f (x) alone. Our multifidelity method depends on learning a GP model that maps the augmented input-fidelity space (x, s) ∈ X × S to output quantities of interestf (x, s).
Remark 2. We mention that by assuming a monotonic cost function c(s) and the fact that fidelity increases as s increases in [0, 1], we assume that a hierarchy of fidelities exists among the set of models. However, such a hierarchy is not a necessary condition for our method to be applicable. In the absence of any known hierarchy among the models, the highest fidelity can be chosen as arg max s∈S c(s). If c(s) is unknown a priori, it can be learned, for example, by placing another GP prior on it and making a few observations.
We assume that we can make evaluations of the functionf at a given x and fidelity level s. That is,ŷ The key idea then is to specify GP prior distributions on the f . That is,f (x, s) ∼ GP (0, k((x, s), (x , s ))). The covariance function (or kernel ) k captures the correlation between the observations in the joint (x, s) space; here we use the product composite form given by k((x, s), (x , s )) = k x (x, x ; γ x ) × k s (s, s ; γ s ), where γ x ∈ R d + and γ s ∈ R d + parameterize the covariance functions for x and s, respectively. We estimate the GP hyperparameters Ω = {γ x , γ s } from data by maximizing the marginal likelihood. We choose an anisotropic Matérn-type kernel for both k x and k s . The posterior predictive distribution of the output Y , conditioned on available observations from the oracle, is given by [40] Y (x, s)|D n , Ω ∼ GP(µ n (x, s), σ 2 n (x, s)), where k n is a vector of covariance between (x, s) and all observed points in D n ; K n is a sample covariance matrix of observed points in D n ; I is the identity matrix; and y n is the vector of all observations in D n . Equation (7) is then used as a surrogate model for f in FPE. Note that µ n and σ 2 n are the posterior mean and variance of the GP, respectively, where the subscript n implies the conditioning based on n past observations. Our surrogate model for the limit state function then is given as g n (x, s) := g(µ n (x, s)).
Remark 3. Note that we assume our models at all fidelities are deterministic and observations are noise-free. Estimating failure probabilities with noisy observations is an interesting problem, but outside the scope of our work.

Remark 4.
For the remainder of this manuscript, we assume that the high fidelity modelf (x, 1) = f (x) is our ground truth. On the other hand, our approach naturally extends to situations involving disparate information sources, where the ground truth could be a physical experiment.
An important component of our proposed method is that the data required to build the surrogate model is adaptively selected. This ensures that the inputs and fidelities are judiciously chosen to be tailored to our specific goals. We call this adaptive approach "sequential experimental design," which we discuss next.

Acquisition-based sequential experiment design
We use principles from Bayesian optimization (BO) [50,12,51] to sequentially update the GP to improve our knowledge about our goal-in this case, estimation of C and the subsequent FPE. Akin to BO, we proceed by defining an acquisition function α(x, s), in terms of the GP posterior that is optimized (typically maximized) to select-given previous n observations-the next point x n+1 at s n+1 , and the process continues until a budget for computation is reached; see Algorithm 1. The choice of α(x, s) is key to the performance of the quality of the sequential point selection, which we discuss in the next section.
Naturally, the acquisition function is designed according to the specific needs of the sequential process. Classical acquisition functions for optimization include the probability of improvement (PI) [25], the expected improvement (EI) [19,31], and the GP upper confidence bound (UCB) [53]; these approaches work with a probabilistic estimate of the optimum via the GP. Other acquisition functions for optimization include entropy-based approaches [61,16,59]. Similarly, acquisition functions for active learning (purely learning the function f (x)) also exist [28,6]. Therefore, acquisition functions to learn a level set (e.g., failure boundary) can also be defined, such as the expected level set improvement of Ranjan et al. [39] and Bichon et al. [5], entropy-based approaches of Cole et al. [7] and Marques et al. [30], stepwise uncertainty reduction by Bect et al. [4], and other methods as in Picheny et al. [34].
In all the aforementioned approaches for failure level set estimation, we point out that application to multifidelity models is not straightforward. Although in [30] the authors show their applicability with multiple information sources, that requires learning individual discrepancies between models at each fidelity and the highest fidelity model. This can be computationally demanding in high-dimensional problems or when some of the lower-fidelity models are not necessarily cheap (compared with the highest-fidelity models). Consequently, such methods are restricted to discrete fidelity spaces with small cardinality. Furthermore, outside of the reliability literature, autoregressive approaches such as that of [22] and [33] suffer from the same limitation. Therefore, in this work we propose a new acquisition approach that applies to both discrete and continuous fidelity spaces. Additionally, our approach can extend any existing acquisition function to the multifidelity setting, as we will show later. We provide more details in the following section.
Algorithm 1: Generic acquisition-based sequential design Update GP hyperparameters Ω

Multifidelity acquisition functions
Drawing from the work of Renganathan et al. [44], where they propose an acquisition function for time-dependent problems with a finite horizon, we define a multifidelity value function framework that we use to develop our multifidelity acquisition function. Given the data from n observations D n , we define the value function υ by where h is a scalar-valued function, ω ∈ R p parametrizes h, Y |D n is the GP posterior given n observations, and y(x, s) is a sample path realized from Y |D n . The reasoning behind defining a value function as in (8) is as follows: 1. we can quantify the utility/value at the highest fidelity (our ground truth) by simply setting s = 1, i.e. υ n (x, 1; ω) : utility at the highest fidelity 2. computing this utility is cheap because it depends only on the posterior multifidelity GP 3. (8) generalizes several of the existing (single-fidelity) acquisition functions in GP-based reliability analysis, and finally 4. as we will show below, (8) allows a multifidelity extension of any given single fidelity acquisition function.
For example, some of the existing single fidelity acquisition functions for GP-based reliability analysis can be expressed as special cases of (8) as shown in Table 1, where [39], see (13). [5], see (15). a is the failure threshold and δ 2 (x) = ησ 2 (x), η ∈ R + defines the width of the uncertainty band around the GP posterior mean prediction. Then, the generalized multifidelity acquisition function is defined as { (x, s), y candidate obs. }   q y (y|D n ) dy, (9) where q y (y|D n ) is the probability density of the GP posterior distribution conditioned on D n . Notice that in the inner maximization we use υ(·, 1); this indicates that we place our value on the highest fidelity. In essence, (9) computes the expected value of the maximum value function at the highest fidelity due to a candidate input-fidelity pair (x, s), where the expectation is with respect to the posterior GP. We write (9) more concisely as In (10), note that υ n,1 = υ|D n {(x, s), y}, which is the value function given the union of D n and the n + 1th observation beinga drawn from the GP posterior Y |D n . In general, (10) is not guaranteed to admit a closed-form expression. Hence we make the MC approximation where υ j n, The inner maximization in (11) can be approximated by evaluating the value function υ n,1 on a discrete setX and the maximization approximated as (forx i ∈X ) As noted above, the value function is a function of the GP only, and hence its evaluation onX is computationally cheap. We now sequentially choose points by maximizing the cost-normalized acquisition function defined as In (12), since we normalize the acquisition function by the cost of evaluations c(s), the optimization penalizes decisions at high s, preventing expensive high-fidelity evaluations unless they offer a significant increase in the maximum value function compared with lower fidelities. For this reason, we consider our approach costaware and thus name it CAMERA: cost-aware, adaptive, multifidelity, efficient reliability analysis. We now discuss the choice of the value function.

Choice of value function
As previously mentioned, we propose a multifidelity acquisition function that provides a framework to extend existing acquisition functions, defined for a singlefidelity setting, toward a multifidelity setting. Therefore, while the choice of the value function can be practically any existing acquisition function, for the experiments in this work we choose the improvement-based acquisition functions for contour estimation by [39] and [5]. Ranjan et al. [39] propose the following improvement function, which we extrapolate to the multifidelity setting as where a is the target level and δ 2 (x, s) = ησ 2 (x, s), η ∈ R + defines the width of the uncertainty band around the GP posterior mean prediction. The intuition behind (13) is that improvement is defined as getting points closer to the target level set or in places where the uncertainty of the GP predictions is large, while still being within a certain confidence region of the GP prediction. Consequently, I r (x) = 0, when (y(x) − a) 2 > 2 (x) and therefore regions far away from the potential contour locations are avoided. Bichon et al. [5] modify (13) as which results in a closed-form expression for the expected improvement , and a ± = a ± δ(x, s). Note that (13) and (14) practically have the same meaning, but (14) leads to a more convenient analytical form upon taking the expectation, as opposed to that in (13). We therefore propose to specify our value function, also extending it to the multifidelity setting, as and define our acquisition function as follows: which essentially quantifies the average value of the maximum EI b at the highest fidelity (s = 1) due to a candidate decision-fidelity pair (x, s). Note that we extend the acquisition function in [5] to general multifidelity cases, where the fidelity space could be continuous or discrete. The acquisition function in (17) is then plugged into Algorithm 1 to obtain a final surrogate modelg q (x, 1) off (x, 1).

Importance sampling
The last step of our approach is the estimation of the failure probability. The multifidelity surrogate model that is adaptively constructed is now used to construct a biasing distribution. This distribution is then used to generate a sample set using which the failure probability is estimated.
Let q (n) x be the density of the biasing distribution, P (n) x , after n rounds of updating our multifidelity surrogate model. Then, the importance sampling estimate is given by where x i , i = 1, . . . , N are iid samples drawn from P (n) x .

Overall CAMERA algorithm
We are now ready to summarize our overall algorithm. We use a Gaussian mixture model (GMM) [18] to learn the biasing density q (n) x from a sample set generated from the GP. For all the cases, we set the number of components of the GMM to 25 and restrict the covariance matrix to be diagonal. A few trials indicated that these settings are sufficient to get a good approximation to the biasing density across all the experiments. The expectation-maximization algorithm [8,62] is used to learn the hyperparameters of the GMM. Once fit, this distribution is sampled to generate a biasing set that is evaluated on the high-fidelity modelf (·, 1) to compute the IS estimate of failure probability. The exact steps are summarized in Algorithm 2.
Note that M in step 6a in Algorithm 2 can be picked arbitrarily large since it does not involve any high-fidelity model queries, only the surrogate model. On the other hand, N must be tractable since it involves high-fidelity queries. Instead of directly setting M , we first uniformly sample X with 10 7 samples and filter out the sample set to obtain X as in step 2 using the surrogate model.

Theoretical Properties
Recall that our surrogate-based IS estimator is unbiased. We now show that our adaptive surrogate model is consistent; that is, with high probability the predicted level set converges to the true level set asymptotically. We first demonstrate that the points selected by maximizing our multifidelity acquisition function (12) asymptotically reduce the GP posterior variance to 0 everywhere; see Lemma 2 and Lemma 3. Then, in Theorem 1, we show that our surrogate model is asymptotically consistent.
We begin by making a smoothness assumption on our unknown function f , followed by an assumption on the limit state function. Assumption 1. We assume that the sample paths drawn from the GP prior on f are twice continuously differentiable. Such a GP is achieved by choosing the covariance kernel to be four times differentiable [13,Theorem 5], for example, the squaredexponential or Matérn-class kernel. Then, the following holds (for L ∈ R + ): Assumption 2. The limit state function g is an affine function of f , defined as g(f (x)) = ρf (x) − a. W.l.o.g., we pick ρ = 1.
First, we show that our proposed two-step lookahead acquisition function will never choose the same point twice. We present this in Lemmas 1 and 2.  1 Input: • Surrogate model: Multifidelity GP posterior of f , GP(µ n (x, s), σ 2 n (x, s)) and observed data D n .
• Cumulative cost C n = n j=1 c(s j ), total budget B. • High-fidelity sample size N .
• Biasing sample size M .
Let q > n be the last iteration of updating the surrogate model.

Fit a biasing distribution P
(a) Obtain a sample set of cardinality M , using the multifidelity surrogate model at s = 1: (b) Fit a GMM with X to obtain the biasing distribution P (q) x with density q (q) x . 7. Obtain a biasing sample set X biasing = {x k ∼P q x , k = 1, . . . , N }.

Compute indicator function (with high fidelity samples) and IS weights
x (x k ) , k = 1, . . . , N, 9. compute the IS failure probability estimate where υ n (·, ·) = υ(·, ·)|D n . Using the definition of the value function, we have where (x i , s i ) and (x j , s j ) are the input-fidelity pairs that maximize α(x, s)|D i−1 and α(x, s)|D j−1 , respectively, and hence are points previously observed by the algorithm.
Proof. For a given (x, s) ∈ X × S, let σ 2 n (x, s) = 0; in other words, a noise-free observation was previously made at (x, s). Then, which follows from the fact that y(x, s)|(σ 2 n (x, s) = 0) = µ n (x, s) and that the value function conditioned on a previously observed point is unchanged. Now consider a candidate point (i.e., an unobserved point) ( where the inequality follows from Jensen's inequality [47,Ch. 3] and we have also used Lemma 1. Furthermore, the equality holds only if (x, s) = (x , s ) or σ 2 n (x, s) = 0 = σ 2 n (x , s ), and therefore (x, s) = (x , s ), since (x , s ) is unobserved and hence σ 2 n (x , s ) > 0. Now, using Lemma 2, we will show that the multifidelity acquisition strategy reduces the GP posterior variance to 0 at s = 1, ∀x ∈ X . Lemma 3. Let {x i }, i = 1, . . . , n be a sequence of points selected via our multifidelity acquisition strategy. Then, lim n→∞ sup x∈X σ 2 n (x, 1) = 0. Proof. W.l.o.g, let k((x, s), (x, s)) = 1. Recall that the posterior variance of the GP is given by σ 2 n (x, s) = k((x, s), (x, s)) − k n K −1 n k n . Note that as n → ∞, and if Lemma 2 holds, then k n → K :,i n , where by K :,i n we mean the ith column of K n , i ∈ [n], and therefore k n K −1 n = e i , where e i is an n-vector with 1 at the ith location and 0 elsewhere. Then We are now ready to show that with high probability, the predicted failure level set via our multifidelity surrogate model converges to the true failure level set. For this purpose we assume a discrete setC ⊂ C and show that the absolute difference between the surrogate prediction atC and the true failure level set is bounded with high probability. Theorem 1 (Asymptotic consistency). Let δ ∈ (0, 1) and Assumptions 1 and 2 hold, with L 2 = 2 log 2 δ . LetC denote a discretization of C with |C| < ∞, where we evaluate the true failure boundary. Let ≥ 0, and a n = 2 log(b n |C|)/δ, where ∞ n=1 1/b n = 1. Then the predicted failure level set converges to the true failure level set with high probability; that is, the following holds: Proof. The proof is based on [53], but we differ in the choice of constants. We first present a simple identity that bounds the probability that a Gaussian random variable is greater than some value. Given r ∼ N (0, 1) and c > 0, consider the following: where the second line is due to the fact that 0 ≤ e −c(r−c) ≤ 1 for r ≥ c and c > 0 and the last line from the fact that the integrand is the density of N (c, 1). Using this result and the prior assumption that the unknown function f has a GP prior, we state, with c = a n and a fixed x i ∈C, Applying the union bound ∀i, we get Pr [|g(f (x i , 1)) − g(µ n−1 (x i , 1))| > a n σ n−1 (x i , 1)] Now we apply the union bound ∀n: We would like to extend the proof to a general x ∈ X . For some x ∈ X , letx denote a point ∈C that is nearest (in Euclidean distance) to x. Then, per Assumption 1 and by setting L 2 = 2 log 2 δ , the following statement holds: Pr |g(f (x, 1)) − g(f (x, 1))| > 2 log Combining (22) and (23), the following statement holds: Pr ∀x ∈ X , ∀n, |g(f (x, 1)) − g(µ n−1 (x, 1))| ≤ a n σ n−1 (x, 1) + 2 log 2 δ x −x > 1−δ.

Numerical Experiments
We first demonstrate our proposed method on synthetic test functions of varying dimensions in Section 5.1 and then on two real-world applications: a gas turbine blade stress analysis test case using finite element analysis in Section 5.2 and a transonic wing aerodynamics test case using a finite volume method in Section 5.3. Our goal is to compare the multifidelity predictions of the failure probability against the predictions with the high-fidelity model only. For the synthetic test cases and the gas turbine blade test case, we estimate the true failure probability using a naive MC with a set of points drawn uniformly at random from X ; note that we use this estimate as the ground truth. Then we compute the failure probability error as . For the transonic wing test case, the true failure probability is unknown and is prohibitive to compute; we therefore report only the comparative predictions between single fidelity and multifidelity.
For all the experiments, we provide a set of seed points D n = {(x i , s i ), y i , i = 1, . . . , n}, selected via a randomized Latin hypercube design, to start the algorithm. Note that we start both the single-fidelity and multifidelity runs with the same budget of information; therefore, for any given experiment, the n for single-fidelity and multifidelity runs are different in order to maintain the same total computational cost with which we start the experiment. Naturally, whereas for the multifidelity experiments the seed points are distributed uniformly in S, for the single-fidelity experiments the seed points are constrained to S = {1}. All the synthetic experiments are repeated a total of 20 times with randomized starting seed points. The turbine stress analysis test case is repeated only 5 times, keeping in consideration the higher computational cost of this experiment. There is no repetition in the transonic wing test case because that would be computationally prohibitive.
For all the experiments, the following cost model is assumed: where c 0 is a normalizing constant that is required to ensure that the acquisition function optimization is not biased by the actual costs of the models; note that c 0 is also the cost of the high-fidelity (s = 1) model, for which we fix c 0 = 500. The c 1 > 0 determines the steepness of the cost model; we set c 1 = 10 for all the experiments. Finally, c 2 > 0 is a constant to ensure that the cost model is bounded away from zero and is roughly the cost of the lowest-fidelity model (f (·, 0)); we fix c 2 = 0.1. We  note that our cost model-by virtue of setting c 1 = 10-is designed to emphasize the disparity between the computational costs of lower-and higher-fidelity model queries; see Figure 1. In Figure A.17 we show the effect of the cost model with varying c 1 and discuss its effect on the predictions.
We evaluate our predictions qualitatively for the 2D cases, in other words, by comparing the failure level sets estimated via the final GP mean µ q (x, 1) against the true level sets in contour plots. For all the synthetic cases (including the 2D cases), we then compute the error in failure probability, as mentioned previously. The true failure probability is estimated via a naive MC with N = 10 6 points sampled uniformly from X . For the turbine case, we first evaluate the high-fidelity model at a total of 150 parameter snapshots and use the observed QoIs to construct a radial basis function (RBF) interpolant. We then perform a naive MC on the RBF interpolant to obtain the true p F . We begin by discussing the synthetic experiments in the next section.

Synthetic experiments
We introduce several new multifideity test functions that can be used to benchmark multifidelity methods beyond the scope of this work. Specifically, we introduce the 2D multifidelity multimodal function [5], the 2D multifidelity four-branches function [49], and the 3D multifidelity Ishigami function [17]. Additionally, we use the 6D multifidelity Hartmann function provided in [55]. The specifics of these test functions are discussed as follows.

2D multimodal function
We extend the test function in [5] to multifidelity setting as shown in (24), a correlated fidelity dependence, which is continuously differentiable in X × S; see Figure 2 for snapshots.

2D multifidelity four-branches function
A system with four distinct component limit-states [49] is a commonly used test problem in reliability analysis to demonstrate the efficacy of the algorithms in the presence of multiple failure regions. Here, in contrast to (24), we introduce fidelity by making a coordinate transformation as x i = x i − a 0 × s, where a 0 ∈ R. The response function can be written as The multifidelity extension we provide offers a translation of the level sets across fidelities, where we set a 0 = −5.0, and is shown in Figure 3.

3D Multifidelity Ishigami function
We modify the 3D Ishigami test function [17], similar to the four-branches function, by introducing a shift in the variables using the fidelity parameter s, as shown in (27). Note that the input x is normalized to be in [0, 1] 3 before introducing the translation: where x = [x 1 , x 2 , x 3 ] and X = [−π, π] 3 . The resulting function is shown for two slices x 3 = −π and x 3 = 0 in Figure 4.

6D Hartmann function
For a higher-dimensional test function, we use the "augmented" Hartmann 6D function introduced in [3] and restated below:   For the 2D test functions previously discussed, the failure threshold was set as a = 0.. For the Ishigami and Hartmann functions, we show the histogram of the output at s = 1 along with the failure thresholds (vertical lines), which are set to a = −9.0 and a = −2.0, respectively, in Figure 5; these thresholds were chosen particularly to make the failure probabilities on the order of 1%.
The results of the 2D test functions are shown in terms of the predicted contours in Figure 6 and Figure 7; in both figures, whereas the true level set is plotted in the red dashed-dotted line, the background (filled contours) shows the final GP posterior mean µ q (x, 1) with the predicted level set as the thick black line. We hope that readers are able to visually compare the predictions against the true failure boundary. Note that in both figures we show only the predictions of one of the twenty random repetitions. Whereas for the multimodal function the predictions from the single-and multifidelity approaches are similar, the multifidelity achieves this prediction with a computational cost that is more than 45% cheaper. For the four-branches test function, the prediction of the multifidelity approach is visually much better than that of the single-fidelity approach, albeit using only half the computational cost.
We use the contours to provide a visual interpretation of the results, since they are feasible for the 2D cases. For a quantitative evaluation, however, we also show the error in computing the failure probability. The absolute errors in computing the failure volume are shown in Figure 8 for all the synthetic test functions. We run both the single-and multifidelity experiments for 100 iterations after providing them with seed samples, as previously mentioned. In all four cases, we observe that the multifidelity approach is able to predict the failure volume with greater accuracy compared with the single-fidelity approach, for a certain computational budget. In order to establish the utility of the adaptive nature of our proposed algorithm, we also compare our predictions with a "non-adaptive" Latin hypercube design. In Figure 8, we show the predictions from a non-adaptive design (horizontal dashed line) with a total budget of the n seed points plus the 100 used in each experiment. This way, we keep the budget consistent with each experiment and repeat this 20 times to randomize the selection of the Latin hypercube points; Figure 8 shows the average of the repetitions. We observe that the non-adaptive design leads to a worse prediction of the failure probability compared to CAMERA for the same budget. This is intuitive because a (space-filling)  Latin hypercube design may provide good coverage of X but can fail to resolve the failure boundary which is critical to fitting a good biasing density and hence to a good estimate of the failure probability. On the other hand, CAMERA attempts to adaptively select points that would best predict the failure boundary at the highest fidelity level.
In addition to the FPE, we are interested in observing how the CAMERA algorithm balances the trade-off between accuracy and cost across the fidelity space. In this regard, we show the distribution of the number of queries to models at each fidelity level, across all the 100 adaptive selections, in Figure 9, where all the 20 repetitions are overlaid on top of each other. We observe that when there is high correlation between the low and high-fidelity models, then the algorithm exploits this to obtain most samples from the lowest fidelity model, e.g., Multimodal, Four branches, and Ishigami functions. For these test functions, only a few queries at the intermediate fidelities were necessary. Most of the evaluations are at s = 0 (lowest fidelity) and s = 1 (highest fidelity). In the case of the Ishigami test function, which also has the rarest failure probability, the algorithm has most queries at s = 1, with queries at intermediate fidelity levels almost uniformly distributed. Due to the lack of sufficient correlation between the high and low fidelity levels, the algorithm is unable to exploit the cheaper information as in the other experiments. However, notice that across all the experiments, the highest fidelity model is queried only about 25% of the time. When sufficient correlations exist, the lowest fidelity model is queried 50-75% of the time. This demonstrates the efficiency of the method in exploiting correlations between models at multiple fidelities. Finally, notice that the distribution of the number of queries is roughly consistent across all the repetitions. This demonstrates the robustness of our algorithm amidst the randomness of the seed point selection.

Turbine reliability analysis with discrete fidelity space
We now demonstrate our method on the structural analysis of a gas turbine blade at steady-state operating conditions. Gas turbines have a radial arrangement of blades downstream of the combustor. Hence, these blades are subject to very high pressures and temperatures that can cause concentrated zones of high stress that can in turn lead to structural deformation and potential catastrophic damage of the engine parts. Therefore, we are interested in the probability that the maximum stress acting on the thermal blade exceeds a certain threshold. The variables are defined as follows.
x 1 : Pressure-side loading x 2 : Suction-side loading s : mesh resolution, In this regard we vary the pressure-side and suction-side pressures as boundary conditions, as shown above, while ignoring thermal effects, to observe the maximum von Mises stress acting on the blade; that is, f (x) = max z u(z), where u(z), z ∈ ∂z is the Von Mises stress distribution on the blade surface ∂z. Note that we negate the output, f (x) × −1, to comply with our definition of failure. The material properties are fixed to the following constants: Young's modulus 227E + 9 Pa, Poisson's ratio 0.27, and coefficient of thermal expansion 12.7E − 6. The governing equations are solved using a finite element method using Matlab's PDE toolbox. The computational mesh density is determined via a minimum element size which also serves as the fidelity parameter. Specifically, we use the minimum mesh element sizes {0.05, 0.025, 0.01}, which correspond to fidelities S = {0.0, 0.5, 1.0}, respectively; note that this case corresponds to a discrete fidelity space. The mesh resolutions for the specified three fidelity levels are shown in Figure 10.   Mises stress is the quantity of interest in this experiment; see Figure 11 for a couple of snapshots of the von Mises stress distribution and blade deformation corresponding to two different pressure loadings.
The suction-side x 1 and pressure-side x 2 pressures are varied in the following ranges: x 1 ∈ [201160, 698390], x 2 ∈ [301460, 598020]. A total of n = 40 solution snapshots are generated and supplied as seed points for the multifidelity run; the equivalent n for the single-fidelity run is determined to match the multifidelity seed budget. Both experiments are then run for 100 iterations and are repeated 5 times with randomized seed points.
The contour plots of maximum von Mises stress are shown in Figure 12 for each of the three fidelities considered. The failure boundary is set as the level set corresponding to 95% of the maximum von Mises stress across a total of 150 high-fidelity simulations that were computed offline. This is shown by the thick black line in Figure 12 (far right) and corresponds to a probability of failure of 0.004236. Notice that for the low-fidelity case (s = 0.0) this level set is outside the domain off (x, 0.0) and for the medium-fidelity case (s = 0.5) this level set is at the top-left corner of the domain. Therefore, this test case poses a unique challenge of the fidelity levels leading to very different function domains; for example, according to the low-fidelity model, the probability of failure is 0.
To evaluate our algorithm, we compute the predicted failure probability via IS at the end of each iteration and take the absolute difference from the true failure probability. This is plotted against the cumulative computational costs for the singleand multifidelity cases in Figure 13. From the figure, notice that for any given cumulative cost, the multifidelity approach results in a failure probability error that is better than or equal to the single-fidelity approach. Because of the availability of only three discrete fidelity levels and the assumed cost model, the discrepancies in the cost between low and high fidelities are very pronounced; the single fidelity curve starts at a cumulative cost of ≈ 550 (which is the cost of a high-fidelity query) and the multifidelity curve starts at a cumulative cost of ≈ 50 (which is the cost of the lowest-fidelity query). Notice that, for the chosen budget of evaluations, CAMERA converges to a slightly higher failure probability error than the high-fidelity; however, this discrepancy (O(0.01%)) can be considered negligible.
On the right of Figure 13, we show the distribution of the number of queries at each fidelity level. Notice that the algorithm is able to exploit the lowest fidelity model while keeping the high-fidelity queries at ∼ 20%, owing to the strong correlations between the fidelity levels for this experiment. The consistency across repetitions is also present, similar to what was observed in the synthetic experiments.

ONERA transonic wing
The ONERA-M6 wing [48] is an analytical (shape prescribed using mathematical equations) wing created in the 1970s by the French Aerospace Laboratory for the purposes of understanding three-dimensional airflow in transonic speeds and high Reynolds numbers. 3 It is a semi-span, swept wing with no twist, and at transonic Mach numbers is characterized by local supersonic flow, shocks, and turbulent boundary layer separation. Therefore, this wing has been frequently used for benchmarking computational software for compressible turbulent flow [29]. Because of its relevance in the commercial and military aircraft flight regimes, which make it an important candidate for reliability analysis, we use it as one of the test cases in this work. The wing, shown in Figure 14a, is surrounded by a fluid domain that is a quadrant of an ellipsoid, as shown in Figure 14b. The two external sides of the fluid domain are designated as "freestream" boundary conditions, with a prescribed ambient temperature, flow speeds, and angle of attack (that is, the incidence angle between the flight direction and the wing mean aerodynamic chord). The fluid domain is discretized with a mesh of approximately 56, 000 hexahedral elements with local refinement to resolve boundary layers, shocks, and shear layers. The three-dimensional compressible RANS equations are solved with a finite-volume method with second-order spatial discretization and Euler implicit discretization. The turbulence closure is provided by the Spalart-Allmaras   For the reliability analysis, we consider the following input parameters: Therefore, x = [x 1 , x 2 , x 3 ], and our f (x) is the wing load factor L/W , where L is the net aerodynamic force acting perpendicular to the flight path and W is the weight of the aircraft. Note that the maximum wing load factor obtainable during flight is limited by the structural strength of the aircraft wing [2,Ch. 4], and hence we want to ensure the L/W does not exceed a certain safety threshold, over a range of operating conditions. The fidelity parameter s is the number of iterations allowed for the solver to converge; fewer iterations result in a faster (and hence cheaper) solution that is considered low fidelity and vice versa. We set the following bounds on the control and fidelity parameters: x 1 ∈ [5 × 10 6 , 11 × 10 6 ], x 2 ∈ [0.7, 0.85], x 3 ∈ [1, 3], and S ∈ [0, 1] (which is normalized and corresponds to the number of iterations ∈ [500, 5000]). Note that we round any fractional values for the number of iterations, recommended via our algorithm, to the nearest integer. We show the pressure coefficient distribution on the wing for the different values of s in Figure 15. As in the turbine test case, we negate the output of the function f and set the failure threshold a = −7.5 and therefore are interested in estimating the probability that g(x) = f (x) − a ≤ 0.
We begin our algorithm with n = 100 seed points obtained via running the simulations at a Latin hypercube design of experiments across the joint X × S space, and we run the algorithm for 100 iterations. For the single-fidelity run, n is set to keep the total cost of seed points identical to that of the multifidelity run. We keep the same cost model as in Figure 1. Note that for this experiment we compute the failure probability only once-at the conclusion of running Algorithm 2 for 100 iterations. Furthermore, we do not repeat the experiment with randomized seed points. We fix N = 400 high-fidelity samples to evaluate p IS F . We show the results in Figure 16, where p IS F and V(p IS F ) are used to construct the confidence intervals. While we do not know the p F (truth) for this case, we see that the predictions from CAMERA are with smaller confidence intervals compared with using the high-fidelity model only. The overall cumulative cost of the high-fidelity run is 44,667 compared with only 38,614 for CAMERA, thereby establishing the cost-effectiveness of the method.

Conclusion
In this work we address the problem of failure probability estimation for complex engineered systems with a computationally expensive model. Additionally, we account for the fact that the expensive model has tunable fidelity parameters that can trade computational cost for predictive accuracy. We propose an approach based on Gaussian process models that can learn the correlations between the different fidelities, given a continuous or discrete fidelity space, and therefore predict the highfidelity function with potentially fewer queries to the expensive high-fidelity model      that are supplemented by cheaper low-fidelity model queries. Further, we propose an acquisition-based sequential experiment design approach to adaptively build the multifidelity surrogate model, so that only the most useful fidelities and inputs are queried. The main highlight of our proposed acquisition function is that it provides a generalized multifidelity framework that can extend existing (single-fidelity) acquisition functions to the multifidelity setting and furthermore applies to both discrete and continuous fidelity spaces. Additionally, we demonstrate that our surrogate modeling approach is asymptotically consistent, meaning that the predicted failure level set approaches the true level set with high probability as the number of iterations approaches infinity. We also show that estimates of the failure probability, predicted by using only the developed multifidelity surrogate model, produces very high accuracy at a lower computational cost, compared with using the high-fidelity model only. Furthermore, using the multifidelity surrogate model with a single batch of (non-adaptive) Latin hypercube design results in poor estimates of the failure probability, compared to the proposed adaptive approach. This confirms the need for adaptive training of surrogate models for failure probability estimation.
With the sequential surrogate model building approach, we observe that the predicted failure level sets are very close to the true level set. Furthermore, we observe that in the presence of multiple fidelity models, our approach performs even better by requiring only ≈ 50% of the computational budget. In terms of predicting the failure probability, our approach shows a decrease in the prediction error with increasing evaluations of the models, while the multifidelity approach always outperforms the single-fidelity case. In addition to the synthetic cases, we demonstrate our method on a real-world turbine reliability case, with discrete fidelities, where the failure probability is ≈ 0.4%, and a transonic wing test case. We observe results consistent with the synthetic cases where we are able to obtain very high accuracy in failure probability error, with the multifidelity approach significantly outperforming using the high-fidelity model only. Overall, our experiments reveal that our approach can lead to cost-effective yet accurate predictions for reliability analysis, for example, using only N = 400 high-fidelity samples for estimating p IS F in the transonic wing example. We anticipate conducting follow-up studies of this work in several directions. First, the CAMERA algorithm can be modified to adaptively provide an estimate for p IS F ; currently, it is estimated at the conclusion of the adaptive surrogate model construction. This is useful for real-world applications, where running expensive highfidelity simulations in multiple batches is less risky and therefore more practical. Along the same lines, we plan on evaluating the performance of our sequential experiment design of adaptive surrogate model construction, in making batch evaluations at every step. Furthermore, whereas we have demonstrated our approach in estimating p F in O(10 −4 ), it would be interesting to evaluate the approach on problems with smaller orders of magnitude in p F .

Appendix A. Impact of the choice of cost model
We vary the choice of the constant c 1 in the cost model to observe its impact on the algorithm performance. Recall that c 1 determines the rate of change of the computational cost of the models with respect to fidelity. In this regard, smaller values of c 1 result in "flatter" cost models, as shown in fig. A.17a. We use the Ishigami test function to evaluate the impact of the choice of cost model on the performance of CAMERA, whose results, with 20 repetitions with randomized seed points, are shown in fig. A.17b. We observe no appreciable difference in the algorithm performance between a choice of c 1 = 10 and c 1 = 5. At c 1 = 1.0, CAMERA still outperforms the high-fidelityonly run, although suffering a higher failure probability error at lower cumulative costs. As c 1 approaches 0, the performance approaches that of the high-fidelity-only run, as expected. Therefore, we conclude that the choice of the cost model, while expected to impact the performance of the algorithm, does not necessarily offer an undue advantage to the algorithm. This experiment justifies the use of c 1 = 10.0 in all of our experiments.   Appendix B. Expected improvement acquisition function of [39] The closed-form expression for the expected improvement acquisition function by Ranjan et al. [39], obtained by setting δ 2 (x, s) − min{(y(x, s) − a) 2 , δ 2 (x, s)} and ω = {a, η} in (8), is given by