# On High-Dimensional Time-Variant Reliability Analysis with the Maximum Entropy Principle

Authors: Fuliang Zhou; Yu Hou; Hong Nie
Journal: International Journal of Aerospace Engineering
 (2022)

---

## Abstract

The structural reliability analysis suffers from the curse of dimensionality if the associated limit state function involves a large number of inputs. This study develops a reliability analysis method that deals with high-dimensional inputs over time. The probability distribution of the structural response is reconstructed by the maximum entropy principle which is achieved by solving an optimization problem derived from the concept of relative entropy. The optimization problem is transformed into a convex one with respect to the orders of fractional moments and the Lagrange multipliers. Additionally, considering the associated computational issues, it is reformulated with side constraints on the parameters of the maximum entropy distribution. Then, a global optimization procedure is performed. The proposed method is successfully applied to the reliability analysis of a linear and a nonlinear structural system, which involves a large number of inputs deriving from the discretization of the input random processes.


---

## Body

## 1. Introduction

Structural reliability analysis has been intensively studied during the past few decades. Based on the consideration of time factor, structural reliability analysis can be roughly divided into two groups: time-invariant reliability analysis and time-variant reliability analysis. For the first group, as its name implies, the time-invariant reliability methods do not consider the effect of time factor, for example, first-/second-order reliability method (FORM/SORM) [1–3], moment methods [4, 5], response surface methods [6–9], and importance sampling [10, 11]. Therefore, they only take random variables as inputs. However, almost all structural systems involve time-dependent parameters during their period of operation and/or service, e.g., stochastic loading and material property degradation. Thus, time-invariant reliability analysis may provide inaccurate reliability assessment and prediction.In the past decade, time-variant reliability analysis has drawn much more attention than its counterpart because of its capacity of dealing with time factor. Thus, the corresponding time-variant reliability methods can take random variables, random processes, and even random fields as inputs. Many time-variant reliability methods have been developed in recent years. They can be roughly classified into three categories: (1) the out-crossing-rate-based methods [12–18], (2) the extreme value methods [19–27], and (3) the composite limit state methods [28–32]. A common feature of these methods is that a time-variant reliability analysis problem is usually transformed into a time-invariant one or broken down into a series of time-invariant ones. Then, the aforementioned time-invariant reliability methods are applied to complete the time-variant reliability analysis. For instance, FORM has been integrated into many time-variant reliability methods [14, 17]. Note that the time-variant reliability methods with FORM inherent the drawbacks of FORM.Many challenges are encountered by time-variant reliability analysis. One main challenge arises from the curse of dimension when a large number of random variables are used to represent the input random processes and/or fields. It should be pointed out that there is also a universal and dimension-free tool for both time-invariant reliability analysis and time-variant ones. It is direct Monte Carlo simulation (MCS). However, MCS suffers from a huge computational effort for estimating a small failure probability when a time-consuming model or experiment is involved, e.g., a large-scale finite element model. Subset simulation [33, 34] is an alternative for MCS for estimating small failure probabilities in both time-invariant reliability analysis and time-variant ones. However, several thousands of samples are still required to obtain a robust and accurate estimation by subset simulation.In this paper, we are interested in developing an efficient and accurate method for high-dimensional time-variant reliability analysis. The high-dimensional issue is solved through the development of a reliability method using statistical moments.Zhao and Ono’s moment method is the most attractive one for structural reliability analysis [4, 5] after the year 2000, where the first three or four moments of structural response are used to fit a parametric distribution. The Pearson system [35] was recommended as the reliability index for the moment method. Although the Pearson system of frequency curves is flexible, they still impose spurious information on the unknown distribution. It is also contentious from both theoretical and practical views [36, 37]. The maximum entropy (MaxEnt) principle was pioneered by Jaynes [38, 39], which minimizes the amount of prior information built into the approximate distribution while maximizing entropy, i.e., a measure of uncertainty. Using the MaxEnt principle, Li and Zhang proposed a combined method for structural reliability analysis [37]. The statistical moments of structural response are calculated by the dimension reduction method [40, 41]. Zhang and Pandey analyzed two challenging issues of using integer moments when approximating the unknown distribution by the MaxEnt principle [36]. They are a large number of required moments and the oscillatory phenomena at the tail of MaxEnt distribution. The usage of fractional moments [42] was proposed to overcome these two issues for time-invariant reliability problems since a fractional moment contains the statistical information of a large number of integer moments. A multiplicative dimensional reduction method was also developed to compute the fractional moments of structural response [36]. Considering correlated input random variables, Li et al. suggested using the unscented transformation technique to calculate fractional moments [43]. Dai et al. employed the MaxEnt principle to construct the optimal importance sampling density for the importance sampling method [44]. However, the above investigations are limited to time-invariant reliability problems.For the time-variant reliability analysis problem, Shi et al. presented the performance of the combination of dimension reduction method and the MaxEnt principle with fractional moments [45]. Xu and his colleagues carried out a series of studies about the MaxEnt principle with fractional moments and their applications for time-variant reliability analysis [46–50]. However, the number of input random variables in the application examples is not beyond 20 [45–50]. Thus, the abovementioned methods cannot be applied to the time-variant reliability analysis problems involving a large number of inputs [47].The objective of this paper is to propose an efficient method for the high-dimensional time-variant reliability problems by using the MaxEnt principle with fractional moments, where the number of input random variables is larger than 100. The contributions of this study are as follows: (1) a new perspective for estimating a failure probability curve from the combination of fractional moments and the maximum entropy principle and (2) a new tool for high-dimensional and nonlinear time-variant reliability analysis with input random variables and processes. High-dimensional issue is very challenging in many science and engineering fields [51–53]. To address this issue, fractional moment and the maximum entropy distribution are used to approximate the probability density function of structural responses of interested, which has a capacity of dealing with zero-shot events [54] for structural reliability analysis. In other words, the proposed method can estimate structural reliability with zero-failure samples.The remainder of the paper is organized as follows: Section2 briefly reviews the problem setup for time-variant reliability analysis. Section 3 describes the procedure and theoretical rationale of the proposed method. Section 4 analyzes the implementation issues and their corresponding possible solutions. Two high-dimensional examples are used to demonstrate the performance of the proposed method in Section 5. Concluding remarks are provided in Section 6.


## 2. Problem Statement

Limit state functions (LSFs) are employed to capture the intended operating condition of structures in the structural reliability community, which divides the operating space into two domains, i.e., the safe domain and the failure domain. Letg⋅ be an LSF of the intended operating condition of a structure, where g⋅>0 indicates the structure is safe while g⋅≤0 means there is a failure in the structure. The general form of LSF is defined as
(1)G=gX,Yt,t,where G is the response quantity, X is the input random variable vector, Yt is the input random process vector, and t is the input time, respectively. Based on the different type of inputs in the LSF g⋅, the LSF may degrade into simpler types [22]. It is well-known that reliability is a function of time snice the LSF of interest is related to time via the time-dependent random process Yt and the explicit operating time t. Suppose that we are interested in the time-variant reliability of the structure over a time interval 0,tf, where tf is the final time instant. The time-variant reliability R0,tf can be expressed as
(2)R0,tf=PrG=g⋅>0,∀t∈0,tf,where the symbol ∀ means for all time instants. In practical estimation, the corresponding time-variant failure probability Pf0,tf is more preferable which is defined as a complementary quantity to the time-variant reliability, i.e.,
(3)Pf0,tf=PrG=g⋅≤0,∃t∈0,tf,where the symbol ∃ has a meaning of “there exists at least one.” To distinguish the failure probability over a time interval from that for a time instant, Pf0,tf in Equation (3) is also called as the cumulative failure probability. Obviously, it is a function of the final time tf.To handle the input random processes, the time interval of interest0,tf is usually divided into a set of discrete time instants with a constant time step Δt. Let ti=iΔt represent the time instant at the ith step, where the index i=1,⋯,n. Spectral decomposition methods are suggested to deal with the input random processes and transform them into a set of random variables [17, 22, 24, 26]. Let Z be the transformed random variables from the input random processes. Then, the LSF becomes G=gX,Z,t. Furtherly, at a time instant ti, an instantaneous LSF is defined as Gi=gX,Z,ti. Using the concept of series-system reliability, the cumulative failure probability in Equation (3) is approximated as
(4)Pf0,tf≈Pr∪i=1ngX,Z,ti≤0=Pr∪i=1nGi≤0=Prmaxi=0,1,⋯,lGi≤0,where W=maxi=0,1,⋯,nGi. The advantage of Equation (4) is that the output random process G=g⋅ is converted into a random variable W.The magnitude of time step sizeΔt affects the approximate accuracy of input random processes with a finite number of random variables firstly. Then, this effect may propagate into Equation (4) through the LSF and furtherly influence the estimation accuracy of the cumulative failure probability. In this study, we assume that the magnitude of time step size is sufficiently small and adequate to capture the uncertainty in the input random processes and failure information in the transformed series system. The current study focuses on the development of a new high-dimensional time-variant reliability analysis method based on Equation (4).


## 3. Proposed Method

### 3.1. Overview

The new high-dimensional time-variant reliability analysis method has three stages. The first stage is to generate random samples for the input random variables and random processes using MCS sampling strategy and furtherly compute the structural responses. The second stage involves the approximation of the MaxEnt distribution. The third stage is related to the calculation of failure probability using the MaxEnt distribution obtained in the second stage. It is clear that the second stage is the key stage of the proposed method. The implementation procedure of the proposed method is shown in Figure1.Figure 1
Procedure of the proposed method.Stage 1 involves four steps. Step 1: discretize the time interval0,tf with a constant time step size Δt and obtain n+1 discrete time instant t0=0,t1,⋯,tn=tf. Step 2: decompose the input random processes Yt with a group of random variables Z. Step 3: generate N random samples for the input random variables and random processes. Step 4: substitute the random samples into the LSF to obtain N trajectories of structural response.Stage 2 is furtherly divided into three steps. Step 1: solve the optimization problem derived from the concept of relative entropy. Step 2: calculate the fractional moments from theN trajectories of structural response. Step 3: determine the MaxEnt distribution based on the solution to the optimization problem in Step 1. In the subsequent sections, we will discuss the details of Step 1 and 3 since these steps are theoretical fundamentals of the proposed method.Stage 3 is straightforward once the probability density function (PDF) of MaxEnt distribution is determined in Stage 2. The estimation of a cumulative failure probability is just a numerical integration over the failure domain.


### 3.2. Maximum Entropy Principle with Fractional Moments

The advantages of usage of fractional moments in the maximum entropy principle have been proven by several studies in the structural reliability community [36, 43, 45–50]. This subsection will give a brief review of the maximum entropy principle with fractional moments.Identifying the probability distribution with statistical moments has been investigated extensively in the literature. However, the use of integer statistical moments suffers from accuracy concern. Zhang and Pandey study the need of fractional moments instead of integer moments, which is evoked by the relationship between them [36].We start with the definition of the fractional binomial coefficientαi that is given by
(5)αi=αα−1⋯α−i+1ii−1⋯1,where α is a real number and i is an integer number. If applying the Taylor series expansion on a function of xα, one has
(6)xα=∑i=0∞αicα−ix−ci,where c is a reference point.Theαth-order fractional moment of a positive random variable X is defined as
(7)EXα=∫XxαfXxdx.Substituting Equation (6) into Equation (7), one has
(8)EXα=∑i=0∞αicα−iEX−ciwhere EX−ci is the ith central moment of the random variable X. If the reference point is put on the mean value of the random variable X, the ith central moment can be rewritten as
(9)EX−ci=∑k=0i−1i−kikci−kEXk,i=1,2,⋯.Equations (8) and (9) indicate that a fractional moment includes the information of a series of integer statistical moments, either central moments or raw moments.If a couple of fractional moments are available at hand, they would provide much more information than the same number of integer moments for identifying an unknown distribution type. In contrast, Gzyl and Tagliani [55] proved that 30 or more integer moments are required by Equations (8) and (9) to guarantee the accuracy of a fractional moment, where the infinite sequence of integer moments is demanded theoretically. This result is extremely important and gives a strong motivation to employ the concept of the fractional moment to reconstruct an unknown distribution type from a group of samples.The information entropy of a random variableX is a measure of the uncertainty containing by the random variable, which is defined as
(10)HX=∫fXxlogfXxdx,where fXx is the probability density function of the random variable X.The core step of the proposed method in this study is to reconstruct an expression for the unknown PDF of the structural responseG, given a sample of experimental data. The MaxEnt principle states that the most unbiased estimation of the unknown PDF is the one that maximizes Equation (10) with given constraints. In this study, the given constraints are the fractional moments determined from the experimental data. Therefore, by the MaxEnt principle, the PDF identification problem is converted into an optimization problem [42].
(11)maxHX,s.t.∫fXxdx=1,∫xαfXXdx=Mxα,where Mxα is the αth fractional moment estimated from a sample
(12)Mxα=1N∑i=1Nxiα≈EXα.Note thatxi are the random samples for the random variable X and N is the number of experimental points in the sample.Solving the optimization problem with the Lagrange multiplier method, one has the following generic form of the MaxEnt distribution [36, 42].
(13)f^Xx=exp−λ0exp−∑i=1mλixαi,where f^Xx denotes the MaxEnt PDF for fXx, λi are the Lagrange multipliers, αi are the orders for the fractional moments, and m is the number of available fractional moments, respectively. The parameter λ0 is a normalizing factor which makes sure that the integration of f^Xx over its support domain is 1 and satisfies the definition of PDF. Thus, from Equation (13), one has
(14)λ0=log∫exp−∑i=1mλixαi.Based on the concept of relative entropy (or the Kullback-Leibler distance), Zhang and Pandey proved that the following optimization problem is equivalent to Equatio (11) [36].
(15)minΙλ,α=log∫Xexp−∑i=1mλixαidx+∑i=1mλiMXαi.The minimization of the relative entropy means the minimization of the generalized distance between the true PDFfXx and the estimated PDF f^Xx, which is quantified by the information entropy. We will further reformulate the optimization problem in Equation (15) and derive an attractive property of it in the next subsection.


### 3.3. Computation of Maximum Entropy Distribution

Recalling the equivalent optimization problem in Equation (15), its objective function is furtherly rewritten as
(16)Ιλ,α=log∫Xexp−∑i=1mλixαidx+∑i=1mλiMXαi=log∫Xexp−∑i=1mλixαi−∑i=1mλiMXαidx.It is obvious that minimizingIλ,α is equivalent to minimize the following integral:
(17)∫Xexp−∑i=1mλixαi−∑i=1mλiMXαidx≡Qλ,α.The integral term is defined as a functionQλ,α for the purpose of simplification. Based on the definition of Qλ,α in Equation (17), the MaxEnt distribution can be rewritten as
(18)f^x=Qλ,α−1exp−∑i=1mλixαi−MXαi.This study examines the properties of the integral functionQλ,α firstly. The derivative of Qλ,α with respect to λi is written
(19)∂Q∂λi=∫X∂∂λiexp−∑i=1mλixαi−∑i=1mλiMXαidx=Qλ,α∫X−xαi−MXαiQλ,α−1exp−∑i=1mλixαi−∑i=1mλiMXαidx=Qλ,α∫X−xαi−MXαifXxdx.Recalling the definition of fractional moment, we have(20)∂Q∂λi=Qλ,α∫X−xαi−MXαifXxdx=Qλ,αMXαi−∫XxαifXxdx=0.The first derivative ofQλ,α with respect to αi is given by
(21)∂Q∂αi=∫X∂∂αiexp−∑i=1mλixαi−∑i=1mλiMXαidx=∫X−λiαixαi−1−∫X∂∂αixαifxdxexp−∑i=1mλixαi−∑i=1mλiMXαidx=Qλ,α∫X−λiαixαi−1−MXαi−1fXxdx=0.From Equations (20) and (21), we know that the solution point λ1,⋯,λm,α1,⋯,αm to Equation (17) is a saddle point of Qλ,α at least.Next, we examine the elements in the Hessian matrix ofQλ,α. There are four types elements in the Hessian matrix, i.e., ∂Q/∂λi∂λj, ∂Q/∂αi∂αj, ∂Q/∂λi∂αj, and ∂Q/∂αi∂λj which consist of the submatrices A, B, C, and D, respectively, as shown in
(22)H=ABCD.For the submatrixA, the elements are
(23)∂Q∂λi∂λj=Qλ,α∫Xxαi−MXαixαj−MXαjfXxdxIt is apparent submatrixA is a covariance-like and symmetric one. In a similar way, we have
(24)∂Q∂αi∂αj=αiαjλiλjQλ,α∫Xxαi−1−MXαi−1xαj−1−MXαj−1fXxdx.Equation (24) indicates that submatrix D is also a covariance-like and symmetric one. For an element in submatrix B, it has an expression
(25)∂Q∂λi∂αj=αjλjQλ,α∫Xxαi−MXαixαj−1−MXαj−1fXxdx.Similarly, an element in submatrixC has an expression
(26)∂Q∂αi∂λj=αiλiQλ,α∫Xxαi−1−MXαi−1xαj−MXαjfXxdx.It can be seen from Equations (25) and (26) that submatrices B and C are also covariance-like ones. An attracting property of submatrices B and C is that the columns of B are the rows of C, i.e., B is the transpose of C.According to above discussion, the Hessian matrix ofQλ,α is proven to be a covariance-like and symmetric one. As long as the investigated fXx is a PDF, the Hessian matrix is positive definite and of full rank. Thus, the solution to Equation (16) is the minimum point of Qλ,α and furtherly of Ιλ,α. Since the optimization problem in Equation(17), (15) is convex and has no constraint; it seems to provide an opportunity to solve it easily.


## 3.1. Overview

The new high-dimensional time-variant reliability analysis method has three stages. The first stage is to generate random samples for the input random variables and random processes using MCS sampling strategy and furtherly compute the structural responses. The second stage involves the approximation of the MaxEnt distribution. The third stage is related to the calculation of failure probability using the MaxEnt distribution obtained in the second stage. It is clear that the second stage is the key stage of the proposed method. The implementation procedure of the proposed method is shown in Figure1.Figure 1
Procedure of the proposed method.Stage 1 involves four steps. Step 1: discretize the time interval0,tf with a constant time step size Δt and obtain n+1 discrete time instant t0=0,t1,⋯,tn=tf. Step 2: decompose the input random processes Yt with a group of random variables Z. Step 3: generate N random samples for the input random variables and random processes. Step 4: substitute the random samples into the LSF to obtain N trajectories of structural response.Stage 2 is furtherly divided into three steps. Step 1: solve the optimization problem derived from the concept of relative entropy. Step 2: calculate the fractional moments from theN trajectories of structural response. Step 3: determine the MaxEnt distribution based on the solution to the optimization problem in Step 1. In the subsequent sections, we will discuss the details of Step 1 and 3 since these steps are theoretical fundamentals of the proposed method.Stage 3 is straightforward once the probability density function (PDF) of MaxEnt distribution is determined in Stage 2. The estimation of a cumulative failure probability is just a numerical integration over the failure domain.


## 3.2. Maximum Entropy Principle with Fractional Moments

The advantages of usage of fractional moments in the maximum entropy principle have been proven by several studies in the structural reliability community [36, 43, 45–50]. This subsection will give a brief review of the maximum entropy principle with fractional moments.Identifying the probability distribution with statistical moments has been investigated extensively in the literature. However, the use of integer statistical moments suffers from accuracy concern. Zhang and Pandey study the need of fractional moments instead of integer moments, which is evoked by the relationship between them [36].We start with the definition of the fractional binomial coefficientαi that is given by
(5)αi=αα−1⋯α−i+1ii−1⋯1,where α is a real number and i is an integer number. If applying the Taylor series expansion on a function of xα, one has
(6)xα=∑i=0∞αicα−ix−ci,where c is a reference point.Theαth-order fractional moment of a positive random variable X is defined as
(7)EXα=∫XxαfXxdx.Substituting Equation (6) into Equation (7), one has
(8)EXα=∑i=0∞αicα−iEX−ciwhere EX−ci is the ith central moment of the random variable X. If the reference point is put on the mean value of the random variable X, the ith central moment can be rewritten as
(9)EX−ci=∑k=0i−1i−kikci−kEXk,i=1,2,⋯.Equations (8) and (9) indicate that a fractional moment includes the information of a series of integer statistical moments, either central moments or raw moments.If a couple of fractional moments are available at hand, they would provide much more information than the same number of integer moments for identifying an unknown distribution type. In contrast, Gzyl and Tagliani [55] proved that 30 or more integer moments are required by Equations (8) and (9) to guarantee the accuracy of a fractional moment, where the infinite sequence of integer moments is demanded theoretically. This result is extremely important and gives a strong motivation to employ the concept of the fractional moment to reconstruct an unknown distribution type from a group of samples.The information entropy of a random variableX is a measure of the uncertainty containing by the random variable, which is defined as
(10)HX=∫fXxlogfXxdx,where fXx is the probability density function of the random variable X.The core step of the proposed method in this study is to reconstruct an expression for the unknown PDF of the structural responseG, given a sample of experimental data. The MaxEnt principle states that the most unbiased estimation of the unknown PDF is the one that maximizes Equation (10) with given constraints. In this study, the given constraints are the fractional moments determined from the experimental data. Therefore, by the MaxEnt principle, the PDF identification problem is converted into an optimization problem [42].
(11)maxHX,s.t.∫fXxdx=1,∫xαfXXdx=Mxα,where Mxα is the αth fractional moment estimated from a sample
(12)Mxα=1N∑i=1Nxiα≈EXα.Note thatxi are the random samples for the random variable X and N is the number of experimental points in the sample.Solving the optimization problem with the Lagrange multiplier method, one has the following generic form of the MaxEnt distribution [36, 42].
(13)f^Xx=exp−λ0exp−∑i=1mλixαi,where f^Xx denotes the MaxEnt PDF for fXx, λi are the Lagrange multipliers, αi are the orders for the fractional moments, and m is the number of available fractional moments, respectively. The parameter λ0 is a normalizing factor which makes sure that the integration of f^Xx over its support domain is 1 and satisfies the definition of PDF. Thus, from Equation (13), one has
(14)λ0=log∫exp−∑i=1mλixαi.Based on the concept of relative entropy (or the Kullback-Leibler distance), Zhang and Pandey proved that the following optimization problem is equivalent to Equatio (11) [36].
(15)minΙλ,α=log∫Xexp−∑i=1mλixαidx+∑i=1mλiMXαi.The minimization of the relative entropy means the minimization of the generalized distance between the true PDFfXx and the estimated PDF f^Xx, which is quantified by the information entropy. We will further reformulate the optimization problem in Equation (15) and derive an attractive property of it in the next subsection.


## 3.3. Computation of Maximum Entropy Distribution

Recalling the equivalent optimization problem in Equation (15), its objective function is furtherly rewritten as
(16)Ιλ,α=log∫Xexp−∑i=1mλixαidx+∑i=1mλiMXαi=log∫Xexp−∑i=1mλixαi−∑i=1mλiMXαidx.It is obvious that minimizingIλ,α is equivalent to minimize the following integral:
(17)∫Xexp−∑i=1mλixαi−∑i=1mλiMXαidx≡Qλ,α.The integral term is defined as a functionQλ,α for the purpose of simplification. Based on the definition of Qλ,α in Equation (17), the MaxEnt distribution can be rewritten as
(18)f^x=Qλ,α−1exp−∑i=1mλixαi−MXαi.This study examines the properties of the integral functionQλ,α firstly. The derivative of Qλ,α with respect to λi is written
(19)∂Q∂λi=∫X∂∂λiexp−∑i=1mλixαi−∑i=1mλiMXαidx=Qλ,α∫X−xαi−MXαiQλ,α−1exp−∑i=1mλixαi−∑i=1mλiMXαidx=Qλ,α∫X−xαi−MXαifXxdx.Recalling the definition of fractional moment, we have(20)∂Q∂λi=Qλ,α∫X−xαi−MXαifXxdx=Qλ,αMXαi−∫XxαifXxdx=0.The first derivative ofQλ,α with respect to αi is given by
(21)∂Q∂αi=∫X∂∂αiexp−∑i=1mλixαi−∑i=1mλiMXαidx=∫X−λiαixαi−1−∫X∂∂αixαifxdxexp−∑i=1mλixαi−∑i=1mλiMXαidx=Qλ,α∫X−λiαixαi−1−MXαi−1fXxdx=0.From Equations (20) and (21), we know that the solution point λ1,⋯,λm,α1,⋯,αm to Equation (17) is a saddle point of Qλ,α at least.Next, we examine the elements in the Hessian matrix ofQλ,α. There are four types elements in the Hessian matrix, i.e., ∂Q/∂λi∂λj, ∂Q/∂αi∂αj, ∂Q/∂λi∂αj, and ∂Q/∂αi∂λj which consist of the submatrices A, B, C, and D, respectively, as shown in
(22)H=ABCD.For the submatrixA, the elements are
(23)∂Q∂λi∂λj=Qλ,α∫Xxαi−MXαixαj−MXαjfXxdxIt is apparent submatrixA is a covariance-like and symmetric one. In a similar way, we have
(24)∂Q∂αi∂αj=αiαjλiλjQλ,α∫Xxαi−1−MXαi−1xαj−1−MXαj−1fXxdx.Equation (24) indicates that submatrix D is also a covariance-like and symmetric one. For an element in submatrix B, it has an expression
(25)∂Q∂λi∂αj=αjλjQλ,α∫Xxαi−MXαixαj−1−MXαj−1fXxdx.Similarly, an element in submatrixC has an expression
(26)∂Q∂αi∂λj=αiλiQλ,α∫Xxαi−1−MXαi−1xαj−MXαjfXxdx.It can be seen from Equations (25) and (26) that submatrices B and C are also covariance-like ones. An attracting property of submatrices B and C is that the columns of B are the rows of C, i.e., B is the transpose of C.According to above discussion, the Hessian matrix ofQλ,α is proven to be a covariance-like and symmetric one. As long as the investigated fXx is a PDF, the Hessian matrix is positive definite and of full rank. Thus, the solution to Equation (16) is the minimum point of Qλ,α and furtherly of Ιλ,α. Since the optimization problem in Equation(17), (15) is convex and has no constraint; it seems to provide an opportunity to solve it easily.


## 4. Computational Issues and Solutions

We have discussed the theoretical rationale of the proposed method in Section3. However, the technical rationale of the proposed method is more important for implementing this algorithm in practical applications.


### 4.1. Problem Reformulation

We have proved that the core optimization problem in Equation (17) is a convex one and there is no constraint on it. That is why most of studies have suggested that the optimization procedure can be carried out by the Nelder-Mead algorithm to search the solution λ1,⋯,λm,α1,⋯,αm [36, 45–50]. However, its technical implementation is not as perfect as the theoretical derivations, based on a lot of numerical implementation experience.The first reason is that the Nelder-Mead algorithm is only suitable for a small number of optimization variables. As the number of optimization variables increases, the performance of the Nelder-Mead algorithm becomes worse. The current study focuses on the high-dimensional time-variant reliability analysis; then, the number of fractional moments is set to be 3 or larger empirically. Thus, the total number of optimization variables is 6 or larger. We also encountered the phenomena of diversity due to different initial values of the optimization variables. This may be caused by the incapability of the Nelder-Mead algorithm.The second reason is that the estimate error of fractional moment increases with the order of moment, i.e.,α [36]. For this reason, the lower-order moments are used in the classical maximum entropy method. To overcome this serious drawback, the values of α should have limits during the optimization procedure. For example, α≤2 was suggested to obtain an unbiased estimation from experimental data [36] and α≤1 was used by Inverardi and Tagliani [42]. In this study, α≤1.5 is used. Additionally, the magnitude of the Lagrange multipliers should also have limits on their values to reduce the searching space and the numerical error associated to the numerical computation. We put the restriction of ±100 on the Lagrange multipliers.The third reason is the numerical error associated with the integral calculation involved in Equation (17). Furthermore, it affects the accuracy of λ0, PDF and CDF of the MaxEnt distribution, and further reliability analysis. The CDF of the MaxEnt distribution may be larger than 1 if the improper fractional exponents and Lagrange multipliers are obtained. A high-order global adaptive quadrature is used to perform the integral calculation in this study, which possesses a relatively high accuracy among the numerical integration algorithms in the literature.Except the above reasons, there is another aspect that is worth to be discussed since it is directly connected with the interval of integration and further the associated accuracy. That is the support of the MaxEnt distribution. Since the samples of structural responses are available when determining the support of the MaxEnt distribution, the maximum, minimum, mean, and standard deviation can be estimated from the samples easily. To avoid arbitrariness, the well-known6σ principle is suggested to determine the support of the MaxEnt distribution. The upper limit is set to be xmax+6s, where xmax and s are the maximum and standard deviation of the sample, respectively. The lower limit is set to be xmin−6s if xmin−6s is larger than 0, where is xmin the minimum observation in the sample. Otherwise, it is set to be 0 as the investigated response quantity has positive support which is required by the definition of fractional moment.Integrating the above discussion, the reformulated optimization problem has side constraints on the fractional exponents and the Lagrange multipliers. Thus, subset simulation optimization [56, 57] is employed to solve the reformulated optimization problem, instead of the Nelder-Mead algorithm in this study.


### 4.2. Parameter Estimation Using Subset Simulation Optimization

Subset simulation optimization is a heuristic algorithm for optimization problems, which is extended from a popular reliability method, i.e., subset simulation [56–58]. An attractive feature of subset simulation optimization is its ability for high-dimensional problems.Based on the discussion in subsection4.1, we have reformulated the optimization problem in the MaxEnt principle as an unconstrained one with side constraints on optimization variables. In subset simulation optimization, a truncated normal distribution is assigned to each optimization variable according to the associated side constraint. Take a moment order αi as an example, the PDF of the corresponding truncated normal distribution is given by [56–58]
(27)fαiαi=ϕαi−μi/σiΦαiu−μi/σi−Φαil−μi/σi,where αiu and αil are the upper and lower limits of the moment order, respectively; μi and σi are the mean value and standard deviation determined by the 3-sigma rule in reliability engineering; and ϕ⋅ and Φ⋅ are the PDF and cumulative distribution function of standard normal distribution, respectively.Since the side constraints are considered in the truncated normal distributions, there is no need to handle them during the optimization procedure. Along with the spirit of subset simulation optimization, subsets are generated to approach the optimal solution progressively. All subsetsFi have a generic form as
(28)Fi=Qλ,α≤bi,where bi are the thresholds for the objective function Qλ,αin Equation (17). If one arranges the thresholds in a descent order, i.e., b1≥b2≥b3⋯, the subsets have a nested structure
(29)F1⊃F2⊃F3⋯.Thus, the final subset in the sequence would be or include the minimum solution theoretically.In subset simulation optimization, the value of the objective functionQλ,α is the driving quantity for optimization, while the optimization procedure is achieved by the evaluation of the level probability pi that is an optimization parameter for subset simulation optimization. The conditional probability between two contiguous subsets is defined as the level probability and simulated by Markov Chain Monte Carlo simulation, i.e.,
(30)pi=PFi+1Fi≈1M∑i=MIQλ,α≥bi+1,where M is the number of candidate designs (or samples) belonging to the subset Fi, and I⋅ is an indicator function which is used to check a candidate design falling in Fi+1 or not. In this study, we adopted a strategy in which the value of pk is reduced as the number of optimization iterations (or simulation level) increases.The standard deviation of the samples in thei+1th simulation level is employed to check the convergence of the optimal searching. The estimator of standard deviation tends to zero while the searching process is approaching the global minimum. At the same time, we also impose a maximum number of simulation levels as another stopping criterion to avoid the overflow of the affordable computational effort.More details of subset simulation optimization can be referred to Refs. [56–58].


## 4.1. Problem Reformulation

We have proved that the core optimization problem in Equation (17) is a convex one and there is no constraint on it. That is why most of studies have suggested that the optimization procedure can be carried out by the Nelder-Mead algorithm to search the solution λ1,⋯,λm,α1,⋯,αm [36, 45–50]. However, its technical implementation is not as perfect as the theoretical derivations, based on a lot of numerical implementation experience.The first reason is that the Nelder-Mead algorithm is only suitable for a small number of optimization variables. As the number of optimization variables increases, the performance of the Nelder-Mead algorithm becomes worse. The current study focuses on the high-dimensional time-variant reliability analysis; then, the number of fractional moments is set to be 3 or larger empirically. Thus, the total number of optimization variables is 6 or larger. We also encountered the phenomena of diversity due to different initial values of the optimization variables. This may be caused by the incapability of the Nelder-Mead algorithm.The second reason is that the estimate error of fractional moment increases with the order of moment, i.e.,α [36]. For this reason, the lower-order moments are used in the classical maximum entropy method. To overcome this serious drawback, the values of α should have limits during the optimization procedure. For example, α≤2 was suggested to obtain an unbiased estimation from experimental data [36] and α≤1 was used by Inverardi and Tagliani [42]. In this study, α≤1.5 is used. Additionally, the magnitude of the Lagrange multipliers should also have limits on their values to reduce the searching space and the numerical error associated to the numerical computation. We put the restriction of ±100 on the Lagrange multipliers.The third reason is the numerical error associated with the integral calculation involved in Equation (17). Furthermore, it affects the accuracy of λ0, PDF and CDF of the MaxEnt distribution, and further reliability analysis. The CDF of the MaxEnt distribution may be larger than 1 if the improper fractional exponents and Lagrange multipliers are obtained. A high-order global adaptive quadrature is used to perform the integral calculation in this study, which possesses a relatively high accuracy among the numerical integration algorithms in the literature.Except the above reasons, there is another aspect that is worth to be discussed since it is directly connected with the interval of integration and further the associated accuracy. That is the support of the MaxEnt distribution. Since the samples of structural responses are available when determining the support of the MaxEnt distribution, the maximum, minimum, mean, and standard deviation can be estimated from the samples easily. To avoid arbitrariness, the well-known6σ principle is suggested to determine the support of the MaxEnt distribution. The upper limit is set to be xmax+6s, where xmax and s are the maximum and standard deviation of the sample, respectively. The lower limit is set to be xmin−6s if xmin−6s is larger than 0, where is xmin the minimum observation in the sample. Otherwise, it is set to be 0 as the investigated response quantity has positive support which is required by the definition of fractional moment.Integrating the above discussion, the reformulated optimization problem has side constraints on the fractional exponents and the Lagrange multipliers. Thus, subset simulation optimization [56, 57] is employed to solve the reformulated optimization problem, instead of the Nelder-Mead algorithm in this study.


## 4.2. Parameter Estimation Using Subset Simulation Optimization

Subset simulation optimization is a heuristic algorithm for optimization problems, which is extended from a popular reliability method, i.e., subset simulation [56–58]. An attractive feature of subset simulation optimization is its ability for high-dimensional problems.Based on the discussion in subsection4.1, we have reformulated the optimization problem in the MaxEnt principle as an unconstrained one with side constraints on optimization variables. In subset simulation optimization, a truncated normal distribution is assigned to each optimization variable according to the associated side constraint. Take a moment order αi as an example, the PDF of the corresponding truncated normal distribution is given by [56–58]
(27)fαiαi=ϕαi−μi/σiΦαiu−μi/σi−Φαil−μi/σi,where αiu and αil are the upper and lower limits of the moment order, respectively; μi and σi are the mean value and standard deviation determined by the 3-sigma rule in reliability engineering; and ϕ⋅ and Φ⋅ are the PDF and cumulative distribution function of standard normal distribution, respectively.Since the side constraints are considered in the truncated normal distributions, there is no need to handle them during the optimization procedure. Along with the spirit of subset simulation optimization, subsets are generated to approach the optimal solution progressively. All subsetsFi have a generic form as
(28)Fi=Qλ,α≤bi,where bi are the thresholds for the objective function Qλ,αin Equation (17). If one arranges the thresholds in a descent order, i.e., b1≥b2≥b3⋯, the subsets have a nested structure
(29)F1⊃F2⊃F3⋯.Thus, the final subset in the sequence would be or include the minimum solution theoretically.In subset simulation optimization, the value of the objective functionQλ,α is the driving quantity for optimization, while the optimization procedure is achieved by the evaluation of the level probability pi that is an optimization parameter for subset simulation optimization. The conditional probability between two contiguous subsets is defined as the level probability and simulated by Markov Chain Monte Carlo simulation, i.e.,
(30)pi=PFi+1Fi≈1M∑i=MIQλ,α≥bi+1,where M is the number of candidate designs (or samples) belonging to the subset Fi, and I⋅ is an indicator function which is used to check a candidate design falling in Fi+1 or not. In this study, we adopted a strategy in which the value of pk is reduced as the number of optimization iterations (or simulation level) increases.The standard deviation of the samples in thei+1th simulation level is employed to check the convergence of the optimal searching. The estimator of standard deviation tends to zero while the searching process is approaching the global minimum. At the same time, we also impose a maximum number of simulation levels as another stopping criterion to avoid the overflow of the affordable computational effort.More details of subset simulation optimization can be referred to Refs. [56–58].


## 5. Test Problems

Two dynamic reliability analysis problems are used to demonstrate the performance of the proposed method. The first problem involves a single-degree-of-freedom (SDOF) oscillator which has a linear behavior, while the second problem is a five-DOF structural system with nonlinear behavior.


### 5.1. A SDOF Oscillator

The investigated SDOF oscillator is shown in Figure2 with a white noise excitation Wt. The governing equation of this structural system is given by
(31)Y¨t+2ςωYt+ω2Yt=Wt,where the natural frequency is =7.85rad/s1.25Hz and the damping ratio is ς=2%, respectively. The spectral intensity of the white noise S takes 1. The following initial conditions are applied to the system:
(32)Y0=0,Ẏ0=0.Figure 2
The SDOF oscillator.The failure event of this system is defined as its maximum displacement exceeding a thresholdb=2.1 over a time interval 0,30 s
(33)F=maxk=1,⋯,n∣Ytk>b,where n is the total number of discrete time instants to represent the Gaussian white noise. Ref. [33] used a constant time step size Δt=0.002s for discretizing the time interval of interest, which leads to the total number of time instants be n=30/Δt+1=1501. Then, 1501 input standard normal random variables φk are employed to represent the input white noise, i.e.,
(34)Wtk=2πSΔtφk:k=1,⋯,n.To meet the definition of the fractional moment, the structural response associated to the failure event in Equation (33) is rewritten as
(35)Gn=gnX,Yt,t=maxk=1,⋯,nYtkb,where Gn is the normalized structural response. Then, the proposed method was applied to solve this problem with a sample size N=500, 1000, 1500, and 2000. For the purpose of comparison, an MCS with 106 samples was carried out to obtain a reference result.Figure3 shows the comparisons of the fractional moments and the relative error for the normalized response in Equation (35). Here, the order of fractional moments varies in the interval of −1.5,1.5 as we discussed in Section 4.1. The fractional moments calculated from different sample sizes almost overlap one another (Figure 3(a)). However, it is clear that the relative error is becoming small with the increase of sample size N (Figure 3(b)). Furthermore, the biggest relative error is 1% when the order of fractional moment takes -1.5, which is close to the used values reported by Xu and Wang [48] using a low-discrepancy sequence. Therefore, the relatively accurate fractional moments are employed in the proposed method by limiting the scale of the order exponents.Figure 3
Fractional moments and relative error of the normalized response (the SDOF oscillator).
(a)
Fractional moment(b)
Relative errorInverardi and Tagliani proposed a penalty termm/N, which is proportional to the order of the MaxEnt distribution m and inversely proportional to the sample size N, to avoid reconstructing a “too elaborate” distribution that cannot be justified by the given experimental data [42]. Table 1 lists the values of the objective function in Equation (17) while Table 2 gives the values of the objective function considering the penalty term. Six digitals must be kept in order to observe the difference among the results for a sample size in Table 1, while four digitals are required in Table 2. The value of the objective function increases as the number of fractional moments used in the MaxEnt distribution and decreases with the number of samples adopted for the estimation of fractional moments. Thus, m=3 is preferable for approximating the MaxEnt distribution of the normalized structural response, and 2000 samples provide the most accurate estimation of fractional moment and further reduce the objective function value.Table 1
The value of the objective function without penalty term.
N=500N=1000N=1500N=2000m=30.4067860.4087040.3982820.397276m=40.4067920.4087090.3982920.397279m=50.4067970.4087080.3982860.397281Table 2
The value of the objective function with penalty term.
N=500N=1000N=1500N=2000m=30.41280.41170.40030.3988m=40.41480.41270.40100.3993m=50.41680.41370.40160.3998Figures4 and 5 show the histograms, PDFs, and probability of exceedances (POEs) [36] obtained by the proposed method and MCS. Since all settings of parameters produce very similarly results, we only report the results for m=3 with N=500 (Figure 4) and 2000 (Figure 5). The term “Exp” refers to the experimental data (sample) used to estimate the fractional moments. The results of the proposed method agree with those of MCS, low to a very small probability (about 10-6). Thus, the proposed method possesses a good accuracy and efficiency for this investigated linear system.Figure 4
Comparison of histogram, PDF, and POE form=3 and N=500.
(a)
PDF(b)
.POE(c)
PDF(d)
POEFigure 5
Comparison of histogram, PDF, and POE form=3 and N=2000.
(a)
PDF(b)
POE(c)
PDF(d)
POEWe completed all calculations on a desktop PC with CPU@3.60 GHz and 16 GB RAM. Table3 reports the average computational time of the SDOF oscillator based on 30 runs of different combinations of N and m. Although it is hard to draw a common conclusion from these data, the average calculational time that are around 20~ 40 s is attractive for engineering optimization design.Table 3
The average computational time of the oscillator example (second).
N=500N=1000N=1500N=2000m=334.041.34.333.4m=423.737.418.230.3m=521.522.540.025.8


### 5.2. A Nonlinear Five-DOF Structural System

The second example is a five-DOF shear building adapted from Refs. [34, 51]. The structural system shown in Figure 6 is subjected to an external excitation. There are 25 input random variables to model uncertainties in the system properties, including masses, stiffness, sliding force stiffness ratios, and damping ratios. Table 4 lists the means of these 25 input random variables which are independently and normally distributed. Note that the damping terms ci in are given by
(36)ci=2ξimiki+ri,where ξi are the damping ratios, mi are the masses, and ki+ri are the stiffness, respectively. Except for the damping variables, the coefficient of variation of the other random property variables is set to be 0.1. The standard deviations of damping variables are set to be 0.01.Figure 6
Schematic of a five-story shear building.Table 4
Means of the random system properties.
StoryiMassmi (kg)Stiffnesski (N/m)Stiffnessri (N/m)Sliding force stiffness ratiosi/ki (m)Damping ratioξi120×10324×1062.4×1068×10−30.06220×10321×1062.1×1068×10−30.06320×10318×1061.8×1068×10−30.06420×10315×1061.5×1068×10−30.06520×10312×1061.2×1068×10−30.06Suppose that the external excitation to the structural system is a general Gaussian excitation acting on the horizontal ground. The Karhunen-Loève expansion is used to represent it with 200 independent standard normal variables. Therefore, there are total 225 input random variables which are beyond the capacity of FORM and surrogate model methods.The failure event is defined by the first-passage of the maximum displacement at the first story with a thresholdb=0.039m For the structural dynamic analysis, the time interval of 0.0s,20.0s is discretized by a time step size of 0.01 s.For this problem, 500 experimental points generated from MCS are employed to estimate the fractional moments of the normalized response. The comparisons of fractional moments and the corresponding relative errors are shown in Figure7. Similar conclusions can be drawn as for the SDOF oscillator, except that the magnitude of relative error grows to 2.75% for α=−1.5.Figure 7
Fractional moments and relative error of the normalized response (the nonlinear structural system).
(a)
Fractional moment(b)
Relative errorAdditionally,m=3, 4, and 5 are adopted in the proposed method. Table 5 lists the values of the objective function with penalty term. As the value of m increases, the values of the objective function also increase, which is consistent with the results of the SDOF oscillator. This indicates that the MaxEnt distribution with m=3 has a good balance of accuracy and model complexity. Table 6 lists the orders of fractional moments and Lagrange multipliers obtained by the proposed method.Table 5
The value of the objective function for the nonlinear structure.
m=3m=4m=5Equation(17)0.30730.30920.3108Table 6
Parameters for the MaxEnt distribution withm=3.
k0123αk0-0.49500.6545-1.1883λk7.1554-14.926512.48753.0281Based on the parameters in Table6, we reconstructed the distribution of the normalized response and compared the proposed method and MCS in Figure 8. The PDF curve of the proposed method matches the histogram of the experimental data closely (Figure 8(a)), while we observed a dispersion between it and the histogram of MCS samples in Figure 8(c). Similarly, the POE curve of the proposed method agrees with that of experimental data (Figure 8(b)), whereas it has a deviation from the curve estimated from MCS samples (Figure 8(d)). For the interval of 10−3,100, the POE of the proposed method has a good agreement with that of MCS. The failure probability is estimated as 1.3400×10−6 by MCS with 107 samples when the threshold level b is 0.039 m (corresponding to the case of the normalized threshold level of 1). The estimation of failure probability yielded by the proposed method is 3.5027×10−5. Therefore, the proposed method overestimated the failure probability for a large threshold level. Compared with the investigated linear system, the overestimations may cause by the nonlinear behavior of the current structural system.Figure 8
Comparison of histogram, PDF, and POE for the nonlinear structure.
(a)
PDF(b)
POE(c)
PDF(d)
POETable7 reports the average computational time of the nonlinear structural system based on 30 runs of different m. Similar conclusion can be drawn from this table as that of the SDOF oscillator example.Table 7
The value of the objective function for the nonlinear structure.
m=3m=4m=544.127.934.5


## 5.1. A SDOF Oscillator

The investigated SDOF oscillator is shown in Figure2 with a white noise excitation Wt. The governing equation of this structural system is given by
(31)Y¨t+2ςωYt+ω2Yt=Wt,where the natural frequency is =7.85rad/s1.25Hz and the damping ratio is ς=2%, respectively. The spectral intensity of the white noise S takes 1. The following initial conditions are applied to the system:
(32)Y0=0,Ẏ0=0.Figure 2
The SDOF oscillator.The failure event of this system is defined as its maximum displacement exceeding a thresholdb=2.1 over a time interval 0,30 s
(33)F=maxk=1,⋯,n∣Ytk>b,where n is the total number of discrete time instants to represent the Gaussian white noise. Ref. [33] used a constant time step size Δt=0.002s for discretizing the time interval of interest, which leads to the total number of time instants be n=30/Δt+1=1501. Then, 1501 input standard normal random variables φk are employed to represent the input white noise, i.e.,
(34)Wtk=2πSΔtφk:k=1,⋯,n.To meet the definition of the fractional moment, the structural response associated to the failure event in Equation (33) is rewritten as
(35)Gn=gnX,Yt,t=maxk=1,⋯,nYtkb,where Gn is the normalized structural response. Then, the proposed method was applied to solve this problem with a sample size N=500, 1000, 1500, and 2000. For the purpose of comparison, an MCS with 106 samples was carried out to obtain a reference result.Figure3 shows the comparisons of the fractional moments and the relative error for the normalized response in Equation (35). Here, the order of fractional moments varies in the interval of −1.5,1.5 as we discussed in Section 4.1. The fractional moments calculated from different sample sizes almost overlap one another (Figure 3(a)). However, it is clear that the relative error is becoming small with the increase of sample size N (Figure 3(b)). Furthermore, the biggest relative error is 1% when the order of fractional moment takes -1.5, which is close to the used values reported by Xu and Wang [48] using a low-discrepancy sequence. Therefore, the relatively accurate fractional moments are employed in the proposed method by limiting the scale of the order exponents.Figure 3
Fractional moments and relative error of the normalized response (the SDOF oscillator).
(a)
Fractional moment(b)
Relative errorInverardi and Tagliani proposed a penalty termm/N, which is proportional to the order of the MaxEnt distribution m and inversely proportional to the sample size N, to avoid reconstructing a “too elaborate” distribution that cannot be justified by the given experimental data [42]. Table 1 lists the values of the objective function in Equation (17) while Table 2 gives the values of the objective function considering the penalty term. Six digitals must be kept in order to observe the difference among the results for a sample size in Table 1, while four digitals are required in Table 2. The value of the objective function increases as the number of fractional moments used in the MaxEnt distribution and decreases with the number of samples adopted for the estimation of fractional moments. Thus, m=3 is preferable for approximating the MaxEnt distribution of the normalized structural response, and 2000 samples provide the most accurate estimation of fractional moment and further reduce the objective function value.Table 1
The value of the objective function without penalty term.
N=500N=1000N=1500N=2000m=30.4067860.4087040.3982820.397276m=40.4067920.4087090.3982920.397279m=50.4067970.4087080.3982860.397281Table 2
The value of the objective function with penalty term.
N=500N=1000N=1500N=2000m=30.41280.41170.40030.3988m=40.41480.41270.40100.3993m=50.41680.41370.40160.3998Figures4 and 5 show the histograms, PDFs, and probability of exceedances (POEs) [36] obtained by the proposed method and MCS. Since all settings of parameters produce very similarly results, we only report the results for m=3 with N=500 (Figure 4) and 2000 (Figure 5). The term “Exp” refers to the experimental data (sample) used to estimate the fractional moments. The results of the proposed method agree with those of MCS, low to a very small probability (about 10-6). Thus, the proposed method possesses a good accuracy and efficiency for this investigated linear system.Figure 4
Comparison of histogram, PDF, and POE form=3 and N=500.
(a)
PDF(b)
.POE(c)
PDF(d)
POEFigure 5
Comparison of histogram, PDF, and POE form=3 and N=2000.
(a)
PDF(b)
POE(c)
PDF(d)
POEWe completed all calculations on a desktop PC with CPU@3.60 GHz and 16 GB RAM. Table3 reports the average computational time of the SDOF oscillator based on 30 runs of different combinations of N and m. Although it is hard to draw a common conclusion from these data, the average calculational time that are around 20~ 40 s is attractive for engineering optimization design.Table 3
The average computational time of the oscillator example (second).
N=500N=1000N=1500N=2000m=334.041.34.333.4m=423.737.418.230.3m=521.522.540.025.8


## 5.2. A Nonlinear Five-DOF Structural System

The second example is a five-DOF shear building adapted from Refs. [34, 51]. The structural system shown in Figure 6 is subjected to an external excitation. There are 25 input random variables to model uncertainties in the system properties, including masses, stiffness, sliding force stiffness ratios, and damping ratios. Table 4 lists the means of these 25 input random variables which are independently and normally distributed. Note that the damping terms ci in are given by
(36)ci=2ξimiki+ri,where ξi are the damping ratios, mi are the masses, and ki+ri are the stiffness, respectively. Except for the damping variables, the coefficient of variation of the other random property variables is set to be 0.1. The standard deviations of damping variables are set to be 0.01.Figure 6
Schematic of a five-story shear building.Table 4
Means of the random system properties.
StoryiMassmi (kg)Stiffnesski (N/m)Stiffnessri (N/m)Sliding force stiffness ratiosi/ki (m)Damping ratioξi120×10324×1062.4×1068×10−30.06220×10321×1062.1×1068×10−30.06320×10318×1061.8×1068×10−30.06420×10315×1061.5×1068×10−30.06520×10312×1061.2×1068×10−30.06Suppose that the external excitation to the structural system is a general Gaussian excitation acting on the horizontal ground. The Karhunen-Loève expansion is used to represent it with 200 independent standard normal variables. Therefore, there are total 225 input random variables which are beyond the capacity of FORM and surrogate model methods.The failure event is defined by the first-passage of the maximum displacement at the first story with a thresholdb=0.039m For the structural dynamic analysis, the time interval of 0.0s,20.0s is discretized by a time step size of 0.01 s.For this problem, 500 experimental points generated from MCS are employed to estimate the fractional moments of the normalized response. The comparisons of fractional moments and the corresponding relative errors are shown in Figure7. Similar conclusions can be drawn as for the SDOF oscillator, except that the magnitude of relative error grows to 2.75% for α=−1.5.Figure 7
Fractional moments and relative error of the normalized response (the nonlinear structural system).
(a)
Fractional moment(b)
Relative errorAdditionally,m=3, 4, and 5 are adopted in the proposed method. Table 5 lists the values of the objective function with penalty term. As the value of m increases, the values of the objective function also increase, which is consistent with the results of the SDOF oscillator. This indicates that the MaxEnt distribution with m=3 has a good balance of accuracy and model complexity. Table 6 lists the orders of fractional moments and Lagrange multipliers obtained by the proposed method.Table 5
The value of the objective function for the nonlinear structure.
m=3m=4m=5Equation(17)0.30730.30920.3108Table 6
Parameters for the MaxEnt distribution withm=3.
k0123αk0-0.49500.6545-1.1883λk7.1554-14.926512.48753.0281Based on the parameters in Table6, we reconstructed the distribution of the normalized response and compared the proposed method and MCS in Figure 8. The PDF curve of the proposed method matches the histogram of the experimental data closely (Figure 8(a)), while we observed a dispersion between it and the histogram of MCS samples in Figure 8(c). Similarly, the POE curve of the proposed method agrees with that of experimental data (Figure 8(b)), whereas it has a deviation from the curve estimated from MCS samples (Figure 8(d)). For the interval of 10−3,100, the POE of the proposed method has a good agreement with that of MCS. The failure probability is estimated as 1.3400×10−6 by MCS with 107 samples when the threshold level b is 0.039 m (corresponding to the case of the normalized threshold level of 1). The estimation of failure probability yielded by the proposed method is 3.5027×10−5. Therefore, the proposed method overestimated the failure probability for a large threshold level. Compared with the investigated linear system, the overestimations may cause by the nonlinear behavior of the current structural system.Figure 8
Comparison of histogram, PDF, and POE for the nonlinear structure.
(a)
PDF(b)
POE(c)
PDF(d)
POETable7 reports the average computational time of the nonlinear structural system based on 30 runs of different m. Similar conclusion can be drawn from this table as that of the SDOF oscillator example.Table 7
The value of the objective function for the nonlinear structure.
m=3m=4m=544.127.934.5


## 6. Conclusions

In this study, a maximum entropy distribution method is proposed to perform structural reliability analysis for high-dimensional problems with random process input. With the development of the maximum entropy distribution with fractional moments, we proved that the solution to a functionQλ,α is the global one to the maximum relative entropy between the true PDF and the approximate one. Thus, the corresponding optimization problem is converted into a convex one with a global solution. To estimate the fractional moments efficiently, a set of experimental points are generated for this purpose, which provides the estimator of a fractional moment for any order. Converting the common optimization problem into an equivalent one can avoid the trapping in a locally optimal solution. Furthermore, we reexamined the associated computational issues of the target optimization problem and reformulated it with side constraints on the optimization variables, i.e., the orders of fractional moments and the corresponding Lagrange multipliers. Subset simulation optimization is used to solve the new optimization problem.Two benchmark problems with high-dimensional inputs show that the proposed method has both good accuracy and efficiency. If the limit state function, however, has a strong nonlinear behavior with respect to inputs, the maximum entropy distribution may not approximate the tail region accurately, and in this case, the proposed method may result in a low efficiency or accuracy, or both.Future studies will involve developing efficiently sampling methods for the estimation of fractional moments for high-dimensional nonlinear problems.



---

*Source: 6612864-2022-03-19.xml*