EEK-SYS: System reliability analysis through estimation error-guided adaptive Kriging approximation of multiple limit state surfaces

Journal Pre-proof

EEK-SYS: System reliability analysis through estimation error-guided adaptive Kriging approximation of multiple limit state surfaces Chen Jiang , Haobo Qiu , Liang Gao , Dapeng Wang , Zan Yang , Liming Chen PII: DOI: Reference:

S0951-8320(19)30772-0 https://doi.org/10.1016/j.ress.2020.106906 RESS 106906

To appear in:

Reliability Engineering and System Safety

Received date: Revised date: Accepted date:

15 June 2019 24 February 2020 25 February 2020

Please cite this article as: Chen Jiang , Haobo Qiu , Liang Gao , Dapeng Wang , Zan Yang , Liming Chen , EEK-SYS: System reliability analysis through estimation error-guided adaptive Kriging approximation of multiple limit state surfaces, Reliability Engineering and System Safety (2020), doi: https://doi.org/10.1016/j.ress.2020.106906

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. 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. © 2020 Published by Elsevier Ltd.

Highlights 

Estimation error-guided adaptive Kriging modeling method is proposed.



The probability of wrong state classification is defined and derived.



The maximum real-time estimation error is derived for series and parallel system.



Both the probability and estimation error assist the refinement of Kriging.

1

EEK-SYS:

System

reliability

analysis

through

estimation

error-guided adaptive Kriging approximation of multiple limit state surfaces Chen Jianga, Haobo Qiu*, a, Liang Gaoa, Dapeng Wanga, Zan Yanga, Liming Chena a

State Key Laboratory of Digital Manufacturing Equipment and Technology, School of Mechanical

Science and Engineering, Huazhong University of Science and Technology, Wuhan, 430074, P. R. China *

Corresponding author. E-mail: [email protected]

Abstract In order to approximate the multiple limit state functions for different failure events, the active learning Kriging model proposed for component reliability analysis has been extended to system reliability analysis. Meanwhile, many efficient sampling strategies have been applied to reduce the high computational burden. However, these strategies meet a challenge in wasting some training points and terminating the training process inappropriately, since they do not directly relate to the estimation error of system failure probability. To address the challenge, this work proposes an estimation error-guided adaptive Kriging method. As Kriging prediction may be inaccurate before being well trained, the predicted system failure probability may deviate from the true result. To quantify this estimation error, the true number of failure points is approximated by adding the number of predicted failure points and the number of wrongly classified points. Since it is impossible to learn the exact number of wrongly classified points, its confidence interval is derived based on the probability of making wrong state classification. Subsequently, the refinement of Kriging is achieved by using the probability to identify new points and using the estimation error to determine the termination, which has been demonstrated by three different cases. Keywords: System reliability analysis; real-time estimation error; wrong state classification; adaptive Kriging; multiple failure modes

2

1 Introduction System reliability analysis (SRA) is to measure the probability that a complex system with multiple failure modes achieves the intended function by taking into account various uncertainties. It will be named as “component reliability analysis (CRA)” if only one single failure mode exists. Generally, it plays a vital role in structural reliability assessment and design optimization [1-5]. In practical engineering, performing CRA or SRA requires to repeatedly evaluate the time-consuming limit state functions (LSFs). Since it takes Nmcs  p function evaluations (p components and Nmcs simulations for each component) to estimate the failure probability by crude Monte Carlo simulation (MCS), many researches are motivated to reduce the number of LSF evaluations as much as possible. To improve the efficiency of MCS , various variance reduction techniques, such as importance sampling [6, 7], line sampling [8, 9] and subset simulation [10], were developed. However, these methods inherit the disadvantages of MCS, as they still need thousands of function evaluations. From the analytical perspective [11-17], the first/second-order reliability methods (FORM/SORM) were proposed to solve the CRA problems and extended to SRA. Although the number of LSF evaluations witnesses a remarkable decrease, the inaccuracy of FORM and inefficiency of SORM are still the unavoidable challenges especially when there are many highly nonlinear LSFs and multiple failure domains. In the last two decades, various surrogate models including response surface model [18, 19], radial basis function [20, 21], support vector regression [22], polynomial chaos expansion [23] and Kriging [24, 25] have been used to approximate the LSFs and combined with simulation methods for reliability analysis. Among them, Kriging is more popular due to its feature in providing the measurement of prediction uncertainty. As for the Kriging-based reliability analysis methods, Efficient global reliability analysis (EGRA) [26] and active learning Kriging-based Monte Carlo simulation (AK-MCS) [27] have made great progress in dealing with CRA. Following the essential ideas of these two methods, many active learning Kriging strategies, such as the LIF function [28], REAK method [29] and FPS framework [30], were proposed for CRA [31-38]. In addition, it is worth noting that adaptive Kriging has also be employed in time-variant (or time-dependent) reliability problem [39-43], which can be regarded as the series system reliability problem since it degrades into CRA in each time instant within 3

the interested time interval. Recently, Hu and Mahadevan [41] developed a single-loop Kriging surrogate modeling method (SILK) by building a global surrogate model for the performance function with the random variables, stochastic processes and time parameter as inputs. Wang and Chen [42] proposed an equivalent stochastic process transformation (eSPT) method to obtain an equivalent time-invariant reliability problem. Jiang et al. [43] proposed to measure the estimation error of the failure probability predicted by current Kriging model in each iteration. These methods achieve good illustrations for applying adaptive Kriging to time-variant reliability problem. Furthermore, the computational efficiency of active learning Kriging is improved by being coupled with variance reduction-based sampling methods [44, 45]. More details of the methods and applications of time-variant reliability can be found in [46-50]. EGRA and AK-MCS were also extended to SRA, which are named as EGRA-SYS [51] and AK-SYS [52] respectively. Both of them pay little attention to the components that have little effect on the system failure domain, and they all have three strategies. The first one is to construct the Kriging model of every LSF by EGRA or AK-MCS. The second one is to build Kriging model of the minimum response function for series system or maximum response function for parallel system. The third one establishes the Kriging model of every LSF like the first one, but it only calibrates one Kriging model that is the most important to the system failure domain at one sample point in each iteration. EGRA-SYS and AK-SYS inherit the characteristics of EGRA and AK-MCS respectively. Compared to EGRA-SYS, AK-SYS requires no optimization program and reduces the sampling region from the whole design space to the areas with a significant probability content [52]. Although both of them can handle the SRA problem, they all encounter a limitation that the false identification of minimum mode or maximum mode may occur by using the last calibrated Kriging and may further mislead the next refinement of Kriging. Therefore, an improved adaptive Kriging model for system reliability analysis (AK-SYSi) [53] was proposed to update Kriging model from the most easily identifiable failure mode among all the failure modes. Additionally, Yang et al. [54, 55] proposed the truncated candidate region for SRA to deal with the situation that large numerical difference arises among the components. The truncated candidate region can be combined with AK-SYSi to enhance the performance. Although AK-SYSi can obtain the accurate results no matter the in-process Kriging models are accurate or not, it does not directly relate the estimation error of system failure probability to the training of surrogate, leading to the limitation that some useless but computationally expensive points may be generated. 4

Moreover, since AK-SYSi uses the learning function U as the stopping criterion, it may tend to terminate prematurely or too late, which has been discussed in [30, 56, 57]. Besides these mentioned methods, more SRA methods and their applications can refer to literatures [58-62]. The main aim of this paper is to make the estimation error of the system failure probability participate in the refinement of Kriging. In this work, an estimation error-guided adaptive Kriging (EEK) modeling method is proposed to approximate the multiple limit state surfaces. Since the estimation error is caused by the wrong classification, it is necessary to quantify the number of the wrongly classified random points after estimating the system failure probability through a large number of MCS random points. To achieve this purpose, the probability of wrong state classification (WSC) is derived based on the probability of making wrong sign estimation (WSE) for every component, which can be calculated by using the local uncertainty information provided by Kriging. In each iteration, the random point with the maximum probability of WSC will be selected as the new training point, and it is used to update the Kriging model for the component LSF that is the most important to the system failure domain. Subsequently, the probability of WSC is used to derive the confidence intervals of the number of the wrongly classified points based on the probabilistic theory, such as the Lindeberg central limit theorem [63]. According to the confidence intervals, the real-time maximum estimation error can be obtained and then be utilized to judge when to stop the training process. The remainder of this paper is organized as follows. Section 2 reviews the background of system reliability analysis. Section 3 gives the details of the proposed EEK-SYS method. Section 4 uses three cases to validate and demonstrate the performance of EEK-SYS. Finally, Section 5 gives the conclusions of this paper.

2 Background of system reliability analysis 2.1 The formulation The system failure probability of a series or parallel system is respectively defined by p  Pfseries  P  gi ( X)  0 ,  i 1 

(1)

p  Pfparallel  P  gi ( X)  0 ,  i 1 

(2)

where p denotes the number of failure modes. To solve Eq. (1) and Eq. (2) by using the crude MCS, the corresponding system failure probability can be given by 5

Pfseries  E ( I FSeries ) 

1 N mcs

Pfparallel  E ( I Fparallel ) 

1 N mcs

N mcs

I j 1

Series F

j (x mcs ),

(3)

N mcs

I j 1

parallel F

j (x mcs ),

(4)

j where N mcs represents the number of random points generated by MCS, x mcs is the jth random input

point. I FSeries and I Fparallel are the failure indicator functions, they can be defined as

 p j gi (x mcs )0 1 min i 1 Series j I F (x mcs )   , p 0 min g (x j )  0  i 1 i mcs

(5)

p  j 1 max gi (x mcs )0  i 1 parallel j IF (x mcs )   . p 0 max g (x j )  0  i 1 i mcs

(6)

2.2 AK-SYS method Since MCS method is not efficient for practical applications, the AK-MCS method is extended to AK-SYS [52] for the SRA problems. AK-SYS consists of the following three different strategies. The first strategy, named as component approach, constructs the Kriging model of every LSF by using AK-MCS method. Each Kriging model is updated with the following U learning function [27] j U i (x mcs )

j  gˆ (x mcs ) i

j  gˆ (x mcs )

, i  1,..., p, j  1..., N mcs ,

(7)

i

j j ) and  gˆi (xmcs ) are the mean prediction and standard deviation respectively. The basic where gˆi (xmcs

theory of Kriging can refer to [64]. This strategy is easy to achieve, but not the most efficient one. The Second strategy, named as composite model approach, only builds an extreme response function for the system limit state surface, i.e. minimum response function for series system or maximum response function for parallel system. However, it is computationally expensive for Kriging to fit the extreme response function which may usually have sharp corners or be highly irregular around the intersections of different LSFs. As a result, many training points far away from the system limit state surface are wasted. The third one makes use of a composite criterion to actively learn the LSF that is the most j important to the system limit sate surface. For each random point x mcs , the corresponding learning

function calculated by 6

U s (x

j mcs

)

j  gˆ (x mcs ) s

j  gˆ (x mcs )

(8)

s

j where s is the index of the LSF that minimizes or maximizes gˆs (xmcs ) for series or parallel system

respectively. Subsequently, only one LSF would be updated in each iteration when the point that minimizes the U function is found. With this composite criterion, no calls will be made for true functions that have little or no influence on the failure domain.

2.3 AK-SYSi method Since AK-SYS finds the most important component based on the last calibrated Kriging that may be not accurate enough, sometimes the insignificant components may be selected and some training points will be wasted. To avoid identifying the minimum or maximum prediction, AK-SYSi [53] was proposed to update the Kriging from the most easily identifiable failure mode among all the multiple failure modes. The learning function of the series system is defined by j   gˆi (x mcs ) j  min  gˆi (x mcs )  0 i  1,..., p j i  1,..., p  ( x )  gˆ i mcs  j U (x mcs )   j  gˆi (x mcs )  j max  gˆi (x mcs )  0 i  1,..., p f i 1,..., p j f  ( x ) gˆ i mcs  

(9)

j where p f (0  p f  p) is the number of failure modes at x mcs . And the learning function of the

parallel system is defined by j   gˆi (x mcs ) j  max  gˆi (x mcs )  0 i  1,..., ps j i  1,..., p s  ( x )  gˆ i mcs  j U (x mcs )   j  gˆi (x mcs )  j min  gˆi (x mcs )  0 i  1,..., p i 1,..., p j  ( x ) gˆ i mcs  

(10)

j where ps (0  ps  p) is the number of safe modes at x mcs .

3 The proposed EEK-SYS method The key of system reliability analysis through surrogate model is to predict the sign and classify the state of random points accurately. However, as we have mentioned, the sign and state of some random points may be wrongly classified by the initial or in-process Kriging model. Therefore, it is necessary to quantify the real-time estimation error produced by the wrong classification and reduce 7

the estimation error as much as possible while constructing Kriging models.

3.1 The estimation error of system failure probability Based on the definition in Eq. (3) and Eq. (4), all the random points (totally N mcs ) are evaluated by Kriging model. According to the predicted signs, the system failure probability is calculated by

pˆ kf  Nˆ kf / Nmcs ,

(11)

where Nˆ kf is the number of predicted failure points in the kth iteration, and the number of predicted safe k points can be represented by Nˆ s . In the next k+1th iteration, the Kriging model is updated with the new k 1 added sample points, and the random points are re-evaluated. Accordingly, Nˆ kf 1 and Nˆ s are

updated. However, the evaluation of random points may be not accurate base on the initial or in-process Kriging model. Some random points that are in fact safe points may be classified into failure points due to the inaccuracy of Kriging model. Thus, we define these points classified into the predicted failure points as the wrongly classified failure points and vice versa. Based on this definition, the number of wrongly classified failure points in Nˆ kf is represented by nˆ kf (0  nˆ kf  Nˆ kf ) , while the number of k k k k wrongly classified safe points in Nˆ s is represented by nˆs (0  nˆs  Nˆ s ) .

On the contrary, the true failure probability can be directly estimated by MCS as follows

p mcs  N Tf / N mcs , f

(12)

where N Tf is the true number of failure points in fact. Then the relative error in the kth iteration can be given as

  k

p mcs  pˆ kf f p mcs f

 100% 

N Tf  Nˆ kf N Tf

 100%.

(13)

Although it is impossible to know the true value of N Tf in practice, we can have access to the following relationship between N Tf and Nˆ kf based on the definition of wrongly classified points

N Tf  Nˆ kf  nˆsk  nˆ kf . Subsequently, the real-time relative error in Eq. (13) can be rewritten as

8

(14)

  1 k

Nˆ kf Nˆ kf  nˆsk  nˆ kf

100%.

(15)

However, nˆ sk and nˆ kf are also unknown. In the next section, we will discuss how to quantify the wrong-classification numbers and obtain the estimation error. Since the relationships among different limit states vary for series system, parallel system and hybrid system, Section 2 and Section 3 will demonstrate the derivation of estimation error of system failure probability for series system and parallel system respectively. As for the hybrid system, the derivation is similar to that of series or parallel system. The only difference is that the relationships among different limit states are more complex for hybrid system which can be decomposed by series system and parallel system.

3.2 The maximum estimation error for series system For CRA with only a single LSF, the probability of wrong sign estimation (WSE) can be defined as j j   gˆ (x mcs ) 0   gˆ (x mcs ) j  F (0)   ( )  1   ( ) if  gˆ ( x mcs )0 j j  ( x )  ( x )  gˆ mcs gˆ mcs j Pwse (x mcs )   . j j  gˆ (x mcs ) 0   gˆ (x mcs )  j )  1  ( ) if  gˆ ( x mcs )  0 1  F (0)  1   ( j j  gˆ (x mcs )  gˆ (x mcs ) 

(16)

While for SRA with multiple LSFs, a probability vector of WSE for all the LSFs at every random point can be then constructed as 1 j p j Pwse   Pwse (xmcs ),..., Pwse (xmcs )  ,

(17)

i j (x mcs ), i  1,..., p the probability of wrong sign estimation for ith where p is the number of LSFs, Pwse

j , j  1,..., N mcs . But different from CRA, the probability vector cannot be LSF at random point x mcs

directly used to search for new training points. Since all the component LSFs are related while determining the ultimate state of a system, a new concept, named as the probability of wrong state classification (WSC), is introduced. More importantly, because the series system and parallel system have different inner principles, we first discuss the series system. In series system, a random point in safe state is expected to be correctly classified, if its signs are j_s correctly estimated for all LSFs. Therefore, for safe point ( xˆ mcs ) predicted by the current Kriging, its

probability of WSC can be calculated by 9

p

i ˆ j_s ˆ j_s Psseries , wsc ( x mcs )  1   (1  Pwse ( x mcs )),

(18)

i 1

where j _ s  1,..., Nˆ sk . Then the number of wrongly classified safe points can be computed by nˆsk 

Nˆ sk

I

j _ s 1

series s, j _ s

j_s (xˆ mcs ).

(19)

series j_s It is defined that indicator function I s , j _ s takes one when xˆ mcs is wrongly classified, and zero when

series correctly classified. Evidently, I s , j _ s follows an independent Bernoulli distribution with the following

mean and variance series ˆ j_s  ˆ j_s E  I sseries , j _ s (xmcs )   Ps , wsc (xmcs ),

(20)

series series ˆ j_s  ˆ j_s ˆ j_s Var  I sseries , j _ s (xmcs )   Ps , wsc (xmcs )(1  Ps , wsc (xmcs )).

(21)

k series Since nˆ sk is the sum of I s , j _ s , it follows a Poisson binomial distribution [57, 65] ( nˆs

PB(nˆk ,  n2ˆk ) ) s

s

with the following mean and variance

nˆ  k s



2 nˆsk



Nˆ sk

P

Nˆ sk

P

j _ s 1

j_s (xˆ mcs ),

(22)

j_s ˆ j_s (xˆ mcs )(1  Psseries , wsc (xmcs )).

(23)

j _ s 1

series s , wsc

series s , wsc

Although the mean value of nˆ sk is obtained, it is still necessary to analytically derive the confidence interval of nˆ sk based on the distribution type and parameter. As a large amount of random points need to be generated by MCS for estimating the system failure probability that is usually small, it is reasonable to assume that the number of safe points Nˆ sk is large enough. According to the Lindeberg central limit theorem (Lindeberg CLT) [63], the distribution of nˆ sk will converge to a normal distribution 2 with mean nˆk and variance  nˆk as Nˆ sk increases to infinite. Hence, the confidence interval can be s

s

approximately given by

nˆsk  nˆsk ,l , nˆsk ,u    Fnˆk1 ( 2), Fnˆk1 (1   2) , s  s 

(24)

1

where Fnˆs is the normal inverse cumulative distribution function,  is the confidence level, which is set to 0.05 in this work. As for the random point classified into failure state, it is assumed that its signs are estimated to be negative for p1 LSFs and positive for p2 LSFs ( p1  p2  p ). Then the failure point is expected to 10

be wrongly classified when the corresponding p1 signs are all wrongly estimated and the p2 signs are all correctly estimated. Accordingly, the probability of WSC is calculated by p1

p2

i1 1

i2 1

i1 i2 ˆ j_ f ˆ j_ f ˆ j_ f Pfseries , wsc (x mcs )   Pwse (x mcs ) (1  Pwse (x mcs )),

(25)

series where i1  1,..., p1 , i2  1,..., p2 and j _ f  1,..., Nˆ kf . Similarly, defining an indicator function I f , j _ f

j_ f which takes one when xˆ mcs is wrongly classified and zero when correctly classified. This indicator

function follows an independent Bernoulli distribution with the following mean and variance series ˆ j_ f  ˆ j_ f E  I series f , j _ f (xmcs )   Pf , wsc (xmcs ),

(26)

series series ˆ j_ f  ˆ j_ f ˆ j_ f Var  I series f , j _ f (xmcs )   Pf , wsc (xmcs )(1  Pf , wsc (xmcs )).

(27)

series The sum of I f , j _ f also follows a Poisson binomial distribution [57, 65] ( nˆ kf

PB( nˆk ,  n2ˆk ) ) with the f

f

following mean and variance

nˆ  

2 nˆ kf



Nˆ kf



j _ f 1

Nˆ kf



ˆ j_ f Pfseries , wsc ( x mcs ),

(28)

series ˆ j_ f ˆ j_ f Pfseries , wsc (x mcs )(1  Pf , wsc (x mcs )).

(29)

k f

j _ f 1

Furthermore, since the number of failure points is not large enough, the approximation for the confidence interval of nˆ sk is not suitable for nˆ kf . But nˆ kf can approximately follow the following Poisson distribution based on the Le cam’s theorem [65] 

e Pr(nˆ  j _ f )  k f

nˆ k f

 j_ f nˆ k f

j_ f !

, j _ f  1,..., Nˆ kf .

(30)

As a result, the confidence interval of nˆ kf can be approximated by

nˆ kf   nˆ kf ,l , nˆ kf ,u    Pnˆk1 ( 2), Pnˆk1 (1   2)  , f  f 

(31)

where Pnˆk1 is the inverse cumulative distribution function of the Poisson distribution. According to the f

confidence intervals of nˆ sk and nˆ kf , the kth maximum real-time estimation error can be calculated by

 Nˆ kf Nˆ kf k ˆmax  max  1  k , 1   Nˆ f  nˆ kf ,u Nˆ kf  nˆsk ,u 

  100%.  

(32)

3.3 The maximum estimation error for parallel system Since the principle of classifying a random point into safe or failure state is different between 11

series system and parallel system, it is necessary to further discuss the maximum estimation error for parallel system. For a random point classified into safe state, it is assumed that its signs are estimated to be negative for p1 LSFs and positive for p2 LSFs ( p1  p2  p ). Then the safe point is expected to be wrongly classified when the corresponding p1 signs are all correctly estimated and the p2 signs are all wrongly estimated. Therefore, the probability of WSC is calculated by p1

p2

i1 1

i2 1

i1 i2 j_s j_s j_s Ps ,parallel (xˆ mcs )   (1  Pwse (xˆ mcs )) Pwse (xˆ mcs ), wsc

(33)

where i1  1,..., p1 , i2  1,..., p2 and j _ s  1,..., Nˆ sk . For a random point classified into failure state, it is expected to be wrongly classified if its signs for the LSFs are all wrongly estimated. In this situation, the probability of WSC is computed by p

i ˆ j_ f ˆ j_ f Pfparallel , wsc ( x mcs )  1   (1  Pwse ( x mcs )).

(34)

i 1

Similar to those mentioned in Section 3.2, the number of wrongly classified safe points in Nˆ sk approximately follow the normal distribution with the following mean and variance

nˆ  k s

 n2ˆ  k s

Nˆ sk

P

(35)

j_s j_s (xˆ mcs )(1  Ps ,parallel (xˆ mcs )), wsc

(36)

j _ s 1

Nˆ sk

P

j _ s 1

j_s (xˆ mcs ),

parallel s , wsc

parallel s , wsc

while the number of wrongly classified failure points in Nˆ f approximately follow the Poisson k

distribution with the following mean

nˆ  k f

Nˆ kf



j _ f 1

ˆ j_ f Pfparallel , wsc (x mcs ).

(37)

Based on Eq. (35)-(37), the confidence interval of nˆ sk and nˆ kf are approximately determined as

nˆsk ,l , nˆsk ,u    Fnˆk1 ( 2), Fnˆk1 (1   2)  s  s 

and

 nˆ kf ,l , nˆ kf ,u    Pnˆk1 ( 2), Pnˆk1 (1   2)   f  f

respectively.

Additionally, the kth maximum estimation error of parallel system can be calculated by Eq. (32).

3.4 The refinement of Kriging As the purpose of adaptive Kriging is to estimate the failure probability with high accuracy and efficiency, the most effective way to achieve this purpose is to reduce the number of the wrongly classified points ( nˆ sk and nˆ kf ) and decrease the estimation error as much as possible when updating 12

Kriging models in each iteration. Therefore, the candidate point with the largest probability of WSC should be added to the training set to increase the accuracy of the local region around it. In each iteration, the new training point can be identified by m ˆ j _ s series ˆ j _ f xnew  xmcs , m  arg max ( Psseries , wsc ( x mcs ), Pf , wsc ( x mcs )) j _ s, j _ f

ˆ j _ s parallel ˆ j _ f or m  arg max ( Ps ,parallel wsc (x mcs ), Pf , wsc (x mcs )).

(38)

j _ s, j _ f

Since some LSFs may have little or no effect on the system failure domain, only one LSF that is the most important to the system failure surface would be refined by this new training point. How to determine the target LSF for series system or parallel system is discussed respectively. For series system, a random point will be classified into the safe state if all the LSFs are in safe state. This safe point is expected to be wrongly classified when at least one LSF’s sign at this point is wrongly estimated. In order to reduce the probability of WSC, the LSF for which the point has the largest WSE should be selected. In contrast, a point will be divided into the failure state if at least one LSF is in failure state. To guarantee the failure state of this point, at least one LSF’s sign among the LSFs in failure state is expected to be correctly estimated. Therefore, the LSF for which the point has the smallest WSE should be selected among the LSFs in failure state. The target LSF is determined by i arg max ( Pwse (xnew )),  gˆi (xnew )  0 i  1,..., p  i 1,..., p N_p i ( Pwse (xnew )),  gˆi (xnew )  0 i  1,..., p1 arg i min 1,..., p1

(39)

where N _ p is the serial number of LSF. In terms of parallel system, it is opposite to series system. Once a point is in safe state for more than one LSF, the parallel system will be in safe state. Hence the refinement can concentrate on increasing the local accuracy for those LSFs in safe state. However, only when a point is in failure state for all the LSFs, the parallel system will be in failure state. Therefore, the refinement should try to guarantee that all the LSFs are correctly estimated. The number of the target LSF is determined by i arg min ( Pwse (xnew )),  gˆi (xnew )  0 i  1,..., p2  i 1,..., p2 N_p i ( Pwse (xnew )),  gˆi (xnew )  0 i  1,..., p arg imax 1,..., p

(40)

Based on Eq. (38)-(40), the N_ pth LSF will be refined by the new training set after the new point x new is added. During the refinement, the stopping criterion is as important as the update of training set. As the 13

added new point tries to decrease the real-time estimation error, it is wise to take this error into consideration while judging whether the reliability result is accurate or not. Therefore, the refinement stops when the maximum estimation error of the failure probability meets the following condition k ˆmax   tar ,

(41)

where  tar is a predefined threshold.

3.5 Implementation procedure of EEK-SYS For the sake of simplicity and clarity, the main implementation procedure is summarized as follows. j Step 1 Generate the random points xmcs , j  1..., N mcs around the mean point of random variables

according to the statistics information. Step 2 Generate the initial training point sets. For every LSF, the training set is defined by DoEi  {(x1 , gi1 ),..., (x h , gih )} , where i  1,..., p and h is the number of initial training points.

Step 3 Construct the Kriging model gˆ i (X) for every LSF. In this paper, Gaussian correlation function is selected for the calculation of the correlation, and the differential evolution algorithm [66] is utilized to optimize the correlation length parameter. After the construction, the system failure probability is estimated. Step 4 Calculate the probability of wrong sign estimation (WSE). For every LSF, the probability of j ) is calculated for all random points based on Kriging prediction. WSE Pwse ( x mcs

Step 5 Compute the probability of wrong state classification (WSC). For series system, the series ˆ j_s ˆ j_ f Psseries , wsc ( x mcs ) and Pf , wsc ( x mcs ) of random points are obtained by Eq. (18) and Eq. (25)

parallel j_s parallel j_ f respectively. For parallel system, the Ps , wsc (xˆ mcs ) and Pf , wsc ( xˆ mcs ) of random points are

obtained by Eq. (33) and Eq. (34) respectively. Step 6 Identify the new training point x new based on WSC and WSE through Eq. (38-40). Step 7 Add the new point x new to the corresponding training set DoE N _ p . k Step 8 Estimate the maximum estimation error ˆmax . Based on the probability of WSC of safe points

and failure points, the corresponding confidence intervals of the wrong-classification numbers 14

k are computed through Eq. (24) and Eq. (31) respectively. Then ˆmax is calculated by Eq. (32).

Step 9 Check the stopping criterion in Eq. (41). If it is satisfied, go to step 10; otherwise, go to step 3. Step 10

Calculate the coefficient of variation of the system failure probability. The formula is given

by

COV ( Pˆf ) 

1  pˆ kf N mcs pˆ kf

.

(42)

If the coefficient of variation is smaller than a predefined threshold that is usually set to 5%, then stop EEK-SYS and output pˆ kf . Otherwise, more random points are generated by MCS and added to the current random point set. Then, go back to Step 4.

4 Examples and discussions In this section, three examples are studied to demonstrate the performance of the proposed EEK-SYS method. The first example is a parallel system with disconnected failure domains. The second one is a roof truss structure with three series failure modes. The third one is an engineering application for the system reliability analysis of a self-balancing vehicle. EEK-SYS is compared with a newly proposed AK-SYSi method. Since it has been proved that AK-SYSi outperforms the classical AK-SYS method [53], the comparison with AK-SYS will not be repeated. To be fair, a dozen of initial sample points [27] are generated for both AK-SYSi and RETK-SYS. To ensure the stability of the results, all the methods are repeated 100 times for example 1 and 2. The average results are listed in the following tables, where Ncall represents the mean value of the number of LSF calls. Moreover, in order to demonstrate the performance of EEK-SYS more specifically, we have tested EEK-SYS with different thresholds. The readers can select the threshold value form interval [0.001, 0.05] according to their requirement about accuracy and efficiency. Additionally, the system failure probability directly estimated by MCS is regarded as the reference. The relative error is calculated by

 mcs

where p f

p mcs  pˆ f f p mcs f

100%

(43)

is the failure probability estimated by MCS, and pˆ f is the failure probability estimated 15

by AK-SYSi or RETK-SYS.

4.1 A parallel system with disconnected failure regions This parallel system [54] has three LSFs, which are defined as

g1 ( X)  (4  2.1X 12  X 14 / 3) X 12  X 1 X 2  (4  4 X 22 ) X 22  0.8 g 2 ( X)  100(4  2.1X 12  X 14 / 3) X 12  100 X 1 X 2  100( 4  4 X 22 ) X 22  60

(44)

g3 ( X)  500(4  2.1X  X / 3) X  500 X 1 X 2  500( 4  4 X ) X  250 2 1

4 1

2 1

2 2

2 2

where both X 1 and X 2 follow the normal distribution with 0 mean and 0.3 standard deviation. Accordingly, the system failure probability is calculated by Pfparallel  P  g1 ( X)  0

g 2 ( X)  0

g 3 ( X)  0 .

(45) 5

This parallel SRA problem is solved by using MCS, AK-SYSi and EEK-SYS. For each LSF, 10

* * random points are generated. The stopping criterion of AK-SYSi is min(U )  U min ( U min is equal to 2

in this case), while EEK-SYS is tested with seven different stopping thresholds. Table 1 lists the results of different methods. All the methods obtain the relatively accurate failure probability according to the true relative error, and the maximum true relative error of EEK-SYS is 1.15%. With the decrease of stopping thresholds, the stopping criterion becomes much stricter and EEK-SYS becomes more accurate. Moreover, Figure 1 indicates one test of EEK-SYS to make the comparison between the real-time true relative error and the real-time maximum estimation error. It is observed that the maximum estimation error of EEK-SYS starts to capture and pursue the true relative error after 10 iterations. Subsequently, the maximum estimation error gradually reduces with the refinement of Kriging until the stopping criterion

ˆmax  0.01 is satisfied. In terms of the number of function evaluations, all the EEK-SYSs with different thresholds are more efficient than AK-SYSi, even the Ncall of EEK-SYS with  tar  0.001 is slightly smaller than that of AK-SYSi. This comparison between AK-SYSi and EEK-SYS is more specific in Figure 3, which demonstrates the variation in Ncall via boxplot with the consideration of different error thresholds. The boxplot shows that EEK-SYS considerably reduces the Ncall compared to AK-SYSi. Generally, as the error threshold increases, the average number of Ncall decreases. Furthermore, it can be observed that the majority of computational cost is paid for approximating the first LSF. Except for the 16

initial points, small training points are generated to build the Kriging model for the second and third LSFs. The main cause of this situation is that the failure domains are completely determined by the first LSF. To illustrate this situation more clearly, the DoEs of one test for AK-SYSi and EEK-SYS (  tar  0.01 ) are shown in Figure 2, where more points are located in the vicinity of the first LSF. Additionally, the difference ( ˆmax   ) between the maximum estimation error and true relative error is also illustrated in Figure 4. The boxplot indicates that the true relative error is smaller than the maximum estimation error in most cases. However, sometimes the former is slightly larger than the latter, and this situation occurs more frequently with the increase of the error threshold. This phenomenon is caused by two facts that ˆmax is derived based on a prescribed confidence level and the relaxed stopping criterion enlarges the possibility of premature convergence. Table 1 Results of the parallel system in Section 4.1 (3.393 3.393 7.994 11.113)

Ncall 2 p parallel ( 10 ) f

Method

Relative error

g1

g2

g3

total

MCS( 3.393 / 7.994 /11.113 10 )

105

105

105

3 105

3.393

-

* AK-SYSi( U min  2)

53.3

15.9

14.7

83.9

3.392

0.0003

EEK-SYS(  tar  0.05 )

24.2

15.1

13.3

52.6

3.354

0.0115

EEK-SYS(  tar  0.04 )

27.0

15.1

13.1

55.2

3.366

0.0080

EEK-SYS(  tar  0.03 )

29.4

15.5

13.3

58.2

3.373

0.0059

EEK-SYS(  tar  0.02 )

30.1

15.7

13.5

59.3

3.376

0.0050

EEK-SYS(  tar  0.01 )

37.8

15.6

13.5

66.9

3.387

0.0018

EEK-SYS(  tar  0.005 )

41.5

16.0

13.7

71.2

3.388

0.0015

EEK-SYS(  tar  0.001 )

48.7

16.6

14.6

79.9

3.394

0.0003

2

Note: 3.393 / 7.994 /11.113 102 represent the component failure probability of g1 , g 2 and g 3 function respectively. It is observed that g1 plays a very important role in this parallel system. That is 17

why the training points are mainly located near the limit state of g1 function shown in Figure 2.

Figure 1 The true relative error and maximum estimation error of EEK-SYS(  tar  0.01 ) in Section 4.1

(a) AK-SYSi

18

(b) EEK-SYS(  tar  0.01 ) Figure 2 DoE of AK-SYSi and EEK-SYS in Section 4.1

Figure 3 Boxplots of Ncall in Section 4.1

19

ˆmax  

 tar Figure 4 Boxplots of error in Section 4.1

4.2 Roof truss structure with three series failure modes As shown in Figure 5, the top chords and compression bars of the roof truss [53] are made of steel reinforced concrete, while the bottom chords and tension bars are made of steel. Assume the uniformly distributed load q that can be transformed into the nodal load P  ql / 4 is applied on the roof truss structure, where l is the length of the truss. This roof truss structure takes into account three failure modes. The first failure mode considers the perpendicular deflection of the truss peak node C, which is derived as C 

ql 2  3.81 1.13     , where AC , AS are the cross sectional areas of the steel 2  AC EC AS ES 

reinforced concrete and the steel bars respectively, and EC , ES are the corresponding elastic modulus. For this mode, the failure event occurs when  C exceeds 0.03 m. The second failure mode is that the internal force of AD bar ( 1.185ql ) exceeds the ultimate stress f C AC , where f C is the compressive strength of AD bar. For the third mode, the failure occurs when the internal force of EC bar ( 0.75ql ) exceeds the ultimate stress f S AS , where f S is the tensile strength of EC bar. Therefore, the LSFs of the three failure modes are given by

20

ql 2  3.81 1.13     2  AC EC AS ES  g 2  fC AC  1.185ql g1  0.03 

(46)

g3  f S AS  0.75ql

As the roof truss structure is a series system, the failure probability is defined as Pfseries  P  g1  0

g2  0

g 3  0 .

(47)

There are eight random variables in this case and they are all assumed to follow the lognormal distribution. Additionally, all the random variables are mutually independent and the statistical information is listed in Table 2. Table 2 The distribution parameters of roof truss structures Variable

Description

Distribution

Mean

coefficient of variation

q(N/m)

Uniform load

Lognormal

20,000

0.07

l(m)

Length

Lognormal

12

0.01

AS(m2)

Cross-sectional area

Lognormal

9.82 104

0.06

AC(m2)

Cross-sectional area

Lognormal

0.04

0.12

ES(N/m2)

Elastic modulus

Lognormal

2 1011

0.06

EC(N/m2)

Elastic modulus

Lognormal

3 1010

0.06

fS(N/m2)

Tensile strength

Lognormal

3.35 108

0.12

fC(N/m2)

Compressive strength

Lognormal

1.34 107

0.18

21

C D

F

A

B

E

G qN/m

P P

C Ac

D

Ac

As

Ac

F

As

Ac

A

B 3As 0.278l

2As

E 0.222l

G 0.25l

l/12 l/12

P

3As 0.25l

Figure 5 The roof truss structure AK-SYSi and EEK-SYS with seven different stopping thresholds are tested in this example. To make the comparison, the result of MCS is regarded as the reference. For each LSF, 5 10 random 5

points are generated by MCS. The following Table 3 lists the results, both AK-SYSi and EEK-SYS can obtain the relatively accurate result compared to MCS. Similar to the first example, the estimated failure probability of EEK-SYS with  tar  0.05 has the largest relative error (equal to 0.77%). Furthermore, Figure 6 illustrates one test of EEK-SYS (  tar  0.01 ) to compare the real-time maximum estimation error with true relative error. During the initial sampling stage, the true relative error is much larger than the maximum estimation error, which means the Kriging model is inaccurate. With the refinement, the maximum estimation error is able to gradually capture and pursue the true relative error after 10 iterations, and the former is always slightly larger than the latter. Gradually, both of them are reduced and the difference between them is narrowed. As for the number of LSF evaluations, there is an interesting thing that both AK-SYSi and EEK-SYS pay no attention to the first LSF and little attention to the third LSF. More training samples are expected to be generated for the second LSF, which is the most significant one for the system failure domain. Nevertheless, AK-SYSi requires more Ncall than EEK-SYS methods except for EEK-SYS 22

* * with  tar  0.001 , although its stopping threshold is modified from U min  2 to U min  1 . This

situation is also illustrated by the boxplot of Ncall shown in Figure 7. Compared to AK-SYSi, Ncall of EEK-SYS with  tar  0.05,0.04,0.03,0.02,0.01,0.005 is reduced by 58.7%, 56.9%, 53.7%, 48.7%, 35.4% and 21.3% respectively. Generally, the average number of calls to the LSFs decreases as the threshold

 tar of EEK-SYS increases for this series system. In addition, Figure 8 demonstrates the average difference between the maximum estimation error and true relative error. Theoretically, the maximum estimation error is expected to be larger than the true relative error. But the opposite situation also occurs for few cases, the corresponding reasons are similar to those mentioned in Section 4.1. Table 3 Results of the series system in Section 4.2

Ncall 3 p series ( 10 ) f

Method

Relative error

g1

g2

g3

total

MCS( 0 / 3.353 / 0.039 103 )

5 105

5 105

5 105

1.5 106

3.391

-

* AK-SYSi( U min 1)

12

91.2

23.2

126.4

3.383

0.0024

EEK-SYS(  tar  0.05 )

12

27.4

12.8

52.2

3.417

0.0077

EEK-SYS(  tar  0.04 )

12

29.3

13.1

54.4

3.410

0.0056

EEK-SYS(  tar  0.03 )

12

33.7

12.8

58.5

3.400

0.0027

EEK-SYS(  tar  0.02 )

12

39.8

13.0

64.8

3.383

0.0024

EEK-SYS(  tar  0.01 )

12

54.7

14.4

81.1

3.376

0.0044

EEK-SYS(  tar  0.005 )

12

71.0

16.5

99.5

3.388

0.0009

EEK-SYS(  tar  0.001 )

12

123.9

35.9

171.8

3.390

0.0003

Note: 0 / 3.353 / 0.039 103 represent the component failure probability of g1 , g 2 and g 3 function respectively. Only g 2 and g 3 determine the final state of this series system and g 2 is more important than g 3 , which can be also observed from the number of function calls of each component function. 23

Figure 6 The comparison between true relative error and maximum estimation error in Section 4.2

Figure 7 Boxplots of Ncall in Section 4.2

24

ˆmax  

 tar Figure 8 Boxplots of error in Section 4.2

4.3 Self-balancing vehicle As a short-distance transport, the self-balancing vehicle [67] has been widely applied in various industries. It intelligently controls the body attitude by the self-balancing technology. As shown in Figure 9, the self-balancing vehicle is comprised of the control unit, driving wheel, lighting module and chassis. It is worth noting that the chassis is one of the most important core parts and its reliability affects the performance and lifetime of self-balancing vehicle. In practice, the self-balancing vehicle is designed to satisfy many requirements in different working conditions, such as acceleration, deceleration, turn, upslope and downhill. As shown in Figure 10, two conditions with different loads are considered in this example. The FEM of the chassis consists of 37948 four-node-shell elements. Control unit

Driving Wheel

Chassis Lighting module

25

(a) Outline drawing

(b) Exploded drawing

Figure 9 The self-balancing vehicle[67]

(a) Working condition 1

(b) Working condition 2 Figure 10 Two conditions with different loads The failure event is defined as the deformation of chassis exceeds the allowable threshold. For the two working conditions, the corresponding limit state functions are respectively defined by g1 =D1allow  D1 ( X, Y1 , E )，

(48)

g 2 =D2allow  D2 ( X, Y2 , E )，

(49)

where D1allow and D2allow are the allowable deformation, Y1 and Y2 represent the random loads under two working conditions. In this example, the length and width of chassis ( X ), elastic modulus ( E ), allowable deformation and random loads are regarded as random variables. The details of statistical information are listed in Table 4. As a series system with two failure modes, the failure probability is 26

defined by Pfseries  P  g1  0

g 2  0 .

(50)

Table 4 Random variables of example 5.3 Variable

Distribution

Mean

Standard deviation

Chassis length X 1 (mm)

Lognormal

450

1

Chassis width X 2 (mm)

Lognormal

250

1

Material elastic modulus E (MPa)

Lognormal

72000

2000

Allowable deformation D1allow (mm)

Normal

5.5

0.11

Allowable deformation D2allow (mm)

Normal

3.5

0.07

load Y1 ( MPa)

Normal

1.5

0.1

load Y2 ( MPa)

Normal

1.5

0.1

Since performing one FEA for the two working conditions is relatively time-consuming compared to the mathematical examples in Section 4.1 and 4.2, it is not appropriate to obtain the reference result by MCS and repeat 100 runs for all the methods. Therefore, only AK-SYSi and EEK-SYS with five different error thresholds are performed for this example, and they are just repeated for five times. Table 5 gives the average results. As for the predicted failure probability, almost all the methods obtain the similar results. Specially, the prediction of EEK-SYS with  tar  0.05 deviates a little compared to others due to the loose stopping criterion. Figure 11 illustrates one test of the real-time maximum estimation error, and the maximum estimation errors of different EEK-SYSs are all smaller than the corresponding error thresholds. In terms of the number of LSF evaluations, it is found that the majority of

Ncall are assigned to the second LSF. Furthermore, all the EEK-SYSs are more efficient than AK-SYSi, and they reduce the number of LSF calls by 62.0%, 56.8%, 52.7%, 39.5% and 20.0% respectively. Table 5 Results of the self-balancing vehicle in Section 4.3

Ncall 2 p series ( 10 ) Maximum estimation error f

Method

g1

g2

total

27

AK-SYSi

25.5

70.3

95.8

1.480

-

EEK-SYS(  tar  0.05 )

16.4

20.0

36.4

1.499

0.0484

EEK-SYS(  tar  0.04 )

16.8

24.6

41.4

1.485

0.0370

EEK-SYS(  tar  0.03 )

17.0

28.3

45.3

1.478

0.0275

EEK-SYS(  tar  0.02 )

19.6

38.4

58.0

1.482

0.0197

EEK-SYS(  tar  0.01 )

22.0

54.6

76.6

1.479

0.0096

Figure 11 The maximum estimation error of EEK-SYS with different thresholds in Section 4.3

5 Conclusions This work aims at solving the system reliability analysis problem through adaptive Kriging approximation with high efficiency and accuracy. Since the estimation error of the system failure probability is expected to be directly related to the refinement of Kriging, this work achieves to quantify the real-time estimation error by deriving the probability of wrong state classification for every point and obtaining the confidence intervals of number of the wrongly classified points. As a result, the proposed EEK-SYS can always pursue the accuracy of the predicted system failure probability. Specifically, the probability of WSE and WSC assist to find the most effective way to refine the Kriging. Furthermore, the real-time maximum estimation error is able to reflect the actual situation of the refinement, and determine when to terminate the training process. 28

The comparison results of the parallel system with disconnected failure regions, roof truss structure with three series failure modes and self-balancing vehicle demonstrate that the proposed method can avoid some useless but computationally expensive samples while training the Kriging model. Moreover, EEK-SYS can estimate the relative error of the predicted system failure probability with high precision, which is very important to the termination of the refinement. In this paper, we do not handle the small failure probability problems since a larger amount of random points (even to 108 ) require to be generated by MCS. This situation can be avoided by combining the proposed EEK-SYS method with the variation reduction approaches, such as importance sampling and subset simulation. In addition, with the help of independent transformation of related variables and sample classification, the proposed method can have a wider range of application. For future work, the proposed method will be extended to the more complex reliability analysis problems with various source of uncertainties (such as dynamic uncertainty) or to the reliability-based design optimization problems with the consideration of multiple systems.

Declaration of interests

☒ √ The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Author Statement Chen Jiang: Conceptualization, Methodology, Writing-Original Draft Haobo Qiu: Writing-Review & Editing, Funding acquisition Liang Gao: Supervision, Funding acquisition Dapeng Wang: Investigation, Validation Zan Yang: Writing-Original Draft Liming Chen: Writing-Original Draft

29

Acknowledgement This research is supported by the National Natural Science Foundation of China under Grant No. 51675198, the National Defense Pre-Research Foundation of China under Grant No. 41423010205, National Science Foundation for Distinguished Young Scholars of China under Grant No. 51825502, and the Program for HUST Academic Frontier Youth Team Grant No. 2017QYTD04.

References [1] Keshtegar B, Kisi O. RM5Tree: Radial basis M5 model tree for accurate structural reliability analysis. Reliability Engineering & System Safety. 2018;180:49-61. [2] Meng Z, Keshtegar B. Adaptive conjugate single-loop method for efficient reliability-based design and topology optimization. Computer Methods in Applied Mechanics and Engineering. 2019;344:95-119. [3] Xu J, Dang C. A novel fractional moments-based maximum entropy method for high-dimensional reliability analysis. Applied Mathematical Modelling. 2019;75:749-68. [4] Yu S, Wang Z. A general decoupling approach for time- and space-variant system reliability-based design optimization. Computer Methods in Applied Mechanics and Engineering. 2019;357:112608. [5] Chen Z, Qiu H, Gao L, Su L, Li P. An adaptive decoupling approach for reliability-based design optimization. Computers & Structures. 2013;117:58-66. [6] Echard B, Gayton N, Lemaire M, Relun N. A combined Importance Sampling and Kriging reliability method for small failure probabilities with time-demanding numerical models. Reliability Engineering & System Safety. 2013;111:232-40. [7] Papaioannou I, Geyer S, Straub D. Improved cross entropy-based importance sampling with a flexible mixture model. Reliability Engineering & System Safety. 2019;191:106564. [8] Depina I, Le TMH, Fenton G, Eiksund G. Reliability analysis with Metamodel Line Sampling. Structural Safety. 2016;60:1-15. [9] Valdebenito MA, Jensen HA, Hernández HB, Mehrez L. Sensitivity estimation of failure probability applying line sampling. Reliability Engineering & System Safety. 2018;171:99-111. [10] Xiao M, Zhang J, Gao L, Lee S, Eshghi AT. An efficient Kriging-based subset simulation method for hybrid reliability analysis under random and interval variables with small failure probability. Structural and Multidisciplinary Optimization. 2019;59:2077-92. [11] Du X, Hu Z. First Order Reliability Method With Truncated Random Variables. Journal of Mechanical Design. 2012;134:091005. [12] Meng Z, Yang D, Zhou H, Wang BP. Convergence control of single loop approach for reliability-based design optimization. Structural and Multidisciplinary Optimization. 2018;57:1079-91. [13] Torii AJ, Lopez RH, Miguel LFF. A second order SAP algorithm for risk and reliability based design optimization. Reliability Engineering & System Safety. 2019;190:106499. [14] Xiaoke L, Fengxiang X, Haobo Q, Zhenzhong C, Wenbin H, Jun M. A MOVING SHIFTING VECTOR

METHOD

FOR

RELIABILITY-BASED

DESIGN

OPTIMIZATION

USING

EFFECTIVENESS CHECKING OF PROBABILISTIC CONSTRAINT. International Journal of 30

Industrial Engineering. 2019;26:34-47. [15] Wu J, Zhang D, Liu J, Han X. A Moment Approach to Positioning Accuracy Reliability Analysis for Industrial Robots. IEEE Transactions on Reliability. 2019. DOI: 10.1109/TR.2019.2919540 [16] Zhang D, Han X. Kinematic Reliability Analysis of Robotic Manipulator. Journal of Mechanical Design. 2020;142: 044502. [17] Jiang C, Qiu H, Li X, Chen Z, Gao L, Li P. Iterative reliable design space approach for efficient reliability-based design optimization. Engineering with Computers. 2020;36:151-69. [18] Zhang D, Han X, Jiang C, Liu J, Li Q. Time-Dependent Reliability Analysis Through Response Surface Method. Journal of Mechanical Design. 2017;139:041404. [19] Dong Y, Teixeira A, Soares CG. Time-variant fatigue reliability assessment of welded joints based on the PHI2 and response surface methods. Reliability Engineering & System Safety. 2018;177:120-30. [20] Jing Z, Chen J, Li X. RBF-GA: An adaptive radial basis function metamodeling with genetic algorithm for structural reliability analysis. Reliability Engineering & System Safety. 2019;189:42-57. [21] Dong H, Sun S, Song B, Wang P. Multi-surrogate-based global optimization using a score-based infill criterion. Structural and Multidisciplinary Optimization. 2019;59:485-506. [22] Roy A, Manna R, Chakraborty S. Support vector regression based metamodeling for structural reliability analysis. Probabilistic Engineering Mechanics. 2019;55:78-89. [23] Hawchar L, El Soueidy C-P, Schoefs F. Principal component analysis and polynomial chaos expansion for time-variant reliability problems. Reliability Engineering & System Safety. 2017;167:406-16. [24] Qian J, Yi J, Cheng Y, Liu J, Zhou Q. A sequential constraints updating approach for Kriging surrogate model-assisted engineering optimization design problem. Engineering with Computers. 2019. DOI: https://doi.org/10.1007/s00366-019-00745-w [25] Shi R, Liu L, Long T, Wu Y, Tang Y. Filter-based adaptive Kriging method for black-box optimization problems with expensive objective and constraints. Computer Methods in Applied Mechanics and Engineering. 2019;347:782-805. [26] Bichon BJ, Eldred MS, Swiler LP, Mahadevan S, McFarland JM. Efficient Global Reliability Analysis for Nonlinear Implicit Performance Functions. AIAA Journal. 2008;46:2459-68. [27] Echard B, Gayton N, Lemaire M. AK-MCS: An active learning reliability method combining Kriging and Monte Carlo Simulation. Structural Safety. 2011;33:145-54. [28] Sun Z, Wang J, Li R, Tong C. LIF: A new Kriging based learning function and its application to structural reliability analysis. Reliability Engineering & System Safety. 2017;157:152-65. [29] Wang Z, Shafieezadeh A. REAK: Reliability analysis through Error rate-based Adaptive Kriging. Reliability Engineering & System Safety. 2019;182:33-45. [30] Jiang C, Qiu H, Yang Z, Chen L, Gao L, Li P. A general failure-pursuing sampling framework for surrogate-based reliability analysis. Reliability Engineering & System Safety. 2019;183:47-59. [31] Meng Z, Zhang Z, Zhou H. A novel experimental data-driven exponential convex model for reliability assessment with uncertain-but-bounded parameters. Applied Mathematical Modelling. 2020;77:773-87. [32] Xiao N-C, Yuan K, Zhou C. Adaptive kriging-based efficient reliability method for structural systems with multiple failure modes and mixed variables. Computer Methods in Applied Mechanics and Engineering. 2019:112649. [33] Xiao N-C, Zuo MJ, Zhou C. A new adaptive sequential sampling method to construct surrogate models for efficient reliability analysis. Reliability Engineering & System Safety. 2018;169:330-8. 31

[34] Zhang J, Xiao M, Gao L, Chu S. A combined projection-outline-based active learning Kriging and adaptive importance sampling method for hybrid reliability analysis with small failure probabilities. Computer Methods in Applied Mechanics and Engineering. 2019;344:13-33. [35] Zhang J, Xiao M, Gao L, Fu J. A novel projection outline based active learning method and its combination with Kriging metamodel for hybrid reliability analysis with random and interval variables. Computer Methods in Applied Mechanics and Engineering. 2018;341:32-52. [36] Zhang Y, Li H, Xiao M, Gao L, Chu S, Zhang J. Concurrent topology optimization for cellular structures with nonuniform microstructures based on the kriging metamodel. Structural and Multidisciplinary Optimization. 2019;59:1273-99. [37] Meng Z, Zhang D, Li G, Yu B. An importance learning method for non-probabilistic reliability analysis and optimization. Structural and Multidisciplinary Optimization. 2019;59:1255-71. [38] Gaspar B, Teixeira AP, Guedes Soares C. Adaptive surrogate model with active refinement combining Kriging and a trust region method. Reliability Engineering & System Safety. 2017;165:277-91. [39] Jiang C, Wang D, Qiu H, Gao L, Chen L, Yang Z. An active failure-pursuing Kriging modeling method for time-dependent reliability analysis. Mechanical Systems and Signal Processing. 2019;129:112-29. [40] Hu Z, Du X. Mixed Efficient Global Optimization for Time-Dependent Reliability Analysis. Journal of Mechanical Design. 2015;137:051401. [41] Hu Z, Mahadevan S. A Single-Loop Kriging Surrogate Modeling for Time-Dependent Reliability Analysis. Journal of Mechanical Design. 2016;138:061406. [42] Wang Z, Chen W. Time-variant reliability assessment through equivalent stochastic process transformation. Reliability Engineering & System Safety. 2016;152:166-75. [43] Jiang C, Qiu H, Gao L, Wang D, Yang Z, Chen L. Real-time estimation error-guided active learning Kriging method for time-dependent reliability analysis. Applied Mathematical Modelling. 2020;77:82-98. [44] Du W, Luo Y, Wang Y. Time-variant reliability analysis using the parallel subset simulation. Reliability Engineering & System Safety. 2019;182:250-7. [45] Ling C, Lu Z, Zhu X. Efficient methods by active learning Kriging coupled with variance reduction based sampling methods for time-dependent failure probability. Reliability Engineering & System Safety. 2019;188:23-35. [46] Cheng K, Lu Z. Time-variant reliability analysis based on high dimensional model representation. Reliability Engineering & System Safety. 2019;188:310-9. [47] Yu S, Wang Z, Zhang K. Sequential time-dependent reliability analysis for the lower extremity exoskeleton under uncertainty. Reliability Engineering & System Safety. 2018;170:45-52. [48] Liu Y, Frangopol DM. Time-dependent reliability assessment of ship structures under progressive and shock deteriorations. Reliability Engineering & System Safety. 2018;173:116-28. [49] Han X, Yang DY, Frangopol DM. Time-variant reliability analysis of steel plates in marine environments considering pit nucleation and propagation. Probabilistic Engineering Mechanics. 2019;57:32-42. [50] Dong Y, Teixeira AP, Guedes Soares C. Application of adaptive surrogate models in time-variant fatigue reliability assessment of welded joints with surface cracks. Reliability Engineering & System Safety. 2020;195:106730. [51] Bichon BJ, McFarland JM, Mahadevan S. Efficient surrogate models for reliability analysis of 32

systems with multiple failure modes. Reliability Engineering & System Safety. 2011;96:1386-95. [52] Fauriat W, Gayton N. AK-SYS: An adaptation of the AK-MCS method for system reliability. Reliability Engineering & System Safety. 2014;123:137-44. [53] Yun W, Lu Z, Zhou Y, Jiang X. 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. 2018;59:263-78. [54] Yang X, Liu Y, Mi C, Tang C. System reliability analysis through active learning Kriging model with truncated candidate region. Reliability Engineering & System Safety. 2018;169:235-41. [55] Yang X, Mi C, Deng D, Liu Y. A system reliability analysis method combining active learning Kriging model with adaptive size of candidate points. Structural and Multidisciplinary Optimization. 2019;60:137-50. [56] Hu Z, Mahadevan S. Global sensitivity analysis-enhanced surrogate (GSAS) modeling for reliability analysis. Structural and Multidisciplinary Optimization. 2015;53:501-21. [57] Wang Z, Shafieezadeh A. ESC: an efficient error-based stopping criterion for kriging-based reliability analysis methods. Structural and Multidisciplinary Optimization. 2019;59:1621-37. [58] Perrin G. Active learning surrogate models for the conception of systems with multiple failure modes. Reliability Engineering & System Safety. 2016;149:130-6. [59] Stern RE, Song J, Work DB. Accelerated Monte Carlo system reliability analysis through machine-learning-based surrogate models of network connectivity. Reliability Engineering & System Safety. 2017;164:1-9. [60] Byun J-E, Noh H-M, Song J. Reliability growth analysis of k-out-of-N systems using matrix-based system reliability method. Reliability Engineering & System Safety. 2017;165:410-21. [61] Sadoughi M, Li M, Hu C. Multivariate system reliability analysis considering highly nonlinear and dependent safety events. Reliability Engineering & System Safety. 2018;180:189-200. [62] Yuan K, Xiao N-C, Wang Z, Shang K. System reliability analysis by combining structure function and active learning kriging model. Reliability Engineering & System Safety. 2020;195:106734. [63] Billingsley P. Probability and measure: John Wiley & Sons; 1986. [64] Sacks J, Welch WJ, Mitchell TJ, Wynn HP. Design and Analysis of Computer Experiments. Statistical Science. 1989;4:409-23. [65] Cam LL. An approximation theorem for the Poisson binomial distribution. Pacific Journal of Mathematics. 1960;10:1181-97. [66] Storn R, Price K. Differential Evolution – A Simple and Efficient Heuristic for global Optimization over Continuous Spaces. Journal of Global Optimization. 1997;11:341-59. [67] Huang ZL, Jiang C, Li XM, Wei XP, Fang T, Han X. A Single-Loop Approach for Time-Variant Reliability-Based Design Optimization. IEEE Transactions on Reliability. 2017;66:651-61.

33