A Multilevel Simulation Method for Time-Variant Reliability Analysis

: Crude Monte Carlo simulation (MCS) is the most robust and easily implemented method for performing time-variant reliability analysis (TRA). However, it is inefﬁcient, especially for high reliability problems. This paper aims to present a random simulation method called the multilevel Monte Carlo (MLMC) method for TRA to enhance the computational efﬁciency of crude MCS while maintaining its accuracy and robustness. The proposed method ﬁrst discretizes the time interval of interest using a geometric sequence of different timesteps. The cumulative probability of failure associated with the ﬁnest level can then be estimated by computing corrections using all levels. To assess the cumulative probability of failure in a way that minimizes the overall computational complexity, the number of random samples at each level is optimized. Moreover, the correction associated with each level is independently computed using crude MCS. Thereby, the proposed method can achieve the accuracy associated with the ﬁnest level at a much lower computational cost than that of crude MCS, and retains the robustness of crude MCS with respect to nonlinearity and dimensions. The effectiveness of the proposed method is validated by numerical examples. efﬁciency of the procedure while retaining its accuracy, we present a sampling method called the multilevel Monte Carlo (MLMC) method. The idea of multilevel Monte Carlo path simulation is introduced into TRA.


Introduction
Reliability is an important measure of structures and systems [1][2][3][4]. It is defined as the probability that a structure or system successfully performs its intended function under specified conditions in a given time interval. Because of material degradation, time-variant loads, etc., a mechanical structure commonly has a time-variant response in performance function and its reliability decreases over service time. Accurately estimating reliability in the time interval of interest is extremely important for design optimization, reliabilitybased maintenance and risk control, which makes time-variant reliability analysis (TRA) methods urgently required in engineering. In recent decades, a number of methods have been developed for assessing time-variant reliability [5]. They can be generally divided into two categories, i.e., analytical methods [6][7][8][9] and sampling-based methods [10][11][12].
The representatives of analytical methods include first-crossing based methods [6,13] and time discretization-based methods [14,15]. After proposing the Rice formula, first-crossing based methods have gained much attention. Renaud et al. [7] presented an approach called PHI2, computing the outcrossing rate by solving a two-component parallel system reliability problem. Sudret [13] provided an analytical expression of the outcrossing rate in order to improve PHI2 in terms of accuracy. Despite the fact that classical first-crossing based methods are popular for TRA, they still experience poor accuracy because of the assumption that outcrossing events are mutually independent. Hu and Du [6] relaxed the assumption using the concept of joint crossing rate, and Yu et al. [16] proposed an approximation of the first-crossing probability density function (PDF) method to avoid directly computing the out-crossing rate. Meanwhile, time discretization-based methods were proposed in recent years, which implicitly accounted for the co-dependence among out-crossing events. Jiang et al. [14,17] transformed the time-variant reliability problem into a time-invariant series system reliability The so called Multilevel Monte Carlo path simulation was first developed by Giles for stochastic differential equations [39]. The method uses a geometric sequence of timesteps to discretize the considered time interval. The discretization error would be small, but the computational cost large if only the finest level were used. The situation would be reversed if only the coarsest level were used. To reduce the computational cost of computing the expectation of a random quantity defined by a stochastic differential equation, Giles [39] solved the problem on the finest level by estimating corrections using all the levels, thereby achieving the accuracy associated with the finest level at a lower cost. The idea is adopted herein, retaining the accuracy associated with the finest level but performing MCS at coarser levels to reduce the variance of the estimate of cumulative failure probability. The total number of function evaluations is minimized by optimizing the number of random trajectories allocated to each level.
The rest of the paper is organized as follows. Section 2 introduces the crude MCS for time-variant reliability analysis, including its main procedure in Section 2.1 and computational complexity in Section 2.2. MLMC for time-variant reliability analysis is elaborated in Section 3 and validated in Section 4. Section 5 provides the conclusions.

Crude MCS
This paper considers the time-variant performance function with the typical expression G(X,Z(t),t), where X is the vector of time-invariant random variables and Z(t) is the vector of stochastic processes. This special case G(X,Z(t),t) appears in many engineering scenarios, e.g., structures under stochastic loads and material degradation. We assume that G(X,Z(t),t) is computationally cheap and a function evaluation can be performed at any time τ. For problems with implicit and time-consuming performance function, one can build an asymptotically explicit expression through employing a surrogate model, such as Kriging and polynomial chaos expansion [34,38].
Without loss of generality, the failure of a structure, whose performance function is fully described by G(X,Z(t),t), is defined as G(X,Z(t),t) ≤ 0. In the time interval of interest [0, T], a failure event occurs if the response value of the performance function G(X,Z(t),t) crosses the failure threshold 0. Accordingly, the cumulative probability of failure can be expressed as, where Pr{·} is the probability operator, and "∃" means "there exists". Herein, the failure indicator function I(x,z(t)) is introduced, whose expression is, where z(t) is a trajectory of Z(t). Thus, P f,c (0,T) can be rewritten as follows: As illustrated in Section 1, crude MCS is the most robust method to estimate P f,c (0, T). It is commonly implemented as follows: Step 1: Discretize the time interval [0, T] into K + 1 uniform instants with time step ∆t. For convenience of the illustration below, we control the number of discretized time instants and the corresponding timestep using two positive integers, i.e., K 0 and L. Then, P f,c (0,T) can be approximated by: Let P (L) f,c (0, T) denote the approximation of the cumulative probability of failure: where I (L) (x,z(t)) is an approximation of the indicator function I(x,z(t)): Step 2: Generate i.i.d. (independent and identically distributed) samples of X and trajectories of Z(t), and estimate P f,c (0,T) according to the law of large numbers: (8) where N MC is the number of random simulations crude MCS performs. Multiple techniques are capable of generating random trajectories of stochastic processes, e.g., the expansion optimal linear estimation method (EOLE) and Karhunen-Loève (KL) expansion are widely used in the literature [40][41][42]. In the section of method validation below, EOLE is adopted. LetP MC f,c (0, T) denote the estimate obtained by crude MCS: It is worth noting that the error in the estimateP MC f,c (0, T) can be asymptotically determined using the following equation: The first term is the bias introduced by the uniform discretization, and the second term corresponds to the error due to the random sampling. Given values of N 0 and L, this paper uses a multilevel method to reduce the error caused by MCS, leaving unchanged the bias due to the discretization.
The coefficient of variation associated withP MC f,c (0, T) is: where

Computational Cost of Crude MCS
Given a realization of X and a trajectory of Z(t) (e.g., x n and z n (t)), K 0 · 2 L + 1 evaluations of the time-variant performance function G(X,Z(t),t) are needed in order to determine the value of I (L) (x n , z n (τ)). This paper ignores the asymptotically negligible cost of generating random realizations of X and trajectories of Z(t), and uses the number of function According to Equations (11) and (12), to achieve a coefficient of variation smaller than δ Pf , i.e., cov Pf < δ Pf (13) the number of trajectories of G(X,Z(t),t), i.e., N MC , needs to satisfy that: Hence, the computational cost of crude MCS corresponding to the accuracy requirement δ Pf , i.e., the total number of function evaluations, is:

Multilevel Monte Carlo for Time-Variant Reliability Analysis
Crude MCS is capable of achieving any required accuracy regardless of nonlinearity and dimensions of the time-variant performance function if K and N MC are sufficiently large. However, its commonly computationally expensive. To overcome this issue, multilevel Monte Carlo is applied herein to deal with TRA.
MLMC discretizes the time interval [0, T] using a geometric sequence of different timesteps ∆t l = T/(K 0 · 2 l ) (l = 0, . . . ,L). The smallest timestep ∆t L corresponds to the discretization bias defined by the first term of Equation (10). The multilevel method optimizes the computational cost on each level in order to reduce the overall computational complexity of estimating the cumulative probability of failure, retaining the accuracy associated with the level L.

Multilevel Monte Carlo
Consider a sequence of approximations P f,c (0, T) can be rewritten as: P Then, it is easily concluded that: where P (0) = E I (0) (X, Z(τ)) and P (l) = E ∆I (l) (X, Z(τ)) ∆I (l) (X, Z(τ)) = I (l) (X, Z(τ)) − I (l−1) (X, Z(τ)) l = 1, 2, . . . , L MLMC independently estimates each expectation on the right side of Equation (17). An unbiased estimate of the cumulative probability of failure is then obtained, which has the following expression:P ML f,c (0, T) =P (0) +P (1) wherê ∆I (l) (x l,n , z l,n (τ))l = 1, 2, . . . , L (20) N l represents the number of random simulations on level l to estimateP (l) (l = 0, . . . ,L). The coefficient of variation of the multilevel estimateP ML f,c (0, T) can be calculated as: where var P ML f, According to Equations (7) and (18), determining the value of ∆I (l) (x, z(τ)) needs K 0 · 2 l + 1 function evaluations. Therefore, the total number of function evaluations MLMC needs is, In order to acquire an estimateP ML f,c (0, T) whose coefficient variation satisfies Equation (13), the combination of [N 0 , . . . ,N L ] must meet with: It is noticed that a great number of combinations of (N 0 , . . . , N L ) are available but computational costs associated with the combinations may be considerably different. For given K 0 and L, the combination of (N 0 , . . . , N L ) is optimal if it requires the least number of function evaluations while maintaining the coefficient variation ofP ML f,c (0, T) smaller than the prescribed threshold δ Pf . That is to say, one can obtain the optimal (N 0 , . . . , N L ) (i.e., [N * 0 , . . . , N * L ]) by solving the optimization problem below: Treating N l (l = 0, . . . , L) as continuous variables, the optimal number of random simulations on level l is easily obtained using the Lagrange Multiplier method: where V l = P (l) · (1 − P (l) ) and C l = K 0 · 2 l + 1 (27) We denote the computational cost associated with [N * 0 , . . . , N * L ] by C * M . According to the optimal allocation of computational resource determined by Equation (26), one can implement the MLMC method and estimate the cumulative probability of failure in a way that minimizes the computational complexity.

Implementation of MLMC
Equation (26) indicates that the knowledge of V l and P (L) f,c (0, T) is essential for calculating the optimal numbers of simulations MLMC needs to perform on each level l and judging whether or not MLMC is more efficient than crude MCS by comparing C S (Equation (15)) with C * M (Equations (23) and (26)). However, values of V l and P (L) f,c (0, T) are unavailable before performing TRA. Therefore, it is difficult to make a decision in practical engineering on employing MLMC or crude MCS to assess the time-variant reliability of a mechanical structure or system, and the optimal allocation of computational cost is hardly obtained in advance if MLMC is adopted.
In order to address the issues above, this subsection proposes a procedure for TRA to adaptively choose the more efficient method between crude MCS and MLMC. The flowchart is given by Figure 1. It consists of eight steps: ability 2021, 13, x FOR PEER REVIEW 8 of 18 Under the assumption made above, the computational cost of crude MCS and MLMC can be pre-estimated according to Equations (15), (23) and (26).
Step 4 If , MLMC is more efficient and adopted to estimate the cumulative probability of failure. Continue the procedure to Step 5. Otherwise, crude MCS is adopted and the procedure goes to Step 7.
and the corresponding f cov P according to Equations (19) and (21), respectively. If T is accurate enough, the procedure stops. The last estimate of the cumulative probability of failure is regarded as the result of the procedure; otherwise, go back to Step 5.
Step 7 Estimate the cumulative probability of failure using crude MCS according to Equations (9) and (14).
Step 8 Output the estimate of cumulative probability of failure.  Step 1 Set values of parameters, including K 0 , L (or K) and δ Pf . This paper sets K 0 = 10. It will be verified in Section 4 that K 0 has little effect on the efficiency of MLMC. One can set L (or K) and δ Pf according to the accuracy requirement of TRA. According to Equation (10), the accuracy of MLMC as well as crude MCS increases with L and decreases with the threshold of the coefficient of variation δ Pf .
Step 2 Discretize the time interval of interest [0,T] into K 0 · 2 l + 1 uniform instants (l = 0, . . . ,L). The relationship between the accuracy level of P (L) f,c (0, T) and the largest number of discretized time instants K 0 · 2 L + 1 is still an open issue, which is not discussed in this paper.
Step 3 Pre-estimate the number of function evaluations of crude MCS and MLMC. According to Equations (15), (23) and (26), P (l) (l = 0, . . . , L) are essential for calculating C S and C * M . However, considerable samples are needed if one independently estimates each P (l) (l = 0, . . . , L). In order to reduce the computational cost of the pre-estimate, this step makes an assumption that P (0) ≈ P (1) ≈ . . . ≈ P (L) . Thus, only P (0) is needed to pre-estimate the computational cost of crude MCS and MLMC.
Generate random sample set S 0 : Then, estimate P (0) as follows: This step sets N 0 meeting with the following accuracy requirement: Under the assumption made above, the computational cost of crude MCS and MLMC can be pre-estimated according to Equations (15), (23) and (26).
Step 4 IfĈ S >Ĉ * M , MLMC is more efficient and adopted to estimate the cumulative probability of failure. Continue the procedure to Step 5. Otherwise, crude MCS is adopted and the procedure goes to Step 7.
Step 5 Implement MLMC as follow: (a) Calculate [N * 0 , . . . , N * L ] according to Equation (26) and the currentP (l) (l = 0, . . . , L). Considering thatP (l) pre-estimated in Step 3 may have a large error, we provisionally set the accuracy requirement equal to 2δ Pf in the first iteration and restore it to δ Pf after the first iteration.
(b) Supplement S l with random samples and increase the number of samples in S l to N * l (l = 0, . . . , L).
Step 6 CalculateP ML f,c (0, T) and the corresponding cov Pf according to Equations (19) and (21), respectively. If cov Pf < δ Pf ,P ML f,c (0, T) is accurate enough, the procedure stops. The last estimate of the cumulative probability of failure is regarded as the result of the procedure; otherwise, go back to Step 5.
Step 7 Estimate the cumulative probability of failure using crude MCS according to Equations (9) and (14).
Step 8 Output the estimate of cumulative probability of failure.

Numerical Validation
This section aims to validate the computational efficiency of the method proposed in Section 3. In the validation, the proposed method is only compared with crude MCS in terms of the total number of function evaluations. As a large number of function evaluations is needed by both methods, examples studied in this section have explicit performance functions.

A Steel Bending Beam
As depicted by Figure 2, the first example studies a steel bending beam whose cross section evolves linearly in time because of isotropic corrosion. The beam is submitted to two loads, i.e., a dead load p = ρ st b 0 h 0 caused by the gravity and a pinpoint load F(t) modeled by a Gaussian process. The time-variant performance function of the beam is given by: where X = [f y ,b 0 ,h 0 ] T . The time interval of interest is (0, 20 year). Table 1 gathers the input parameters of the time-variant performance function.
Section 3. In the validation, the proposed method is only compared with crude MCS i terms of the total number of function evaluations. As a large number of function evalua tions is needed by both methods, examples studied in this section have explicit perfor mance functions.

A Steel Bending Beam
As depicted by Figure 2, the first example studies a steel bending beam whose cros section evolves linearly in time because of isotropic corrosion. The beam is submitted t two loads, i.e., a dead load p = ρstb0h0 caused by the gravity and a pinpoint load F(t) mod eled by a Gaussian process. The time-variant performance function of the beam is give by: (31 where X = [fy,b0,h0] T . The time interval of interest is (0, 20 year). Table 1 gathers the inpu parameters of the time-variant performance function.    Table 2. According t the presented results, the proposed method remarkably reduces the total number of func tion evaluations while both methods achieve the level of accuracy satisfying the threshol of the coefficient of variation δPf. It is worth mentioning that all results coming from th proposed procedure are achieved using MLMC. Integrating crude MCS into the propose procedure is to account for special cases in which MLMC needs more computational cos even though we have not yet encountered such cases.  This subsection regards crude MCS with 10 9 random simulations as the reference to compute P (l) and P (l) f,c (0, T) (K 0 = 5 and L = 8), as shown by Figure 3. It can be found that P (l) decreases dramatically with the increase of l. To investigate the performance of the proposed method, it is implemented with several combinations of (K, δ Pf , K 0 ). Results coming from the proposed MLMC and crude MCS are summarized in Table 2. According to the presented results, the proposed method remarkably reduces the total number of function evaluations while both methods achieve the level of accuracy satisfying the threshold of the coefficient of variation δ Pf . It is worth mentioning that all results coming from the proposed procedure are achieved using MLMC. Integrating crude MCS into the proposed procedure is to account for special cases in which MLMC needs more computational cost even though we have not yet encountered such cases.
According to the theory of MLMC, information about P (l) (or P (l) f,c (0, T), l = 0, . . . , L) is needed to optimize the number of random simulations allocated to each level of MLMC. However, the information is unavailable in advance. Thus, we propose a method in Section 3.2 to pre-estimate P (l) . Based on the pre-estimate, the computational cost of crude MCS and MLMC is compared and the size of random population MLMC allocates to each level l is optimized. Therefore, in order to validate the effectiveness of the proposed procedure, crude MCS, the proposed MLMC and the optimal MLMC in terms of total number of function evaluations are also compared by Figure 4, and the numbers of random simulations on each level determined according to the proposed procedure and the optimal MLMC are depicted by Figure 5. Figure 4 indicates that the optimal MLMC and the proposed MLMC are remarkably more efficient than crude MCS. The total number of function evaluations of the proposed method does not significantly increase with K, whereas that of crude MCS linearly increases with K. According to Figure 5, the computational cost allocated according to the proposed procedure is very close to that of the optimal MLMC, i.e., the true optimal allocation.    lations on each level determined according to the proposed procedure and the optimal MLMC are depicted by Figure 5. Figure 4 indicates that the optimal MLMC and the proposed MLMC are remarkably more efficient than crude MCS. The total number of function evaluations of the proposed method does not significantly increase with K, whereas that of crude MCS linearly increases with K. According to Figure 5, the computational cost allocated according to the proposed procedure is very close to that of the optimal MLMC, i.e., the true optimal allocation.

A Cantilever Tube Structure
The second example studies the cantilever tube structure shown by Figure 6. Four time-variant loads, i.e., F1(t), F2(t), P(t) and T(t), are applied to the structure. The yield strength of the structure linearly decreases with time caused by material degradation. The performance function is defined as: Figure 5. The allocation of computational cost determined by the optimal MLMC and the proposed method ((K, δ Pf , K 0 ) = (1280, 0.03,10)).

A Cantilever Tube Structure
The second example studies the cantilever tube structure shown by Figure 6. Four time-variant loads, i.e., F 1 (t), F 2 (t), P(t) and T(t), are applied to the structure. The yield strength of the structure linearly decreases with time caused by material degradation. The performance function is defined as: where σ x (t) = F 1 (t) sin(θ 1 )+F 2 (t) sin(θ 2 )+P(t) Sustainability 2021, 13, x FOR PEER REVIEW 13 of 18 1 θ 2 θ Figure 6. The cantilever tube structure (example 2). Table 3. The input parameters (example 2).  Table 4. According to Figure 3 in Section 4.1 and Figure 7, P (l) with smaller l (e.g., l < 5) contributes most of the cumulative probability of failure ( ) f,c (0, ) L P T . Considering that the computational cost on lower levels is much less than that on the highest level, optimizing the allocation of the  Table 3.

Autocorrelation Coefficient Function
The benchmark values of P (l) and P  cost on lower levels is much less than that on the highest level, optimizing the allocation of the computational cost is necessary in order to reduce the variance of MCS in a way that minimizes the over computational cost, which is exactly the purpose of the proposed procedure. The performance of the proposed method considering different combinations of (K, δPf, K0) is investigated and compared to crude MCS. Results are summarized in Table 4. It can be noticed from Table 4 that achieving similar levels of accuracy, the proposed method needs significantly less evaluations of performance function than crude MCS. The saving of the proposed method increases with the accuracy requirement (small δPf means high accuracy requirement). Figure 8 compares crude MCS, the proposed MLMC and the optimal MLMC. It shows that MLMCs hold considerable advantage over crude MCS in terms of computational cost. The allocation of computational cost determined according to the proposed method is near to that of the optimal MLMC, which is concluded from Figure 9.   The performance of the proposed method considering different combinations of (K, δ Pf , K 0 ) is investigated and compared to crude MCS. Results are summarized in Table 4. It can be noticed from Table 4 that achieving similar levels of accuracy, the proposed method needs significantly less evaluations of performance function than crude MCS. The saving of the proposed method increases with the accuracy requirement (small δ Pf means high accuracy requirement). Figure 8 compares crude MCS, the proposed MLMC and the optimal MLMC. It shows that MLMCs hold considerable advantage over crude MCS in terms of computational cost. The allocation of computational cost determined according to the proposed method is near to that of the optimal MLMC, which is concluded from Figure 9.

Conclusions
This paper presents a random simulation method called MLMC for TRA to enhance the efficiency of MCS. The proposed method discretizes the time interval of interest using a geometric sequence of timesteps and then estimates the cumulative probability of failure  (17) and (19)) coming from all levels. Considering the computational complexity of each level (Cl, Equation (27)) and the contribution that each P (l) makes to ( ) f,c (0, ) L P T , the allocation of computational cost is optimized (Equation (26)).

Conclusions
This paper presents a random simulation method called MLMC for TRA to enhance the efficiency of MCS. The proposed method discretizes the time interval of interest using a geometric sequence of timesteps and then estimates the cumulative probability of failure  (17) and (19)) coming from all levels. Considering the computational complexity of each level (Cl, Equation (27)) and the contribution that each P (l) makes to ( ) f,c (0, ) L P T , the allocation of computational cost is optimized (Equation (26)).

Conclusions
This paper presents a random simulation method called MLMC for TRA to enhance the efficiency of MCS. The proposed method discretizes the time interval of interest using a geometric sequence of timesteps and then estimates the cumulative probability of failure P (L) f,c (0, T) (Equation (17)) associated with the finest level by computing corrections P (l) (Equations (17) and (19)) coming from all levels. Considering the computational complexity of each level (C l , Equation (27)) and the contribution that each P (l) makes to P (L) f,c (0, T), the allocation of computational cost is optimized (Equation (26)).
In the proposed procedure of MLMC, the switch between MLMC and crude MCS is proposed, which is merely to circumvent the complex deduction of conditions that MLMC is more efficient than crude MCS. The switch performs practically no function in Section 4 because all pre-estimates of the computational cost of crude MCS and MLMC show that MLMC needs less function evaluations, which indicates that MLMC is commonly more efficient. As the random trajectories of the time-variant performance function are generated in the same way with crude MCS, the proposed method retains most of the characteristics of robustness with respect to nonlinearity and dimensions of performance function. However, since a geometric sequence of different timesteps is essential for the proposed methods, the number of discretized time instants is limited by K 0 and L (Equation (4)).
The effectiveness of the proposed method is validated by two examples. Results demonstrate that the proposed method remarkably reduces the total number of function evaluations while maintaining the accuracy, meeting the prescribed requirement. The saving increases with the accuracy requirement. Moreover, results indicate that the total number of function evaluations of the proposed method increases slowly with K. This characteristic means that one can use large K to reduce the discretization error (Equation (10)) without a significant increase in the computational cost, like crude MCS. The proposed method can be combined, in principle, with surrogate models (e.g., the Kriging model and polynomial chaos expansion) and variance reduction techniques (e.g., IS and SS for TRA) to obtain great savings.

Conflicts of Interest:
The authors declare no conflict of interest.