Current surrogate modeling methods for time-dependent reliability analysis implement a double-loop procedure, with the computation of extreme value response in the outer loop and optimization in the inner loop. The computational effort of the double-loop procedure is quite high even though improvements have been made to improve the efficiency of the inner loop. This paper proposes a single-loop Kriging (SILK) surrogate modeling method for time-dependent reliability analysis. The optimization loop used in current methods is completely removed in the proposed method. A single surrogate model is built for the purpose of time-dependent reliability assessment. Training points of random variables and over time are generated at the same level instead of at two separate levels. The surrogate model is refined adaptively based on a learning function modified from time-independent reliability analysis and a newly developed convergence criterion. Strategies for building the surrogate model are investigated for problems with and without stochastic processes. Results of three numerical examples show that the proposed single-loop procedure significantly increases the efficiency of time-dependent reliability analysis without sacrificing the accuracy.

## Introduction

Time-dependent reliability analysis computes the probability that a system or component fulfills its intended function over a time interval of interest without failures. Since it is related to system performance degradation [1,2], lifetime cost estimation [1], maintenance [3], lifetime testing [4], and system resilience [5], time-dependent reliability analysis has gained a lot of attention during the past decades. The biggest challenge in time-dependent reliability analysis is how to efficiently and accurately perform reliability analysis and account for the correlation of the response over time during the analysis.

This paper focuses on component reliability, which has been extensively studied in recent years. For instance, a PHI2 method has been developed to compute the time-variant reliability by estimating the upcrossing rate at each time instant [6]; Singh and Mourelatos [7] proposed a composite limit-state function method and investigated importance sampling [8], Markov Chain Monte Carlo simulation (MCS) [9], subset simulation [10], and total probability theory [11] for this estimation; Hagen and Tvedt [12,13] proposed a parallel system approach to solve time-dependent problems with binomial distributions; Li and Chen developed a reliability analysis method for dynamic response using a new probability density evolution approach [14]; Along with these methods, Zhang and Du [15] derived analytical expressions for the upcrossing rate of function generator mechanisms; Du [16] proposed an envelope function method based on first-order approximation; Hu and Du investigated the upcrossing rate method [17], joint upcrossing rate method [18], and the first-order sampling approach [19]; Wang and Wang [20,21] presented a nested extreme value response method to estimate the time-dependent reliability; and Jiang et al. [22] studied the strategy of time-dependent reliability assessment through stochastic process discretization.

From the above literature review, it is found that most of the current methods rely on a first-order approximation of the limit state, which is inaccurate for problems with nonlinear responses and multimodal statistical properties. In that situation, surrogate modeling-based methods become more promising than the methods based on linearization at the most probable point (MPP). In terms of surrogate model-based time-dependent reliability analysis, a representative method is the nested extreme value response method proposed by Wang and Wang [20,21]. In the nested method, the efficient global optimization (EGO) method [23] is embedded in an extreme value surrogate model to identify the extreme values. There are two loops in the nested method. The inner loop is the identification of global extremes of the response and the outer loop is used to build the extreme value surrogate model. This method is then improved by Hu and Du [24] by developing a mixed EGO method in the inner loop to improve the efficiency of identifying the extreme values and implementing an adaptive sampling approach in the outer loop to reduce the number of training points required to build the extreme value surrogate model. Using adaptive sampling in the outer loop for the extreme value, surrogate modeling has also been investigated by Wang and Wang [25]. The surrogate modeling method for time-dependent reliability analysis has also been studied by Schoefs et al. [26] and Wang et al. [27] using polynomial chaos expansion.

Analyzing the current surrogate modeling methods for time-dependent reliability analysis, it is found that current methods implement a double-loop procedure. As discussed above, the outer loop constructs the extreme value surrogate model for the latent extreme value response of the limit state function as a function of random variables and the inner loop performs the global optimization for the limit state function over the time duration of interest under given realization of the random variables. The double-loop procedure has two main drawbacks: (1) The accuracy of global optimization in the inner loop will affect the accuracy of the extreme value surrogate model in the outer loop. (2) Identifying the extreme value in the inner loop for problems with stochastic processes over a long time period is computationally expensive since the realization of the stochastic process may have numerous peaks.

This paper develops a SILK surrogate modeling method for time-dependent reliability analysis, in order to overcome the above drawbacks. The global optimizations in the inner loop of current methods are completely removed in the proposed method. It is shown that global optimization in surrogate model-based time-dependent reliability analysis is not indispensable. The surrogate modeling strategy is investigated for problems with and without stochastic processes, which makes the proposed method applicable to general time-dependent reliability analysis problems instead of a special group of problems. Thus, the contributions of this paper are summarized as: (1) a novel perspective for surrogate modeling-based time-dependent reliability analysis; (2) a strategy for surrogate modeling for time-dependent reliability analysis with stochastic processes; (3) a modified learning function for training the surrogate model.

The rest of the paper is organized as follows. Section 2 provides background concepts in time-dependent reliability analysis and related surrogate modeling approaches. Section 3 presents the proposed SILK surrogate method. Three numerical examples are used to demonstrate the proposed method in Sec. 4. Concluding remarks are provided in Sec. 5.

## Background

### Time-Dependent Reliability Analysis.

Let $X=[X1, X2, ⋯, Xn]$ be a vector of random variables, $Y(t)=[Y1(t), Y2(t), ⋯, Ym(t)]$ be a vector of stochastic processes, and $G(t)=g(X, Y(t), t)$ be the limit state function. The time-dependent probability of failure is given by
$pf(t0, te)=Pr{g(X, Y(t), t)>0, ∃t∈[t0, te]}$
(1)

where $[t0, te]$ is the time interval of interest, $t0$ is the initial time instant, $te$ is the last time instant, “ $Pr{⋅}$ ” is probability, and “ $∃$ ” means “there exists.”

The most commonly used method for solving Eq. (1) uses Rice's formula for upcrossing rate [28]. $pf(t0, te)$ is estimated using the upcrossing rate as [17]
$pf(t0, te)=1−R(t0)exp {∫t0tev+(t)dt}$
(2)

in which $R(t0)$ is the reliability at the initial time instant and $v+(t)$ is the upcrossing rate at time instant $t$. $v+(t)$ is often estimated using the first-order reliability method (FORM) based on Rice's formula [17].

The upcrossing rate method based on Rice's formula is accurate when the failure probability is low and FORM is accurate for time-instantaneous reliability analysis. When FORM is not accurate due to the nonlinear response of the system, surrogate modeling-based method is a promising way.

### Surrogate Modeling for Time-Dependent Reliability Analysis.

Current surrogate modeling methods for time-dependent reliability analysis are mainly developed for problems with limit-state functions $G(t)=g(X, t)$. These methods are based on the following probability equivalency [24]:
$pf(t0, te)=Pr{g(X, t)>0, ∃t∈[t0, te]}=Pr{gmax(X)>0}$
(3)
in which $gmax(X)$ is the extreme value response on $[t0, te]$, for any given $X=x$, where $x$ is a realization of random variables $X$, $gmax(x)$ is given by
$gmax(x)=maxt∈[t0, te]{g(x, t)}, ∀X=x$
(4)

Since the extreme value response is an unknown function, a surrogate model $ĝmax(X)$ is built and the time-dependent reliability is estimated based on that. The key issue is how to build $ĝmax(X)$. Due to the requirement of the extreme value responses, the time-dependent reliability estimate based on Eqs. (3) and (4) is a double-loop procedure, which is summarized as below.

• Outer loop: Build a surrogate model $ĝmax(X)$.

• Inner loop: Identify $gmax(x)$ for given $X=x$ (as shown in Fig. 1
Fig. 1

Identification of extreme value response ($gmax(x)$) for given $X=x$

Fig. 1

Identification of extreme value response ($gmax(x)$) for given $X=x$

).

Next, we will briefly review the surrogate modeling methods, which implement the double-loop procedure.

#### Double-Loop Procedure With Independent Maxima.

The double-loop procedure with independent maxima refers to the method presented in Ref. [20], which identifies the maxima independently in the inner loop using the EGO method [23]. In the rest of this paper, we call it the independent EGO method for the sake of illustration. The basic idea of the independent EGO method is summarized as follows:

• Outer loop: Generate training points $xs=[x(1), x(2), ⋯, x(nt)]$ for $X$ and build a surrogate model $ĝmax(X)$ based on $x(i)$ and $ĝmax(x(i))$, $i=1, 2, ⋯, nt$.

• Inner loop: For each training point $x(i)$, build a surrogate model $ĝi(t)$ to identify $ĝmax(x(i))$ independently.

#### Double-Loop Procedure With Simultaneous Maxima.

The double-loop procedure with simultaneous maxima, which is also called the mixed EGO method in Ref. [24], refers to the method that identifies the maxima simultaneously in the inner loop. In the mixed EGO method, an adaptive sampling approach is employed to reduce the number of training points required to build the surrogate model, $ĝmax(X)$ in the outer loop. In the inner loop, a surrogate model $ĝ(X, t)$ is built to identify the extreme values of all the training points simultaneously. Every time a new training point is added in the outer loop, the surrogate model $ĝ(X, t)$ is updated in the inner loop to identify the corresponding extreme value response of the new training point. The basic framework of the mixed EGO method is summarized as

• Outer loop: Build a surrogate model $ĝmax(X)$ using the adaptive sampling approach.

• Inner loop: Build a surrogate model $ĝ(X, t)$ to identify $ĝmax(x(i)), ∀i=1, 2, ⋯, nt$, simultaneously.

#### Drawback Analysis of the Double-Loop Procedure.

Even if both the independent EGO and mixed EGO methods are more accurate and efficient than the upcrossing rate methods, as discussed in Sec. 1, there are two main drawbacks (explained as below) for all the double-loop procedure methods.

1. (1)

Since the extreme values identified in the inner loop are used to build the extreme value response $ĝmax(X)$ in the outer loop, the accuracy of the identified extreme values will affect the accuracy of predictions in the outer loop at the untrained sample points, which will in turn affect the accuracy of time-dependent reliability analysis.

2. (2)

Current surrogate modeling method for time-dependent reliability analysis is mainly developed for problems with response function $G(t)=g(X, t)$. These methods can be extended to problems with stochastic processes by using the Karhunen–Loeve (KL) expansion [29,30]. Although the extension is easy and possible, the number of function (NOF) evaluations will increase significantly for problems with stochastic processes over a long time interval of interest. The reason is that there may be numerous peaks (extremes) in the stochastic process, which creates a large computation demand on the global optimization in the inner loop.

In Sec. 3, we propose a single-loop procedure for time-dependent reliability analysis, which completely removes the optimization inner loop in the double-loop procedure.

## Proposed SILK Method

### Basic Principle of Surrogate Modeling for Time-Dependent Reliability Analysis.

The basic idea of SILK is to build a single surrogate model $ĝ(X, t)$ to perform time-dependent reliability analysis, instead of using a double-loop procedure. The main principle of SILK is explained as below.

The time-dependent probability of failure given in Eq. (3) can be rewritten as follows based on the principle of MCS:
$pf(t0, te)=Pr{g(X, t)>0, ∃t∈[t0, te]}=∑i=1NIt(x(i))/N$
(5)
where $It(x(i))$ is the time-dependent failure indicator for given $X=x(i)$
$It(x(i))={1, if g(x(i), t)>0, ∃t∈[t0, te]0, otherwise, ∀i=1, 2, ⋯, N$
(6)

where “ $∀$ ” means “for all.”

After discretizing $[t0, te]$ into $Nt$ time instants, the time-dependent failure indicator can be rewritten as
$It(x(i))={1, if any I(g(x(i), t(j)))=1, ∀j=1, 2, ⋯, Nt0, otherwise$
(7)
in which $Nt$ is problem-dependent, typically the larger the better, and $I(g(x(i), t(j)))$ is given by
$I(g(x(i), t(j)))={1, if g(x(i), t(j))>00, otherwise$
(8)

The above equations imply that the accuracy of time-dependent reliability analysis is mainly affected by $It(x(i))$, which is actually a one-dimensional classification problem for given $X=x(i)$ (as shown in Eqs. (7) and (8)). As shown in Fig. 2, the accuracy of the one-dimensional classification problem is mainly affected by the accuracy near the crossing points [31,32]. Or in other words, we may not really need the extreme value as shown in Fig. 1.

Further analysis shows that Eq. (6) sometimes does not require high accuracy near the crossing points. This is different from the one-dimensional classification problem which requires accurate classification at each sample point. The time-dependent failure indicator function (Eq. (6)) can be accurately determined if there exists a failed point on the trajectory and the sign of the failed point is accurately classified no matter the signs of crossing points are accurately classified or not. Equation (6) can be further divided into the following two cases:

• Case 1: $g(x(i), t)≤0$ for all the time instant over $[t0, te]$. In this case, the sign of surrogate model prediction $ĝ(x(i), t(j))$ for all $t(j)$, $j=1, 2, ⋯, Nt$ need to be accurately classified to get an accurate estimate of $It(x(i))$.

• Case 2: There exists a time instant $t*∈[t0, te]$ such that $g(x(i), t*)>0$. In this case, if the sign of any $ĝ(x(i), t*)>0$ is correctly classified, Eq. (6) will be accurately estimated.

Based on the above analysis, the following function can be defined:
$It(x(i))={1, if ĝ(x(i), t(j))>0 and the sign of ĝ(x(i), t(j)) is correctly classified, ∃j=1, 2, ⋯, Nt0, if ĝ(x(i), t(j))≤0 and the sign of ĝ(x(i), t(j)) is correctly classified, ∀j=1, 2, ⋯, Nt$
(9)
For all the sample points $x(i)$, if the prediction from the surrogate model $ĝ(X, t)$ can satisfy the conditions given in Eq. (9), the value of $It(x(i))$ can be accurately determined. Otherwise, $ĝ(X, t)$ needs to be refined until one of the conditions presented in Eq. (9) is satisfied. Next, we will discuss how to determine whether the conditions given in Eq. (9) are satisfied or not, and how to refine $ĝ(X, t)$ if none of the conditions is met.

### Surrogate Modeling of ĝ(X, t)

#### Initial Surrogate Modeling of ĝ(X, t).

According to the basic principle of SILK as discussed in Sec. 3.1, an initial surrogate model $ĝ(X, t)$ needs to be constructed. Suppose $nin$ initial training points are generated (the samples of $X$ and $t$ are generated together using Latin Hypercube Sampling approach [33] or Hammersley sampling approach [34]), we then have the matrix of training points as
$xs=[x1(1)⋯xn(1)t(1)x1(2)⋱xn(2)t(2)⋮⋮⋱⋮x1(nin)⋯xn(nin)t(nin)]$
(10)

where $xi(j)$ is the jth training point of the ith random variable.

The response function $G(t)=g(X, t)$ is then evaluated at these training points. Based on the training points and the associated responses, a surrogate model $Ĝ(t)=ĝ(X, t)$ is built using the Kriging method reviewed in the appendix. For any untrained point $[x, t]$, the surrogate model prediction is
$Ĝ(t)∼N(ĝ(x, t), σĝ2(x, t))$
(11)

in which $N(⋅, ⋅)$ stands for normal distribution, and $ĝ(x, t)$ and $σĝ2(x, t)$ are the expected value and variance of the prediction, which are obtained from Eqs. (A4) and (A5), respectively.

After discretizing $[t0, te]$ into $Nt$ time instants, if $Ĝ(t)=ĝ(X, t)$ is well trained, $pf(t0, te)$ is obtained by
$pf(t0, te)=∑i=1NI(maxj=1, 2, ⋯, Nt{ĝ(x(i), t(j))})/N$
(12)

In Secs. 3.2.2 and 3.2.3, we will investigate how to determine whether the surrogate model is well trained (stopping criterion of training) and how to refine the surrogate model if it is not well trained.

#### Stopping Criterion for Training.

In order to determine whether the surrogate model is well trained or not, for any given $X=x(i)$, we first need to determine whether the conditions given in Eq. (9) are satisfied.

There are two types of conditions presented in Eq. (9). The first condition checks the sign, whether $ĝ(x(i), t(j))>0$ or $ĝ(x(i), t(j))≤0$. The second condition checks whether the sign of $ĝ(x(i), t(j))>0$ or $ĝ(x(i), t(j))≤0$ is accurately determined or not. In other words, the first condition performs the classification and the second condition checks the accuracy of the classification. The first condition can be easily checked based on the surrogate model prediction. For the second condition, the widely used U function defined in the adaptive Kriging Monte Carlo simulation (AK-MCS) method [32] for time-independent reliability analysis is employed.

If $ĝ(x(i), t(j))>0$, the probability of making a mistake on the sign of $g(x(i), t(j))$ is given by
$P1=Φ(0−|ĝ(x(i), t(j))|σĝ(x(i), t(j)))$
(13)

in which $Φ(⋅)$ is the cumulative density function of a standard normal variable.

If $ĝ(x(i), t(j))≤0$, the probability of making a mistake on the sign of $g(x(i), t(j))$ is given by
$P2=1−Φ(0+|ĝ(x(i), t(j))|σĝ(x(i), t(j)))=Φ(−|ĝ(x(i), t(j))|σĝ(x(i), t(j)))$
(14)
The above equations indicate that no matter $ĝ(x(i), t(j))>0$ or $ĝ(x(i), t(j))≤0$, the probability of making a mistake on the sign of $g(x(i), t(j))$ is given by [32]
$Perror=Φ(−U(x(i), t(j)))$
(15)
in which
$U(x(i), t(j))=|ĝ(x(i), t(j))|σĝ(x(i), t(j))$
(16)
It can be seen that the smaller the value of $U(x(i), t(j))$, the higher the probability of making a mistake on the sign of $g(x(i), t(j))$ (i.e., the higher the probability that the point is close to the crossing point). It is recommended that we can assume that the sign of a response is correctly determined based on the surrogate model prediction if $U(x(i), t(j))≥2$, which corresponds to a probability of making a mistake on the sign of 0.023 [32]. Based on this, Eq. (9) is rewritten as
$It(x(i))={1, if ĝ(x(i), t(j))>0 and U(x(i), t(j))≥2, ∃j=1, 2, ⋯, Nt0, if ĝ(x(i), t(j))≤0 and U(x(i), t(j))≥2, ∀j=1, 2, ⋯, Nt$
(17)
Alternatively, if we define $U′(x(i), t(j))=ĝ(x(i), t(j))/σĝ(x(i), t(j))$, Eq. (9) can be rewritten as
$It(x(i))={1, if U′(x(i), t(j))≥2, ∃j=1, 2, ⋯, Nt0, if U′(x(i), t(j))≤−2, ∀j=1, 2, ⋯, Nt$
(18)
In this paper, Eq. (17) is used. For any given $X=x(i)$, we can then determine whether there is a training point required on the trajectory $g(x(i), t)$, $t∈[t0, te]$ to refine the surrogate model $ĝ(X, t)$ based on another indicator defined as follows:
$Umin(x(i))={ue, if ĝ(x(i), t(j))>0 and U(x(i), t(j))≥2, ∃j=1, 2, ⋯, Ntminj=1, 2, ⋯, Nt{U(x(i), t(j))}, otherwise$
(19)
where $ue$ is any number so that $ue>2$.

If $Umin(x(i))>2$, it means that the trajectory $g(x(i), t)$, $t∈[t0, te]$ does not need to be trained. Otherwise, a new training point is needed to refine the trajectory. The $Umin(x(i))$ indicator is just for one trajectory under given $X=x(i)$. Next, we will discuss how to connect this indicator with the stopping criterion of the surrogate modeling.

In order to connect $Umin(x(i))$ with the overall surrogate model training, in this work, a recently developed stopping criterion for surrogate modeling-based time-independent reliability analysis [35] is extended to the time-dependent reliability analysis. In order to define the stopping criterion, we first partition the samples $x(i), i=1, 2, ⋯, N$ into two groups, namely, group-one samples $xg1MCS$ and group-two samples $xg2MCS$. The group-one samples $xg1MCS$ correspond to samples with $Umin(x(i))≥2$ and $xg2MCS$ are the rest of samples. Based on this partition, Eq. (5) is rewritten as
$pf(t0, te)=(Nf1+Nf2)/N$
(20)

where $Nf1=∑It(xg1MCS)$, $Nf2=∑It(xg2MCS)$, and $It(xg1MCS)$ and $It(xg2MCS)$ are obtained using Eq. (6) based on the mean prediction of the surrogate model $ĝ(X, t)$.

Since the group-one samples correspond to $Umin(x(i))≥2$, it can be assumed that $Nf1$ is accurate. Due to the uncertainty in the $It(xg2MCS)$, there is uncertainty in $Nf2$. Let the true value of $Nf2$ be $Nf2*$, the true time-dependent probability of failure is given by
$pf*(t0, te)=(Nf1+Nf2*)/N$
(21)
The percentage error of the current time-dependent failure probability estimate given in Eq. (20) is computed by
$εr=|pf(t0, te)−pf*(t0, te)|pf*(t0, te)×100%=|Nf2−Nf2*|Nf1+Nf2*×100%$
(22)
Let the number of samples in group two be $N2$. The true value of $Nf2$ (i.e., $Nf2*$) can be any number in the range $[0, N2]$. The maximum percentage error of the failure probability estimate given in Eq. (20) is estimated as
$εrmax=maxNf2*∈[0, N2]{|Nf2−Nf2*|Nf1+Nf2*×100%}$
(23)

The training of $ĝ(X, t)$ stops if $εrmax<5%$. The above procedure determines whether a surrogate model is well trained or not (i.e., stopping criterion). Next, we will investigate how to identify a new training point to refine the surrogate model if the surrogate model is not well trained.

#### Identification of New Training Point.

If the stopping criterion given in Eq. (23) is not satisfied, a potential new training point needs to be identified to refine the surrogate model $ĝ(X, t)$. Since $Umin(x(i))$ corresponds to the point with the highest probability of making a mistake on the sign of $g(x(i), t)$, we identify the new training point by minimizing the $Umin(x(i))$. Similar strategy has been implemented in surrogate model-based time-independent reliability analysis in the AK-MCS method [32].

For each trajectory, we first identify the associated time instant that corresponds to $Umin(x(i))$ as follows:
$tmin(i)=argminj=1, 2, ⋯, Nt{U(x(i), t(j))}$
(24)
The new training point is then identified as
$xnewt=[x(imin), tmin(imin)]$
(25)
where $imin$ is obtained by
$imin=argmini=1, 2, ⋯, N{Umin(x(i))}$
(26)

After the new training point is identified, the response function is evaluated and $ĝ(X, t)$ is trained again and the stopping criterion given in Eq. (23) is checked for all iterations until the stopping criterion is satisfied.

• (a)

Additional criterion for the selection of new training point

In addition to the above criterion for the selection of new training points, in this paper, we also introduce another criterion based on the correlation analysis between the potential new training points and the current training points.

Considering that the U function given in Eq. (16) is used to check the conditions given in Eq. (17) and the new training point is selected based on that, there may be clustering of the training points when the surrogate model is not well trained (i.e., $σĝ(x(i), t(j))$ are large and close to each other for different samples) [35]. The clustering of training points refers to the phenomenon that the selected training points are very close to each other. When this phenomenon happens, the correlation matrix used in the prediction of Kriging model will be ill-conditioned. In that case, some of the clustered training points will have negligible effect on the accuracy improvement of the surrogate model (i.e., some of the training points are not useful). An example of the clustering of training points is given in example 2 in Sec. 4. In order to avoid the clustering issue, we define a correlation condition for the refinement of the surrogate model.

Let the current training points be $xs$ (Eq. (10)). The correlation condition is defined as follows:
$max{ρ(xnewt, xs)}<0.95$
(27)

where $ρ(xnew, xs)$ is computed by substituting the normalized $xnewt$ and $xs$ into Eq. (A4). The above correlation constraint avoids the situation that the new training point is too close to current training points. Note that the above correlation condition is only applied to new training points to make sure the new training points are not clustered with current training points. An empirical range is suggested for the correlation threshold as [0.95, 1] based on our numerical examples. A high value of the threshold (i.e., 1) implies a loose correlation condition. On the other hand, a low value of threshold means a strict correlation condition.

Until now, we have discussed the basic idea for the refinement of $ĝ(X, t)$. Table 1 gives a detailed procedure for the identification of a new training point based on the conditions given in Eqs. (24)(27).

Note: In step 2 (c), $Unew$ is replaced with a larger number due to the application of correlation condition (Eq. (27)) to control the distance between current training points and new training points.

The above training of the surrogate model and stopping criterion are for problems without stochastic processes. The basic idea and main procedures are applicable to problems with and without stochastic processes. In Sec. 3.3, we will investigate the extension of the proposed method to general problems with both random variables and stochastic processes.

### Extension to Problems With Stochastic Processes.

As discussed in Sec. 2.1, the response function for problems with stochastic processes is given by $G(t)=g(X, Y(t), t)$. In both the double-loop procedure and SILK, stochastic processes $Y(t)$ first need to be represented as independent random variables using the KL expansion [29]. For a stochastic process, $Yi(t)$, the KL expansion is given by
$Yi(t)=μYi(t)+σYi(t)∑j=1neλj ξj fj(t)$
(28)

in which $μYi(t)$ and $σYi(t)$ are the mean and standard deviation of the stochastic process, $ξj$, $j=1, 2, ⋯, ne$ are independent random variables, $λj$ and $fj(t)$ are the eigenvalues and eigenvectors of the covariance function of $Yi(t)$, and $ne$ is the number of eigenvectors used to represent the stochastic process.

After the KL expansion, the response function for problems with stochastic processes becomes $G(t)=g(X, ξ, t)$, where $ξ$ is a vector of random variables used to represent the stochastic processes. Based on the transformation of the response function, if the double-loop procedure is applied to time-dependent reliability analysis, an extreme value surrogate model $ĝmax(X, ξ)$ needs to be constructed. Since the number of random variables in $ξ$ is usually very large (i.e., the terms in KL expansion will increase with the time duration of interest), the dimension of $ĝmax(X, ξ)$ is high. Constructing a high-dimensional surrogate model usually requires a large number of training points and is difficult for current surrogate modeling methods. In the proposed method, a surrogate model is directly built for $G(t)=g(X, Y(t), t)$.

Similar to the procedure for problems with response function $ĝ(X, t)$, $nin$ initial training points are first generated for $X$, $ξ$, and $t$. The training points of $X$, $ξ$, and $t$ are then transformed into training points of $X$, $Y$, and $t$ as follows:
$[x(1)ξ(1)t(1)⋮⋮⋮x(nin)ξ(nin)t(nin)]→Eq. (28)→ [x(1)y(1)t(1)⋮⋮⋮x(nin)y(nin)t(nin)]$
(29)

Based on the initial training points, an initial surrogate model $ĝ(X, Y, t)$ is built. In order to perform time-dependent reliability analysis based on $ĝ(X, Y, t)$, $N$ samples are first generated for $X$ and $ξ$ using MCS, and $[t0, te]$ is discretized into $Nt$ time instants. The MCS samples of $ξ$ are then converted into MCS samples of $Y$ based on Eq. (28) and the discretized time instants $t(j),$$j=1, 2, ⋯,Nt$. In order to determine the new training points to refine $ĝ(X, Y, t)$, the procedure presented in Table 1 is modified for problems with stochastic processes. Table 2 gives the procedure of identifying a new training point for $ĝ(X, Y, t)$.

Similar to the problems without stochastic processes, the stopping criterion in each iteration is computed using Eq. (22) and $Umin(x(i), ξ(i))$, $i=1, 2, ⋯, N$. Since in the proposed method, the surrogate model is constructed directly for $ĝ(X, Y, t)$ instead of $ĝmax(X, ξ)$, the dimensionality of $ĝ(X, Y, t)$ is much lower than that of $ĝmax(X, ξ)$. For example, suppose there are $nX$ random variables ($X$) and $nY$ stochastic processes ($Y(t)$), the dimension of $ĝ(X, Y, t)$ is $nX+nY+1$. If $nξ$ standard normal variables ($ξ$) are used in the KL expansion of each stochastic process, the dimension of $ĝmax(X, ξ)$ will be $nX+nYnξ+1$, where $nξ$ is problem dependent and will increase with time duration of interest. The value of $nξ$ is usually larger than 5 and can go up to a very large number (i.e., 1000). For example 3 presented in Sec. 4, the dimensions of $ĝ(X, Y, t)$ and $ĝmax(X, ξ)$ are 5 and 10, respectively. This property makes the proposed method applicable to problems with and without stochastic processes.

### Implementation Procedure.

We now summarize the implementation procedure of SILK. Table 3 gives the main steps of SILK by integrating Tables 1 and 2 into the implementation framework.

Note that the $Covpf$ criterion is chosen as 0.02 in this paper. Other threshold values, such as 0.05, can be used as well according to the requirement of the decision maker. This threshold is used to account for the statistical uncertainty in MCS.

## Numerical Examples

In this section, three numerical examples are used to demonstrate the proposed method. Examples 1 and 2 do not have stochastic process in the limit state function while example 3 has both random variables and stochastic process. The NOF evaluations and the percentage error of time-dependent failure probability estimate of the proposed method (i.e., SILK) are compared with the following four methods in the three examples.

• Rice: the upcrossing rate method based on the Rice's formula as reviewed in Sec. 2.1 and presented in Refs. [15] and [17]

• independent EGO: the double-loop procedure with independent EGO (Sec. 2.2.1) [20,21]

• mixed EGO: the double-loop procedure with mixed EGO (Sec. 2.2.2) [24]

• MCS: the brute force MCS performed on the original limit-state function

The percentage error of the time-dependent failure probability estimate is computed by
$ε%=|pf(t0, te)−pfMCS(t0, te)|/pfMCS(t0, te)×100%$
(30)

where $pfMCS(t0, te)$ is the time-dependent failure probability estimate obtained from MCS.

### A Mathematical Example.

A mathematical example given in Eq. (31), which is adopted from Ref. [24], is used as our first example
$g(X, t)=[sin(2.5X)cos (t+0.4)2]/(X2+4)$
(31)

where $X$ is a random variable following a normal distribution $X∼N(10, 0.52)$.

The time-dependent probability of failure over time interval $[1, 2.5]$ is to be estimated and is defined as
$pf(t0, ts)=Pr{g(X, τ)>0.014, ∃τ∈[1, 2.5]}$
(32)

We first solve this example using SILK. The time interval $[1, 2.5]$ is discretized into 100 time instants. Ten initial training points are generated for $X$ and $t$ using the Hammersley sampling approach [34]. The initial training points of $X$ are generated in the interval $[7.5, 12.5]$. A Kriging model with a constant trend function is used to build the initial surrogate model. The surrogate model $ĝ(X, t)$ is then refined using SILK. Table 4 gives the results comparison between SILK, Rice formula-based upcrossing rate method, double-loop procedure with independent EGO, double-loop procedure with mixed EGO, and MCS. Note that the NOF of the proposed method is not an integer because SILK is run for 20 times and the average results are reported to reduce the uncertainty in the MCS samples. The results of Rice, independent EGO, and mixed EGO are taken from Ref. [24]. It shows that the proposed SILK method is much more efficient than the other methods. However, the accuracy of SILK is a little bit lower than the mixed EGO method even though it satisfies the accuracy requirement defined in Table 3 ($εrmax≤0.05$). Further analysis shows that SILK requires 22 NOF to achieve the same accuracy level as mixed EGO.

### A Function Generator Mechanism.

A function generator mechanism as shown in Fig. 3 is used as our second example. This example is taken from Ref. [15].

The time-dependent probability of failure is given by [15]
$pf(θ0, θs)=Pr{ψ(X, θ)−ψd(θ)>θeπ/180, ∃θ∈[97, 217]}$
(33)
where $X=[L1, L2, L3, L4]$, $θe$ is failure threshold, $ψd(θ)=60π/180+60π/180 sin(0.75(θ−97)π/180)$, and $ψ(X, θ)$ is given by
$ψ(X, θ)=2arctan((−K1±K12+K22−K32)/(K3−K2))$
(34)

in which $K1=−2L2L4 sin(θπ/180)$, $K2=2L4(L1−L2 cos(θπ/180))$, and $K3=L12+L22+L42−L32−2L1L2 cos(θπ/180)$.

Table 5 gives the random variables of example 2.

We then perform time-dependent reliability analysis for this example using SILK, rice, independent EGO, mixed EGO, and MCS. Table 6 gives the results comparison of different methods with two different failure thresholds (i.e., $θe=0.9 deg$ and $θe=0.95 deg$). Similar to example 1, the SILK method is run for 20 times and the average results are reported. The results indicate that SILK is more efficient and accurate than the other methods.

• (a)

Parameter study

We also investigated the effects of the correlation threshold given in Eq. (27) on the results of time-dependent reliability analysis. Figure 4 shows the training points identified by SILK with and without the correlation condition. It shows that the correlation condition successfully avoids the clustering of training points in SILK.

We also performed SILK with different correlation threshold to study the effect of correlation condition on the accuracy and efficiency of time-dependent reliability analysis. Table 7 gives the NOF and percentage error of SILK (average results of 20 runs) with different values of correlation threshold for the $θe=0.9 deg$ case. It shows that by controlling the correlation between new training point and current training points, the efficiency can be improved while the effect on the accuracy is minor if the threshold is not set too low. Since the correlation is related to the distance between the training points in Kriging surrogate model, it also means controlling the distance between new and current training points helps to avoid the clustering issue.

In addition to the correlation threshold, we studied the effects of the number of initial training points and the number of discretization points of the time interval on the accuracy and efficiency of time-dependent reliability analysis. Tables 8 and 9 give the average results of 20 runs of SILK with different values of $nin$ and $Nt$. It shows that the effect of $nin$ and $Nt$ on the results of time-dependent reliability analysis is very small. However, very small values of $nin$ and $Nt$ may result in large error of reliability analysis. A recommended value of $nin$ is given in Refs. [31] and [32] as $nin=(nr+1)(nr+2)/2$, where $nr$ is the number of random variables. For the value of $Nt$, the results indicate that the larger the better.

A corroded beam (Fig. 5) subjected to stochastic load is adopted from Ref. [24] as our third example.

The time-dependent probability of failure is given by
$pf=Pr{g(X, Y(t), t)>0, ∃t∈[0, 35]}$
(35)
where
$g(X, Y(t), t)=(F(t)L/4+ρsta0b0L2/8)−(a0−2kt)(b0−2kt)2σu/4$
(36)
in which $ρst=7.85×104 N$, $k=5×10−5 m/year$, $L=5 m$, and $F(t)$ is a stochastic process modeled by [24]
$F(t)=6500+∑i=17ξi(∑j=17(aij sin(bijt+cij)))$
(37)

where $ξi$, $i=1, 2, ⋯, 7$ are independent random variables, $aij$, $bij$, $cij$, ∀i, j = 1,2,…,7 are coeffiecients of the sine wave basis functions. Details of Eq. (37) can be found in Ref. [24].

Table 10 gives the random variables of this example and Table 11 gives the results comparison of different methods. Similar conclusions can be obtained as examples 1 and 2. Thus SILK achieves efficiency and accuracy in both problems with and without stochastic processes.

## Conclusion

For problems with nonlinear response, the MPP-based time-dependent reliability analysis methods may have large error in the reliability estimate due to linearization of response function at the MPP. In that situation, a surrogate model-based method is more promising. Current surrogate model methods for time-dependent reliability analysis implement a double-loop procedure, which is computationally expensive and unaffordable for problems with stochastic processes over a long time period.

This paper develops a SILK surrogate model method for time-dependent reliability analysis with and without stochastic processes. The proposed method generates training points and builds the surrogate model for random variables, stochastic processes, and time at the same level instead of at two levels. The global optimization in current methods is completely removed in the proposed method. The single-loop surrogate model is refined adaptively based on criteria developed according to the properties of the time-dependent problem. Three numerical examples demonstrated the effectiveness of the proposed method.

There are several advantages for the removal of the global optimization loop (inner loop) in the time-dependent reliability analysis. First, in the independent EGO and mixed EGO methods, the NOF evaluations will increase with the number of initial samples of random variables in the outer loop. In the SILK method, a single surrogate model is trained and no global optimization is required. The NOF of SILK method therefore will not increase with the number of initial random variable samples. Second, the removal of global optimization enables the surrogate model-based time-dependent reliability analysis method to be applied to problems with stochastic processes and over a long time period.

Since the proposed method makes use of the Kriging surrogate modeling method, it also has the limitations of the Kriging surrogate model. When the Kriging surrogate model cannot effectively model the time-dependent problem or the limit-state functions are highly nonlinear, it will result in the ineffectiveness of the proposed method. The basic idea of the proposed method is general and applicable with different types of surrogate models. Extension of the proposed method to other kinds of surrogate models for time-dependent reliability analysis may be investigated in future.

## Acknowledgment

The research reported in this paper was supported by the Air Force Office of Scientific Research (Grant No. FA9550-15-1-0018, Technical Monitor: Dr. David Stargel). The support is gratefully acknowledged.

### Appendix: Kriging Surrogate Model

For an unknown function $g(x)$ with inputs of $x$, the Kriging model is given by [30,36]
$ĝ(x)=h(x)Tυ+ε(x)$
(A1)

where $υ=[υ1, υ2, ⋯, υp]T$ is a vector of unknown coefficients, $h(x)=[h1(x), h2(x), ⋯, hp(x)]T$ is a vector of regression functions, $h(x)Tυ$ is the trend of prediction, and $ε(x)$ is assumed to be a Gaussian process with zero mean and covariance $Cov[ε(xi), ε(xj)]$.

The covariance between two points $xi$ and $xj$ is given by
$Cov[ε(xi), ε(xj)]=σε2R(xi, xj)$
(A2)

in which $σε2$ is the process variance and $R(⋅, ⋅)$ is the correlation function.

With $ns$ training points, $[xi, g(xi)]i=1, 2, ⋯, ns$, the coefficients $υ$ are obtained by
$υ=(HTR−1H)−1HTR−1g$
(A3)

Where $R$ is the correlation matrix with element, $R(xi, xj)$, $i, j=1, 2, ⋯, ns$, $H=[h(x1)T, h(x2)T, ⋯, h(xns)T]T$, and $g=[g(x1), g(x2), ⋯, g(xns)]T$.

For a new point $x$, the expected value of the prediction is given by
$ĝ(x)=h(x)Tυ+r(x)TR−1(g−Hυ)$
(A4)

where $r(x)=[R(x, x1), R(x, x2), ⋯, R(x, xns)]$

The mean square error of the prediction is given by [37]
$MSE(x)=σε2{1−r(x)TR−1r(x)+[HTR−1r(x)−h(x)]T(HTR−1H)−1[HTR−1r(x)−h(x)]$
(A5)

## References

References
1.
Stewart
,
M. G.
,
2001
, “
Reliability-Based Assessment of Ageing Bridges Using Risk Ranking and Life Cycle Cost Decision Analyses
,”
Reliab. Eng. Syst. Saf.
,
74
(
3
), pp.
263
273
.
2.
Hu
,
Z.
, and
,
S.
,
2015
, “
Time-Dependent System Reliability Analysis Using Random Field Discretization
,”
ASME J. Mech. Des.
,
137
(
10
), p.
101404
.
3.
Hu
,
Z.
, and
Du
,
X.
,
2014
, “
Lifetime Cost Optimization With Time-Dependent Reliability
,”
Eng. Optim.
,
46
(
10
), pp.
1389
1410
.
4.
Hu
,
Z.
, and
,
S.
,
2015
, “
Accelerated Life Testing (ALT) Design Based on Computational Reliability Analysis
,”
Qual. Reliab. Eng. Int.
, (online).
5.
Youn
,
B. D.
,
Hu
,
C.
, and
Wang
,
P.
,
2011
, “
Resilience-Driven System Design of Complex Engineered Systems
,”
ASME J. Mech. Des.
,
133
(
10
), p.
101011
.
6.
Andrieu-Renaud
,
C.
,
Sudret
,
B.
, and
Lemaire
,
M.
,
2004
, “
The PHI2 Method: A Way to Compute Time-Variant Reliability
,”
Reliab. Eng. Syst. Saf.
,
84
(
1
), pp.
75
86
.
7.
Singh
,
A.
, and
Mourelatos
,
Z. P.
,
2010
, “
On the Time-Dependent Reliability of Non-Monotonic, Non-Repairable Systems
,”
SAE
Paper No. 2010-01-0696.
8.
Singh
,
A.
,
Mourelatos
,
Z.
, and
Nikolaidis
,
E.
,
2011
, “
Time-Dependent Reliability of Random Dynamic Systems Using Time-Series Modeling and Importance Sampling
,”
SAE Int. J. Mater. Man.
,
4
(
1
), pp.
929
946
.
9.
Wang
,
Z.
,
Mourelatos
,
Z. P.
,
Li
,
J.
,
Baseski
,
I.
, and
Singh
,
A.
,
2014
, “
Time-Dependent Reliability of Dynamic Systems Using Subset Simulation With Splitting Over a Series of Correlated Time Intervals
,”
ASME J. Mech. Des.
,
136
(
6
), p.
061008
.
10.
Wang
,
Z.
,
Mourelatos
,
Z. P.
,
Li
,
J.
,
Singh
,
A.
, and
Baseski
,
I.
,
2014
, “
Time-Dependent Reliability of Dynamic Systems Using Subset Simulation With Splitting Over a Series of Correlated Time Intervals
,”
ASME J. Mech. Des.
,
136
(
6
), p.
061008
.
11.
Mourelatos
,
Z. P.
,
Majcher
,
M.
,
Pandey
,
V.
, and
Baseski
,
I.
,
2015
, “
Time-Dependent Reliability Analysis Using the Total Probability Theorem
,”
ASME J. Mech. Des.
,
137
(
3
), p.
031405
.
12.
Hagen
,
Ø.
, and
Tvedt
,
L.
,
1991
, “
Vector Process Out-Crossing as Parallel System Sensitivity Measure
,”
J. Eng. Mech.
,
117
(
10
), pp.
2201
2220
.
13.
Hagen
,
O.
, and
Tvedt
,
L.
,
1992
, “
Parallel System Approach for Vector Out-Crossing
,”
ASME J. Offshore Mech. Arct. Eng.
,
114
(
2
), pp.
122
128
.
14.
Chen
,
J.-B.
, and
Li
,
J.
,
2005
, “
Dynamic Response and Reliability Analysis of Non-Linear Stochastic Structures
,”
Probab. Eng. Mech.
,
20
(
1
), pp.
33
44
.
15.
Zhang
,
J.
, and
Du
,
X.
,
2011
, “
Time-Dependent Reliability Analysis for Function Generator Mechanisms
,”
ASME J. Mech. Des.
,
133
(
3
), p.
031005
.
16.
Du
,
X.
,
2014
, “
Time-Dependent Mechanism Reliability Analysis With Envelope Functions and First-Order Approximation
,”
ASME J. Mech. Des.
,
136
(
8
), p.
081010
.
17.
Hu
,
Z.
, and
Du
,
X.
,
2012
, “
Reliability Analysis for Hydrokinetic Turbine Blades
,”
Renewable Energy
,
48
(
1
), pp.
251
262
.
18.
Hu
,
Z.
, and
Du
,
X.
,
2013
, “
Time-Dependent Reliability Analysis With Joint Upcrossing Rates
,”
Struct. Multidiscip. Optim.
,
48
(
5
), pp.
893
907
.
19.
Hu
,
Z.
, and
Du
,
X.
,
2015
, “
First Order Reliability Method for Time-Variant Problems Using Series Expansions
,”
Struct. Multidiscip. Optim.
,
51
(
1
), pp.
1
21
.
20.
Wang
,
Z.
, and
Wang
,
P.
,
2012
, “
A Nested Extreme Response Surface Approach for Time-Dependent Reliability-Based Design Optimization
,”
ASME J. Mech. Des.
,
134
(
12
), p.
121007
.
21.
Wang
,
Z.
, and
Wang
,
P.
,
2013
, “
A New Approach for Reliability Analysis With Time-Variant Performance Characteristics
,”
Reliab. Eng. Syst. Saf.
,
115
(
1
), pp.
70
81
.
22.
Jiang
,
C.
,
Huang
,
X.
,
Han
,
X.
, and
Zhang
,
D.
,
2014
, “
A Time-Variant Reliability Analysis Method Based on Stochastic Process Discretization
,”
ASME J. Mech. Des.
,
136
(
9
), p.
091009
.
23.
Jones
,
D. R.
,
Schonlau
,
M.
, and
Welch
,
W. J.
,
1998
, “
Efficient Global Optimization of Expensive Black-Box Functions
,”
J. Global Optim.
,
13
(
4
), pp.
455
492
.
24.
Hu
,
Z.
, and
Du
,
X.
,
2015
, “
Mixed Efficient Global Optimization for Time-Dependent Reliability Analysis
,”
ASME J. Mech. Des.
,
137
(
5
), p.
051401
.
25.
Wang
,
Z.
, and
Wang
,
P.
,
2015
, “
A Double-Loop Adaptive Sampling Approach for Sensitivity-Free Dynamic Reliability Analysis
,”
Reliab. Eng. Syst. Saf.
,
142
(
1
), pp.
346
356
.
26.
Schoefs
,
F.
,
El Soueidy
,
C.-P.
, and
Hawchar
,
L.
,
2015
, “
Time-Variant Reliability Analysis Using Polynomial Chaos Expansion
,”
ICASP12
12th International Conference on Applications of Statistics and Probability in Civil Engineering
27.
Wang
,
Y.
,
Zeng
,
S.
, and
Guo
,
J.
,
2013
, “
Time-Dependent Reliability-Based Design Optimization Utilizing Nonintrusive Polynomial Chaos
,”
J. Appl. Math.
,
2013
, p.
513261
.
28.
Rice
,
S. O.
,
1944
, “
Mathematical Analysis of Random Noise
,”
Bell Syst. Tech. J.
,
23
(
3
), pp.
282
332
.
29.
Huang
,
S.
,
,
S.
, and
Rebba
,
R.
,
2007
, “
Collocation-Based Stochastic Finite Element Analysis for Random Field Problems
,”
Probab. Eng. Mech.
,
22
(
2
), pp.
194
205
.
30.
Li
,
C.
, and
,
S.
,
2016
, “
Relative Contributions of Aleatory and Epistemic Uncertainty Sources in Time Series Prediction
,”
Int. J. Fatigue
,
82
(
Part 3
), pp.
474
486
.
31.
Bichon
,
B. J.
,
Eldred
,
M. S.
,
Swiler
,
L. P.
,
,
S.
, and
McFarland
,
J. M.
,
2008
, “
Efficient Global Reliability Analysis for Nonlinear Implicit Performance Functions
,”
AIAA J.
,
46
(
10
), pp.
2459
2468
.
32.
Echard
,
B.
,
Gayton
,
N.
, and
Lemaire
,
M.
,
2011
, “
AK-MCS: An Active Learning Reliability Method Combining Kriging and Monte Carlo Simulation
,”
Struct. Saf.
,
33
(
2
), pp.
145
154
.
33.
Stein
,
M.
,
1987
, “
Large Sample Properties of Simulations Using Latin Hypercube Sampling
,”
Technometrics
,
29
(
2
), pp.
143
151
.
34.
Hammersley
,
A.
,
Svensson
,
S.
,
Hanfland
,
M.
,
Fitch
,
A.
, and
Hausermann
,
D.
,
1996
, “
Two-Dimensional Detector Software: From Real Detector to Idealised Image or Two-Theta Scan
,”
Int. J. High Pressure Res., 14(4–6), pp.
235
248
.
35.
Hu
,
Z.
, and
,
S.
,
2016
, “
Global Sensitivity Analysis Enhanced Surrogate (GSAS) Modeling for Reliability Analysis
,”
Struct. Multidiscip. Optim.
,
53
(
3
), pp.
501
521
.
36.
Rasmussen
,
C. E.
,
2006
,
Gaussian Processes for Machine Learning
,
The MIT Press
,
Cambridge, MA.
37.
Lophaven
,
S. N.
,
Nielsen
,
H. B.
, and
Søndergaard
,
J.
,
2002
, “
DACE-A Matlab Kriging Toolbox, Version 2.0
,” Technical University of Denmark, Technical Report No. IMM-TR-2002-12.