Kriging Model for Time-Dependent Reliability: Accuracy Measure and Efficient Time-Dependent Reliability Analysis Method

As the performance function of a mechanical structure is usually based on time-consuming computer codes, predicting time-dependent reliability analysis requires a large number of costly simulations in engineering. To reduce the number of evaluations of time-consuming models and enhance efficiency of time-dependent reliability analysis, Kriging is employed as a surrogate of original performance function. Firstly, a quantitative measure of the error of Kriging-based estimation of time-dependent failure probability is obtained by derivation. Dividing the quantitative measure by the Kriging-based estimation, the associated relative error is approximately estimated. Secondly, to construct an accurate Kriging model using fewer samples, a DoE (the design of experiments) strategy is developed. The idea is to adaptively refresh Kriging model with the best next sample that could enhance Kriging model the most with regard to expectation. Finally, a Kriging-based time-dependent reliability analysis method is constructed. In the method, Kriging model is adaptively refreshed according to the proposed DoE strategy, until the relative error of estimated probability of failure below a given threshold. The threshold could be quantitatively adjusted to accuracy requirement of reliability analysis. The proposed method can predict the evolution of failure probability over time. Its advantage is validated by three numerical examples.


I. INTRODUCTION
A great many of mechanical structures are subject to time dependent uncertainties, such as time-dependent loadings and decay of material properties. Their probabilities of failure generally increase over service time [1]. Herein, the probability of failure refers to the probability that a structure will fail according to its performance function. Time-dependent structural reliability analysis aims at computing the evolution of failure probability in time accounting for uncertainties involved in performance function. Owing to the time dimension, time-dependent problems are much more complicated than time-independent ones.
In time-dependent structural reliability analysis, the most general performance function is defined as G(X, Z(t), t). X and {Z(t)} are time-independent random variables The associate editor coordinating the review of this manuscript and approving it for publication was Cristian Zambelli .
(X = [X 1 , X 2 , . . . , X M 1 ] T ) and stochastic processes (Z(t) = [Z 1 (t), Z 2 (t), . . . , Z M 2 (t)] T ), respectively [2]. Without loss of generality, this study defines the failure of a structure as G(X, Z(t), t) ≤ 0. In engineering, G(X, Z(t), t) is usually defined based on time-consuming computer codes (such as finite element models), and the number of calls to time-consuming performance function is limited during predicting the evolution of failure probability. Therefore, the most important metric measuring the efficiency of a time-dependent reliability analysis method is the number of calls to performance function G(X, Z(t), t). Two categories of methods have been proposed in the past decades, i.e. outcrossing-based methods and extreme-value based methods.
As subject to stochastic processes {Z(t)} and time t, G(X, Z(t), t) is a stochastic process. Consequentially, the distribution of the first time to failure is related to the mean number of outcrossings that G(X, Z(t), t) reaches the limit state surface G(X, Z(t), t) = 0 from the safe domain G(X, Z(t), t) > 0.
Outcrossing-based methods have been developed to predict the evolution of failure probability. Based on the Rice's formula, studies derived the analytical outcrossing rates for stationary Gaussian processes [3], [4], non-stationary Gaussian processes [5], [6], non-stationary lognormal processes [7], etc.. However, stochastic process {G(X, Z(t), t)} of an engineering structure is always more general, and the associated outcrossing rate is hardly obtained. Asymptotic methods (AsM) approximately compute outcrossing rates without any assumption on the stochastic process {G(X, Z(t), t)} [8], [9]. Tvedt et al. [10] formulated the mean outcrossing rate with a sensitivity measure of the reliability of a parallel system. Inspired by [10], Andrieu-Renaud et al. [11] and Sudret [12] calculated the outcrossing rate by solving a two-component parallel-system-reliability problem, which was the widely referred PHI2. Compared with AsM, PHI2 was apparently more accurate. Jiang et al. [13] defined an outcrossing rate for time-dependent systems that had multiple failure modes and provided a method to calculate it. Most of outcrossing-based methods were built based on the assumption that outcrossings were mutually independent, which may result in overestimating failure probability [14]. To relax the assumption, Hu and Du [15] employed the joint outcrossing rate firstly proposed by Madsen and Krenk [16] and remarkably enhanced the accuracy of the prediction of time-dependent failure probability.
Within a predefined time interval, failure occurs if the extreme value of G(X, Z(t), t) exceeds a given threshold. Extreme value methods [17], [18] were to research the distribution of the extreme value. Hu and Du [2] addressed the reliability analysis of structures subject to only one stochastic process and proposed a sampling approach to determine the distribution of extreme value of G(X, Z(t), t). To reduce the number of calls to G(X, Z(t), t), surrogate models, including Hu and Mahadevan [19], Drignei et al. [20], polynomial response surface [21], and polynomial chaos expansion (PCE) [22], were usually used. Singh et al. [23] introduced the so-called ''composite'' limit state, and combined it with the niching genetic algorithm and lazy learning meta-model. By employing the total probability theorem and composite limit state, Mourelatos et al. [24] dealt with the timedependent problems involving time-independent variables and time t. Wang and Wang [25], [26] used Kriging model and efficient global optimization (EGO) algorithm [27] to construct the nested extreme response surface of time corresponding to the extreme value of performance function. Hu and Du [28] improved the method and proposed a mixed EGO-based method, in which the time-independent method AK-MCS (the adaptive Kriging and Monte Carlo simulation) [29] was used. Hu and Mahadevan [19] developed a single-loop Kriging (SILK) surrogate modeling method to further reduce the number of calls to G(X, Z(t), t). Wang and Chen [30], [31] adapted the time-independent reliability analysis method based on cumulative confidence level (CCL) [32] to time-dependent reliability analysis.
Besides outcrossing and extreme-value based methods mentioned above, there are some methods hard to sort. Jiang et al. [33], [34] predicted the time-dependent probability of failure of structures using stochastic process discretization. Jiang et al. [35] also proposed the non-probabilistic convex model process which was available for time-dependent reliability analysis. Hawchar et al. [22] combined principal component analysis and PCE, and the evolution of failure probability could be obtained with an unique reliability analysis. Wang et al. [36] employed the convex process to deal with problems with insufficient timedependent uncertainties. Importance sampling (IS) [37] and subset simulation (SS) [38] were available for some special cases.
As reviewed above, outcrossing-based methods predict the evolution of failure probability, while extreme value-based methods only estimate the cumulative failure probability over a given time interval. In terms of the number of calls to performance function, extreme value-based methods [25], [26], [28], [39] are much more efficient than outcrossing-based ones [11], [12], [22], [40]. To obtain the evolution of failure probability, existing methods always need evaluate performance functions thousands of times [11], [12], [22]; this is computationally prohibitive when time-consuming performance functions are involved. The present study aims to predict the evolution of failure probability with much less calls to performance function. Herein, Kriging model constructed based on limit training points is employed as the surrogate of performance function G(X, Z(t), t). Reasons why we adopt Kriging are summarized as follows: 1) Kriging is an interpolation model; 2) A well trained Kriging model can not only predict the performance value at any untrained point, but also provide the mean square error (MSE) of its prediction; 3) It has already been widely applied to time-independent reliability analysis of structures whose performance functions are time-consuming, which is reviewed below.
As introduced above, the computational demand associated with reliability analysis depends on how many times performance function is evaluated when performance function is based on time-consuming codes. To release computational burden, researches focus on adaptive DoE (the design of experiment) strategy, i.e., constructing accurate Kriging model using fewer DoE points. For time-independent problems, learning function-based adaptive Kriging models have gained much popular. The basic idea is adaptively enhancing the accuracy of Kriging model with the maximum (or minimum) point of a learning function until meeting with a predefined stopping criterion. The maximum (or minimum) point, also called as the best next point, is expected to maximally enhance the Kriging. Several learning functions have been proposed from different perspectives, including U function [41], the expected feasibility function (EFF) [29], cumulative confidence level (CCL) [32], the expected risk function (ERF) [42], information entropy-based function H [43], least improvement function [44], the reliability-based VOLUME 8, 2020 expected improvement function [45], and error rate-based adaptive Wang and Shafieezadeh [46]. In spite of their differences, all learning function intends to enhance Kriging with points in the vicinity of limit state surface. Li and Xiu [47], Li et al. [48] proposed a hybrid approach by sampling both the surrogate model in a ''large'' portion of the probability space and the original system in a ''small'' portion. By combining the idea with IS, their method was capable of capturing failure probability as small as 10 −6 [48]. Besides adaptive DoE strategies, suitable stopping criteria are also helpful to reduce the number of calls to performance function. Gaspar et al. [49], [50] exploited the stabilization of the Kriging-based estimation of failure probability during the active refinement procedure, and provided a compromise between the accuracy and computational cost. Balesdent et al. [51] proposed a Kriging-based adaptive IS algorithm and used the statistical properties of Kriging to compute a confidence measure for failure probability estimation. Authors researched accuracy measures of Kriging model in terms of time-independent reliability analysis [52], [53], and proposed the stepwise accuracyimprovement strategy (SAIS) to adaptively enhance Kriging model [54].
Kriging model applied in the field of time-dependent reliability analysis still needs to be perfected, even though it is widely adopted in time-independent reliability analysis and performs outstandingly. To efficiently and accurately predict the evolution of structural failure probability over time, this study systematically researches of the Kriging-based method for time-dependent reliability analysis based on the previous works [44], [52], [54]. Herein, Kriging is employed to surrogate G(X, Z(t), t) which is usually time-consuming in engineering. A quantitative accuracy measure of Krigingbased estimation of failure probability is derived. Based on the derivation, a qualitative accuracy measure of Kriging model in the context of time-dependent reliability is proposed. Referring to SAIS [54], a DoE strategy, named as timedependent SAIS (tSAIS), is built. tSAIS iteratively refreshes Kriging model with the best next sample. Finally, a Krigingbased time-dependent reliability analysis method, which can predict the evolution of cumulative failure probability over time, is developed. In the method, to construct a well-trained Kriging model with fewer DoE samples, tSAIS is employed. The stopping criterion of the proposed method considers the relative error of the Kriging-based estimation of failure probability.
The remainder of this paper is organized as follows. Section II introduces the background concepts of time-dependent reliability analysis, Monte Carlo simulation (MCS), and Kriging model. The quantitative accuracy measure of Kriging-based predictor of cumulative failure probability is derived in Section III, which is followed by the proposed tSAIS strategy in Section IV. A time-dependent reliability-analysis method is presented in Section V, and the advantage of the method is validated in Section VI. Section VII provides conclusions.

A. PROBLEM STATEMENT
As illustrated in Section I, this study defines the failure of a structure as G(X, Z(t), t) ≤ 0, where G(X, Z(t), t) is the timedependent performance function. The cumulative probability of failure within [t 0 , t] (t ≥ t 0 ) is expressed as, where t 0 is the starting time instant of interest. Apparently, cumulative failure probability P f,c (t 0 , t) considers the evolution uncertainty of {G (X, Z(·), ·)} within [t 0 , t]. This study is to propose a Kriging-based method to predict the change of P f,c (t 0 , t) over time interval MCS is the most robust method for structural reliability analysis. However, it needs a great many of independent and identically distributed (i.i.d.) samples of random vector X and discretized trajectories of {Z(t)} on {t k } (k = 0, . . . , K ). The spectral decomposition and Karhunen-Loeve expansion-based techniques are available to discretize {Z(t)} (see [30], [55]- [57] for details). Discretization techniques are not discussed here. This study assumes that the discretized trajectories of {Z(t)} are generated by an existing technique and they are shown along with i.i.d. realizations of X as follows: where N MC is the number of MCS samples. {z n,0 , z n,1 , . . . , z n,K } is a discretized trajectory of {Z(t)}. Each sample [x MC,n , z n,k , t k ] (n = 1, . . . , n MC ; k = 0, . . . , K ) in (3) is associated with a performance value G(x MC,n , z n,k , t k ), as (4), shown at the bottom of the next page. Then, by using MCS, (2) can be further simplified as where I f (n, t k ) is the failure indicator function, presented as The coefficient of variation ofP f,c (t 0 , t k ) is According to (6) and (7), I f (n, t k ) andP f,c (t 0 , t k ) are apparently non-decreasing with the increase of t k . Therefore, δ MC (t 0 , t k ) is monotonically non-increasing with respect to t k , and satisfies As this study concerns the cumulative failure probabilities

C. KRIGING
According to Section II A, MCS needs to evaluate performance function (K + 1)N MC times to estimate the evolution For an implicit G(X, Z(t), t) based on time-consuming codes, hours even days are required to evaluate G(X, Z(t), t) once, thus, MCS is computationally prohibitive in engineering. The present study employs Kriging model to overcome the aforementioned situation. As illustrated in Section I, Kriging is an interpolation model. It can estimate the performance value at any point, as well as provide the associated MSE which is called the Kriging variance. In Kriging model, the performance value at any point is a variable with epistemic uncertainty and is subject to normal distribution. The mean and variance are the Kriging predictor of performance value and the Kriging variance, respectively.
Given a DoE S DoE N including N samples and the corre- The Kriging predictor of the performance value at [x, z, t] is calculated as follows:  [58]- [60] for details). R and r x are written as R(·, ·; θ) is the correlation function. This study employs the Gaussian correlation function, expressed as The MSE ofĜ(x, z, t) is Equation (17) measures the local accuracy of the Kriging model. σ 2 is estimated using generalized least squares.
According to the theory of Kriging model, the epistemic uncertainty of G(x, z, t) is normally distributed as To avoid evaluating performance model too many times, this study predicts performance values in (4) with Kriging model, as shown by (19), at the bottom of the next page.
By replacing G(·, ·, ·) in (5)- (7) with Kriging predictor G(·, ·, ·), Kriging-based estimations of cumulative probabilities of failure {P f,c (t 0 , t k )} (k = 0, 1, . . . , K ) are obtained. where, As shown by (20), the error ofP f,c (t 0 , t k ) (k = 0, 1, . . . , K ) consists of two parts: These result from: x parameters of MCS, i.e., N MC and K ; y the employed discretization technique of stochastic processes; and z the Kriging model. MCS can achieve any prescribed accuracy by increasing N MC and K . The available discretization technique of stochastic processes can handle Gaussian processes and some non-Gaussian processes. The former two factors directly influ- | caused by the Kriging model, which is barely explored in literatures.

A. A SIMPLIFICATION
According to (7), Based on [52], [54] and (18), the probability that Kriging model inaccurately predicts the sign of G(x MC,n , z n,k , t k ) is where (·) is the cumulative distribution function of the standard normal distribution, and Thus, it is easily concluded that, The probability that the signs of G(x MC,n , z n,k , t k ) and G(x MC,n , z n,k , t k ) are contradictory is negligible as U N tends to infinity.
To simplify the derivation, a threshold value [U ] is introduced. If U N (x MC,n , z n,k , t k ) meets (26), this study approximates (22) with (27).

B. ACCURACY MEASURE
By considering (7) and (20), the error ofP f,c (t 0 , t k ) caused by the inaccuracy of Kriging model is expressed below, Given S DoE and Y DoE , G(x MC,n , z n,k , t k ) is normally distributed according to (18). Therefore, failure indicator I f (n, t k ) (Equation (6)) is subject to Bernoulli distribution. By implementing expectation operator on both sides of (28), an upper bound of the expectation of |P f,c (t 0 , t k )−P f,c (t 0 , t k )| is obtained.
As E |P f,c (t 0 , t k ) −P f,c (t 0 , t k )| is computationally complex, we propose the right side of (29) as a quantificational measure of |P f,c (t 0 , t k ) −P f,c (t 0 , t k )| and denote it by e R (t 0 , t k ). where, According to (6), (21) and (24), Therefore, where, As authors have derived in [52], [G(x MC,n , z n,1 , t 1 ), . . . , G(x MC,n , z n,k , t k )] T is k-variate normally distributed. P (I f (n, t k ) = 0) and P (I f (n, t k ) = 1) are easily calculated with MATLAB. However, the expression of e R (t 0 , t k ) consists of N MC terms, and the evaluation time needed may not be satisfying. As far as authors know, for a given Kriging model, most of [x MC,n , z n,k , t k ] (n = 1, . . . , n MC ; k = 0, . . . , K ) meets (26), i.e. most of signs of their performance values are almost confirmed. Thus, we further simplify the calculation of E |I f (n, t k ) −Î f (n, t k )| . According to the simplification in Section III A, we rewrite (31) and discuss it in four cases, Case 1 (Fig. 1) andÎ f (n, t k ) = 0 According to (26) and (27), Case 2 (Fig. 2 However, at least one of them is located in the relative safety domain According to (31) and (32), The last approximation of (35) is based on (26) and (27). In this case, there exists at least one sample of [x MC,n , z n,i , Kriging model may improperly predict the sign(s) of the associated performance value(s). Case 3 (Fig. 3 According to (31) and (32), Case 4 (Fig. 4): min andÎ f (n, t k ) = 1 According to (31) and (32), The last approximation of (37) is based on (26)  Based on (33)-(37), all terms E |I f (n, t k ) −Î f (n, t k )| (n = 1, . . . , N MC ) in (30) can be approximately calculated. Then, e R (t 0 , t k ), i.e. the proposed quantificational measure of |P f,c (t 0 , t k ) −P f,c (t 0 , t k )|, is obtained.

IV. THE PROPOSED DoE STRATEGY tSAIS
The derivation in Section III indicates that only a part of samples in (3) contribute to e R (t 0 , t k ), and e R (t 0 , t k ) and the set of the samples that make the contribution vary with t k . Consequentially, none of {e R (t 0 , t k )} (k = 0, . . . , K ) is suitable for measuring the accuracy of Kriging model. Firstly, this section proposes an accuracy measure of Kriging model in the context of time-dependent reliability analysis. Then, the best next sample is defined as the sample that can enhance Kriging model the most in terms of expectation. Finally, the algorithm for determining the best next sample is developed.

A. ACCURACY OF KRIGING MODEL
To appropriately measure how accurate a Kriging model is in the context of time-dependent reliability analysis, two factors are considered: x the samples in (3) that may disturb values of {P f,c (t 0 , t k )} (k = 0, . . . , K ) and y the probability that signs of their performance values are misclassified.
According to the simplification in Section III A and the derivation in Section III B, the sample [x MC,n , z n,i , t i ] that may disturb the accuracy of {P f,c (t 0 , t k )} (k = 0, . . . , K ) must simultaneously satisfy (38) and (39).
The sign of G(x MC,n , z n,i , t i ) is uncertain only if it satisfies (38). Equations (39) and (26)  Let S U represent the sample set that collects the samples in (3) satisfying both (38) and (39).
In the framework of Kriging, G(x, z, t) is subject to normal distribution (Equation (18)). The probability that Kriging model wrongly predicts the sign of the performance value at [x MC,n , z n,i , t i ] ∈ S U is −U N (x MC,i , z n,i , t i ) . By considering both S U and −U N (x MC,i , z n,i , t i ) , this study proposes e N shown by (41) as the accuracy measure of Kriging model in the context of time-dependent reliability analysis.
It is not difficult to identify that e N is in positive correlation with {e R (t 0 , t k )}. e N qualitatively reflects the accuracy degree of the Kriging model as well as that of {P f,c (t 0 , t k )} (k = 0, . . . , K ).

B. THE PROPOSED DoE STRATEGY tSAIS
It has been well demonstrated that an excellent DoE strategy remarkably improves the efficiency of structural reliability analysis [19], [25], [29], [54]. Referring to SAIS for timeindependent reliability analysis [54], this subsection first derives the degree that an untried sample could enhance Kriging model with respect to the proposed accuracy measure in Section IV A. Then, the best next sample is defined as the sample that maximizes the degree. Finally, the best next sample is searched.

1) THE BEST NEXT SAMPLE
By adding an untried sample [x , z , t ] and G(x , z , t ), the current DoE S DoE N (Equation (10)) and Y DoE N (Equation (11)) are refreshed.
A new Kriging model is then constructed based on S DoE  N +1 ([x , z , t ]) and Y DoE  N +1 G(x , z , t ) . The accuracy measure proposed in Section IV A is employed to measure the quality of the new Kriging model, where S U is defined by (40). For the Kriging model based x refreshing the Kriging model with a new sample would generally shrink S U , i.e., (44) tends to be conservative, and y the reconstruction of S U would result in massive computation of the proposed tSAIS.
The enhancement of the Kriging model due to [x , z , t ] and G(x , z , t ) is given as Apparently, e(·, ·) depends on both [x , z , t ] and G(x , z , t ). Without evaluating G(x , z , t ), the enhancement e(·, ·) is unknown; this is unsatisfactory. A measure depending only on [x , z , t ] is essential to measure the degree that [x , z , t ] could improve the original Kriging model.
According to (18), G(x , z , t ) is normally distributed as Thus, for given [x , z , t ], e [x , z , t ], G(x , z , t ) is a function of normal random variable G(x , z , t ). This study uses the expectation of e [x , z , t ], G(x , z , t ) (Equation (47)) to measure the enhancement of the original Kriging (47) where f G(x ,z ,t ) (y) is the probability density function (PDF) of G(x , z , t ).
The maximum point of (47) tends to maximize the enhancement defined by (45). Therefore, this study defines VOLUME 8, 2020 the best next sample as follows: (47) can be rewritten as . The Gauss-Hermite quadrature is employed to calculate the infinite integral in (50).
where n G is the number of quadrature points, and v j and w j (j = 1, . . . , n G ) denote the quadrature points and associated weights, respectively. It is foreseeable that searching for the best next sample in the whole feasible region will be time-consuming. According to the simplification in (26) and (26), Kriging model is very accurate in the domain {U N (X, Z(t), t) > [U ]}. Therefore, referring to [54], this study treats the sample set S U (Equation (40)) as the candidate pool of the best next sample [x DoE N +1 , z DoE N +1 , t DoE N +1 ], instead of searching in the whole feasible region

1) THE PSEUDOCODE OF THE PROPOSED tSAIS
In summary, tSAIS proposed herein for time-dependent reliability analysis is to search for the best next sample defined by (52). The expectation [x , z , t ] quantifies how much [x , z , t ] enhances the Kriging model in terms of the accuracy measure e N (Equation (41)); the sample set S U , collecting all the samples in (3) which may affect the accuracy of {P f,c (t 0 , t k )} (k = 0, . . . , K ), is used to construct the accuracy measure of the Kriging model (e N ) and performs as the candidate pool of the best next sample.
Given i.i.d. realizations of the input variables (3), the pseudocode of tSAIS is summarized as follows: x Construct the Kriging model based on S DoE N and Y DoE N .
y Determine the candidate pool of the best next sample S U according to (40).
{ Identify the approximate maximum point of [x , z , t ] , which is the best next sample.

V. PROPOSED TIME-DEPENDENT RELIABILITY ANALYSIS METHOD
This section presents a time-dependent reliability analysis method. The evolution of cumulative probability of failure over a time interval is evaluated by a single time-dependent reliability analysis. In the proposed method, Kriging model surrogates the performance function G (X, Z(t), t); DoEs of the initial Kriging model is generated by Latin hypercube sampling (LHS); then, tSAIS proposed in Section IV is used to iteratively refresh the DoEs and Kriging model until the stopping criterion is satisfied. The adopted stopping criterion considers the relative error of Kriging-based estimations of cumulative failure probabilities {P f,c (t 0 , t k )} (k = 0, 1, . . . , K ). The main steps are listed as follows: Step 1 Discretize the time interval [t 0 , t E ] into K equal sub-intervals [t k−1 , t k ] (k = 1, . . . , K , and t K = t E ).
Step 2 Generate N 0 initial DoE samples with LHS, calculate their corresponding performance values by evaluating G (X, Z(t), t), and construct the initial Kriging model. PDF for generating {z DoE n } (n = 1, . . . , N 0 ) is calculated as (53) [31].f where f Z(t k ) (·) represents the PDF of the random vector Z(t k ).
Step 3 Generate N MC i.i.d. realizations of X and discretized trajectories of {Z(t)}, as shown by (3). N MC should meet (54).
Step 4 Refresh the current DoEs according to the pseudocode of tSAIS.
Step 6 If (55) is satisfied, the Kriging model is converged.
This study sets [e] = 0.05. The stopping criterion (55) implies that {P f,c (t 0 , t k )} and Kriging model are accurate enough if all of the relative errors of {P f,c (t 0 , t k )} (k = 0, 1, . . . , K ) are lower than [e]. It must be emphasized that one can quantificationally adjust the threshold value [e] according to accuracy requirement.
The flowchart of the reliability analysis method is depicted in Fig. 5.

VI. VALIDATION
This section presents three numerical examples to validate the contribution of this study.
x Demonstrate the availability of the proposed accuracy measure ofP f,c (t 0 , t k ) by comparing e R (t 0 , t k ) (Equation (30)) with |P f,c (t 0 , t k ) −P f,c (t 0 , t k )|. As analyzed in Section III, the error ofP f,c (t 0 , t k ) consists of |P f,c (t 0 , t k ) − P f,c (t 0 , t k )| and |P f,c (t 0 , t k ) −P f,c (t 0 , t k )|. The former results from MCS and discretization techniques of stochastic processes. The latter depends on the accuracy of Kriging model. This study focuses on efficiently reducing|P f,c (t 0 , t k ) −P f,c (t 0 , t k )|. Section III proposes e R (t 0 , t k ) as a quantitative measure of |P f,c (t 0 , t k ) −P f,c (t 0 , t k )|, and it is also the basic of the stopping criterion of the proposed time-dependent reliability analysis method. Therefore, this section systematically investigates the availability of e R (t 0 , t k ).
y Demonstrate the advantage of the time-dependent reliability analysis method proposed. The method is constructed mainly based on tSAIS, and can predict the evolution of cumulative failure probability with only a single reliability analysis. tSAIS is designed to adaptively refresh Kriging model with the best next sample. This section validates its efficiency in terms of the number of calls to performance function by comparing with several existing methods.

A. A CORRODED BEAM SUBJECTED TO TIME-DEPENDENT STOCHASTIC LOAD
The example was first used in [11], and it was then modified by [28]. Herein, the modified version is adopted. As shown in Fig. 6, the length of the steel bending beam is L = 5m, its cross-section is rectangular, and it isotropically corrodes. It is assumed that the cross-section evolves linearly in time, where a 0 and b 0 are the initial sizes of the cross-section, ω represents the corrosion rate (ω = 5 × 10 −5 m/year). The beam is subject to two loads, i.e. gravity (the force density ρ st = 7.85 × 10 4 N/m 3 ) and pinpoint time-dependent stochastic load F(t). Reference [28] expressed F(t) using the spectral representation method, where, The steel beam fails as soon as plastic deformation occurs. Therefore, the performance function is where X = [σ u , a 0 , b 0 ] T , and σ u represents the steel yield strength. Table 1 lists all the random variables involved. They are normally distributed. This study focuses on the evolution of the cumulative failure probability over [0,35], i.e. P f,c (0, t) (0 ≤ t ≤ 35 years).
We perform the reliability analysis of the steel beam using the proposed method. The time interval [0, 35] is equally discretized into 100 subintervals (K = 100). Ten initial samples are generated with LHS (N 0 = 10). Then, the Kriging model is enhanced by tSAIS. It meets with the stopping criterion after 17 iterations (N call = 10 + 17), where N call represents the number of calls to the performance function (58). In the second column of Fig. 7, the solid lines represent   Fig. 7.
The rapid decrease of max{e R (t 0 , t k )/P f,c (t 0 , t k )} shown in Fig. 8 indicates that the best next samples obtained with tSAIS indeed enhance Kriging model and the accuracy of  To demonstrate the efficiency of the proposed method, Table 2 compares it with several existing methods. The proposed method is randomly run seven times. The results show that the proposed method is considerably more efficient than the Rice-based method, independent EGO-based method, and the mixed-EGO based method, and compatible with SILK in terms of N call . It can accurately predict the evolution of the target cumulative failure probability {P f,c (t 0 , t)} by calling to (58) about 27 times, while SILK, independent EGO-based method, and mixed-EGO based method only predict the cumulative failure probability P f,c (0, 35) rather than the evolution of {P f,c (t 0 , t)}. Further, PHI2 has already analyzed the original example of the steel beam with N call = 1774 [11], which indirectly exhibits the advantage of the proposed method.  [28], and that of SILK are from [19]. For a fair comparison, ε represents the relative error of these methods with regard to the corresponding reference values, i.e., the reference value of Rice, independent EGO, mixed EGO and SILK is 3.03 × 10 −2 [19], [28], and that of the proposed method is 3.04 × 10 −2 .  Fig. 9 shows a cantilever tube structure. It was studied by [22], [33]. The cantilever tube structure is subject to two time-independent stochastic loads (F 2 and P) and two time-dependent stochastic loads (F 1 (t) and T (t)). The material yield strength σ u is assumed to linearly decrease over time.

B. A CANTILEVER TUBE STRUCTURE
where σ 0 is the initial yield strength.
The performance function of the structure is where X = [d, h, σ 0 , F 2 , P] T and Z(t) = [F 1 (t), T (t)] T . All parameters and variables are listed in Table 3. The time interval analyzed in the example is [0, 5] (unit: year). For further details, please see [22], [33]. First, F 1 (t) and T (t) are approximately discretized using expansion optimal linear-estimation method (EOLE) (see [56] for further details). Then, we use the proposed method to predict the evolution of P f,c (0, t) (0 ≤ t ≤ 5) with K = 100 and N 0 = 15. The stopping criterion is satisfied until N call = 26 (15 + 11), as shown in Fig. 10.
The subfigures in the first column of Fig. 11 depict the convergence of  Results of PHI2 and t -PCE (polynomial chaos expansion for time-dependent reliability analysis) are from [22]. Their reference value of P f,c (0,5) is [9.46 × 10 −3 , 9.84 × 10 −3 ], while that of the proposed method is 1.01 × 10 −2 . The difference may result from the discretization techniques of stochastic processes employed. Reference [22] uses Karhunen-Loè ve (KL) expansion-based method to represent F 1 (t ) and T (t ), while the present study employs EOLE. Here, for a fair comparison, for each method, ε(%) is calculated with its corresponding reference value of P f,c (0,5).
show the decrease in {e R (t 0 , t k )/P f,c (t 0 , t k )} with the increase of N call . Table 4 summarizes results obtained using the methods from [11], [22] and the proposed method. All methods in Table 4 can predict the evolution of P f,c (0, t) over the interval [0,5]. In [33], another version of this structure was analyzed using the stochastic process discretization-based method with N call ≥ 270. Apparently, the proposed method outperforms methods mentioned above in terms of the number of function evaluations.

C. A TRUSS STRUCTURE
The truss structure shown in Fig. 12 has been analyzed in various studies to validate time-independent reliability analysis methods [44], [54], [61], [62]. Reference [22] modifies it to fit the time-dependent reliability analysis, which is employed herein. The example involves four time-independent random variables (Young's modulus E 1 and E 2 , cross-sections A 1 and A 2 ) and six time-dependent loads P 1 (t), . . . , P 6 (t) modelized by Gaussian processes. E 1 and A 1 are for horizontal bars, while E 2 and A 2 are for diagonal ones.
The performance function of the example is where X = [A 1 , A 2 , E 1 , E 2 ] T and Z(t) = [P 1 (t), . . . , P 6 (t)] T . The failure of the structure is defined as its maximum deflection δ max (t) greater than the threshold value [δ] = 0.113m. In the example, the time interval analyzed is [0, 15] (unit: year). The distribution of random variables is listed in Table 5. See [22] for further details about the truss structure. The time-dependent loads P i (t) (i = 1, . . . , 6) are approximately discretized with EOLE [56]. The proposed method is applied to perform the time-dependent reliability analysis    with K = 100 and N 0 = 15. As shown in Fig. 13 and  (N call = 15 + 34). Fig. 13 and Fig. 14 also provide details of the convergence of the proposed method. The availability of e R (t 0 , t k ), the quantitative accuracy measure ofP f,c (t 0 , t k ), is demonstrated in the subfigures in the second column of Fig. 14. Table 6 summarizes results of t-PCE [22], MCS, and the proposed method. Table 6 validates the efficiency of the proposed method in terms of the number of performance   [22]. In [22], the time-dependent loads P i (t ) (i = 1, . . . , 6) are discretized using the KL-based method, while this study uses EOLE. The reference value of P f,c (0,5) in [22] is [2.72 × 10 −2 , 2.78 × 10 −2 ], and that of this study is 2.64 × 10 −2 . function evaluations. Table 6, Fig. 13, and Fig. 14 also indicate the feasibility of the proposed stopping criterion.

VII. CONCLUSION
This study systematically researches the application of Kriging model to time-dependent reliability analysis of mechanical structures. Kriging model is used as the surrogate of performance function to reduce the number of evaluations of time-consuming models and release the computational burden of time-dependent reliability analysis of mechanical structures. The contributions of the study include: 1) deriving a quantitative accuracy measure of Kriging-based predictor of cumulative failure probability; 2) proposing an adaptive DoE strategy (tSAIS) of Kriging by adapting the idea of SAIS to time-dependent reliability analysis and 3) proposing a Kriging-based time-dependent reliability analysis method by combining tSAIS and MCS.
Three benchmark examples, including moderately highdimensional, (strongly) nonlinear and implicit performance functions, validate the contributions. According to the results presented in Section VI, it can be concluded as follows: (1) The derived accuracy measure, e R (t 0 , t k ), is highly competent to measure the error of Kriging-based estimations of time-dependent failure probability. It is actually an upper bound of the expectation of |P f,c (t 0 , t k )−P f,c (t 0 , t k )|. As visualized by Fig. 7, Fig. 10, and Fig. 13, e R (t 0 , t k ) is greater than the real error |P f,c (t 0 , t k ) −P f,c (t 0 , t k )|, and evolves consistently with |P f,c (t 0 , t k ) −P f,c (t 0 , t k )|. Therefore, it fulfills its intended function. e R (t 0 , t k ) is used to construct the stopping criterion of the proposed reliability analysis method in Section V. Fig. 7, Fig. 10 and Fig. 13 show evolutions of failure probability the proposed method predicts are quite accurate relative to the reference evolution obtained with MCS when Kriging model meets with the stopping criterion. Results summarized in Table 2, Table 4 and Table 6 indicate that relative errors ofP f,c (t 0 , t k ) the proposed method provides are all below [e] (=0.05). Results mentioned above validate the good quality of the stopping criterion based on e R (t 0 , t k ).
(2) tSAIS and the proposed reliability analysis method do well in reducing the number of calls to performance function. According to comparisons performed by Table 2, Table 4 and Table 6, the proposed method needs much less evaluations of performance function than t-PCE, PHI2, EGO-based method, and only SILK is compatible with the proposed method in terms of N call . However, the proposed method can predict the evolution of failure probability over time as shown by Fig. 7,  Fig. 10, and Fig. 13, while SILK is an extreme value-based method and only predicts the cumulative failure probability over a given time interval rather than the evolution. Therefore, the advantage of the proposed method over other methods is noticeable.
In the proposed method, the estimated failure probabilities and candidate pool of the best next sample are obtained with MCS. Therefore, its calculation efficiency is still potential to improve when dealing with small probability of failure. Adapting IS and SS to time-dependent reliability analysis and combining them with Kriging will be a part of our future work.