A Time-Variant Reliability Analysis Method Based on the Stochastic Process Discretization under Random and Interval Variables

: In practical engineering, it is a cost-consuming problem to consider the time-variant reliability of both random variables and interval variables, which usually requires a lot of calculation. Therefore, a time-variant reliability analysis approach with hybrid uncertain variables is proposed in this paper. In the design period, the stochastic process is discretized into random variables. Simultaneously, the original random variables and the discrete random variables are converted into independent normal variables, and the interval variables are changed into standard variables. Then it is transformed into a hybrid reliability problem of static series system. At different times, the limited state functions are linearized at the most probable point (MPP) and at the most unfavorable point (MUP). The transformed static system reliability problem with hybrid uncertain variables can be solved effectively by introducing random variables. To solve the double-loop nested optimization in the hybrid reliability calculation, an effective iterative method is proposed. Two numerical examples and an engineering example demonstrate the validity of the present approach.


Introduction
Due to structural material performance degradation, changing working environment, time-variant load effects, etc., the reliability of the structure exhibits time-variant properties [1][2][3][4][5][6]. Over the past decades, different time-variant reliability analysis methods have been developed including first-passage approaches, numerical simulation approaches, extreme value approaches and quasi-static approaches. The first-passage methods based on out-crossing events [7] have been developed. The representative methods are the PHI2 method [8], improved PHI2 method [9] and joint out-crossing rate method [10,11]. Although many first-passage methods [10,12,13] aimed at accuracy and efficiency have been developed for two decades, it is generally hard to achieve the first-time out-crossing rate due to complicated mathematical characteristics. Different numerical simulation methods have been developed including Monte Carlo methods [14] and their improved versions, namely important sampling methods [15,16] and subset methods [17][18][19]. Although the numerical simulation method is accurate, it demands huge computational cost. The extremum method [20][21][22][23] mainly focuses on the worst scenario over the time scale, in which surrogate model or probability distribution are used to describe the uncertainty of response extremum for time-variant problem. In some cases, the extreme value distribution may be subject to multi-modal or highly non-linear distribution, and the realistic application of the extreme value method is hindered. Recently, the quasi-static methods have been developed to improve efficiency, such as the stochastic process discretization approach [24] and the envelope method [25]. These methods translate the estimation of time-dependent failure

Problem of General Time-Variant Reliability Model
Time-variant reliability of structures refers to the possibility of completing predetermined functions within a specified time and under specified conditions for structures subjected to dynamic uncertainties. For a specific structure whose limit-state function is g(X(t), Y, t), the probability of structural reliability within the time period [0, T] can be defined as: P s (T) = Prob{∀t ∈ [0, T]|g(X(t), Y, t) > 0 } (1) where Prob stands for probability operation, t is the time, T is the design lifetime, X(t) = [X 1 (t), X 2 (t), . . . , X n (t)] is the n-dimensional stochastic process vector, and Y = [Y 1 , Y 2 , . . . , Y h ] is the h-dimensional random vector. The out-crossing rate method is the most commonly used method to solve time-variant reliability problems. But sometimes it is difficult to calculate the crossing rate accurately and effectively. The time-variant reliability analysis method based on stochastic process discretization avoids the calculation of out-crossing rate and simplifies the solving process of timevariant reliability. According to the method in literature [26,27], the design lifetime T Symmetry 2021, 13, 568 3 of 15 is discretized into m equal periods with each time step size ∇t = T m . According to the reliability calculation theory of series system, Equation (1) can be changed into: . . , m, (X i,1 , X i,2 , . . . , X i,n ), i= 1, 2, . . . , m, and (X T 1 , X T 2 , . . . , X T m ) is m × n dimensional random variable obtained by the discrete process X(t). Obviously, Equation (2) is a time-variant reliability problem with only random variables. If all random variables can provide accurate probability distribution, the reliability calculation can be carried out by using the conventional probability model.

Problem of General Time-Variant Reliability Model with Random and Interval Variables
Due to the lack of accurate probability distribution of uncertain parameters in the time-variant reliability problem, the interval ranges of uncertain parameters are relatively easy to be given. Therefore, a time-variant reliability analysis model for random variables and interval variables is constructed. The limit-state function of hybrid uncertain variables can be formulated as where q = [q 1 , q 2 , . . . , q l ] is the l-dimensional interval vector. The corresponding interval vector can be defined as where superscript I, L and R represent interval, upper and lower ranges of interval, respectively.

Normalization of Random Variables
In the process of time-variant reliability analysis, the Nataf transformation [37] is first adopted to transform the random variable (X, Y) at time t i into the standard normal space: where Nataf (·) represents the Nataf transformation. ρ X and ρ Y are the correlation coefficient matrices of the random variable X i and Y, respectively. ρ U and ρ V are the correlation coefficient matrices of the standard normal variable U and V, respectively.
Let the covariance of the random variable U= (U 1 , where i = 1, 2, . . . , m, j = 1, 2, . . . , m. The diagonal element is the variance σ 2 U i of the variable U i , the non-diagonal element is the covariance C U = Cov(U i , U j ) of the variable U i and U j , and the C U = Cov(U i , U j ) is a symmetric positive definite matrix of order n, which can be expressed as follows: The above m-order symmetric positive definite matrix C U has m linearly independent eigenvectors α 1 , α 2 , . . . , α m . Let the linear transformation matrix be A = [α 1 , α 2 , . . . , α m ], then A −1 C U A = Λ. Λ denotes an N-order diagonal matrix whose diagonal element is the eigenvalue λ i , i = 1, . . . , m of the orthogonal matrix A. Relevant normal random variable U can be changed into independent normal random variable P by the method of orthogonal transformation [37].
The random vector U is changed into a linearly independent random vector P. It is known from the matrix theory that A −1 = A T , and Equation (9) can be written as Superscript T indicates the transpose of the matrix. Through the orthogonal transformation of Equation (10), the covariance matrix C P can be transformed into the following form Similarly, the correlated random variable V can be transformed into the uncorrelated random variable Q by orthogonal transformation.
where B is a linear transformation matrix.

Normalization of Interval Variables
The interval variables q can be expressed as where superscript C and W are the midpoint and radius of the interval, respectively.
The interval is normalized The uncertainty domain . . , l} becomes a standard multi-dimensional cube, and the center coincides with the origin.

Solution of Time-Invariant Reliability Problem
After the normalized transformation of the above hybrid uncertain variables, at the time t i , the limit-state function g(X i , Y, q, t i ) in the original uncertainty variable space is mapped to its standardized form G(P i , Q, δ, t i ).
Because of involving random variables and interval variables, the limit-state equation of a structure G(P i , Q, δ, t i ) = 0 forms a critical region of a strip. Thus, the entire standard space is divided into three parts: the safety zone, the critical zone and the failure zone. Schematic diagram for two random variables (P, Q) in two-dimensional standard space is shown in Figure 1. Therefore, we define the reliability of the structure in the mixed model: the probability that the structure can perform at least the predetermined function of the Symmetry 2021, 13, 568 5 of 15 structure for any possible implementation of the interval parameter at a certain time. The mathematical formula is expressed as where subscript h represents the mixed model, G = min δ (P i , Q, δ, t i ) and G= 0 represents the critical failure of the structure.
mapped to its standardized form ( ) Because of involving random variables and interval variables, the limit-state equation of a structure ( ) forms a critical region of a strip. Thus, the entire standard space is divided into three parts: the safety zone, the critical zone and the failure zone. Schematic diagram for two random variables (P, Q) in two-dimensional standard space is shown in Figure 1. Therefore, we define the reliability of the structure in the mixed model: the probability that the structure can perform at least the predetermined function of the structure for any possible implementation of the interval parameter at a certain time. The mathematical formula is expressed as where subscript h represents the mixed model, the critical failure of the structure. As an extension of the probabilistic reliability index, the mixed reliability index h β can be defined as the shortest distance from the origin to the most probable failure sur- As an extension of the probabilistic reliability index, the mixed reliability index β h can be defined as the shortest distance from the origin to the most probable failure surface G(P, Q) = 0 in the standard space. Mathematically, it is the solution to the following optimization problem: where δ * is the most unfavorable point (MUP), which can be obtained through the following optimization: min Equations (17) and (18) that define the hybrid reliability index are nested optimization problems. The most probable point (MPP) P i , Q is obtained in the outer loop optimization, while the most unfavorable point (MUP) δ * is achieved in the inner loop optimization. The nested optimization problem in Equations (17) and (18) are solved by decoupling method. In turn, the inner and outer optimization problems can be solved according to the idea of the two-layer method.
In each iteration, MPP can be expressed as following: Symmetry 2021, 13, 568 6 of 15 In the above algorithm, δ k+1 is a fixed value. Once MPP points are obtained, interval variables δ k+1 can be obtained through the following optimization.
In each iteration of MPP search, standardized interval variable δ is fixed. MPP algorithm updates the random variable, and optimizes the minimum limit-state function when random variable (P, Q) is fixed. The MUP δ * can be obtained by Equation (20).
For the convenience of analysis, the linear expansion of the static limit-state equations is carried out at the MPP and MUP δ * k .
Equation (21) can be converted to the following form: A new random vector ξ = (ξ 1 , . . . , ξ m ) T is introduced, which can be expressed as: Since P i,j represents a standard normal random variable, the random vector ξ is the m-dimensional normal distribution. According to the properties of mean and covariance matrix of random vectors, we can define mean vector µ and covariance matrix C as: where µ P i,j is the mean of µ P i,j , and (i − 1)n + j and (v − 1)n + k in parentheses represent the numbers of rows and columns in the matrix, respectively. In order to understand the construction of C P , this paper gives the specific form of the matrix in four simple cases (See reference [24] for details). Integrating the m-dimensional normal distribution function, we change Equation (22) into the following expression: φ m is the probability distribution function of m-dimensional normal distribution. f Q (Q) denotes the joint probability density function of the random vector Q. Since Equation (26) is a multi-dimensional integration problem, the numerical solution requires a large amount of calculation. Therefore, this paper applies a new random variable E and converts the formula Equation (26) into the following form.
where Ω denotes the integral region, which can be expressed as: F E (e) and f E (e) denote the probability distribution function and the probability density function of the random variable E, respectively; F −1 E (·) is the inverse function of F E (·). Equation (28) is a static reliability analysis model with new limit-state equations.
From Equation (29), we can see that the ultimate expression has nothing to do with the specific distribution of E. For the convenience of calculation, we can choose a normal distribution.
Detailed steps of obtaining MPP and MUP are as follows: (1) Input initial start points P (3) Interval analysis is carried out and the MUP δ * can be obtained through the optimization problem of Equation (20). The value of the random variable is achieved from the above probability analysis.

Numerical Examples
The accuracy and effectiveness of this method are verified by Monte Carlo method. The calculation steps of Monte Carlo method are as follows: (1) In outer loop, random variable sample Y is generated, while using the EOLE model [38] to generate stochastic process samples X(t). (2) In inner loop, at each sample, the optimal solution δ* resulting in the minimum response of the limit-state function is obtained by the optimization method.
(3) The obtained sample and δ* are used to calculate value of the limit-state function. If min G(X(t),Y, δ * , t) ≤ 0, t ∈ [0, T], the number of failures n f = n f + 1. (4) Repeat the above steps until the total sample size n s is reached. Thus, the cumulative failure probability P f = n f n s is obtained.

Steel Beam
The simply supported steel beam structure [23,24,39] has a span of L = 5m and a rectangular section of b 0 × h 0 as shown in Figure 2. It bears uniform load P, and its middle point is affected by a concentrated dynamic load Q(t). The uniform load P can be expressed as P = ρ st b 0 h 0 , where ρ st = 78,500 N/m 3 is the steel force density. Assuming that the corroded part will lose mechanical strength in time, the change rule of the remaining section area A(t) can be expressed as: 03mm/year. Thus, on the basis of the strength failure criterion, the limited state function is established as follows: where σ represents the material yield stress. The initial dimensions of the beam section b 0 and h 0 and the material yield stress σ are regarded as random variables, and the dynamic load F(t) is treated as a stationary Gaussian process. Tables 1 and 2 respectively list the distribution of specific random parameters and the range of interval parameters.

Autocorrelation Coefficient
Yield stress σ (MPa)   Figure 3 shows that the time-variant reliability index of the steel beam structure continuously decayed as the design period increased. The algorithm proposed calls the limit-state equations 540 times, 1055 times and 2170 times respectively, and obtains stable calculation results. The Monte Carlo sample size is set to 100,000. The limit-state equation is called 74,441,300 times by using the Monte Carlo method. The proposed approach shows high calculation efficiency. The calculation error can be expressed as the difference between the current method and the Monte Carlo method and divided by the Monte Carlo method. As can be seen from Table 3, with the decrease of the time step, the accuracy of calculation becomes higher, and the reliability index gradually approaches the accurate value. When Δt = 2 years, the maximum error is 24.5%. When Δt =1 year, the maximum error is 15.6%. When Δt =0.5 year, the maximum error is 6.67%. Thus, the calculation result is accurate and meets the engineering needs.    Figure 3 shows that the time-variant reliability index of the steel beam structure continuously decayed as the design period increased. The algorithm proposed calls the limit-state equations 540 times, 1055 times and 2170 times respectively, and obtains stable calculation results. The Monte Carlo sample size is set to 100,000. The limit-state equation is called 74,441,300 times by using the Monte Carlo method. The proposed approach shows high calculation efficiency. The calculation error can be expressed as the difference between the current method and the Monte Carlo method and divided by the Monte Carlo method. As can be seen from Table 3, with the decrease of the time step, the accuracy of calculation becomes higher, and the reliability index gradually approaches the accurate value. When ∆t = 2 years, the maximum error is 24.5%. When ∆t =1 year, the maximum error is 15.6%. When ∆t =0.5 year, the maximum error is 6.67%. Thus, the calculation result is accurate and meets the engineering needs.

A Cantilever Tube Structure
The structure of a tubular cantilever beam [24] is shown in Figure 4, which is subjected to external forces Q(t), F, P and torque U(t). F and P are permanent loads; Q(t) and U(t) are dynamic loads. According to [24], material strength R decays with time due to material degradation, and its decay rule is assumed to be ( )

A Cantilever Tube Structure
The structure of a tubular cantilever beam [24] is shown in Figure 4, which is subjected to external forces Q(t), F, P and torque U(t). F and P are permanent loads; Q(t) and U(t) are dynamic loads. According to [24], material strength R decays with time due to material degradation, and its decay rule is assumed to be R(t) = R 0 (1 − 0.01t), where R 0 is the initial yield strength.
In the formula, σ max (t) is calculated as follows: M(t) = FL 1 cos θ 1 + Q(t)L 2 cos θ 2 (35) In this problem, the initial yield strength R 0 , the size parameters h and d, and the loads F and P are seen as random variables. The dynamic load Q(t) and U(t) are treated as Gaussian random processes. Table 4 lists the parameters' distributions. As shown in Table 5, the L 1 and L 2 are treated as interval variables.
In this problem, the initial yield strength R0, the size parameters h and d, and the loads F and P are seen as random variables. The dynamic load Q(t) and U(t) are treated as Gaussian random processes. Table 4 lists the parameters' distributions. As shown in Table 5, the L1 and L2 are treated as interval variables.  Figure 4. A cantilever tube [24].  Three cases are considered, that is, ∆t = 1 year, ∆t = 0.5 year and ∆t = 1/3 year. We adopt this method and Monte Carlo to solve the reliability. The random sample number of Monte Carlo is set to 100,000. Table 6 and Figure 5 show the calculation results. It can be observed that when ∆t = 1 year, ∆t = 0.5 year and ∆t = 1/3 year, the maximum errors are 25.4%, 9.6% and 5.3% respectively. It can be found that with the decrease of the time step, the accuracy of reliability solution is gradually improved. In terms of efficiency, Monte Carlo needs to call the limit-state function 89,175,543 times, while this method calls the limit-state function 485, 960 and 1340 times, respectively. of Monte Carlo is set to 100,000. Table 6 and Figure 5 show the calculation results. It can be observed that when Δt = 1 year, Δt = 0.5 year and Δt = 1/3 year, the maximum errors are 25.4%, 9.6% and 5.3% respectively. It can be found that with the decrease of the time step, the accuracy of reliability solution is gradually improved. In terms of efficiency, Monte Carlo needs to call the limit-state function 89,175,543 times, while this method calls the limit-state function 485, 960 and 1340 times, respectively.

A Vehicle Frame
Consider the frame structure of a commercial vehicle as shown in Figure 6, which consists of two longitudinal beams and multiple transverse beams. The frame is the base of the whole car. Most parts and assemblies of the car are fixed by the frame. These connecting parts produce loads on the frame. The static finite element model of the frame is obtained by simplifying various constraints and structural loads. The maximum displacement in Y direction after frame deformation can represent its stiffness, which is an evaluation standard of vehicle performance. Therefore, the limit-state function can be expressed as follows: where d m (t) = d 0 e −0.01t denotes maximum allowable vertical displacement of structure, d 0 represents the maximum allowable initial displacement. d(th, E, ρ, Q 1 (t), t) is the maximum displacement calculated by finite element software. The thickness th 1 -th 5 of the key components of the frame are random variables, Q 1 (t) is time-variant load, and the density ρ and elastic modulus E of the material are interval variables. Tables 7 and 8       The finite element model of the frame consists of 558,178 shell elements and 297,484 nodes in Figure 6. The limit-state function is constructed through the Kriging model to realize parameterization and improve the calculation efficiency of reliability analysis. Table 9 and Figures 7 and 8 show the calculation results. It can be found that when considering the stiffness failure, the reliability of the frame is on the decline. When T = 1, the reliability index is 3.21; when T = 10, the reliability index is 2.33. That means the probability of failure increases. The failure probability is 0.066% in the first year and 0.99% in the tenth year. realize parameterization and improve the calculation efficiency of reliability analysis. Table 9 and Figures 7 and 8 show the calculation results. It can be found that when considering the stiffness failure, the reliability of the frame is on the decline. When T = 1, the reliability index is 3.21; when T = 10, the reliability index is 2.33. That means the probability of failure increases. The failure probability is 0.066% in the first year and 0.99% in the tenth year.

Conclusions
A time-variant reliability analysis method considering interval variables and stochastic processes is proposed in this paper. This method can handle time-variant reliability calculation when some structural parameters cannot obtain accurate probability distribution due to insufficient samples. By discretizing the stochastic process in time, the time-variant reliability problem is changed into a static hybrid reliability problem. The limit-state function at each time point is expanded at MPP and MUP points. This method avoids the difficulty in solving the crossing rate problem and is easy to understand conceptually. The efficient solution format greatly simplifies the multi-layer nested solution process of time-variant reliability.
The results of numerical example analysis show that the structure's reliability will not be a fixed value due to the degradation of materials and dynamic uncertain loads but will gradually decrease with the increase of design reference period. Therefore, to ensure

Conclusions
A time-variant reliability analysis method considering interval variables and stochastic processes is proposed in this paper. This method can handle time-variant reliability calculation when some structural parameters cannot obtain accurate probability distribution due to insufficient samples. By discretizing the stochastic process in time, the time-variant reliability problem is changed into a static hybrid reliability problem. The limit-state function at each time point is expanded at MPP and MUP points. This method avoids the difficulty in solving the crossing rate problem and is easy to understand conceptually. The efficient solution format greatly simplifies the multi-layer nested solution process of time-variant reliability.
The results of numerical example analysis show that the structure's reliability will not be a fixed value due to the degradation of materials and dynamic uncertain loads but will gradually decrease with the increase of design reference period. Therefore, to ensure the performance of the structure in the whole service period, it is necessary to carry out time-variant reliability analysis of the structure. As the size of time step decreases, the analysis results of the present method approach those of Monte Carlo method. Besides, the approach has high computational efficiency and can satisfy the needs of complex engineering problems. However, this method is mainly applied to the case that the load changes slowly in the time-variant stochastic process. For the dynamic reliability problem that the stochastic process changes violently, the crossing rate method can be used for reliability analysis.
In the future, it is necessary to develop the mixed time-variant reliability of random variables and other uncertain variables, such as fuzzy variables and evidence variables. The hybrid time-varying reliability optimization method is further extended, such as multi-disciplinary reliability design optimization, multi-objective reliability optimization and so on.