An adaptive multiple-Kriging-surrogate method for time-dependent reliability analysis

Accepted Manuscript

An adaptive multiple-Kriging-surrogate method for time-dependent reliability analysis Yan Shi , Zhenzhou Lu , Liyang Xu , Siyu Chen PII: DOI: Reference:

S0307-904X(19)30070-8 https://doi.org/10.1016/j.apm.2019.01.040 APM 12646

To appear in:

Applied Mathematical Modelling

Received date: Revised date: Accepted date:

26 September 2018 26 December 2018 29 January 2019

Please cite this article as: Yan Shi , Zhenzhou Lu , Liyang Xu , Siyu Chen , An adaptive multipleKriging-surrogate method for time-dependent reliability analysis, Applied Mathematical Modelling (2019), doi: https://doi.org/10.1016/j.apm.2019.01.040

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT Highlights Multiple-Kriging-surrogate is simultaneously constructed for time-dependent structure. Adaptive strategy is established to identify surrogates and new training samples. Man-made subjectivity for regression trend is avoided with adaptive selecting procedure.



Time-dependent failure probability can be accurately estimated by the most accurate local model.

AC

CE

PT

ED

M

AN US

CR IP T

  

ACCEPTED MANUSCRIPT

An adaptive multiple-Kriging-surrogate method for time-dependent reliability analysis Yan Shi a, Zhenzhou Lu a, Liyang Xu a, Siyu Chen a a School of Aeronautics, Northwestern Polytechnical University, Xi’an, 710072, P. R. China

1.

AN US

CR IP T

Abstract: For accurately and efficiently estimating the time-dependent failure probability (TDFP) of the structure, a novel adaptive multiple-Kriging-surrogate method is proposed. In the proposed method, the multiple Kriging models with different regression trends (i.e., constant, linear and quadratic) are simultaneously constructed with the highest accuracy, on which the TDFP can be obtained. The multiple regression trends are adaptively selected based on the size of sample base, the maximum differences of multiple models and the global accuracy of multiple models. After that, the most suitable multiple regression trends are identified. The proposed method can avoid man-made subjectivity for regression trend in general Kriging surrogate method. Furthermore, better accuracy and efficiency will be obtained by the proposed multiple surrogates than just using a fixed regression model for some engineering applications. Five examples involving four applications with explicit performance function and one tone arch bridge under hurricane load example with implicit performance function are introduced to illustrate the effectiveness of the proposed method for estimating TDFP. Keywords: Multiple-Kriging-surrogate; Regression trend; Time-dependent failure probability; Single-loop procedure; Kullback-Leibler divergence. Introduction

AC

CE

PT

ED

M

Nowadays, time-dependent reliability analysis has gained more and more attentions due to its important significance in measuring the safety degree of the time-dependent structure. Comparing with the time-independent reliability analysis, the time-dependent reliability is generally more complicated because the time-dependent output varies not only with the input variables but also with the time parameter [1-2]. Therefore, it is vital to give some accurate and efficient methods to analyze the time-dependent reliability. The purpose of the time-dependent reliability analysis is to estimate the time-dependent failure probability (TDFP). Up to now, many methods have been proposed to estimate the TDFP. The first passage based methods [3-9] are widely used for estimating the TDFP. Andrieu-Renaud et al. [3] proposed the PHI2 approach to estimate the TDFP and Sudret [4] further developed this approach by deriving the analytical expression of the outcrossing rate. Zhang, Du [5] and Jiang et al. [6] proposed the outcrossing rate models to estimate the TDFP for the function generator mechanisms and time-dependent system respectively. Hu and Du [7] gave the joint outcrossing rates model to estimate the TDFP, which considers the dependency of the outcrossing events at different time instants. For the time-dependent structure with temporal and spatial multiple parameters, Shi et al. [8] proposed the effective first-crossing point method to estimate the TDFP. Beck and Melchers [9] proposed a simulation-based technique to estimate the error functions in ensemble first-crossing rate to obtain the TDFP. The extreme value based methods [10-15] are other efficient ones to estimate the TDFP. Li et al. [10] employed the equivalent extreme value event to estimate the TDFP, and Chen and Li [11] used the extreme value distribution to estimate the TDFP. Du [12] proposed the envelope function method by linearized the performance function of the time-dependent structure, and the envelope function method 

Corresponding author E-mail: [email protected]

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

was further developed by Shi et al. [13] for the time-dependent structure with temporal and spatial multiple parameters. For the time-dependent structure with only one input stochastic process and the performance function is implicit with respect to time parameter case, Hu and Du [14] employed a sampling method to calculate the extreme value distribution of the input stochastic process and further estimate the TDFP with the extreme value distribution. By computing the extreme value moments of the performance function, Shi et al. [15] proposed some efficient methods for estimating the TDFP. For obtaining an accurate result of the TDFP, the first passage based methods [3-9] and the extreme value based methods [9-15] may need large computational cost especially for the time-dependent structure with highly nonlinear performance function. In view of this, the surrogate model based methods [16-28] have gained researchers’ attention because of their good trade-off between accuracy and efficiency. Zhang et al. [16] proposed a time-dependent reliability analysis method by introducing the response surface model, in which a gradient projection technique is employed to generate new training samples to update the surrogate. By employing a double-loop surrogate procedure to construct an extreme value surrogate model, Wang and Wang [17] proposed the nested extreme value response surface method and employed it to time-dependent reliability based design optimization [18]. Hu and Du [19] improved the extreme value response surface method [17] by using a mixed efficient global optimization technique to obtain the extreme value. Wang and Chen [20] constructed the Gaussian process models to predict the outputs at each time instants to approximate the extreme response surface. These surrogate model based methods [16-20] implement a double-loop procedure to construct the surrogate model, and one shortcoming of this double-loop procedure is the estimation accuracy of the extreme value in the inner loop will affect the accuracy of the constructed surrogate model in the outer loop. The single-loop surrogate model method [21-22] directly constructs the surrogate model for the performance function and overcomes the shortcoming of the double-loop surrogate procedure, thus it is more efficient than those of the double-loop surrogate model methods. Wang and Chen [21] presented an equivalent stochastic process transformation method which firstly constructs the surrogate of the instantaneous failure surface, and then updates the surrogate by using a maximum confidence enhancement approach. Hu and Mahadevan [22] proposed to directly construct the surrogate model of original limit state function and then adaptively update the surrogate by using a learning function. Hawchar et al. [23] established a time-dependent reliability analysis method by combining the principle component analysis and polynomial chaos expansion, it uses the principle component analysis to reduce the number of components and then constructs surrogate model by the reduced solutions. Other surrogate model based methods have also been studied in Refs. [24-29]. All above surrogate model based methods [16-29] employ only one surrogate model to approximate the performance function or the extreme value distribution. It has been showed that the accuracy of the surrogate model predictions can be increased by combining individual surrogate model in the form of an ensemble [30]. Therefore, multiple surrogates [31-33] have gained more and more attention recently. Goel et al. [31] approximated the function relationship by using the best surrogate or a weighted average surrogate model instead of individual surrogate model, in which the weights associated with each surrogate model are identified by the errors in surrogates. Zhang et al. [32] proposed a high fidelity surrogate modeling technique to exploit the advantages of each component surrogate. Vahedi et al. [33] proposed an adaptive divergence-based multiple surrogates method to the time-independent reliability analysis. It should be noted that there are generally three types of regression trends in the Kriging model, i.e., constant regression, linear regression and quadratic regression. In the single-loop Kriging surrogate method [22], the choice of the regression trend in the Kriging model relies on specific problems and generally there is no prior information for identifying a suitable regression trend and it is always man-made. The aim of this work is to propose an adaptive multiple-Kriging-surrogate method which combines different regression trends to the time-dependent reliability analysis, and thus it can avoid

ACCEPTED MANUSCRIPT

2.

CR IP T

man-made subjectivity for regression trend in general Kriging surrogate method. In the proposed method, multiple Kriging models are directly constructed in the single-loop procedure. Furthermore, the Kullback-Leibler (KL) divergence [33-34] is employed to measure the differences among multiple surrogate models and the new training samples are identified by the differences among multiple surrogate models. The multiple Kriging surrogate models are used simultaneously to reduce the sensitivity of the Kriging surrogate to the type of regression trend. To the authors’ knowledge, multiple Kriging surrogates have not been employed for time-dependent reliability analysis of structure. This work is organized as follows. The definition of the time-dependent reliability is given in Section 2. The proposed adaptive multiple-Kriging-surrogate method is showed in Section 3. The estimation procedure of the proposed method is provided in Section 4. Several examples are introduced in Section 5. Conclusions are drawn in Section 6. Definition of the time-dependent reliability

It is well known that the performance function of the time-dependent structure can be expressed by

G  g  X, Y (t ), t  , where X =  X1 , X 2 ,..., X n 

T

random variables. Y(t )  Y1 (t ), Y2 (t ),..., Ym (t )

is the n -dimensional vector of independent input

with

mean

vector

is the m -dimensional input stochastic processes

AN US

T

 Y (t )   Y1 (t ), Y2 (t ),..., Ym (t )  T

T

and

standard

deviation

vector

 Y (t )   Y1 (t ),  Y2 (t ),...,  Ym (t )  . The auto-covariance function of the input stochastic process

Y j (t )( j  1, 2,..., m) is described by CY j (t , )   Y j (t ) Y j ( ) Y j (t , ) , in which Y j (t , ) represents the autocorrelation coefficient of the input stochastic process Y j (t )( j  1, 2,..., m) at time instants t

M

and  . t  t0 , te  is the time parameter. The failure corresponds to the output less than or equal to zero, thus the TDFP Pf  t0 , te  of the time-dependent structure is defined as follows:

ED

Pf  t0 , te   P G  g  X, Y(t ), t   0, t  t0 , te 

(1)

PT

in which  means there exists a value. Because the values of the stochastic process Y j (t )( j  1, 2,..., m) at different time instants are

CE

dependent, thus it is difficulty to directly deal with the stochastic process in estimating the TDFP. Fortunately, several discretization techniques such as the Expansion Optimal Linear Estimation (EOLE) [23] and the Orthogonal Series Expansion (OSE) [35] have been proposed to simulate the stochastic process by a set of independent random variables. In this work, the OSE method is used to describe the

AC

input stochastic process Y j (t )( j  1, 2,..., m) as follows: Mj

Y j (t )  Y j (t )    jk  jk (t ) jk

(2)

k 1

Mj

k where  jk (t )    kji hi (t ) ,  jk is a set of independent standard normal variable,  jk and  ji are i 1

the eigenvalue and the corresponding eigenvector of the auto-covariance function CY j (t , ) , hi (t ) means a set of orthogonal function and M j is the number of the truncated terms. More details about the OSE method are showed in Ref. [35] and Appendix A. Therefore, the input stochastic process

Y j (t )( j  1, 2,..., m)



is



expressed

by

a

set

of

independent

 j   j1 ,  j 2 ,...,  jM j ( j  1, 2,..., m) . Let us define   1 ,  2 ,...,  m 

T

standard

normal

variable

to represent the vector of input

ACCEPTED MANUSCRIPT stochastic processes Y(t ) .

3.

CR IP T

Time-dependent reliability analysis is more complicated than the time-independent one, and directly solving Eq. (1) by Monte Carlo Simulation method will introduce huge computational cost especially for the engineering application. As discussed in the Introduction, the surrogate model based methods can give good trade-off between accuracy and efficiency, and thus it is necessary to be used for estimating TDFP. The double-loop surrogate method treats the input variables and time parameter at different levels and constructs the surrogate model based on the extreme values of output, but the estimation accuracy of the extreme values in the inner loop will affect the accuracy of the constructed surrogate model in the outer loop. The single-loop surrogate method directly constructs the surrogate model of original limit state function and shows good role in accessing the safety of time-dependent structure. Therefore, studying the single-loop surrogate method for estimating the TDFP is significant. In the next section, a new multiple-Kriging-surrogate single-loop method is proposed to efficiently estimate the TDFP. The proposed adaptive multiple-Kriging-surrogate method

AN US

3.1 Basic principle of constructing multiple-Kriging-surrogate

theory is expressed as follows:

M

The Kriging model [36] is a nonparametric interpolation model that hypothesizes the true function is the realization of a Gaussian process. This model employs the training samples and the corresponding outputs to construct the surrogate model for forecasting the output. Up to now, the Kriging model has been widely used in the field of reliability assessment. The advantage of the Kriging model is that it not only provides the mean value of prediction but also gives the standard deviation of prediction, they can usually be used to construct the learning function to identify new training samples to update model. For the performance function G  g  X, Y (t ), t  , the Kriging model Gˆ  gˆ  X, Y, t  based on single-loop T Gˆ  f  x, y, t      x, y, t 

ED

(3)

in which f  x, y,t   is the regression part of the Kriging model,   x , y ,t  is stationary zero mean T

f  x, y, t    f1  x, y, t  , f 2  x, y, t  ,..., f p  x, y, t 

PT

value Gaussian process.

basic function vector, and    1 ,  2 ,...,  p 

T

T

represents the p -size

is the corresponding regression coefficient vector.

CE

More details about the Kriging can be found in Ref. [36]. Actually, there are three types of regression trends in the Kriging model, i.e., constant regression, linear regression and quadratic regression [36]. Let us denote the inputs as z   x , y ,t  . Because the

AC

dimensional size of the inputs is n  m  1 , then the three types of regression trends with unknown coefficients i (i  1, 2,..., p) can be respectively expressed as follows: The constant regression corresponds to p  1 :

f  z    0 T

(4)

The linear regression corresponds to p  n  m  2 : f  z    0  T

The quadratic regression corresponds to p 

n  m 1

 i 1

i zi

 n  m  2 n  m  3 2

(5) :

ACCEPTED MANUSCRIPT

f  z    0  T

n  m 1

 i 1

i zi 

n  m 1 n  m 1

  i 1

j 1

ij zi z j

(6)

That’s to say, we will obtain different surrogates by using different types of regression trends with the same group training samples and corresponding outputs. For brief, we firstly define the an ensemble of models as  which contains one or several these models, and let us denote the particular model member in  as  . Then the mean value and standard deviation of prediction for the  -th Kriging model can be expressed by Gˆ  x , y , t  and  Gˆ  x , y , t  respectively. Thus the U-learning function [22] for the time-dependent structure is defined by

Gˆ  x, y, t 

 Gˆ  x, y, t 

(7)

CR IP T

U  x , y , t  

The U-learning function U  x , y , t  is equivalent to the index of the misclassification probability of the sign of the performance function by the  -th surrogate model. Generally, the corresponding prediction can be considered as accurate value when U   x , y , t   2 , and the corresponding prediction is considered as inaccurate enough one when U   x , y , t   2 .

AN US

The single-loop Kriging surrogate method [22] firstly directly constructs the surrogate model of original performance function with a few training samples by treating the time parameter as uniform distribution in the time interval, and then updates the surrogate model by using the U-learning function

U  x , y , t  showed in Eq. (7). One thing to be mentioned is that this method uses a fixed regression

ED

M

trend to construct the surrogate model. The regression trend must be given as one of these three trends, i.e., constant regression, linear regression or quadratic regression, before constructing the surrogate model. But the choice of regression trends will affects the accuracy of the final surrogate model, and it is unknown which one is the best before constructing the surrogate model. For dealing with this problem, the multiple-Kriging-surrogate method is proposed in this work. It should be noted that when constructing the Kriging model with different regression trends, the number of the training samples should not be less than the corresponding p because that the number of the regression coefficients i (i  1, 2,..., p ) equals to p . Therefore, multiple Kriging models are constructed by considering the size relationship between the number of training samples and the value

PT

of p . Let us denote the number of training samples as N S , then the strategy for constructing multiple Kriging models is showed below.

CE

(a) If 1  N S  n  m  2 , one Kriging model with constant regression trend is constructed. (b) If n  m  2  N S 

 n  m  2 n  m  3

AC

2 regression trends are constructed.

 n  m  2 n  m  3

, two Kriging models with constant and linear

 N S , three Kriging models with constant, linear and quadratic 2 regression trends are constructed. It should be mentioned that the multiple Kriging models are constructed by the same training samples and the corresponding performance responses, thus it doesn’t introduce extra computational cost of calling the original performance function when constructing Kriging models. The basic principle for constructing multiple Kriging surrogates can be briefly summarized as follows. Firstly, construct the multiple Kriging models with the initial training samples by above strategy. Secondly, select the new training samples by the method proposed in subsection 3.3 to update the multiple Kriging models and identify whether it is necessary to add the Kriging model or delete the Kriging model in the ensemble of models. Thirdly, employ the final multiple Kriging surrogate models (c) If

ACCEPTED MANUSCRIPT to estimate the TDFP. 3.2 Global accuracy and local accuracy Based on the definition of the TDFP Pf  t0 , te  showed in Eq. (1), Pf  t0 , te  can be estimated by the minimum value of performance response with respect to time parameter as follows: Pf  t0 , te   P Gmin  gmin  X, Y   0

(8)

where g min  X, Y   min  g  X, Y(t ), t  . Therefore, for each Kriging model, the minimum value of tt0 , te 

TDFP, where tmin   arg min Gˆ  x , y , t  tt0 , te 

The corresponding U-learning function is showed below

U  x, y, tmin   

Gˆ  x, y, tmin  

 Gˆ  x, y, tmin  

CR IP T

Gˆ  x , y , t  with respect to time parameter, i.e., Gˆ  x , y , tmin   , is the key point to estimate the (9)

(10)

AN US

where  Gˆ  x , y , tmin   represents the standard deviation of prediction at time instant tmin  . Therefore, the following measure  can be used to compare the global accuracy of multiple surrogate models.

 

NU  x , y ,tmin  2 N MCS

(11)

where N MCS is the size of input variable samples in the database, NU  x , y ,tmin  2 represents the

M

number of samples corresponding to U  x, y, tmin    2 .  can be regarded as the error of prediction for the  -th surrogate model. Thus a surrogate model with the smallest value of  can

ED

be regarded as possessing the best global accuracy. Furthermore, we standardize this measure  to

PT

the following expression, which represents the ration of measure for the  -th surrogate model to the average measure.

 

 N   K    NK  

 

 

in which N K is the number of the surrogate models in  ,

CE

(12)

 

represents the sum of the measure

 

 in the ensemble of models  . It is easy to understand that the smaller  is, the more accurate of the  -th surrogate model is. Therefore, if  is larger than a specified value  0 , the  -th

AC

surrogate model may be inaccurate enough. Referring to the study in Ref. [33],  0 is set to be 1.5. If

 is larger than 1.5, then the  -th surrogate model is regarded as an inaccurate surrogate, and this

surrogate should be removed from  . After defining the global accuracy of multiple surrogates, the local accuracy of multiple surrogates is given below, which is relative to the estimation of the TDFP. In the proposed method, there are multiple Kriging surrogate models and each Kriging model has the predicted mean value and standard deviation. In other words, there are multiple group predictions with the multiple surrogates, but only one group prediction is needed to estimate the TDFP. Considering the discussion of the U-learning function U  x , y , t  in subsection 3.1, a surrogate model corresponding to the maximum value of U-learning function for the input

 x , y ,t 

is

ACCEPTED MANUSCRIPT considered as the most accurate local model. Note that the most accurate local model may be different from one sample to another. Thus the most accurate local model denoted by local can be obtained by

local  arg max U  x, y, t 

(13)

 

Let us denote the mean value and standard deviation of the prediction corresponding to the most accurate local model as Gˆ  x , y , t  and  Gˆ  x , y , t  respectively, and let us denote the U-learning function as U  x , y , t  corresponding to the local model at

 x, y, t 

(14)

 Gˆ  x, y, t    Gˆ

 x, y, t 

(15)

Gˆ

local

 Gˆ

local

 x, y, t   x, y, t 

local

local

U  x, y, t  

CR IP T

Gˆ  x, y, t   Gˆ

 x , y ,t  , then they can be obtained by

(16)

With the most accurate local model local , the multiple surrogates can be handled as one. That’s to say, we can obtain one prediction mean value Gˆ  x , y , t  and one prediction standard deviation

 x , y ,t 

with the most accurate local model local .

AN US

 Gˆ  x , y , t  for the input

For illustrating the process of integrating the multiple surrogates to obtain the most accurate local model local , a specific example showed in Table 1 is employed. Suppose there are three surrogates in the ensemble of model, and we employ these surrogates to estimate the values of Gˆ  x , y , t  ,

 Gˆ  x , y , t 

and U  x , y , t 

for two samples

x

M

respectively. One can see that for the first sample

x

(1)

(1)

, y (1)  t (1)  , t (1)



and

x

(2)

, y (2)  t (2)  , t (2)





, y (1)  t (1)  , t (1) , the surrogate with the constant

regression trend has the biggest value of U  x , y , t  , thus the surrogate with the constant regression trend is the most accurate local model local (2)



, y (2)  t (2)  , t (2) , the surrogate with the quadratic regression trend has the biggest value of

ED

x

for this sample. For the second sample

U  x , y , t  , thus the surrogate with the quadratic regression trend is the most accurate local model

PT

local for the second sample. Therefore, if there are N MCS samples, then we can obtain N MCS most accurate local models and the corresponding N MCS prediction mean values Gˆ  x , y , t  and standard

CE

deviations  Gˆ  x , y , t  , and the TDFP can be estimated by these solutions. Table 1 The illustration of identifying the most accurate local model local

AC

 x , y ,t 

x

x

(1)

(2)

, y (1)  t (1)  , t (1)

, y (2)  t (2)  , t (2)





Items

Constant

Linear

Quadratic

local

U  x , y , t 

24.0

20.8

18.6

24.0

Gˆ  x , y , t 

2.40

2.50

2.60

2.40

 Gˆ  x , y , t 

0.10

0.12

0.14

0.10

U  x , y , t 

18.3

16.3

23.3

23.3

Gˆ  x , y , t 

2.38

2.45

2.56

2.56

 Gˆ  x , y , t 

0.13

0.15

0.11

0.11

ACCEPTED MANUSCRIPT 3.3 Selecting new training samples In this paper, the Kullback-Leibler (KL) divergence [33-34] is used to identify the new training sample in each iterative step, which has been employed for the time-independent structure [33]. The KL divergence measures the differences among multiple surrogate models, and the new training sample corresponds to the maximum average KL divergence in each iteration. The average KL divergence

 KL  x , y , t  can be expressed as follows: 1 NK

 P Gˆ  0  x, y, t    ˆ    P G  0 x , y , t log        ˆ      P G  0  x, y, t    

(17)

CR IP T

 KL  x, y, t  

P Gˆ  0  x, y, t     P Gˆ  0  x, y, t   log     ˆ P G  0  x, y, t      where

NK

is the number

of the surrogate

models in

 ,

P Gˆ  0  x, y, t   

AN US

P Gˆ  0  x, y, t  represent the probabilities of Gˆ  0 and Gˆ  0 respectively at   based on the  -th surrogate model.

and

 x , y ,t 

 0  Gˆ  x, y, t   P Gˆ  0  x, y, t          Gˆ  x, y, t  

(18)

 0  Gˆ  x, y, t   P Gˆ  0  x, y, t    1         Gˆ  x, y, t  

(19)

M

P Gˆ  0  x, y, t  and P Gˆ  0  x, y, t  mean the average probabilities of Gˆ  0 and Gˆ  0     respectively with    , and they are defined by

 P Gˆ  0  x, y, t 

(20)

 

1 P Gˆ  0  x, y, t      NK

P Gˆ  0  x, y, t     

(21)

PT

ED

1 P Gˆ  0  x, y, t      NK



Then, the new training sample can be obtained by maximizing the average KL divergence  KL  x , y , t 

CE

across all the members in  .

x

new

, y new , t new   arg max KL  x, y, t 

(22)

AC

However, it does not require high estimation accuracy for all the inputs in the time-dependent reliability analysis, which is different from the classification problem in the time-independent reliability analysis. For the time-dependent structure, the time-dependent failure event can be accurately determined if there exists a failure point (time instant) on the trajectory and the failure point is accurately classified no matter the signs of other points (time instants) are accurately classified or not. Therefore, for any given input variable

 x, y  , whether the trajectory

g  x , y , t  , t  t0 , te  needs to

be refined or not can be classified into two cases. The first case is that the performance satisfies

Gˆ  x , y , t   0 for all the time instants in t0 , te  , where Gˆ  x , y , t  is estimated by the most accurate local model given in subsection 3.2. In this case, the signs of surrogate models for all the time instants in t0 , te  need to be accurately classified, and new training samples are needed to refine this trajectory until the minimized U  x , y , t  satisfies the threshold requirement, where U  x , y , t  is

ACCEPTED MANUSCRIPT estimated by the most accurate local model with Gˆ  x , y , t  and  Gˆ  x , y , t  . The second case is that  there exists a time instant t  t0 , te  which satisfies that Gˆ  x, y, t    0 and U  x, y, t    2 . It

illustrates that the time-dependent failure event has been accurately identified and the trajectory

g  x , y , t  , t  t0 , te  does not need to be refined. Therefore, we should check if there exists a time instant

t    t0 , te 

Gˆ  x, y, t    0

meeting

g  x , y , t  , t  t0 , te  for the input variable

 x, y 

U  x, y, t    2

and

on

the

trajectory

in each iterative step.

 if Gˆ  x, y, t    0 and U  x, y, t    2 check :    subject to : t  t0 , te 

CR IP T

(23)

new Based on Eq. (23), a new average KL divergence  KL  x , y , t  can be defined by

    0 Gˆ  x, y, t   0 and U  x, y, t   2, t  t0 , te 

new KL  x, y, t   

 KL  x, y, t  otherwise

Suppose the database of the input samples is

x

(i )

(24)



, y (i )  t ( j )  , t ( j ) (i  1, 2,..., N MCS , j  1, 2,..., N t ) ,

Then, the new training sample

x

new

x

new

AN US

where N MCS is the number of random inputs and N t is the number of the discretized time instants. , y new ,t new  can be identified by

 

new , y new , t new   arg max KL x (i ) , y (i )  t ( j )  , t ( j )



i  1, 2,..., NMCS , j  1, 2,..., Nt



new x (i ) , y (i )  t ( j )  , t ( j ) It is impossible that all the new average KL divergences  KL





(25)

for the samples in

ED

M

database equal to zero. The reason is that if all the new average KL divergences equal to zero, then all the time-dependent failure events have been accurately identified and the iteration should be stopped with the stopping criterion showed in next subsection 3.4. Thus there should exist one new training sample being selected in each iteration. It should be noted that if there is only one surrogate model in  , then the process for selecting the new training samples will degrade into that procedure in the single-loop Kriging surrogate model method [22], which employs the U-learning function to identify the new training samples.

PT

In subsection 3.2, the metric  is proposed to measure the global accuracy of multiple surrogate models, and the surrogate model with   1.5 should be removed from the ensemble of models. Next, another strategy by considering the differences between multiple surrogates is established.

CE

The average KL divergence  KL  x , y , t  measures the differences among multiple surrogates, and

AC

the larger  KL  x , y , t  , the more differences among multiple surrogates. Because  KL  x , y , t  is varying with different input

 x , y ,t  , thus it is inconvenient to measure the differences among multiple

surrogates by using  KL  x , y , t  directly. In this work, the maximum average KL divergence

KL max = max KL  x, y, t  of all the sample base is a value which can be conveniently employed to

measure the differences among multiple surrogates. If  KL max is less than a specified value of KL 0 , which means the differences among multiple surrogates are small and the multiple surrogates are not diverse enough. This illustrates that it is no necessary to employ multiple surrogates, and only one surrogate should be used. In view of this, the surrogate model with the best global accuracy is reserved, and other surrogates should be removed from the ensemble of models. By referring to the study in Ref. [33], KL 0 is set to be 0.15 in this work.

ACCEPTED MANUSCRIPT 3.4 Stopping criterion In this work, the following stopping criterion [37] is employed.

P 

Pf  t0 , te   Pf  t0 , te 

(26)

Pf  t0 , te 

f

where Pf  t0 , te   P Gˆ min  x, y   0

(27)

Pf  t0 , te   P Gˆ min  x, y   Gˆ min  x, y   0

(28)

in which Gˆ min  x , y   min Gˆ  x , y , t  , tmin  arg min Gˆ  x , y , t  ,  Gˆ min  x , y    Gˆ  x , y , tmin  , tt0 , te 

CR IP T

tt0 , te 

Gˆ  x , y , t  and  Gˆ  x , y , t  are estimated by the most accurate local model defined in subsection 3.2.  is the confidence level typically equals to  1  97.5%   1.96 . The acceptable threshold value is set to be  0  0.05 , thus the adaptive process will be stop when  Pf   0 . Estimation procedure

AN US

4.

In the adaptive process, if one surrogate model is removed from the ensemble of models, then it will never be added in the ensemble although the number of training samples N S satisfies the condition of

M

adding new surrogate model given in subsection 3.1. If there is only one surrogate model in the ensemble of models, then the process for selecting the new training samples will degrade into that procedure in the single-loop Kriging surrogate model method [22]. The estimation procedure of the proposed adaptive multiple-Kriging-surrogate models method for estimating the TDFP can be briefly summarized as follows: Step 1: Set the acceptable threshold values of  0 , KL 0 and  0 . Uniformly discretize the time

the sample base

x

(i )

ED

(N ) (1) (2) parameter into N t time instants in the interval t0 , te  , i.e., t0  t  t  ...  t t  te . Generate

, (i )  (i  1, 2,..., NMCS ) of input variables, and construct the sample matrices

PT

A N MCS  Nt with the vector quantity element

x

(i )

, y (i )  t ( j )  , t ( j )



in the i(i  1, 2,..., N MCS ) row and

j ( j  1, 2,..., Nt ) column, where y (i )  t ( j )  is converted by  (i ) at t ( j ) based on Eq. (2).

CE

Step 2: Generate N S initial training samples

x

( i )

, (i ) , t (i )  (i  1, 2,..., NS ) , where the samples of

time parameter are generated by treating the time parameter as the uniform distribution in the interval

 t0 , te 



, y   t (i )  , t (i ) (i  1, 2,..., N S ) . Estimate the corresponding performances for these samples by the

AC

x

. Convert   ( i ) to y  t (i )  by Eq. (2), then we can obtain the training samples

( i )

real performance function. Step 3: Construct the multiple-Kriging-surrogate models with the existing training samples and the corresponding performances by the strategy proposed in subsection 3.1.



Step 4: Estimate the prediction mean value Gˆ x (i ) , y ( i )  t ( j )  , t ( j )



 Gˆ x (i ) , y ( i )  t ( j )  , t ( j )



multiple-Kriging-surrogate



(i )

U x , y

(i )

t  , t  . ( j)

( j)

for

the

models,

sample and

matrices

calculate

the



and prediction standard deviation

A N MCS  Nt values

by of

the U-learning

constructed function

ACCEPTED MANUSCRIPT



Step 5: Compute the average KL divergence  KL x (i ) , y ( i )  t ( j )  , t ( j )



Gˆ x (i ) , y ( i )  t ( j )  , t ( j )







with the solutions of



and  Gˆ x (i ) , y ( i )  t ( j )  , t ( j ) . Obtain the maximum average KL divergence

 





KL max = max KL x (i ) , y (i )  t ( j )  , t ( j ) , i  1, 2,..., NMCS , j  1, 2,..., Nt . Step 6: Estimate the global accuracy  of the multiple-Kriging-surrogate models. If KL max  KL 0 is satisfied, only the surrogate with the best global accuracy is reserved and removed other surrogates from the ensemble of models. If KL max  KL 0 , check if there are some models satisfy   0 in two consecutive iterations, and removed these models from the ensemble of models.



 Gˆ x (i ) , y (i )  t ( j )  , t ( j )





and prediction standard deviation

CR IP T



Step 7: Obtain the prediction mean value Gˆ x (i ) , y (i )  t ( j )  , t ( j )

with the most accurate local model, and estimate the solution of  Pf . If

 Pf   0 , stop the adaptive process and go to Step 9; otherwise, go to Step 8.





new x (i ) , y (i )  t ( j )  , t ( j ) , and identifying the new training sample by Step 8: Estimate the value of  KL

new

 





new , y new , t new   arg max KL x (i ) , y (i )  t ( j )  , t ( j ) , i  1, 2,..., N MCS , j  1, 2,..., Nt . Adding this new

AN US

x

training sample into the training sample base and estimating the output for this new training sample by the real performance function. Let N S  N S  1 and go to Step 3. Step 9: Obtain the final TDFP as Pf  t0 , te  . Examples

M

5.

5.1 Numerical example

ED

Considering the time-dependent structure with the following performance function,

g  X, t   X 12 X 2  5 X 1t   X 2  1 t 2  13

(29)

where X 1 and X 2 are the input normal random variables with same mean value 3.5 and standard

PT

deviation 0.3. The considered time interval is

 0, 5 .

In numerical example, the number of sample base of the input variables is set to be N MCS  105 , and

CE

the time interval is uniformly discretized into Nt  100 time instants. Because the single-loop Kriging

AC

surrogate model method [22] constructs the surrogate model of the performance function by a fixed regression trend, thus we provide the TDFP solutions estimated by the single-loop Kriging surrogate model method with different regression trends. The TDFP solutions estimated by different methods are listed in Table 2. The maximum average KL divergence and the global accuracy of the proposed method are showed in Fig. 1. The convergence history of the TDFP of the proposed method is showed in Fig. 2. The solutions of the finial limit state estimated by the proposed method in the standard normal space are showed in Fig. 3. From Table 1, one can see that comparing with the single-loop Kriging surrogate method with different regression trends, i.e., constant, linear and quadratic, the proposed multiple-Kriging-surrogate method can give more accurate estimation for the TDFP by less computational cost. This is because that the proposed multiple-Kriging-surrogate method can adaptively select the multiple suitable regression trends to construct the surrogate models of the performance function, and the final TDFP is estimated by the most accurate local model which has most confidence to give accurate prediction. Furthermore, samples which have the most contribution to the TDFP are identified to update the

ACCEPTED MANUSCRIPT multiple surrogates. It can be seen from Fig. 1 that three surrogates with constant, linear and quadratic regression trends are selected to construct the surrogate models of the performance function in the

 2  0  2  2  0  3

 12 in the beginning of 2 the update process, and three regression trends are chosen to the ensemble of models. The model with

initial phase of the adaptive process. This is because

linear regression trend satisfies   0 in two consecutive iterations, thus this model is removed

AN US

CR IP T

from the ensemble of models after the second iteration. The models with constant and quadratic regression trends are reserved in the ensemble of models to the end. Fig. 1 shows that the final ensemble of models contains two models with constant and quadratic regression trends, and the TDFP is estimated by the most accurate local model between these two models (the most accurate local model may be different for different samples). Fig. 2 shows that the TDFP gradually converges to a convergent solution as the increasing of training sample number, which illustrates the good convergence of the proposed method. Fig. 3 shows that not all the approximate instantaneous limit states match well with the real instantaneous limit states. Generally, it is unnecessary to make all the approximate instantaneous limit states have good accuracy on real instantaneous limit states. For time-dependent reliability analysis, only the time-dependent limit state decides the failure of the time-dependent structure. The solutions show that the approximate time-dependent limit state estimated by the proposed method can match well with that of the real time-dependent limit state in the 5-sigma domain of the input space. Table 2 The solutions of the TDFP by different methods MCS

SLK-Constant

SLK-Linear

SLK-Quadratic

PR

Pf  t0 , te 

1.3310 102

1.3020 102

1.3450 102

1.3280 102

1.3290 102

Ncall

107

21

22

22

20

 P (%)

---

1.04

0.23

0.15

f

M

Methods

2.18

ED

Note: MCS represents Monte Carlo Simulation method. SLK-Constant, SLK-Linear and SLK-Quadratic mean the single-loop Kriging surrogate model method [22] with constant, linear and quadratic regression trends respectively. PR is the proposed method, and  Pf represents the relative

PT

error of the TDFP estimated by these methods comparing with the MCS method.

AC

CE

 KL max

N call

(a) The maximum average KL divergence

ACCEPTED MANUSCRIPT

Three Surrogates

Two Surrogates

CR IP T



N call

AN US

(b) The global accuracy Fig. 1 The results of the maximum average KL divergence and the global accuracy

M

Pf  t0 , te 

ED

PfMCS  t0 , te 

PT

N call

AC

CE

Fig. 2 The convergence history of the TDFP of the proposed method

Time-dependent limit state Time=5

Time=0

Fig. 3 The results of the multiple surrogates

ACCEPTED MANUSCRIPT

5.2 The lower extremity exoskeleton The bionic lower extremity exoskeleton is a wearable robotic derive by integrating human brain and mechanical force to increase speed, strength and endurance [38-39]. The wearer and derive under the uncertainty should be consistency so to achieve the comfort and safety of the wearer. Fig. 4 shows the simplified configuration of a powered anthropomorphic leg for the lower extremity exoskeleton, which involves the hip joint O1 , the knee joint O2 , the ankle joint O3 and the corresponding linear hydraulic actuators Lh , Lk , La .  h  t  ,  k  t  and  a  t  mean the motion angle of the hip,

6

  t    ai sin  bi t  ci  i 1

CR IP T

knee and ankle joint respectively. The ideal values of the motion of joints are from clinical gait analysis data and expressed as follows: (30)

where the coefficients ai , bi and ci for the three joints are showed in Table 3. The three joints configurations are similar, thus the hip joint is taken for a general illustration, which is showed in Fig. 5.

AN US

T Considering the errors from manufacturing and assembly, the practical motion angle  h  t  for the

hip joint is showed below:

 h4  J 2   h1  J1     arctan   h2  O1   h3 

 hT  t     arctan 

  2 2 2 2 2   h1  J1    h2  O1    h4  J 2   h3  Lh  t    arccos    2  h1  J1 2   h2  O1 2    h4  J 2 2  h32       

(31)

M

in which hi (i  1, 2,3, 4) are the length of the links, J1 , J 2 and O1 mean the clearance of the joint, Lh  t  is the length of the hydraulic cylinder. hi (i  1, 2,3, 4) , J1 , J 2 and O1 are input

ED

normal random variables, and Lh  t  is the input Gaussian process. The distribution parameters of these inputs are showed in Table 4. The mean value Lh  t  of Lh  t  is represented as follows:

PT

 4 Lh  t    hi 2  2  i 1 

h

2 1

 

 h22  h32  h42

1



2   h  h    cos arctan  1   arctan  4   h  t    h2   h3     

(32)

CE

where hi (i  1, 2,3, 4) are the mean values of hi (i  1, 2,3, 4) , and  h  t  is the ideal value of motion

AC

of the hip joint. Therefore, the performance function for the hip joint of the lower extremity exoskeleton can be expressed as follows:

Gh   h  hT  t   h  t 

(33)

where the allowable motion angle is set to be  h  8 by analyzing the practical case, and let us h denote the TDFP for the hip joint as Pf  t0 , te  . The corresponding distribution parameters of the

inputs for the knee and ankle joint are provided in Table 5 and Table 6 respectively. Let us set the allowable motion angles for the knee and ankle joint be  k  6 and  a  9

respectively, then the

k a TDFPs for the knee and ankle joint can be estimated as Pf  t0 , te  and Pf  t0 , te  respectively by the

same procedure. The TDFP Pf  t0 , te  of the lower extremity exoskeleton can be estimated by

Pf  t0 , te   1  1  Pfh  t0 , te   1  Pfk t0 , te   1  Pfa t0 , te 

(34)

ACCEPTED MANUSCRIPT The considered time interval is t0 , te    0,10 .

O1

Lh

h Lk

O2

CR IP T

k

La O3

a

Fig. 4 The simplified configuration of a powered anthropomorphic leg

Lh  t 

h1

AN US

J1

O1

h t 

M

J2

h4

h2

h3

O2

ED

Fig. 5 The configuration for the hip joint

Table 3 The solutions of the coefficients ai , bi and ci for the three joints

a3

a4

a5

a6

27.46

34.79

3.309

2.004

6.395

6.136

33.35

19.05

23.4

12.81

34.08

34.34

a t 

42.54

24.41

5.12

8.571

1.708

0.5637

Coefficients

b1

b2

b3

b4

b5

b6

h t 

0.03314

0.06207

0.1739

0.2142

0.2861

0.2906

k t 

0.02532

0.1541

0.07117

0.1776

0.2856

0.2835

a t 

0.01413

0.03969

0.2125

0.1322

0.3177

0.4185

Coefficients

c1

c2

c3

c4

c5

c6

h t 

-0.1087

1.713

-4.579

-2.774

-2.98

-0.006893

k t 

0.2035

-2.459

2.255

-0.1569

0.6679

-2.384

a t 

4.982

0.9175

-2.604

-3.661

-2.531

-8.677

h t 

AC

CE

k t 

a1

PT

a2

Coefficients

ACCEPTED MANUSCRIPT Table 4 The distribution parameters of the inputs for the hip joint Distribution

Mean (mm)

Standard deviation (mm)

h1

Normal

65.17

0.5

h2

Normal

55.08

0.5

h3

Normal

400.59

0.5

h4

Normal

56.85

0.5

J1

Normal

0.2

0.1

J 2

Normal

0.2

0.1

O1

Normal

0.2

0.1

Lh  t 

Gaussian process

Lh  t 

CR IP T

Inputs

3.5

Table 5 The distribution parameters of the inputs for the knee joint Distribution

Mean (mm)

h1

Normal

45.14

h2

Normal

h3

Normal

h4

Normal

J1

Normal

J 2

Normal

O2

Normal

Lk  t 

Gaussian process

M

ED

Standard deviation (mm) 0.5

AN US

Inputs

360.53

0.5

50.07

0.5

65.17

0.5

0.2

0.1

0.2

0.1

0.2

0.1

Lk  t 

3

Table 6 The distribution parameters of the inputs for the ankle joint Mean (mm)

Standard deviation (mm)

Normal

75.23

0.5

Normal

363.87

0.5

Normal

40.06

0.5

h4

Normal

85.27

0.5

J1

Normal

0.2

0.1

J 2

Normal

0.2

0.1

O3

Normal

0.2

0.1

La  t 

Gaussian process

La  t 

3

h1 h2

AC

CE

h3

Distribution

PT

Inputs

The number of the sample base for the input variables is set to be N MCS  5 105 and the time parameter is discretized into Nt  100 time instants for this example. The solutions of the TDFP estimated by different methods are showed in Table 7. The results of the maximum average KL divergence (or minimal U-learning function) and the global accuracy for the hip joint are provided in Fig. 6. The convergence history of the TDFP of the proposed method for hip joint is showed in Fig. 7.

ACCEPTED MANUSCRIPT From Table 7, one can see that the TDFP solutions of the hip joint, knee joint and ankle joint estimated by the proposed method can match well with the corresponding solutions of the MCS method, which illustrate the accuracy of the proposed method. The TDFP Pf  t0 , te  of the lower extremity exoskeleton is obtained by Eq. (34) based on the TDFP solutions of the hip joint, knee joint and ankle joint. Fig. 6 shows the adaptive process of constructing multiple surrogates of the proposed method as the increasing of the training sample number. It can been drawn that the adaptive process can be classified into four stages. In the first stage, two surrogates involving the Kriging models with constant and linear regression trends are adaptively constructed because that the adding sample condition (b) in subsection 3.1 is set up, i.e., 7  1  2  12 

 7  1  2 7  1  3

CR IP T

. Because the Kriging model with 2 linear regression trend satisfies   0 in two consecutive iterations, which means this model is

AN US

inaccurate and removed from the ensemble of models. In the second stage, only the Kriging model with constant regression trend is reserved to construct the surrogate model of performance function. In this stage, the U-learning function in the single-loop Kriging method [22] is used to identify the new training samples, and Fig. 6(a) lists the variation of minimum U-learning function for this stage. As the increase of the size of the training samples, the Kriging model with quadratic regression trend is added in the ensemble of models and the adaptive procedure goes to the third stage. Because the same reason in the first stage, the Kriging model with quadratic regression trend is removed from the ensemble of models and only the surrogate with constant regression trend is reserved in the fourth stage. At the end of the adaptive procedure, one can see that only the Kriging model with constant regression trend exists h in the ensemble of models, thus the TDFP Pf  t0 , te  of the hip joint can be obtained by this model.

AC

CE

PT

ED

M

Fig. 7 displays that the TDFP converges to the convergence solution as the training sample number increasing, which demonstrates the good convergence of the proposed method. It also can been found from Fig 7 that there is a cusp point at the moment that the Kriging model with quadratic regression trend is added in the ensemble of models. Because the Kriging model with quadratic regression trend is not the suitable surrogate for this example, thus the adding of this surrogate will cause large error. That’s why there is a cusp point in the Fig. 7. Actually, this model is quickly deleted from the ensemble of models by our proposed strategy. Because the performance functions of the hip joint, knee joint and ankle joint are similar to each other, thus the adaptive procedures of the knee joint and ankle joint are similar to that of the hip joint. It should be noted that the TDFP solutions estimated by the proposed method and the single-loop Kriging method with constant regression trend have the similar estimation accuracy, which are higher than that of the single-loop Kriging method with linear and quadratic regression trends. This is because that the finial surrogate of the proposed method only has the Kriging model with constant regression trend and this model is the most suitable model in these three models. However, there is general no prior information for identifying a regression trend in the single-loop Kriging method and the regression trend is always man-made. Therefore, we cannot guarantee that the man-made regression trend in the single-loop Kriging method is always the best for constructing the surrogate model of the performance function. Comparing with the single-loop Kriging method, the proposed method can adaptively identify the suitable regression trends so to construct the multiple-Kriging-surrogate, which can avoid man-made subjectivity for regression trend in single-loop Kriging method. Furthermore, the computational cost of the proposed method is different from that of the single-loop Kriging method with constant regression trend, which is mostly caused by the different stopping criterions in these two methods.

ACCEPTED MANUSCRIPT Table 7 The solutions of the TDFP by different methods SLK-Constant

SLK-Linear

SLK-Quadratic

PR

P  t0 , te 

1.2120 103

1.2180 103

1.2260 103

1.2340 103

1.2160 103

Ncall

5 107

90

120

116

86

 P (%)

---

0.50

1.16

1.82

0.33

Pfk  t0 , te 

1.9048 102

1.9100 102

2.0142 102

1.9114 102

1.9090 102

Ncall

5 107

77

79

82

80

 P (%)

---

0.27

5.74

Pfa  t0 , te 

1.8320 103

1.8300 103

Ncall

5 107

99

125

 P (%)

---

0.11

0.44

Pf  t0 , te 

2.2032 102

2.2088 102

Ncall

1.5 108

266

 P (%)

---

0.25

Hip

f

Knee

f

Ankle

f

System

f

 KL max

0.35

0.22

1.8520 103

1.8340 103

103

98

1.09

0.11

U min

2.3144 102

2.2139 102

2.2080 102

324

301

264

5.05

0.49

0.22

 KL max

U min

AC

CE

PT

ED

 KL max / U min

1.8400 103

CR IP T

MCS

h f

AN US

Methods

M

Parts

N call

(a) The maximum average KL divergence (or minimal U-learning function)

ACCEPTED MANUSCRIPT

 Two Surrogates

Two Surrogates

One Surrogates

CR IP T

One Surrogates

N call

AN US

(b) The global accuracy Fig. 6 The results of the maximum average KL divergence (or minimal U-learning function) and the global accuracy for hip joint

PT

ED

M

Pf  t0 , te 

PfMCS  t0 , te 

N call

CE

Fig. 7 The convergence history of the TDFP of the proposed method for hip joint

AC

5.3 A beam under stochastic loads A beam under stochastic loads [40] showed in Fig. 8 is employed to illustrate the effectiveness of the

proposed method. The concentrated load F (t ) acting at the midpoint of the beam and the load q (t ) distributed on the top surface of the beam uniformly. The cross-section A-A is rectangular with the width a0 and height b0 . The limit state function for this beam is as follows:

g  X, Y(t ), t  

 F (t ) L q(t ) L2  st a0b0 L2 1 (a0  2kt )(b0  2kt )2  u      4 4 8 8 

  

(35)

3 where  st  78.5kN / m is the density, L is the length of the beam,  u represents the ultimate

strength and k  1  104 m / yr is the decreasing rate of the length and the width of the beam. The

ACCEPTED MANUSCRIPT

time parameter t is considered in

 0, 5

years. The stochastic processes F (t ) and q(t ) are

supposed to be Gaussian processes, and the autocorrelation coefficients are showed below:

 F (t , )  exp (  t ) 2   4(  t ) 2   9  

(36)

 q (t , )  exp 

The distribution parameters of these input variables are showed in Table 8. L

CR IP T

AA

F (t )

L 2 A

kt

kt

q (t )

b0

a0

A

AN US

Fig. 8 Beam under stochastic loads

Table 8 The distribution parameters of the inputs Inputs

Distribution

a0 (m)

Lognormal

b0 (m)

Normal

L (m)  u (Pa)

Normal

F (t ) (N)

q (t ) (N/m)

Mean

Standard deviation

0.2

0.01

0.04

4 103

5

0.08

1.5 10

2.4 105

Gaussian process

3000

400

Gaussian process

450

50

ED

M

Normal

8

The number of the sample base for the input variables is set to be N MCS  105 and the time

PT

parameter is discretized into Nt  100 time instants. The TDFP solutions estimated by different

AC

CE

methods are showed in Table 9. The results of the maximum average KL divergence and the global accuracy are provided in Fig. 9. The convergence history of the TDFP of the proposed method is showed in Fig. 10. From Table 9, one can see that the TDFP estimated by the proposed method matches well with that of the MCS method. Comparing these methods with each other, it is easy to know that the proposed multiple-Kriging-surrogate method is the most efficient method for estimating TDFP. Furthermore, when comparing the TDFP solutions estimated by the single-loop Kriging method with different regression trends, we can draw that the TDFP estimated by the surrogate with the constant regression trend is more accurate than that of the surrogate with linear and quadratic regression trends. Therefore, we cannot say that the surrogate with complicated regression trend (i.e., quadratic regression) in single-loop Kriging method can give more accurate estimation result, and we should choose the model with suitable regression trend for different problem. It can been seen from Fig. 9 that the proposed multiple-Kriging-surrogate method adaptively identified two regression trends involving both constant and linear to construct the final surrogate of the performance function, and the solution illustrates that these two regression trends are the suitable models for constructing the surrogates. Fig. 10 shows that the proposed method has good convergence in estimating the TDFP.

ACCEPTED MANUSCRIPT Table 9 The solutions of the TDFP by different methods Methods

MCS

SLK-Constant

SLK-Linear

SLK-Quadratic

PR

Pf  t0 , te 

2.8260 102

2.7840 102

2.7760 102

2.7520 102

2.8090 102

Ncall

107

33

33

34

29

 P (%)

---

1.49

1.77

2.62

0.60

f

AN US

CR IP T

 KL max

N call

(a) The maximum average KL divergence

CE

PT

ED

M



N call

AC

(b) The global accuracy Fig. 9 The results of the maximum average KL divergence and the global accuracy

ACCEPTED MANUSCRIPT

Pf  t0 , te 

CR IP T

PfMCS  t0 , te 

N call

Fig. 10 The convergence history of the TDFP of the proposed method

AN US

5.4 A roof truss structure

A roof truss [41, 42] showed in Fig. 11 is modified to show the effectiveness of the proposed method for estimating TDFP. The top boom and the compression bars of this roof truss are made of concrete, and the bottom boom and the tension bars are reinforced by steel. The uniformly distributed load q  t  is applied on the roof truss structure, and the load q  t  can be transformed into nodal load

q t  l 4

, where l is the length of bottom boom. The limit state function of this structure is showed

M

P

below:

ED

g  X, Y(t )   0.03 

q  t  l 2  3.81 1.13     2  AC EC AS ES 

(37)

which means that the perpendicular deflection of node C exceeds 0.03m corresponds to the structure

PT

failure state. AC , AS , EC , ES represent sectional area and elastic modulus respectively. The distribution parameters of these variables are showed in Table 10, where the load q  t  is considered to be the stochastic process with the following autocorrelation coefficient. (38)

C

F

D

A

B

q t  N / m

E

G

P

Ac

D

0.278l

C

5 A c As 0.7

A 3 As

Ac

P

Ac As

F 0.7 5A c

2 As

E

0.222l

G

0.25l

l 12

P

Ac

B

l 12

AC

CE

q (t , )  exp (  t )2 

3 As

0.25l

Fig. 11 The schematic diagram of a roof truss

ACCEPTED MANUSCRIPT Table 10 The distribution parameters of the inputs Inputs

Distribution

Mean

Standard deviation

l (m)

Normal

12

0.12

AS ( m )

Normal

9.82 104

5.892 105

ES ( Pa )

Normal

11011

6 109

AC ( m2 )

Normal

4 102

4.8 103

EC ( Pa )

Normal

2 1010

1.2 109

q (t ) (N/m)

Gaussian process

2 104

1.4 103

CR IP T

2

M

AN US

Table 11 shows the TDFP solutions estimated by these methods. It can been seen that the SLK-Constant, SLK-Linear, SLK-Quadratic and the proposed method all can give acceptable TDFP solutions comparing with the MCS method. However, when comparing the proposed method with the single-loop Kriging method with different regression trends, one can see that the proposed method has the most accuracy and efficiency. From Fig. 12, it is easy to know that there are two surrogates involving the Kriging models with linear and quadratic regression trends exist in the ensemble of models after the adaptive process, and the TDFP is estimated by the most accurate local model based on the two models. The reason of why the proposed method is more accurate and efficient than the single-loop Kriging method is that the proposed method can adaptively identify the most suitable models, i.e., Kriging models with linear and quadratic regression trends, to approximate the time-dependent limit state of original performance function. Therefore, fewer training samples are needed and higher accuracy is obtained than the single-loop Kriging method. Fig. 13 shows the convergence history of TDFP of the proposed method, one can see that the TDFP converges to the accurate solution as the training sample number increasing, which illustrates the good convergence of the proposed method.

ED

Table 11 The solutions of the TDFP by different methods MCS

SLK-Constant

SLK-Linear

SLK-Quadratic

PR

Pf  t0 , te 

1.7980 102

1.7600 102

1.7800 102

1.8180 102

1.8060 102

Ncall

107

76

67

72

53

---

2.11

1.00

1.11

0.44

 P (%)

CE

f

PT

Methods

AC

 KL max / U min

 KL max

U min

 KL max

N call

(a) The maximum average KL divergence (or minimal U-learning function)

ACCEPTED MANUSCRIPT

One Surrogate

Two Surrogates

CR IP T



Two Surrogates

N call

AN US

(b) The global accuracy Fig. 12 The results of the maximum average KL divergence (or minimal U-learning function) and the global accuracy

M

Pf  t0 , te 

PT

ED

PfMCS  t0 , te 

N call

CE

Fig. 13 The convergence history of the TDFP of the proposed method

5.5 A stone arch bridge under hurricane load

AC

Hurricane is a catastrophic natural disaster, which may cause serious human damages and huge economic losses. Considering a stone arch bridge under hurricane load example. The finite element model of a stone arch bridge showed in Fig. 14 is built in ABAQUS 6.14-1. The width of this stone arch bridge is L  6m . This bridge is symmetrical structure and the geometric size of the profile is showed in Fig. 15. A hurricane load F  t  uniformly loads on the profile of this bridge along the negative direction of z-axis. The hurricane load F  t  is expressed as follows: t  0.5s   F0 F t      F0  0.4 F0 cos  2  t  0.5   0.6 F0 cos  4  t  0.5  t  0.5s

(39)

which means that the bridge sustains equally load F0 before t  0.5s and bears time-variant load

ACCEPTED MANUSCRIPT after t  0.5s . F0 is considered to be normal distributed variable with mean value 1.5 103 Pa and standard deviation 75Pa. The hurricane load F  t  varies with the time parameter when F0 is fixed at the mean value is showed in Fig. 16. It can be found that the middle part of the bridge deck general has the maximum deformation under the hurricane load. Thus the displacement  Z of node P (showed in Fig. 14) along the direction of z-axis is employed to construct the limit state function of this bridge. Considering the absolute displacement  Z

exceeds 0.011m corresponding to the failure state

of this stone bridge, thus the limit state function can be expressed as follows: g  X, Y(t )   0.011   Z

(40)

CR IP T

The considered time interval is t   0, 5 s .The input random variables involve the density  , the Elastic modulus E , the Poisson’s ratio  , the Rayleigh damping coefficient  and the load F0 .

AN US

The distribution parameters of these input variables are showed in Table 12. When all these inputs are fixed at their mean values, the first order model of the bridge and the deformation of the bridge at time instant t  5 are obtained and showed in Fig. 17 and Fig. 18 respectively, and the corresponding displacement of node P varies with the time parameter is displayed in Fig. 19.

M

P

PT

ED

F t 

AC

CE

Fig. 14 The Finite element model of the stone arch bridge

1m

3.5m 1.75m

1m

9m

0.5m

Fig. 15 The geometric size of the stone arch bridge

1m

0.5m

9m

1.5m

ACCEPTED MANUSCRIPT

CR IP T

F t 

t

ED

M

AN US

Fig. 16 The hurricane load at the mean value of F0

AC

CE

PT

Fig. 17 The first order model of the stone arch bridge

Fig. 18 The deformation of the stone arch bridge at time instant t  5

ACCEPTED MANUSCRIPT

CR IP T

Z

t

Fig. 19 The displacement of the node P at the inputs mean values Table 12 The distribution parameters of the inputs Distribution

Mean

 ( Kg / m )

Normal

2.43 103

E (Pa)

Normal

 

Normal Normal

F0 (Pa)

Normal

3

AN US

Inputs

Standard deviation

1.215 102

5.23 107 0.3 2

2.615 106 0.015 0.1

1.5 103

75

M

The sample size of MCS method is set to be 5 103 , and the displacement  Z of node P varies with the time parameter at each input random sample can be obtained by once call of finite element

ED

model. Thus the total computational cost of the MCS method is 5 103 . For the single-loop Kriging method and the proposed method, we set the number of the sample base to be 1104 , and the time parameter is discretized into Nt  100 time instants.

AC

CE

PT

Table 13 shows the TDFP solutions estimated by these methods. One can see that the TDFP solution estimated by the proposed method can match well with that of the MCS method, which illustrates the accuracy of the proposed method. Comparing the proposed method with the single-loop Kriging method with different regression trends, we can draw that the proposed method has the highest accuracy among these methods while by the least computational cost. Simultaneously, the single-loop Kriging method with constant regression trend has large error in this example. Fig. 20 illustrates that there are two models involving the Kriging models with linear and quadratic regression trends exist in the final ensemble of models, and the TDFP is estimated by the most accurate local model which is determined by these two models. With the suitable models, the time-dependent limit state can be better identified than only fixed the model. That’s why the proposed method is more accurate and efficient than the single-loop Kriging method. Fig. 21 shows the convergence history of TDFP of the proposed method. It can be found that the TDFP converges to a determined solution as the training sample size increasing. Table 14 shows the comparison of the computational times of the single-loop Kriging methods with different regression trends and the proposed method. It can been seen that the proposed method is the most efficient method among these methods. This is because that the proposed method has the least numbers of calling the finite element model. Calling the finite element model is general time-consuming in engineering application. Simultaneously, the single-loop Kriging method with the

ACCEPTED MANUSCRIPT quadratic regression trend needs the most computational time among these methods because it needs the most numbers of calling the finite element model. Table 13 The solutions of the TDFP by different methods MCS

SLK-Constant

SLK-Linear

SLK-Quadratic

PR

Pf  t0 , te 

2.9800 102

3.2400 102

3.1200 102

2.8600 102

2.9400 102

Ncall

5 103

81

85

89

66

 P (%)

---

8.72

4.70

4.03

1.34

f

CR IP T

Methods

Table 14 The comparison of computational time Methods

Time (s)

SLK-Constant SLK-Linear SLK-Quadratic PR

48199 53912 60222 43928

 KL max / U min

ED

M

 KL max

AN US

Note: Computations are conducted on the computer with Intel(R) Xeon(R) CPU W3550 @ 3.07GHz 3.06GHz, 12.0GB of RAM and 64bit Windows 7 system.

 KL max N call

PT

U min

AC

CE

(a) The maximum average KL divergence (or minimal U-learning function) 

Two Surrogates

One Surrogate

Two Surrogates

N call

(b) The global accuracy Fig. 20 The results of the maximum average KL divergence (or minimal U-learning function) and the global accuracy

ACCEPTED MANUSCRIPT

Pf  t0 , te 

CR IP T

PfMCS  t0 , te 

N call

Fig. 21 The convergence history of the TDFP of the proposed method Conclusions

AN US

6.

A multiple-Kriging-surrogate method is established to estimate the TDFP for the time-dependent structure. It is well known that the single-loop Kriging method constructs the surrogate model of performance function with a fixed regression trend. However, there is general no prior information for identifying a suitable regression trend and it is always man-made. Therefore, it is hard to say that the

M

man-made regression trend is the best for constructing the surrogate model of the performance function. For dealing with this issue, we proposed the multiple-Kriging-surrogate method which can adaptively identify the suitable regression trend to construct multiple-Kriging-surrogate, thus the man-made

ED

subjectivity for regression trend can be avoided in single-loop Kriging method. Generally, two cases of surrogate solutions may be drawn in the final ensemble of models by the proposed method, and they are shown as follows.

The first case is that there exists only one model with a regression trend in the final ensemble of

PT

1.

models just like the example 5.2. This case illustrates that the performance function can be

CE

approximated well by only a regression trend, i.e., constant regression, linear regression or quadratic regression. The single-loop Kriging method cannot guarantee to select the most suitable regression trend before constructing surrogate model. But the proposed method can

AC

adaptively identify the most suitable regression trend to construct the final surrogate model. In this case, the accuracy and efficiency of the proposed method are similar to those of the single-loop Kriging method with the best regression trend.

2.

The second case is that there exist multiple models with different regression trends in the final ensemble of models just like the example 5.1, 5.3, 5.4 and 5.5. This case shows that the performance function can be approximated more accurately and efficiently with multiple models rather than one model. Fortunately, the multiple suitable models can be adaptively identified by the established procedure in our contribution. In this case, the accuracy and efficiency of the proposed method are both superior than those of the single-loop Kriging method with only one regression trend.

ACCEPTED MANUSCRIPT Comparing with the single-loop Kriging method, the proposed method needs to construct multiple surrogate models, thus it needs more numbers of calling the surrogate models to identify the new training samples in the adaptive process, which may be a limitation of the proposed method. Perhaps several optimization techniques can be combined to determine the new training samples to avoid repeatedly calling the surrogate models. Studying several approaches by combining these techniques will be our future work. Furthermore, we will try to extend the proposed method to the sensitivity analysis [43, 44] of structure and the structure with epistemic uncertainty [45, 46].

CR IP T

Acknowledgement This work was supported by the Natural Science Foundation of China (Grants 51475370 and 51775439).

Appendix A: The OSE method

AN US

With the OSE method, the orthogonal function hi (t )(i  1, 2,..., ) is expressed as follows:

2i  1 Lei (tL ) 2(te  t0 )

hi (t ) 

(A-1)

where Lei () is the i -th Legendre polynomial and t L is given by

tL  2

t 1 te  t0

(A-2)

M

Generally, the stochastic process Y j (t )( j  1, 2,..., m) can be expressed as follows: 

Y j (t )  Y j (t )    ji hi (t )

(A-3)

ED

i 1

where  ji (i  1, 2,..., ) is a set of correlated zero-mean normal random variables. With the theory of

PT

the OSE, the set of random variables  ji (i  1, 2,..., ) should be transformed to be uncorrelated standard normal variables. Based on the orthogonal property of the basic orthogonal function

CE

hi (t )(i  1, 2,..., ) , the covariance   ji  jk between  ji and  jk can be estimated as follows:

 ji  jk  

te

t0



te

t0

CYj (t , )hi (t )hk ( )dtd

(A-4)

AC

  ji  jk can be expressed as follows by the spectral decomposition.   ji  jk  j   j  j

(A-5)

in which  j is a diagonal matrix with the diagonal elements  jk being the eigenvalues of   ji  jk and  j is a matrix with each column  ji being the corresponding eigenvectors of   ji  jk . The random variable  ji (i  1, 2,..., ) can be expressed by the independent standard normal vector  j as follows: 

 ji    kji  jk  jk

(A-6)

k 1

k where  ji (i  1, 2,..., ) means the coordinate of the k -th eigenvector. Then the stochastic process

Y j (t )( j  1, 2,..., m) can be simulated by a set of independent normal variable by substituting Eq. (A-6)

ACCEPTED MANUSCRIPT into Eq. (A-3). 



i 1

k 1

Y j (t )  Y j (t )   (  kji  jk  jk )hi (t ) 

(A-7)



 Y j (t )    jk  jk   h (t ) k 1

i 1

k ji i

It is unrealistic to estimate the infinite eigenvalues and eigenvectors of   ji  jk . Given the

M j ( j  1, 2,..., m) terms for the approximation, the stochastic process Y j (t )( j  1, 2,..., m) can be

Mj

Mj

k 1

i 1

Y j (t )  Y j (t )    jk  jk   kji hi (t ) Mj

 Y j (t )    jk  jk (t ) jk k 1

Mj

in which  jk (t )    kji hi (t ) .

AN US

i 1

CR IP T

expressed as follows:

References

(A-8)

AC

CE

PT

ED

M

[1] SI Ammar, YF Alharbi. Time-dependent analysis for a two-processor heterogeneous system with time-varying arrival and service rates. Applied Mathematical Modelling, 2018, 54: 743-751. [2] Y Shi, ZZ Lu, SY Chen, LY Xu. A reliability analysis method based on analytical expressions of the first four moments of the surrogate model of the performance function. Mechanical Systems and Signal Processing, 2018, 111: 47-67. [3] C Andrieu-Renaud, B Sudret, M Lemaire. The PHI2 method: a way to compute time-variant reliability. Reliability Engineering and System Safety, 2004, 84: 75-86. [4] B Sudret. Analytical derivation of the outcrossing rate in time-variant reliability problems. Structure and Infrastructure Engineering, 2008, 4(5): 353-362. [5] JF Zhang, XP Du. Time-dependent reliability analysis for function generator mechanisms. Journal of Mechanical Design, 2011, 133: 031005. [6] C Jiang, XP Wei, ZL Huang. An outcrossing rate model and its efficient calculation for time-dependent system reliability analysis. Journal of Mechanical Design, 2017, 139: 041402. [7] Z Hu, XP Du. Time-dependent reliability analysis with joint upcrossing rates. Structural and Multidisciplinary Optimization, 2013, 48: 893-907. [8] Y Shi, ZZ Lu, KC Zhang. Reliability analysis for structure with multiple temporal and spatial parameters based on the effective first-crossing point. Journal of Mechanical Design, 2017, 139(12): 121403. [9] AT Beck, RE Melchers. On the ensemble crossing rate approach to time variant reliability analysis of uncertain structures. Probabilistic Engineering Mechanics, 2004, 19: 9-19. [10] J Li, JB Chen, WL Fan. The equivalent extreme-value event and evaluation of the structural system reliability. Structural Safety, 2007, 29: 112-131. [11] JB Chen, J Li. The extreme value distribution and dynamic reliability analysis of nonlinear structures with uncertain parameters. Structural Safety, 2007, 29: 77-93. [12] XP Du. Time-dependent mechanism reliability analysis with envelope functions and first-order approximation. Journal of Mechanical Design, 2014, 136: 081010. [13] Y Shi, ZZ Lu, YC Zhou. Temporal and spatial reliability and global sensitivity analysis with envelope functions. Journal of Northwestern Polytechnical University, 2017, 35(4): 591-598.

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

[14] Z Hu, XP Du. A Sampling approach to extreme value distribution for time-dependent reliability analysis. Journal of Mechanical Design, 2013, 135: 071003. [15] Y Shi, ZZ Lu, K Cheng. Temporal and spatial multi-parameter dynamic reliability and global reliability sensitivity analysis based on the extreme value moments. Structural and Multidisciplinary Optimization, 2017, 56(1): 117-129. [16] DQ Zhang, X Han, C Jiang, J Liu, Q Li. Time-dependent reliability analysis through response surface method. Journal of Mechanical Design, 2017, 139: 041404. [17] ZQ Wang, PF Wang. A new approach analysis with time-variant performance characteristics. Reliability Engineering and System Safety, 2013, 115: 70-81. [18] ZQ Wang, PF Wang. A nested extreme response surface approach for time-dependent reliability-based design optimization. Journal of Mechanical Design, 2012, 134: 121007. [19] Z Hu, XP Du. Mixed efficient global optimization for time-dependent reliability analysis. Journal of Mechanical Design, 2015, 137: 051401. [20] ZQ Wang, W Chen. Confidence-based adaptive extreme response surface for time-variant reliability analysis under random excitation. Structural Safety, 2017, 64: 76-86. [21] ZQ Wang, W Chen. Time-variant reliability assessment through equivalent stochastic process transformation. Reliability Engineering and System Safety, 2016, 152: 166-175. [22] Z Hu, S Mahadevan. A single-loop kriging surrogate modeling for time-dependent reliability analysis. Journal of Mechanical Design, 2016, 138: 061406. [23] L Hawchar, CPE Soueidy, F Schoefs. Principal component analysis and polynomial chaos expansion for time-variant reliability problems. Reliability Engineering and System Safety, 2017, 167: 406-416. [24] D Drignei, I Baseski, ZP Mourelatos, E Kosova. A random process metamodel approach for time-dependent reliability. Journal of Mechanical Design, 2016, 138: 011403. [25] HP Lal, S Sarkar, S Gupta. Stochastic model order reduction in randomly parametered linear dynamical systems. Applied Mathematical Modelling, 2017, 51: 744-763. [26] HB Zhao, SJ Li, ZL Ru. Adaptive reliability analysis based on a support vector machine and its application to rock engineering. Applied Mathematical Modelling, 2017, 44: 508-522. [27] N Pedroni, E Zio. An adaptive metamoel-based subset importance sampling approach for the assessment of the functional failure probability of a thermal-hydraulic passive system. Applied Mathematical Modelling, 2017, 48: 269-288. [28] WY Yun, ZZ Lu, X Jiang. AK-SYSi: An improved adaptive Kriging model for system reliability analysis with multiple failure modes by a refined U learning function. Structural and Multidisciplinary Optimization. In press. [29] YC Zhou, ZZ Lu, K Cheng, WY Yun. A Bayesian Monte Carlo-based method for efficient computation of global sensitivity indices. Mechanical Systems and Signal Processing, 2019, 117(15), 498-516. [30] E Acar. Various approaches for constructing an ensemble of metamodels using local measures. Structural and Multidisciplinary Optimization, 2010, 42: 879-896. [31] T Goel, RT Haftka, W Shyy, NV Queipo. Ensemble of surrogates. Structural and Multidisciplinary Optimization, 2007, 33: 199-216. [32] J Zhang, S Chowdhury, A Messac. An adaptive hybrid surrogate model. Structural and Multidisciplinary Optimization, 2012, 46: 223-238. [33] J Vahedi, MR Ghasemi, M Miri. An adaptive divergence-based method for structural reliability analysis via multiple kriging models. Applied Mathematical Modelling, 2018, 63: 542-561. [34] B Settles. Active learning literature survey. University of Wisconsin at, Madison, 2010.

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

[35] B Sudret, A Der Kiureghian. Stochastic finite element methods and reliability: a state-of-the-art report. Department of Civil and Environmental Engineering, University of California, Berkeley, 2000. [36] SN Lophaven, HB Nielsen, J Søndergaard. DACE-A MATLAB Kriging toolbox. Technical University of Denmark, 2002. [37] R Schobi, B Sudret, S Marelli. Rare event estimation using polynomial-chaos Kriging. ASCE-ASME Journal of Risk and Uncertainty in Engineering Systems, Part A: Civil Engineering, 2017, 3(2): D4016002. [38] S Yu, ZL Wang, KW Zhang. Sequential time-dependent reliability analysis for the lower extremity exoskeleton under uncertainty. Reliability Engineering and System Safety, 2018, 170: 45-52. [39] S Yu, ZL Wang. A novel time-variant reliability analysis method based on failure processes decomposition for dynamic uncertain structures. Journal of Mechanical Design, 2018, 140: 051401. [40] Z Hu, XP Du. First order reliability method for time-variant problems using series expansions. Structural and Multidisciplinary Optimization, 2015, 51: 1-21. [41] QJ Pan, D Dias. Sliced inverse regression-based sparse polynomial chaos expansions for reliability analysis in high dimensions. Reliability Engineering and System Safety, 2017, 167: 484-493. [42] Y Shi, ZZ Lu, YC Zhou. Time-dependent safety and sensitivity analysis for structure involving both random and fuzzy inputs. Structural and Multidisciplinary Optimization, 2018, 58(6): 2655-2675. [43] WY Yun, ZZ Lu, X Jiang. An efficient method for moment-independent global sensitivity analysis by dimensional reduction technique and principle of maximum entropy. Reliability Engineering and System Safety. In press. [44] LY Xu, ZZ Lu, SN Xiao. Generalized sensitivity indices based on vector projection with multivariate outputs. Applied Mathematical Modelling, 2019, 66, 592-610. [45] CQ Fan, ZZ Lu, Y Shi. Time-dependent failure possibility analysis under consideration of fuzzy uncertainty. Fuzzy Sets and Systems. In press. [46] KX Feng, ZZ Lu, WY Yun. Aircraft icing severity analysis with hybrid parameters under considering epistemic uncertainty. AIAA journal. In press.